Adding Material Models

This guide describes how to add a new material constitutive model to Sentinel. Each model is a singleton type that drives dispatch for element integration, gradient computation, and assembly.

AbstractMaterialModel Interface

All material models subtype AbstractMaterialModel. Models are identified by an integer ID in runfiles (the Fortran mtrlmodel flag) and mapped to Julia types via material_model.

To register a new model, define a singleton struct and add entries to the dispatch functions:

struct MyNewModel <: AbstractMaterialModel end

# Integer mapping (used by runfile parser)
material_model(8) = MyNewModel()
model_index(::MyNewModel) = 8

# DOF structure
has_pressure_dofs(::MyNewModel) = true          # or false
pressure_dof_type(::MyNewModel) = :elemental    # :elemental, :nodal, or :none
uses_fiber_directions(::MyNewModel) = false      # true if DTI data needed

See the Material Models API page for full docstrings of AbstractMaterialModel, material_model, model_index, has_pressure_dofs, pressure_dof_type, and uses_fiber_directions.

Element Kernel

Implement element_stiffness! for the new model. This is the core physics routine that computes the complex element stiffness matrix Ke for a single element.

function element_stiffness!(Ke::Matrix{ComplexF64},
                             cv_disp::CellValues,
                             cv_press,              # CellValues or nothing
                             ::MyNewModel,
                             props_at_gp::Vector{GaussPointMaterial},
                             ω::Float64)
    # Loop over quadrature points, accumulate into Ke
    for q in 1:getnquadpoints(cv_disp)
        dΩ = getdetJdV(cv_disp, q)
        μ  = props_at_gp[q].μ
        # ... integrate B^T C B dΩ ...
    end
    return Ke
end

Requirements:

• Ke must be complex symmetric (Ke = Ke^T, NOT Hermitian). This is a fundamental property of the viscoelastic formulation with complex moduli.

• For hex27 elements with elemental pressure: Ke is 85 x 85 (81 displacement DOFs + 4 pressure DOFs). The displacement block occupies rows/columns 1:81, and the pressure block occupies rows/columns 82:85.

• For hex27 pure-displacement models: Ke is 81 x 81.

• Ke is zeroed before entry by the assembly loop; accumulate with +=.

• Material properties at each Gauss point are provided via GaussPointMaterial, which carries the complex shear modulus, bulk modulus, full stiffness tensor, fiber direction, and density.

See the Element Kernels API page for full docstrings of element_stiffness! and GaussPointMaterial.

Gradient Computation

Implement compute_gradient! for the new model. The adjoint method gives the gradient of the objective function with respect to material parameter theta:

$$\frac{\partial \mathrm{\Phi }}{\partial \theta }=\text{Re}\left[{\bar{\lambda }}^T\frac{\partial K}{\partial \theta }u\right]$$

where lambda is the adjoint displacement (from the adjoint solve) and u is the forward displacement.

The gradient routine loops over elements and Gauss points, accumulating contributions into the flat gradient vector grad.rvalue. For each material property, implement a private _accumulate_*_gradient! method:

function compute_gradient!(grad::MaterialGradient,
                            dh::DofHandler, cv_disp::CellValues,
                            calc::Displacement, adjnt::Displacement,
                            model::MyNewModel,
                            material::Material,
                            gp2mtrs::Vector{GP2Mtr},
                            meshes::Vector{MaterialMesh},
                            ω::Float64; numdispsets::Int=1)
    # For each element, Gauss point, and displacement set:
    #   1. Evaluate material properties at the Gauss point
    #   2. Gather local forward (u) and adjoint (lambda) displacements
    #   3. Call _accumulate_*_gradient! for each active property
end

The _accumulate_*_gradient! methods compute Re(conj(lambda)^T * dK/d(theta) * u) for a specific material parameter and scatter the result into grad.rvalue weighted by the material mesh basis function at the current Gauss point.

Assembly Integration

The global assembly function assemble_stiffness! dispatches to the element kernel based on the model type. No changes to the assembly loop are needed for new models – the dispatch is automatic as long as:

1. _n_press_per_el(::MyNewModel) returns the correct count (0 or 4)

2. _n_dofs_local(::MyNewModel) returns the correct element DOF count (81 or 85)

3. element_stiffness! is implemented for the model

The assembly loop handles the global DOF mapping, including the elemental pressure DOF indexing scheme (ndofs(dh) + 4*(el-1) + p for pressure DOF p of element el).

Testing

Every new model should include three levels of testing:

1. Forward solve patch test: Apply a linear displacement field as Dirichlet boundary conditions. The forward solver must reproduce this field exactly (to machine precision) in the interior. This validates the element kernel integration.

# Example: patch test for a unit cube with linear u = [x, 0, 0]
@test maximum(abs.(u_calc - u_exact)) < 1e-12

2. Finite-difference gradient check: Compare the adjoint gradient against central finite differences for each material property. Use a tolerance of approximately 5% to account for the finite-difference approximation error.

# Central differences: (Phi(theta+h) - Phi(theta-h)) / (2h)
fd_grad = (obj_plus - obj_minus) / (2 * h)
@test abs(adjoint_grad - fd_grad) / abs(fd_grad) < 0.05

3. Integration test: Run a full conjugate gradient optimization on a synthetic problem (known ground truth). Verify that the optimizer converges and recovers the true material properties within acceptable tolerance.

Existing Models

ID	Julia Type	Description	Pressure DOFs	Fiber Dirs
1	IsotropicIncompressible	Rayleigh-damped isotropic, incompressible	4 elemental	No
2	Orthotropic	Rayleigh-damped orthotropic	4 elemental	Yes
3	IsotropicCompressible	Rayleigh-damped isotropic, compressible	None	No
4	Poroelastic	Poroelastic with fluid pressure (Tet4)	4 nodal	No
5	TransverseIsotropicV1	Transverse isotropic, formulation 1	4 elemental	Yes
6	TransverseIsotropicV2	Transverse isotropic, formulation 2	4 elemental	Yes
7	GeneralizedAnisotropic	Generalized anisotropic	4 elemental	Yes

Models 5, 6, and 7 share the same assembly kernel (genanisoforwardassemble in the original Fortran) with different parameterizations of the 6x6 stiffness tensor.