Forward Solver

The forward solver assembles the global stiffness matrix, applies boundary conditions, and solves the linear system $\mathbf{K}\mathbf{u}=\mathbf{f}$.

DOF Handler Setup

Sentinel.setup_hex27_dofhandler Function

setup_hex27_dofhandler(grid) -> DofHandler

Create a Ferrite DofHandler with a 3-component vector displacement field on a hex27 mesh. Each of the 27 nodes per element gets 3 DOFs (u_x, u_y, u_z).

Returns a closed DofHandler ready for DOF queries and assembly iteration.

source

Stiffness Assembly

The global system $\mathbf{K}\mathbf{u}=\mathbf{f}$ is assembled from element contributions. For mixed formulations (Models 1, 2, 5, 6, 7), the stiffness matrix includes both displacement and pressure DOFs.

Sentinel.allocate_stiffness Function

allocate_stiffness(dh::DofHandler, model::AbstractMaterialModel) -> SparseMatrixCSC{ComplexF64}

Preallocate the global complex stiffness matrix with the correct sparsity pattern for the given DofHandler and material model.

For mixed (incompressible) models (1, 2, 5, 6, 7), elemental pressure DOFs are appended after the Ferrite displacement DOFs. The global index for element el and local pressure DOF p (1-based) is ndofs(dh) + 4*(el-1) + p.

For pure-displacement models (IsotropicCompressible), uses Ferrite's native allocate_matrix for optimal sparsity structure.

Returned matrix size

Model	Size
Incompressible (1,2,5,6,7)	(ndofs(dh)+4*ne) × same
Compressible (3)	ndofs(dh) × ndofs(dh)

source

Sentinel.assemble_stiffness! Function

assemble_stiffness!(K, dh, cv_disp, cv_press, model, props::GaussPointMaterial, ω) -> K

Assemble the global complex stiffness matrix using uniform material properties (the same GaussPointMaterial at every Gauss point of every element).

Arguments

• K: preallocated SparseMatrixCSC{ComplexF64} from allocate_stiffness(dh, model)

• dh: DofHandler with vector field Lagrange{RefHexahedron, 2}()^3

• cv_disp: CellValues for scalar Lagrange{RefHexahedron, 2}(), used for shape function evaluation inside the element kernels

• cv_press: CellValues for pressure field (or nothing; currently unused)

• model: AbstractMaterialModel singleton for kernel dispatch

• props: GaussPointMaterial applied uniformly at all Gauss points

• ω: angular frequency [rad/s]

Notes

• K is zeroed on entry; existing values are discarded.

• For incompressible models, pressure DOFs are at global indices ndofs(dh) + 4*(el-1) + p for element el and local DOF p (1..4).

• The assembled matrix is complex symmetric (K = Kᵀ), not Hermitian.

Modifies

K in-place; also returns K.

source

assemble_stiffness!(K, dh, cv_disp, cv_press, model, material, gp2mtrs, meshes, ω;
                    ρ=1000.0, fiber_field=nothing) -> K

Assemble the global complex stiffness matrix from a Material struct with spatially varying properties interpolated via the GP2MTR mapping.

For each element and each Gauss point:

1. Maps the Ferrite flat quadrature index to the GP2Mtr (ig, jg, kg) indices

2. Evaluates each material property at the Gauss point location using the material mesh interpolation (evaluate_property_at_gp)

3. Constructs a GaussPointMaterial via model-specific build_gp_material

4. Calls element_stiffness! with the per-GP materials

5. Scatters into the global K

Arguments

• K: preallocated SparseMatrixCSC{ComplexF64} from allocate_stiffness

• dh: DofHandler with vector field Lagrange{RefHexahedron, 2}()^3

• cv_disp: CellValues for scalar Lagrange{RefHexahedron, 2}()

• cv_press: CellValues for pressure field (or nothing)

• model: AbstractMaterialModel singleton for kernel dispatch

• material: Material struct with per-property SingleProperty arrays

• gp2mtrs: vector of GP2Mtr (one per material mesh)

• meshes: vector of MaterialMesh (one per material mesh)

• ω: angular frequency [rad/s]

Keyword Arguments

• ρ: density kg/m³

• fiber_field: optional per-element fiber direction vectors (for TI models)

Notes

• K is zeroed on entry; existing values are discarded.

• Currently requires prop.rmesh == prop.imesh for all properties (i.e., real and imaginary parts on the same material mesh).

source

Boundary Conditions

Dirichlet boundary conditions are applied via symmetric elimination: the constrained row and column are zeroed, and the diagonal is set to 1.

