Merge pull request #179 from kaipartmann/params

Enhancements to point parameters
kaipartmann · Oct 14, 2024 · 766cfc1 · 766cfc1
2 parents a309e1c + 2ab0148
commit 766cfc1
Show file tree

Hide file tree

Showing 23 changed files with 252 additions and 200 deletions.
diff --git a/src/Peridynamics.jl b/src/Peridynamics.jl
@@ -71,6 +71,7 @@ include("auxiliary/function_arguments.jl")
 include("auxiliary/io.jl")
 include("auxiliary/logs.jl")
 include("auxiliary/mpi.jl")
+include("auxiliary/errors.jl")
 
 include("physics/boundary_conditions.jl")
 include("physics/initial_conditions.jl")
@@ -92,9 +93,9 @@ include("discretization/body_chunk.jl")
 
 include("core/job.jl")
 include("core/submit.jl")
-include("core/parameter_handler.jl")
 include("core/systems.jl")
-include("core/materials.jl")
+include("core/parameters.jl")
+include("core/parameter_handler.jl")
 include("core/storages.jl")
 include("core/time_solvers.jl")
 include("core/halo_exchange.jl")

diff --git a/src/auxiliary/errors.jl b/src/auxiliary/errors.jl
@@ -0,0 +1,19 @@
+struct InterfaceError <: Exception
+    type::DataType
+    func::String
+    function InterfaceError(_type::T, _func::F) where {T,F}
+        func = string(_func)
+        type = isa(_type, DataType) ? _type : T
+        return new(type, func)
+    end
+end
+
+function Base.showerror(io::IO, e::InterfaceError)
+    print(io, "interface method not correctly defined!")
+    print(io, "\n  type:    ")
+    printstyled(io, string(e.type); bold=true, color=:red)
+    print(io, "\n  method:  ")
+    printstyled(io, string(e.func); bold=true, color=:red)
+    println(io)
+    return nothing
+end
diff --git a/src/core/materials.jl → src/core/parameters.jl b/src/core/materials.jl → src/core/parameters.jl
@@ -1,22 +1,29 @@
 function point_param_type(mat::AbstractMaterial)
-    throw(MethodError(point_param_type, mat))
+    return throw(InterfaceError(mat, "point_param_type"))
 end
 
 function material_type(params::AbstractPointParameters)
-    throw(MethodError(material_type, params))
+    return throw(InterfaceError(params, "material_type"))
 end
 
 function get_point_params(mat::AbstractMaterial, ::Dict{Symbol,Any})
-    throw(MethodError(get_point_params, mat))
+    return throw(InterfaceError(mat, "get_point_params"))
+end
+
+function required_point_parameters(mat::Type{<:AbstractMaterial})
+    return throw(InterfaceError(mat, "required_point_parameters"))
+end
+
+function allowed_material_kwargs(mat::AbstractMaterial)
+    return throw(InterfaceError(mat, "allowed_material_kwargs"))
 end
 
 macro params(material, params)
     macrocheck_input_material(material)
     macrocheck_input_params(params)
-    local _checks = quote
+    local _pre_checks = quote
         Peridynamics.typecheck_material($(esc(material)))
-        #TODO: make this check material dependent -> overloadable for own types!
-        Peridynamics.typecheck_params($(esc(params)))
+        Peridynamics.typecheck_params($(esc(material)), $(esc(params)))
     end
     local _pointparam_type = quote
         Peridynamics.point_param_type(::$(esc(material))) = $(esc(params))
@@ -29,7 +36,12 @@ macro params(material, params)
             return $(esc(params))(m, p)
         end
     end
-    return Expr(:block, _checks, _pointparam_type, _material_type, _get_pointparams)
+    local _post_checks = quote
+        Peridynamics.constructor_check($(esc(material)), $(esc(params)))
+    end
+    exprs = Expr(:block, _pre_checks, _pointparam_type, _material_type, _get_pointparams,
+                 _post_checks)
+    return exprs
 end
 
 function macrocheck_input_material(material)
@@ -53,17 +65,24 @@ function typecheck_material(::Type{Material}) where {Material}
     return nothing
 end
 
-function typecheck_params(::Type{Param}) where {Param}
+function typecheck_params(::Type{Material}, ::Type{Param}) where {Material,Param}
     if !(Param <: AbstractPointParameters)
