Merge branch 'main' into l1-gs-smoother

Author: Abdelrahman912

.github/workflows/ci.yml

@@ -11,10 +11,9 @@ jobs:
     name: Julia ${{ matrix.julia-version }} - ${{ matrix.os }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
     timeout-minutes: 40
-    continue-on-error: ${{ matrix.julia-version == 'nightly' }}
     strategy:
       matrix:
-        julia-version: ['1', 'nightly']
+        julia-version: ['1', '1.10']
         os: ['ubuntu-latest']
         include:
           # - os: windows-latest

docs/src/api-reference/models.md

@@ -38,7 +38,7 @@ PressureFieldBC
 ## Solid Mechanics
 
 ```@docs
-StructuralModel
+QuasiStaticModel
 ExtendedHillModel
 GeneralizedHillModel
 ActiveStressModel

docs/src/literate-tutorials/cm01_simple-active-stress.jl

@@ -94,7 +94,7 @@ active_stress_model = ActiveStressModel(
 weak_boundary_conditions = (NormalSpringBC(1.0, "Epicardium"),)
 
 # We finalize the mechanical model by assigning a symbol to identify the unknown solution field and connect the active stress model with the weak boundary conditions.
-mechanical_model = StructuralModel(:displacement, active_stress_model, weak_boundary_conditions)
+mechanical_model = QuasiStaticModel(:displacement, active_stress_model, weak_boundary_conditions)
 
 # !!! tip
 #     A full list of all models can be found in the [API reference](@ref models-api).

docs/src/literate-tutorials/cm03_3d0d-coupling.jl

@@ -286,7 +286,7 @@ active_stress_model = ActiveStressModel(
     microstructure,
 )
 weak_boundary_conditions = (NormalSpringBC(1.0, "Epicardium"),)
-solid_model = StructuralModel(:displacement, active_stress_model, weak_boundary_conditions);
+solid_model = QuasiStaticModel(:displacement, active_stress_model, weak_boundary_conditions);
 
 # The solid model is now couple with the circuit model by adding a Lagrange multipliers constraining the 3D chamber volume to match the chamber volume in the 0D model.
 p3D = LVc.p3D

src/Thunderbolt.jl

@@ -21,7 +21,7 @@ include("solver/operator_splitting.jl")
 solution_size(f::GenericSplitFunction) = OS.function_size(f)
 
 @reexport using Ferrite
-import Ferrite: AbstractDofHandler, AbstractGrid, AbstractRefShape, AbstractCell, get_grid
+import Ferrite: AbstractDofHandler, AbstractGrid, AbstractRefShape, AbstractCell, get_grid, get_coordinate_eltype
 import Ferrite: vertices, edges, faces, sortedge, sortface
 import Ferrite: get_coordinate_type, getspatialdim
 import Ferrite: reference_shape_value
@@ -125,14 +125,14 @@ export
     # Generic models
     ODEProblem,
     TransientDiffusionModel,
-    TransientDiffusionFunction,
+    AffineODEFunction,
     # Local API
     PointwiseODEProblem,
     PointwiseODEFunction,
     # Mechanics
-    StructuralModel,
+    QuasiStaticModel,
     QuasiStaticProblem,
-    QuasiStaticNonlinearFunction,
+    QuasiStaticFunction,
     PK1Model,
     PrestressedMechanicalModel,
     # Passive material models

src/discretization/fem.jl

@@ -1,10 +1,13 @@
 """
 Descriptor for a finite element discretization of a part of a PDE over some subdomain.
+    
+!!! note
+    The current implementation is restricted to Bubnov-Galerkin methods. Petrov-Galerkin support will come in the future.
 """
 struct FiniteElementDiscretization
     """
     """
-    interpolations::Dict{Symbol, InterpolationCollection}
+    interpolations::Dict{Symbol}#, Union{<:InterpolationCollection, Pair{<:InterpolationCollection, <:QuadratureRuleCollection}}}
     """
     """
     dbcs::Vector{Dirichlet}
@@ -13,17 +16,42 @@ struct FiniteElementDiscretization
     subdomains::Vector{String}
     """
     """
-    function FiniteElementDiscretization(ips::Dict{Symbol, <: InterpolationCollection}, dbcs::Vector{Dirichlet} = Dirichlet[], subdomains::Vector{String} = [""])
-        new(ips, dbcs, subdomains)
+    mass_qrc::Union{<:QuadratureRuleCollection,Nothing} # TODO maybe an "extras" field should be used instead :)
+    """
+    """
+    function FiniteElementDiscretization(ips::Dict{Symbol}, dbcs::Vector{Dirichlet} = Dirichlet[], subdomains::Vector{String} = [""], mass_qrc = nothing)
+        new(ips, dbcs, subdomains, mass_qrc)
     end
 end
 
+_extract_ipc(ipc::InterpolationCollection) = ipc
+_extract_ipc(p::Pair{<:InterpolationCollection, <:QuadratureRuleCollection}) = first(p)
+
+function _extract_qrc(ipc::InterpolationCollection)
+    ansatzorder = getorder(ipc)
+    return QuadratureRuleCollection(max(2ansatzorder-1,2))
+end
+_extract_qrc(p::Pair{<:InterpolationCollection, <:QuadratureRuleCollection}) = last(p)
+
 # Internal utility with proper error message
 function _get_interpolation_from_discretization(disc::FiniteElementDiscretization, sym::Symbol)
     if !haskey(disc.interpolations, sym)
         error("Finite element discretization does not have an interpolation for $sym. Available symbols: $(collect(keys(disc.interpolations))).")
     end
-    return disc.interpolations[sym]
+    return _extract_ipc(disc.interpolations[sym])
+end
+function _get_quadrature_from_discretization(disc::FiniteElementDiscretization, sym::Symbol)
+    if !haskey(disc.interpolations, sym)
+        error("Finite element discretization does not have an interpolation for $sym. Available symbols: $(collect(keys(disc.interpolations))).")
+    end
+    return _extract_qrc(disc.interpolations[sym])
+end
+function _get_facet_quadrature_from_discretization(disc::FiniteElementDiscretization, sym::Symbol)
+    if !haskey(disc.interpolations, sym)
+        error("Finite element discretization does not have an interpolation for $sym. Available symbols: $(collect(keys(disc.interpolations))).")
+    end
+    intorder = getorder(_extract_ipc(disc.interpolations[sym]))
+    return FacetQuadratureRuleCollection(intorder)
 end
 
 semidiscretize(::CoupledModel, discretization, mesh::AbstractGrid) = @error "No implementation for the generic discretization of coupled problems available yet."
@@ -33,22 +61,34 @@ function semidiscretize(model::TransientDiffusionModel, discretization::FiniteEl
 
     sym = model.solution_variable_symbol
     ipc = _get_interpolation_from_discretization(discretization, sym)
+    qrc = _get_quadrature_from_discretization(discretization, sym)
     dh = DofHandler(mesh)
     for name in discretization.subdomains
         add_subdomain!(dh, name, [ApproximationDescriptor(sym, ipc)])
     end
     close!(dh)
 
-    return TransientDiffusionFunction(
-        model.κ,
-        model.source,
-        dh
+    T = get_coordinate_eltype(get_grid(dh))
+    return AffineODEFunction(
+        BilinearMassIntegrator(
+            ConstantCoefficient(T(1.0)),
+            discretization.mass_qrc === nothing ? qrc : mass_qrc, # Allow e.g. mass lumping for explicit integrators.
+            sym,
+        ),
+        BilinearDiffusionIntegrator(
+            model.κ,
+            qrc,
+            sym,
+        ),
+        model.source, # TODO qrc for source term
+        dh,
     )
 end
 
 function semidiscretize(model::SteadyDiffusionModel, discretization::FiniteElementDiscretization, mesh::AbstractGrid)
     sym = model.solution_variable_symbol
     ipc = _get_interpolation_from_discretization(discretization, sym)
+    qrc = _get_quadrature_from_discretization(discretization, sym)
     dh = DofHandler(mesh)
     for name in discretization.subdomains
         add_subdomain!(dh, name, [ApproximationDescriptor(sym, ipc)])
@@ -61,11 +101,15 @@ function semidiscretize(model::SteadyDiffusionModel, discretization::FiniteEleme
     end
     close!(ch)
 
-    return SteadyDiffusionFunction(
-        model.κ,
-        model.source,
+    return AffineSteadyStateFunction(
+        BilinearDiffusionIntegrator(
+            model.κ,
+            qrc,
+            sym,
+        ),
+        model.source, # TODO qrc for source term
         dh,
-        ch
+        ch,
     )
 end
 
@@ -109,9 +153,10 @@ function semidiscretize(split::ReactionDiffusionSplit{<:MonodomainModel}, discre
     return semidiscrete_ode
 end
 
-function semidiscretize(model::StructuralModel{<:QuasiStaticModel}, discretization::FiniteElementDiscretization, mesh::AbstractGrid)
+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)
     for name in discretization.subdomains
         add_subdomain!(dh, name, [ApproximationDescriptor(sym, ipc)])
@@ -124,11 +169,16 @@ function semidiscretize(model::StructuralModel{<:QuasiStaticModel}, discretizati
     end
     close!(ch)
 
-    semidiscrete_problem = QuasiStaticNonlinearFunction(
+    semidiscrete_problem = QuasiStaticFunction(
         dh,
         ch,
-        model.mechanical_model,
-        model.face_models
+        NonlinearIntegrator(
+            model,
+            model.face_models,
+            [sym],
+            qrc,
+            _get_facet_quadrature_from_discretization(discretization, sym),
+        ),
     )
 
     return semidiscrete_problem