Sentinel.apply_dirichlet! Function

apply_dirichlet!(K, f, constrained_dofs, prescribed_values)

Apply Dirichlet boundary conditions by modifying K and f in-place.

Uses the symmetric elimination method:

1. For each constrained DOF k with prescribed value v:

• Update f[i] -= K[i,k] * v for all rows i (read column k before zeroing)

• Zero column k of K (efficient in CSC format)

• Zero row k of K (scan all columns)

• Set K[k,k] = 1, f[k] = v

This preserves complex symmetry of K (K = Kᵀ) after BC application.

Arguments

• K: preallocated SparseMatrixCSC{ComplexF64} to modify in-place

• f: right-hand side vector, length size(K, 1)

• constrained_dofs: 1-based global DOF indices to constrain

• prescribed_values: corresponding prescribed complex values

source

apply_dirichlet!(K, f, bcs::BoundaryConditions, dispset::Int)

Apply Dirichlet boundary conditions from a BoundaryConditions struct for the given displacement set index.

Iterates bcs.bcind[:, dispset] to find all constrained DOFs with BC_DIRICHLET type, then delegates to the DOF-vector overload.

source

Linear Solvers

Sentinel.AbstractLinearSolver Type

AbstractLinearSolver

Abstract type for sparse linear solver backends. Subtypes must implement linear_solve!(x, K, f, solver).

source

Sentinel.DirectSolver Type

DirectSolver <: AbstractLinearSolver

Julia's built-in sparse direct solver (SuiteSparse/UMFPACK). Works for any sparse complex system without additional dependencies.

source

Sentinel.MUMPSSolver Type

MUMPSSolver <: AbstractLinearSolver

MUMPS sparse direct solver. Requires the MUMPS.jl package extension.

Fields

• sym::Int: MUMPS symmetry flag. 2 for complex symmetric (default, used for viscoelastic models 1,2,3,5,6,7), 0 for general unsymmetric (poroelastic model 4).

source

Sentinel.linear_solve! Function

linear_solve!(x, K, f, solver::AbstractLinearSolver) -> x

Solve the sparse linear system K·x = f in-place, storing the result in x.

Arguments

• x::Vector{ComplexF64}: preallocated solution vector (overwritten)

• K::SparseMatrixCSC{ComplexF64}: assembled stiffness matrix (with BCs applied)

• f::Vector{ComplexF64}: right-hand side vector (with BCs applied)

• solver: solver backend

source

Forward Solve

High-level forward solve that handles DOF scattering:

Sentinel.forward_solve! Function

forward_solve!(disp::Displacement, K, f, model::AbstractMaterialModel,
               dispset::Int; solver=DirectSolver()) -> Displacement

Solve the forward problem K·u = f and store the solution in disp.

K and f must already have boundary conditions applied (via apply_dirichlet!). The solution vector is scattered into disp for the given dispset using scatter_solution!.

Arguments

• disp: preallocated Displacement struct (from init_displacement!)

• K: assembled stiffness matrix with BCs applied

• f: right-hand side vector with BCs applied

• model: material model (used for dispatch/documentation only)

• dispset: displacement set index (1-based)

• solver: linear solver backend (default: DirectSolver())

source

Adjoint Solve

The adjoint problem $\bar{\mathbf{K}}\mathit{\lambda }=-({\mathbf{u}}_{\mathrm{calc}}-{\mathbf{u}}_{\mathrm{meas}})$ is solved for gradient computation. Since $\mathbf{K}$ is complex symmetric, $\bar{\mathbf{K}}={\mathbf{K}}^H$.

Sentinel.adjoint_solve! Function

adjoint_solve!(adjnt::Displacement, K, calc::Displacement, meas::Displacement,
               bcs::BoundaryConditions, dispset::Int;
               solver::AbstractLinearSolver = DirectSolver()) -> Displacement

Solve the adjoint problem conj(K)·λ = -(u_calc - u_meas) and store the result in adjnt for the given displacement set.

K should be the stiffness matrix with forward Dirichlet BCs already applied (rows/columns zeroed, diagonal = 1 at constrained DOFs). The adjoint uses homogeneous BCs on the same DOFs by zeroing the RHS there.

Arguments

• adjnt: preallocated Displacement struct (from init_displacement!)

• K: assembled stiffness matrix with forward BCs applied

• calc: calculated displacement (from forward solve)

• meas: measured displacement (experimental data)

• bcs: BoundaryConditions identifying constrained DOFs

• dispset: displacement set index (1-based)

• solver: linear solver backend (default: DirectSolver())

source