-        msg = "$Param is not a valid point parameter type!\n"
+        msg = "$Param is not a subtype of AbstractPointParameters!\n"
         throw(ArgumentError(msg))
     end
     parameters = fieldnames(Param)
-    for req_param in required_point_parameters()
+    for req_param in required_point_parameters(Material)
         if !in(req_param, parameters)
             msg = "required parameter $req_param not found in $(Param)!\n"
             error(msg)
         end
     end
     return nothing
 end
+
+function constructor_check(::Type{Material}, ::Type{Param}) where {Material,Param}
+    if !hasmethod(Param, Tuple{Material,Dict{Symbol,Any}})
+        throw(InterfaceError(Param, "$(Param)(::$(Material), ::Dict{Symbol,Any})"))
+    end
+    return nothing
+end
diff --git a/src/core/systems.jl b/src/core/systems.jl
@@ -1,17 +1,7 @@
 function system_type(mat::AbstractMaterial)
-    throw(MethodError(system_type, mat))
+    return throw(InterfaceError(mat, "system_type"))
 end
 
-# function get_system_type(mat::AbstractMaterial)
-#     Sys = system_type(mat)
-#     FullSys = full_system_type(Sys, mat)
-#     return FullSys
-# end
-
-# function full_system_type(::Type{System}, mat::AbstractMaterial) where {System}
-#     return system_type(mat)
-# end
-
 function get_system(::AbstractBody{M}, ::PointDecomposition, ::Int) where {M}
     msg = "system for material $M not specified!\n"
     return error(msg)
@@ -21,36 +11,3 @@ function log_system(options::AbstractJobOptions, dh::AbstractDataHandler)
     log_system(system_type(dh), options, dh)
     return nothing
 end
-
-# macro system(material, system)
-#     macrocheck_input_material(material)
-#     macrocheck_input_system(system)
-#     local _checks = quote
-#         Peridynamics.typecheck_material($(esc(material)))
-#         Peridynamics.typecheck_system($(esc(system)))
-#     end
-#     local _system_type = quote
-#         Peridynamics.system_type(::$(esc(material))) = $(esc(system))
-#     end
-#     local _get_system = quote
-#         function Peridynamics.get_system(body::Peridynamics.AbstractBody{$(esc(material))},
-#                                          args...)
-#             return $(esc(system))(body, args...)
-#         end
-#     end
-#     return Expr(:block, _checks, _system_type, _get_system)
-# end
-
-# function macrocheck_input_system(system)
-#     system isa Symbol && return nothing
-#     (system isa Expr && system.head === :.) && return nothing
-#     return throw(ArgumentError("argument `$system` is not a valid system input!\n"))
-# end
-
-# function typecheck_system(::Type{System}) where {System}
-#     if !(System <: AbstractSystem)
-#         msg = "$System is not a valid system type!\n"
-#         throw(ArgumentError(msg))
-#     end
-#     return nothing
-# end
diff --git a/src/discretization/bond_system.jl b/src/discretization/bond_system.jl
@@ -251,3 +251,15 @@ function calc_n_bonds(dh::AbstractMPIBodyDataHandler)
     n_bonds = MPI.Reduce(length(dh.chunk.system.bonds), MPI.SUM, mpi_comm())
     return n_bonds
 end
+
+function required_point_parameters(::Type{<:AbstractBondSystemMaterial})
+    return (:δ, :rho, elasticity_parameters()...)
+end
+
+function get_required_point_parameters(::AbstractBondSystemMaterial, p::Dict{Symbol,Any})
+    return (; get_horizon(p)..., get_density(p)..., get_elastic_params(p)...)
+end
+
+function allowed_material_kwargs(::AbstractBondSystemMaterial)
+    return (discretization_kwargs()..., elasticity_kwargs()..., fracture_kwargs()...)
+end
diff --git a/src/discretization/interaction_system.jl b/src/discretization/interaction_system.jl
@@ -81,6 +81,16 @@ function check_interaction_system_compat(::AbstractInteractionSystemMaterial)
     return nothing
 end
 
+function get_c2(params::AbstractPointParameters)
+    hasproperty(params, :C2) || return 0.0
+    return params.C2
+end
+
+function get_c3(params::AbstractPointParameters)
+    hasproperty(params, :C3) || return 0.0
+    return params.C3
+end
+
 function has_two_nis(body::AbstractBody)
     for params in body.point_params
         get_c2(params) ≈ 0 || return true
