Multi-Material support for mechanical problems (#203)

* Solid multi-material type
* Add some pretty printing

---------

Co-authored-by: termi-official <9196588+termi-official@users.noreply.github.com>

Author: termi-official

src/discretization/fem.jl

@@ -153,16 +153,44 @@ function semidiscretize(split::ReactionDiffusionSplit{<:MonodomainModel}, discre
     return semidiscrete_ode
 end
 
-function semidiscretize(model::QuasiStaticModel, discretization::FiniteElementDiscretization, mesh::AbstractGrid)
+# Solid mechanics semidiscretize interface
+function semidiscretize_register_subdomains!(dh, lvh, model, material_model::AbstractMaterialModel, discretization::FiniteElementDiscretization)
     sym = model.displacement_symbol
     ipc = _get_interpolation_from_discretization(discretization, sym)
     qrc = _get_quadrature_from_discretization(discretization, sym)
-    dh = DofHandler(mesh)
-    lvh = InternalVariableHandler(mesh)
     for name in discretization.subdomains
         add_subdomain!(dh, name, [ApproximationDescriptor(sym, ipc)])
-        add_subdomain!(lvh, name, gather_internal_variable_infos(model.material_model), qrc, dh)
+        add_subdomain!(lvh, name, gather_internal_variable_infos(material_model), qrc, dh)
+    end
+end
+
+function semidiscretize_register_subdomains!(dh, lvh, model, material_models::MultiMaterialModel, discretization::FiniteElementDiscretization)
+    semidiscretize_register_subdomains_multi!(dh, lvh, model, material_models.materials, material_models.domains, material_models.domain_names, discretization)
+end
+@unroll function semidiscretize_register_subdomains_multi!(dh, lvh, model, material_models, domains, domain_names, discretization)
+    sym = model.displacement_symbol
+    ipc = _get_interpolation_from_discretization(discretization, sym)
+    qrc = _get_quadrature_from_discretization(discretization, sym)
+    idx = 1
+    @unroll for material_model in material_models
+        add_subdomain!(dh,  domain_names[idx], [ApproximationDescriptor(sym, ipc)])
+        add_subdomain!(lvh, domain_names[idx], gather_internal_variable_infos(material_model), qrc, dh)
+        idx += 1
     end
+end
+
+
+function semidiscretize(model::QuasiStaticModel, discretization::FiniteElementDiscretization, mesh::AbstractGrid)
+    sym = model.displacement_symbol
+    ipc = _get_interpolation_from_discretization(discretization, sym)
+    qrc = _get_quadrature_from_discretization(discretization, sym)
+    dh = DofHandler(mesh)
+    lvh = InternalVariableHandler(mesh)
+    # for name in discretization.subdomains
+    #     add_subdomain!(dh, name, [ApproximationDescriptor(sym, ipc)])
+    #     add_subdomain!(lvh, name, gather_internal_variable_infos(model.material_model), qrc, dh)
+    # end
+    semidiscretize_register_subdomains!(dh, lvh, model, model.material_model, discretization)
     close!(dh)
     close!(lvh)
 

src/discretization/operator.jl

@@ -161,6 +161,12 @@ function AssembledNonlinearOperator(integrator::NonlinearIntegrator, dh::Abstrac
     )
 end
 
+function Base.show(io::IO, cache::AssembledNonlinearOperator)
+    println(io, "AssembledNonlinearOperator:")
+    Base.show(io, typeof(cache.integrator))
+    Base.show(io, MIME"text/plain"(), cache.dh)
+end
+
 getJ(op::AssembledNonlinearOperator) = op.J
 
 function update_linearization!(op::AssembledNonlinearOperator, u::AbstractVector, time)

src/mesh/simple_meshes.jl

@@ -208,6 +208,7 @@ end
 @inline Ferrite.getcells(mesh::SimpleMesh, setname::String) = Ferrite.getcells(mesh.grid, setname)
 @inline Ferrite.getcelltype(mesh::SimpleMesh) = Ferrite.getcelltype(mesh.grid)
 @inline Ferrite.getcelltype(mesh::SimpleMesh, i::Int) = Ferrite.getcelltype(mesh.grid, i)
+@inline Ferrite.getcellset(mesh::SimpleMesh, name::String) = Ferrite.getcellset(mesh.grid, name)
 @inline Ferrite.getfacetset(mesh::SimpleMesh, name::String) = Ferrite.getfacetset(mesh.grid, name)
 @inline Ferrite.getnodeset(mesh::SimpleMesh, name::String) = Ferrite.getnodeset(mesh.grid, name)
 

src/modeling/core/nonlinear.jl

@@ -1,5 +1,5 @@
 @doc raw"""
-    NonlinearIntegrator{CoefficientType}
+    NonlinearIntegrator
 
 Represents the integrand a the nonlinear form over some function space.
 """

src/modeling/functions.jl

@@ -93,8 +93,7 @@ end
 
 solution_size(f::QuasiStaticFunction) = ndofs(f.dh)+ndofs(f.lvh)
 function local_function_size(f::QuasiStaticFunction)
-    length(f.lvh.dh.subdofhandlers) == 0 && return 0
-    return sum(Ferrite.n_components.(f.lvh.dh.subdofhandlers[1].field_interpolations); init=0)
+    error("Local function size of QuasiStaticFunction can vary!")
 end
 function default_initial_condition!(u::AbstractVector, f::QuasiStaticFunction)
     fill!(u, 0.0)
@@ -105,10 +104,12 @@ function default_initial_condition!(u::AbstractVector, f::QuasiStaticFunction)
         qr = getquadraturerule(f.integrator.qrc, sdh)
         # ivsize_per_qp = sum(sdh.field_n_components; init=0) # FIXME broken...
         ivsize_per_qp = sum(Ferrite.n_components.(sdh.field_interpolations); init=0)
+        ivsize_per_qp == 0 && continue
+        material_model = get_material_model(f, sdh)
         for cell in CellIterator(sdh)
             for qp in QuadratureIterator(qr)
                 q = @view uq[offset:(offset+ivsize_per_qp-1)]
-                default_initial_state!(q, f.integrator.volume_model.material_model)
+                default_initial_state!(q, material_model)
                 offset += ivsize_per_qp
             end
         end
@@ -117,3 +118,17 @@ end
 
 gather_internal_variable_infos(model::QuasiStaticModel) = gather_internal_variable_infos(model.material_model)
 gather_internal_variable_infos(model::AbstractMaterialModel) = InternalVariableInfo[]
+
+@unroll function __get_material_model_multi(materials, domains, sdh)
+    idx = 1
+    @unroll for material ∈ materials
+        if first(domains[idx]) ∈ sdh.cellset
+            return material
+        end
+        idx += 1
+    end
+    error("MultiDomainIntegrator is broken: Requested to construct an internal cache for a SubDofHandler which is not associated with the integrator.")
+end
+__get_material_model(model::MultiMaterialModel, sdh) = __get_material_model_multi(model.materials, model.domains, sdh)
+__get_material_model(model::AbstractMaterialModel, sdh) = model
+get_material_model(f::QuasiStaticFunction, sdh) = __get_material_model(f.integrator.volume_model.material_model, sdh)

src/modeling/solid/elements.jl

@@ -117,16 +117,60 @@ function assemble_element!(residualₑ::AbstractVector, uₑ::AbstractVector, ge
     end
 end
 
+struct MultiMaterialModel{MaterialTuple} <: AbstractMaterialModel
+    materials::MaterialTuple
+    domains::Vector{OrderedSet{Int}} # These must match the subdofhandler sets and hence be disjoint
+    domain_names::Vector{String}     # These must match the subdofhandler sets and hence be disjoint
+    function MultiMaterialModel(materials::MaterialTuple, subdomain_names::Vector{String}, mesh::AbstractGrid) where MaterialTuple
+        return new{MaterialTuple}(materials, [getcellset(mesh, subdomain_name) for subdomain_name in subdomain_names], subdomain_names)
+    end
+end
+
+function setup_quasistatic_element_cache(material_model::MultiMaterialModel, qr::QuadratureRule, sdh::SubDofHandler, cv::CellValues)
+    return setup_quasistatic_element_cache_multi(material_model.materials, material_model.domains, qr, sdh, cv)
+end
+@unroll function setup_quasistatic_element_cache_multi(materials::Tuple, domains::Vector, qr::QuadratureRule, sdh::SubDofHandler, cv::CellValues)
+    idx = 1
+    @unroll for material ∈ materials
+        if first(domains[idx]) ∈ sdh.cellset
+            return QuasiStaticElementCache(
+                material,
+                setup_coefficient_cache(material, qr, sdh),
+                setup_internal_cache(material, qr, sdh),
+                cv
+            )
+        end
+        idx += 1
+    end
+    error("MultiDomainIntegrator is broken: Requested to construct an element cache for a SubDofHandler which is not associated with the integrator.")
+end
+function setup_quasistatic_element_cache(material_model::AbstractMaterialModel, qr::QuadratureRule, sdh::SubDofHandler, cv::CellValues)
+    return QuasiStaticElementCache(
+        material_model,
+        setup_coefficient_cache(material_model, qr, sdh),
+        setup_internal_cache(material_model, qr, sdh),
+        cv
+    )
+end
 function setup_element_cache(model::QuasiStaticModel, qr::QuadratureRule, sdh::SubDofHandler)
     @assert length(sdh.dh.field_names) == 1 "Support for multiple fields not yet implemented."
     field_name = first(sdh.dh.field_names)
     ip         = Ferrite.getfieldinterpolation(sdh, field_name)
     ip_geo     = geometric_subdomain_interpolation(sdh)
     cv = CellValues(qr, ip, ip_geo)
-    return QuasiStaticElementCache(
-        model.material_model,
-        setup_coefficient_cache(model.material_model, qr, sdh),
-        setup_internal_cache(model.material_model, qr, sdh),
-        cv
-    )
+    return setup_quasistatic_element_cache(model.material_model, qr, sdh, cv)
+end
+
+@unroll function setup_internal_cache_multi(materials::Tuple, domains::Vector, qr::QuadratureRule, sdh::SubDofHandler)
+    idx = 1
+    @unroll for material ∈ materials
+        if first(domains[idx]) ∈ sdh.cellset
+            return setup_internal_cache(material, qr, sdh)
+        end
+        idx += 1
+    end
+    error("MultiDomainIntegrator is broken: Requested to construct an internal cache for a SubDofHandler which is not associated with the integrator.")
+end
+function setup_internal_cache(model::MultiMaterialModel, qr::QuadratureRule, sdh::SubDofHandler)
+    return setup_internal_cache_multi(model.materials, model.domains, qr, sdh)
 end

src/modeling/solid/materials.jl

@@ -251,9 +251,9 @@ function _solve_local_sarcomere_dQdF(dQdλ, dλdF, λ, F, coefficients, active_t
 end
 
 function solve_local_constraint(F::Tensor{2,dim}, coefficients, material_model::ActiveStressModel, state_cache::GenericFirstOrderRateIndependentCondensationMaterialStateCache, geometry_cache, qp, time) where dim
+    Qflat, Qprevflat = _query_local_state(state_cache, geometry_cache, qp)
     # Early out if any of the previous local solves failed
     if state_cache.local_solver_cache.retcode ∉ (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success)
-        Qflat, _ = _query_local_state(state_cache, geometry_cache, qp)
         return Qflat, zero(Tensor{4,dim,Float64,4^dim})
     end
 
@@ -310,11 +310,9 @@ function solve_local_constraint(F::Tensor{2,dim}, coefficients, material_model::
         return Q, J
     end
 
-    Qflat, Qprevflat = _query_local_state(state_cache, geometry_cache, qp)
     Q, J = solve_internal_timestep(material_model, state_cache, λ, dλdt, Qflat, Qprevflat)
     # Abort if local solve failed
     if state_cache.local_solver_cache.retcode ∉ (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success)
-        Qflat, _ = _query_local_state(state_cache, geometry_cache, qp)
         return Qflat, zero(Tensor{4,dim,Float64,4^dim})
     end
     Qflat .= Q
@@ -336,13 +334,14 @@ function solve_local_constraint(F::Tensor{2,dim}, coefficients, material_model::
 end
 
 # Some debug materials
-Base.@kwdef struct LinearMaxwellMaterial{T} <: AbstractMaterialModel
+Base.@kwdef struct LinearMaxwellMaterial{T, sdim} <: AbstractMaterialModel
     E₀::T
     E₁::T
     μ::T
     η₁::T
     ν::T
 end
+LinearMaxwellMaterial(E₀::T, Eₗ::T, μ::T, η₁::T, ν::T) where T = LinearMaxwellMaterial{T,3}(E₀, Eₗ, μ, η₁, ν)
 
 local_function_size(model::QuasiStaticModel) = local_function_size(model.material_model)
 function local_function_size(model::AbstractMaterialModel)
@@ -452,7 +451,7 @@ function stress_function(material::LinearMaxwellMaterial, ε, coefficients, ε
     return E₀ * ℂ ⊡ ε + E₁ * ℂ ⊡ (ε - εᵛ)
 end
 
-function stress_and_tangent(material_model::LinearMaxwellMaterial, F::Tensor{2}, coefficients, εᵛ)
+function stress_and_tangent(material_model::LinearMaxwellMaterial, F::Tensor{2}, coefficients, εᵛ::SymmetricTensor{2})
     ε = symmetric(F - one(F))
     ∂σ∂ε, σ = Tensors.gradient(ε->stress_function(material_model, ε, coefficients, εᵛ), ε, :all)
     return σ, ∂σ∂ε
@@ -463,9 +462,13 @@ function setup_coefficient_cache(m::LinearMaxwellMaterial, qr::QuadratureRule, s
 end
 
 function setup_internal_cache(material_model::LinearMaxwellMaterial, qr::QuadratureRule, sdh::SubDofHandler)
-    return nothing # FIXME what should we do here? :)
+    return EmptyRateIndependentCondensationMaterialStateCache()
 end
 
-function gather_internal_variable_infos(model::LinearMaxwellMaterial)
-    return InternalVariableInfo(:εᵛ, 6) # TODO dimension info
+function gather_internal_variable_infos(model::LinearMaxwellMaterial{T, sdim}) where {T, sdim}
+    if sdim == 3
+        return InternalVariableInfo(:εᵛ, 6)
+    else
+        return InternalVariableInfo(:εᵛ, 4)
+    end
 end

src/solver/nonlinear/multilevel_newton_raphson.jl

@@ -12,6 +12,30 @@ Base.@kwdef mutable struct GenericLocalNonlinearSolverCache{JacobianType, Residu
     retcode::SciMLBase.ReturnCode.T = SciMLBase.ReturnCode.Default
 end
 
+function Base.show(io::IO, cache::GenericLocalNonlinearSolverCache)
+    println(io, "NewtonRaphsonSolverCache:")
+    Base.show(io, cache.params)
+    println(io, "J=$(typeof(cache.J)) with size $(size(cache.J))")
+    println(io, "R=$(typeof(cache.residual)) with size $(size(cache.residual))")
+    println(io, "C=$(typeof(cache.rhs_corrector)) with size $(size(cache.rhs_corrector))")
+    println(io, "outer_tol=$(cache.outer_tol)")
+    println(io, "status=$(cache.retcode)")
+end
+
+function check_local_solve_covergence(local_solver_cache::GenericLocalNonlinearSolverCache)
+    return local_solver_cache.retcode ∉ (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success)
+end
+function check_local_solve_covergence(local_solver_cache::Tuple)
+    return any(check_local_solve_covergence.(local_solver_cache))
+end
+
+function set_local_solver_tol(local_solver_cache::GenericLocalNonlinearSolverCache, tol)
+    local_solver_cache.outer_tol = tol
+end
+function set_local_solver_tol(local_solver_cache::Tuple, tol)
+    set_local_solver_tol.(local_solver_cache, tol)
+end
+
 """
     MultilevelNewtonRaphsonSolver{T}
 
@@ -29,6 +53,18 @@ struct MultiLevelNewtonRaphsonSolverCache{gCacheType, lCacheType} <: AbstractNon
     local_solver_cache::lCacheType
 end
 
+function Base.show(io::IO, cache::MultiLevelNewtonRaphsonSolverCache)
+    println(io, "MultiLevelNewtonRaphsonSolverCache:")
+    Base.show(io, cache.global_solver_cache)
+    if cache.local_solver_cache isa Tuple
+        for local_solver_cache in cache.local_solver_cache
+            Base.show(io, local_solver_cache)
+        end
+    else
+        Base.show(io, cache.local_solver_cache)
+    end
+end
+
 function nlsolve!(u::AbstractVector, f::AbstractSemidiscreteFunction, mlcache::MultiLevelNewtonRaphsonSolverCache, t)
     cache = mlcache.global_solver_cache
 
@@ -39,21 +75,21 @@ function nlsolve!(u::AbstractVector, f::AbstractSemidiscreteFunction, mlcache::M
     residualnormprev = 0.0
     Θ1prev = length(Θks) > 0 ? first(Θks) : 0.0
     resize!(Θks, 0)
-    mlcache.local_solver_cache.outer_tol = Inf
+    set_local_solver_tol(mlcache.local_solver_cache, Inf)
     while true
         cache.iter += 1
         residual .= 0.0
         @timeit_debug "update operator" update_linearization!(op, residual, u, t)
         # Check if local solve failed
-        if mlcache.local_solver_cache.retcode ∉ (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success)
+        if check_local_solve_covergence(mlcache.local_solver_cache)
             @debug "Some local newton did not converge. Aborting. ||r|| = $residualnorm" _group=:nlsolve
             return false
         end
         @timeit_debug "elimination" eliminate_constraints_from_linearization!(cache, f)
         linear_solver_cache.isfresh = true # Notify linear solver that we touched the system matrix
 
         residualnorm = residual_norm(cache, f)
-        mlcache.local_solver_cache.outer_tol = residualnorm
+        set_local_solver_tol(mlcache.local_solver_cache, residualnorm^2)
         if residualnorm < cache.parameters.tol && cache.iter > 1 # Do at least two iterations to get a sane convergence estimate
             break
         elseif cache.iter > cache.parameters.max_iter

src/solver/nonlinear/newton_raphson.jl

@@ -27,6 +27,12 @@ mutable struct NewtonRaphsonSolverCache{OpType, ResidualType, T, NewtonType <: N
     iter::Int
 end
 
+function Base.show(io::IO, cache::NewtonRaphsonSolverCache)
+    println(io, "NewtonRaphsonSolverCache:")
+    Base.show(io, cache.parameters)
+    Base.show(io, cache.op)
+end
+
 function setup_solver_cache(f::AbstractSemidiscreteFunction, solver::NewtonRaphsonSolver{T}) where {T}
     @unpack inner_solver = solver
     op = setup_operator(f, solver)