@@ -415,3 +425,54 @@ function calc_n_interactions(dh::AbstractMPIBodyDataHandler)
     n_three_nis = MPI.Reduce(length(dh.chunk.system.three_nis), MPI.SUM, mpi_comm())
     return n_one_nis, n_two_nis, n_three_nis
 end
+
+function required_point_parameters(::Type{<:AbstractInteractionSystemMaterial})
+    return (:δ, :rho, elasticity_parameters()..., :C1, :C2, :C3)
+end
+
+function get_required_point_parameters(mat::AbstractInteractionSystemMaterial,
+                                       p::Dict{Symbol,Any})
+    params = (; get_horizon(p)..., get_density(p)..., get_elastic_params(p)...)
+    params = (; params..., get_interaction_parameters(mat, p, params)...)
+    return params
+end
+
+function get_interaction_parameters(mat::AbstractInteractionSystemMaterial,
+                                    p::Dict{Symbol,Any}, params)
+    (; δ, μ, λ) = params
+    if haskey(p, :C1)
+        C1::Float64 = float(p[:C1])
+    else
+        C1 = 0.0
+    end
+
+    if haskey(p, :C2)
+        C2::Float64 = float(p[:C2])
+    else
+        C2 = 0.0
+    end
+
+    if haskey(p, :C3)
+        C3::Float64 = float(p[:C3])
+    else
+        C3 = 0.0
+    end
+
+    if C1 ≈ 0 && C2 ≈ 0 && C3 ≈ 0
+        C1 = 30 / π * μ / δ^4
+        C2 = 0.0
+        C3 = 32 / π^4 * (λ - μ) / δ^12
+    else
+        msg = "interaction parameters for $(typeof(mat)) specified manually!\n"
+        msg *= "Be careful when adjusting these parameters to avoid unexpected outcomes!"
+        @mpiroot @warn msg
+    end
+
+    return (; C1, C2, C3)
+end
+
+function allowed_material_kwargs(::AbstractInteractionSystemMaterial)
+    kwargs = (discretization_kwargs()..., elasticity_kwargs()..., fracture_kwargs()...,
+              :C1, :C2, :C3)
+    return kwargs
+end
diff --git a/src/physics/bond_based.jl b/src/physics/bond_based.jl
@@ -76,19 +76,15 @@ struct BBPointParameters <: AbstractPointParameters
     bc::Float64
 end
 
-function BBPointParameters(::BBMaterial, p::Dict{Symbol,Any})
-    δ = get_horizon(p)
-    rho = get_density(p)
-    if haskey(p, :nu)
-        msg = "keyword `nu` is not allowed for BBMaterial!\n"
-        msg *= "Bond-based peridynamics has a limitation on the possion ratio.\n"
-        msg *= "Therefore, when using BBMaterial, `nu` is hardcoded as 1/4.\n"
+function BBPointParameters(mat::BBMaterial, p::Dict{Symbol,Any})
+    if haskey(p, :nu) && !isapprox(p[:nu], 0.25)
+        msg = "Bond-based peridynamics has a limitation on the Poisson's ratio!\n"
+        msg *= "With BBMaterial, no other values than nu=0.25 are allowed!\n"
         throw(ArgumentError(msg))
-    else
-        p[:nu] = 0.25
     end
-    E, nu, G, K, λ, μ = get_elastic_params(p)
-    Gc, εc = get_frac_params(p, δ, K)
+    p[:nu] = 0.25
+    (; δ, rho, E, nu, G, K, λ, μ) = get_required_point_parameters(mat, p)
+    (; Gc, εc) = get_frac_params(p, δ, K)
     bc = 18 * K / (π * δ^4) # bond constant
     return BBPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc)
 end

diff --git a/src/physics/continuum_kinematics_inspired.jl b/src/physics/continuum_kinematics_inspired.jl
@@ -72,58 +72,14 @@ struct CKIPointParameters <: AbstractPointParameters
     C3::Float64
 end
 
-function CKIPointParameters(::CKIMaterial, p::Dict{Symbol,Any})
-    δ = get_horizon(p)
-    rho = get_density(p)
-    E, nu, G, K, λ, μ = get_elastic_params(p)
-    Gc, εc = get_frac_params(p, δ, K)
-    C1, C2, C3 = cki_parameters(p, δ, λ, μ)
+function CKIPointParameters(mat::CKIMaterial, p::Dict{Symbol,Any})
+    (; δ, rho, E, nu, G, K, λ, μ, C1, C2, C3) = get_required_point_parameters(mat, p)
+    (; Gc, εc) = get_frac_params(p, δ, K)
     return CKIPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, C1, C2, C3)
 end
 
-function cki_parameters(p::Dict{Symbol,Any}, δ, λ, μ)
-    if haskey(p, :C1)
-        C1::Float64 = float(p[:C1])
-    else
-        C1 = 0.0
-    end
-
-    if haskey(p, :C2)
-        C2::Float64 = float(p[:C2])
-    else
-        C2 = 0.0
-    end
-
-    if haskey(p, :C3)
-        C3::Float64 = float(p[:C3])
-    else
-        C3 = 0.0
-    end
-
-    if C1 ≈ 0 && C2 ≈ 0 && C3 ≈ 0
-        C1 = 30 / π * μ / δ^4
-        C2 = 0.0
-        C3 = 32 / π^4 * (λ - μ) / δ^12
-    else
-        msg = "parameters for CKIMaterial specified manually!\n"
-        msg *= "Be careful when adjusting these parameters to avoid unexpected outcomes!"
-        @mpiroot @warn msg
-    end
-
-    return C1, C2, C3
-end
-
-function allowed_material_kwargs(::CKIMaterial)
-    return (DEFAULT_POINT_KWARGS..., :C1, :C2, :C3)
-end
-
 @params CKIMaterial CKIPointParameters
 
-@inline get_c2(params::CKIPointParameters) = params.C2
-@inline get_c2(params) = 0.0
-@inline get_c3(params::CKIPointParameters) = params.C3
-@inline get_c3(params) = 0.0
-
 struct CKIVerletStorage <: AbstractStorage
     position::Matrix{Float64}
     displacement::Matrix{Float64}
@@ -297,9 +253,8 @@ function force_density_point_two_ni!(storage::CKIStorage, system::InteractionSys
         aijky = Δxijz * Δxikx - Δxijx * Δxikz
         aijkz = Δxijx * Δxiky - Δxijy * Δxikx
         surface = sqrt(aijkx * aijkx + aijky * aijky + aijkz * aijkz)
-        if surface == 0 # to avoid divide by zero error for failed interactions
-            surface = 1e-40
-        end
+        # avoid to divide by zero error for failed interactions
+        isapprox(surface, 0; atol=eps()) && (surface = 1e-40)
         failure = storage.one_ni_active[oni_j_id] * storage.one_ni_active[oni_k_id]
         params_j = get_params(paramhandler, j)
         params_k = get_params(paramhandler, k)
@@ -345,9 +300,8 @@ function force_density_point_three_ni!(storage::CKIStorage, system::InteractionS
         aijky = Δxijz * Δxikx - Δxijx * Δxikz
         aijkz = Δxijx * Δxiky - Δxijy * Δxikx
         volume = aijkx * Δxilx + aijky * Δxily + aijkz * Δxilz
-        if volume == 0 # avoid to divide by zero error for failed interactions
-            volume = 1e-40
-        end
+        # avoid to divide by zero error for failed interactions
+        isapprox(volume, 0; atol=eps()) && (volume = 1e-40)
         abs_volume = abs(volume)
         failure = storage.one_ni_active[oni_j_id] * storage.one_ni_active[oni_k_id] *
                   storage.one_ni_active[oni_l_id]

diff --git a/src/physics/correspondence.jl b/src/physics/correspondence.jl
@@ -94,11 +94,9 @@ struct NOSBPointParameters <: AbstractPointParameters
     bc::Float64
 end
 
-function NOSBPointParameters(::NOSBMaterial, p::Dict{Symbol,Any})
-    δ = get_horizon(p)
-    rho = get_density(p)
-    E, nu, G, K, λ, μ = get_elastic_params(p)
-    Gc, εc = get_frac_params(p, δ, K)
+function NOSBPointParameters(mat::NOSBMaterial, p::Dict{Symbol,Any})
+    (; δ, rho, E, nu, G, K, λ, μ) = get_required_point_parameters(mat, p)
+    (; Gc, εc) = get_frac_params(p, δ, K)
     bc = 18 * K / (π * δ^4) # bond constant
     return NOSBPointParameters(δ, rho, E, nu, G, K, λ, μ, Gc, εc, bc)
 end