Mathematical Reference

This page provides the complete mathematical formulations used in Sentinel.jl, including the governing equations, all seven material model element kernels, the inverse problem, regularization functionals, and optimizer update rules.

General Framework

Governing Equation

Sentinel solves the frequency-domain viscoelastic wave equation. Given harmonic excitation at angular frequency $\omega$, the time-harmonic displacement $\mathbf{u}(\mathbf{x})e^{i\omega t}$ satisfies:

$$-\mathrm{\nabla }\cdot \mathit{\sigma }(\mathbf{u})-\rho {\omega }^2\mathbf{u}=\mathbf{0}$$

where $\mathit{\sigma }$ is the Cauchy stress tensor and $\rho$ is the tissue density.

Mixed Finite Element Formulation

For incompressible and nearly-incompressible models (Models 1, 2, 5, 6, 7), a mixed $(\mathbf{u},p)$ formulation avoids volumetric locking. The displacement field uses triquadratic (Hex27) basis functions $N_i$, while the pressure field uses an element-local basis:

$${\psi }_k(\xi ,\eta ,\zeta )=\{1,\xi ,\eta ,\zeta \},k=1,\dots ,4$$

This gives 81 displacement DOFs (27 nodes $\times$ 3 components) and 4 pressure DOFs per element, for a total of 85 local DOFs.

Element Stiffness Block Structure

The element stiffness matrix has the $2\times 2$ block form:

$${\mathbf{K}}_e=\left[\begin{array}{cc}{\mathbf{K}}_{uu}& {\mathbf{K}}_{up}\\ {\mathbf{K}}_{pu}& {\mathbf{K}}_{pp}\end{array}\right]$$

where ${\mathbf{K}}_{uu}$ is $81\times 81$, ${\mathbf{K}}_{up}$ is $81\times 4$, ${\mathbf{K}}_{pu}={\mathbf{K}}_{up}^T$, and ${\mathbf{K}}_{pp}$ is $4\times 4$.

Complex Symmetry

All stiffness matrices are complex symmetric: $\mathbf{K}={\mathbf{K}}^T$ (transpose, NOT Hermitian conjugate). This arises from the viscoelastic constitutive law with complex-valued moduli ${\mu }^{\ast }={\mu }^{\prime }+i{\mu }^{″}$.

Static Pressure Condensation

The pressure basis ${\psi }_k$ is element-local and discontinuous: a pressure DOF belongs to a single element and is shared with no neighbor. Globally this makes the pressure–pressure block ${\mathbf{K}}_{pp}$ block-diagonal, with one independent $4\times 4$ block per element. For the penalty-regularized incompressible models (1, 2, 5, 6) this block is an invertible $\kappa$-scaled penalty — for Model 1 it is ${\mathbf{K}}_{pp}^{kl}=-\frac{1}{{\kappa }^{\ast }}\int {\psi }_k{\psi }_{l}d\mathrm{\Omega }$ (see the per-model kernels below) — so the pressure DOFs can be eliminated element by element before the global solve. (Model 7 uses a pure Lagrange multiplier with ${\mathbf{K}}_{pp}=\mathbf{0}$, which is singular and cannot be condensed this way.)

Writing the global mixed system as a $2\times 2$ saddle-point system,

$$\left[\begin{array}{cc}{\mathbf{K}}_{uu}& {\mathbf{K}}_{up}\\ {\mathbf{K}}_{pu}& {\mathbf{K}}_{pp}\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \mathbf{p}\end{array}\right]=\left[\begin{array}{c}{\mathbf{f}}_u\\ {\mathbf{f}}_p\end{array}\right],$$

eliminating $\mathbf{p}={\mathbf{K}}_{pp}^{-1}({\mathbf{f}}_p-{\mathbf{K}}_{pu}\mathbf{u})$ yields the Schur complement system in the displacement unknowns alone:

$$\mathbf{S}\mathbf{u}={\mathbf{f}}_u^{\prime },\mathbf{S}={\mathbf{K}}_{uu}-{\mathbf{K}}_{up}{\mathbf{K}}_{pp}^{-1}{\mathbf{K}}_{pu},{\mathbf{f}}_u^{\prime }={\mathbf{f}}_u-{\mathbf{K}}_{up}{\mathbf{K}}_{pp}^{-1}{\mathbf{f}}_p.$$

Because ${\mathbf{K}}_{pp}$ is block-diagonal, the condensation is performed entirely at the element level — each ${\mathbf{S}}_e={\mathbf{K}}_{uu}^e-{\mathbf{K}}_{up}^e({\mathbf{K}}_{pp}^e{)}^{-1}{\mathbf{K}}_{pu}^e$ is formed from the $85\times 85$ element matrix via one $4\times 4$ solve, then scattered like an ordinary $81\times 81$ displacement stiffness.

The reasoning behind condensation — why the elimination is exact and free of fill on this discretization, and why it improves accuracy rather than merely speed — is in About Static Pressure Condensation.

Applicability

The forward solve, the objective, and the ${\mu }^{\ast }$- and $\rho$-gradients depend only on the displacement solution $\mathbf{u}$; the pressure $\mathbf{p}$ need not be recovered for these. The only gradient terms that read $\mathbf{p}$ are the bulk-modulus gradient $\partial \mathrm{\Phi }/\partial \kappa$ (which differentiates ${\mathbf{K}}_{pp}$) and the anisotropic-model pressure-coupling gradients. Static condensation therefore applies to the Model 1 reconstruction of the complex shear modulus ${\mu }^{\ast }={\mu }^{\prime }+i{\mu }^{″}$ (both the storage modulus ${\mu }^{\prime }$ and loss modulus ${\mu }^{″}$), optionally with density $\rho$, while the bulk modulus ${\kappa }^{\ast }$ is held fixed — which is the standard viscoelastic MRE reconstruction, and is the default on the batched solver path. See How to Tune Performance for usage and benchmarks.

Model 1: Isotropic Incompressible

Parameters: complex shear modulus ${\mu }^{\ast }$, complex bulk modulus ${\kappa }^{\ast }$, density $\rho$

Constitutive law: $\mathit{\sigma }=2{\mu }^{\ast }{\mathit{\epsilon }}^{\mathrm{dev}}-p\mathbf{I}$

Displacement-Displacement Block

Diagonal terms ($\alpha =\beta$):

$$K_{ij}^{u_x{u}_x}={\int }_{{\mathrm{\Omega }}_e}{\mu }^{\ast }(\frac{4}{3}{N}_{i,x}{N}_{j,x}+N_{i,y}{N}_{j,y}+N_{i,z}{N}_{j,z})-{\omega }^2\rho N_i{N}_{j}d\mathrm{\Omega }$$

Off-diagonal terms ($\alpha \ne \beta$):

$$K_{ij}^{u_x{u}_y}={\int }_{{\mathrm{\Omega }}_e}{\mu }^{\ast }(N_{i,y}{N}_{j,x}-\frac{2}{3}{N}_{i,x}{N}_{j,y})d\mathrm{\Omega }$$

with cyclic permutations for other component pairs.

Displacement-Pressure Block

$$K_i^{u_{\alpha },p_k}=-{\int }_{{\mathrm{\Omega }}_e}{\psi }_k{N}_{i,\alpha }d\mathrm{\Omega }$$

Pressure-Pressure Block

$$K_{kl}^{pp}=-\frac{1}{{\kappa }^{\ast }}{\int }_{{\mathrm{\Omega }}_e}{\psi }_k{\psi }_{l}d\mathrm{\Omega }$$

Model 2: Orthotropic

Parameters: 3 Young's moduli $E_1^{\ast },E_2^{\ast },E_3^{\ast }$ (complex), density $\rho$

Compliance Matrix

The $6\times 6$ compliance matrix $\mathbf{S}$ in Voigt notation ($\{{\epsilon }_{xx},{\epsilon }_{yy},{\epsilon }_{zz},{\gamma }_{yz},{\gamma }_{xz},{\gamma }_{xy}\}$):

$$\mathbf{S}=\left[\begin{array}{cccccc}1/E_1& -{\nu }_{12}/E_1& -{\nu }_{13}/E_1& 0& 0& 0\\ -{\nu }_{21}/E_2& 1/E_2& -{\nu }_{23}/E_2& 0& 0& 0\\ -{\nu }_{31}/E_3& -{\nu }_{32}/E_3& 1/E_3& 0& 0& 0\\ 0& 0& 0& 1/G_{23}& 0& 0\\ 0& 0& 0& 0& 1/G_{13}& 0\\ 0& 0& 0& 0& 0& 1/G_{12}\end{array}\right]$$

Derived Poisson Ratios

$${\nu }_{12}=\frac{1}{2}(\frac{E_1}{E_2}-\frac{E_1}{E_3}+1),{\nu }_{13}=\frac{1}{2}(\frac{E_1}{E_3}-\frac{E_1}{E_2}+1)$$$${\nu }_{23}=\frac{1}{2}(\frac{E_2}{E_3}-\frac{E_2}{E_1}+1)$$

Derived Stiffness Coefficients

Define the compliance sub-terms:

$$S_{11}=\frac{1}{E_1},S_{22}=\frac{1}{E_2},S_{33}=\frac{1}{E_3}$$$$S_{12}=\frac{{\nu }_{12}}{E_1},S_{13}=\frac{{\nu }_{13}}{E_1},S_{23}=\frac{{\nu }_{23}}{E_2}$$$$D_1=S_{11}{S}_{23}-S_{12}{S}_{13}$$

The $3\times 3$ normal stiffness sub-block is obtained from the compliance inverse. The stiffness assembly then follows the general anisotropic kernel (Section "General Anisotropic Framework" below), with the $6\times 6$ stiffness matrix $\mathbf{C}={\mathbf{S}}^{-1}$.

Pressure Coupling

Compliance-weighted pressure terms:

$$S{P}_x=\frac{p_{\mathrm{ref}}}{D_1},K_i^{u_x,p_k}=-S{P}_x{\int }_{{\mathrm{\Omega }}_e}{\psi }_k{N}_{i,x}d\mathrm{\Omega }$$

with analogous terms for $y$ and $z$ directions, scaled by the corresponding compliance column sums.

Model 3: Isotropic Compressible

Parameters: complex shear modulus ${\mu }^{\ast }$, complex bulk modulus ${\kappa }^{\ast }$, density $\rho$

No pressure DOFs — element stiffness is $81\times 81$.

Lame Parameters

$${\lambda }^{\ast }={\kappa }^{\ast }-\frac{2}{3}{\mu }^{\ast }$$

Stiffness Matrix

Diagonal ($\alpha =\beta$):

$$K_{ij}^{u_x{u}_x}={\int }_{{\mathrm{\Omega }}_e}(2{\mu }^{\ast }+{\lambda }^{\ast })N_{i,x}{N}_{j,x}+{\mu }^{\ast }(N_{i,y}{N}_{j,y}+N_{i,z}{N}_{j,z})-{\omega }^2\rho N_i{N}_{j}d\mathrm{\Omega }$$

Off-diagonal ($\alpha \ne \beta$):

$$K_{ij}^{u_x{u}_y}={\int }_{{\mathrm{\Omega }}_e}{\lambda }^{\ast }{N}_{i,x}{N}_{j,y}+{\mu }^{\ast }{N}_{i,y}{N}_{j,x}d\mathrm{\Omega }$$

Model 4: Poroelastic (Tet4)

Parameters: solid shear ${\mu }_s$, solid bulk ${\kappa }_s$, fluid bulk ${\kappa }_f$, permeability $k$, porosity $\varphi$, fluid viscosity $\eta$, solid density ${\rho }_s$, fluid density ${\rho }_f$, added mass ${\rho }_a$

Uses linear tetrahedral (Tet4) elements with 4 displacement nodes (12 DOFs) and 4 nodal pressure DOFs, giving a $16\times 16$ element matrix.

Biot Coupling Coefficient

$$\beta (\omega )=\frac{\omega {\eta }^2{\rho }_f\kappa }{i{\eta }^2+\omega \kappa ({\rho }_a+\eta {\rho }_f)}$$

where $\kappa$ is the permeability.

Non-Symmetry

Warning

The poroelastic element matrix is NOT symmetric: ${\mathbf{K}}_{pu}\ne {\mathbf{K}}_{up}^T$. The ${\mathbf{K}}_{pu}$ block contains an additional $i\omega$ factor absent from ${\mathbf{K}}_{up}$.

Analytical Integration

Tet4 uses analytical integration (no Gauss quadrature) based on the constant-gradient property of linear tetrahedra:

$${\int }_{{\mathrm{\Omega }}_e}{N}_i{N}_{j}d\mathrm{\Omega }=\frac{V}{20}(1+{\delta }_{ij}),{\int }_{{\mathrm{\Omega }}_e}{N}_{i,\alpha }{N}_{j,\beta }d\mathrm{\Omega }=\frac{b_{i\alpha }{b}_{j\beta }}{36V}$$

where $V$ is the tetrahedral volume and $b_{i\alpha }$ are the shape function gradient coefficients.

General Anisotropic Framework (Models 2, 5, 6, 7)

Models 2, 5, 6, and 7 all use the general $6\times 6$ stiffness tensor contraction. Given the Voigt stiffness matrix $\mathbf{S}$ (in global coordinates), the displacement-displacement block is:

$$K_{ij}^{u_{\alpha }{u}_{\beta }}={\int }_{{\mathrm{\Omega }}_e}{\mathbf{B}}_{\alpha i}^T\mathbf{S}{\mathbf{B}}_{\beta j}d\mathrm{\Omega }-{\delta }_{\alpha \beta }{\omega }^2\rho N_i{N}_j$$

Expanding the $9$ block contributions from the $6\times 6$ stiffness $S_{pq}$:

Block	Formula
$K^{u_x{u}_x}$	$S_{11}{N}_{i,x}{N}_{j,x}+S_{66}{N}_{i,y}{N}_{j,y}+S_{55}{N}_{i,z}{N}_{j,z}$
$K^{u_x{u}_y}$	$S_{12}{N}_{i,x}{N}_{j,y}+S_{66}{N}_{i,y}{N}_{j,x}$
$K^{u_x{u}_z}$	$S_{13}{N}_{i,x}{N}_{j,z}+S_{55}{N}_{i,z}{N}_{j,x}$
$K^{u_y{u}_x}$	$S_{12}{N}_{i,y}{N}_{j,x}+S_{66}{N}_{i,x}{N}_{j,y}$
$K^{u_y{u}_y}$	$S_{22}{N}_{i,y}{N}_{j,y}+S_{66}{N}_{i,x}{N}_{j,x}+S_{44}{N}_{i,z}{N}_{j,z}$
$K^{u_y{u}_z}$	$S_{23}{N}_{i,y}{N}_{j,z}+S_{44}{N}_{i,z}{N}_{j,y}$
$K^{u_z{u}_x}$	$S_{13}{N}_{i,z}{N}_{j,x}+S_{55}{N}_{i,x}{N}_{j,z}$
$K^{u_z{u}_y}$	$S_{23}{N}_{i,z}{N}_{j,y}+S_{44}{N}_{i,y}{N}_{j,z}$
$K^{u_z{u}_z}$	$S_{33}{N}_{i,z}{N}_{j,z}+S_{44}{N}_{i,y}{N}_{j,y}+S_{55}{N}_{i,x}{N}_{j,x}$

Each integral includes the inertial mass term $-{\omega }^2\rho N_i{N}_j$ on the diagonal blocks only.

Model 5: Transverse Isotropic V1

Parameters: isotropic shear ${\mu }^{\ast }$, axial shear ${\mu }_s^{\ast }$, tensile anisotropy $E_t^{\ast }$, bulk ${\kappa }^{\ast }$, density $\rho$, fiber direction $\mathbf{n}$

Local Stiffness (Fiber-Aligned Frame)

$$\begin{array}{rlrlrl}{S}_{11}^{\lg }& =\frac{4}{9}{E}_t+\frac{4}{3}\mu ,& S_{12}^{\lg }& =\frac{2}{9}{E}_t-\frac{2}{3}\mu ,& S_{13}^{\lg }& =-\frac{2}{9}{E}_t\\ S_{22}^{\lg }& =\mu +\frac{1}{9}{E}_t,& S_{23}^{\lg }& =-\mu +\frac{1}{9}{E}_t,& S_{33}^{\lg }& =\mu +\frac{1}{9}{E}_t\\ S_{44}^{\lg }& ={\mu }_s,& S_{55}^{\lg }& ={\mu }_s,& S_{66}^{\lg }& =\mu \end{array}$$

The global stiffness is obtained via Bond rotation: ${\mathbf{S}}_{\mathrm{global}}={\mathbf{B}}^T{\mathbf{S}}_{\mathrm{local}}\mathbf{B}$ (see Bond Rotation).

Model 6: Transverse Isotropic V2

Parameters: shear ${\mu }^{\ast }$, anisotropy ratio $\varphi$, tensile ratio $\zeta$, bulk ${\kappa }^{\ast }$, density $\rho$, fiber direction $\mathbf{n}$

Local Stiffness (Fiber-Aligned Frame)

$$\begin{array}{rlrlrl}{S}_{11}^{\lg }& =\frac{4}{3}\mu (1+\frac{4}{3}\zeta ),& S_{12}^{\lg }& =-\frac{2}{3}\mu (1-\frac{2}{3}\zeta ),& S_{13}^{\lg }& =-\frac{2}{3}\mu (1+\frac{1}{3}\zeta )\\ S_{22}^{\lg }& =\mu (1+\frac{1}{3}\zeta )+\frac{1}{3}\mu ,& S_{23}^{\lg }& =-\mu (1-\frac{1}{3}\zeta )+\frac{1}{3}\mu \zeta ,& S_{33}^{\lg }& =S_{22}^{\lg }\\ S_{44}^{\lg }& =\varphi \mu ,& S_{55}^{\lg }& =\varphi \mu ,& S_{66}^{\lg }& =\mu \end{array}$$

Model 7: Generalized Anisotropic

Parameters: 2 Young's moduli $E_1^{\ast },E_2^{\ast }$, shear ${\mu }^{\ast }$, bulk ${\kappa }^{\ast }$, density $\rho$, fiber direction $\mathbf{n}$

Based on the Itzkov-Askel incompressible transversely isotropic formulation.

Local Stiffness (Fiber-Aligned Frame)

Define $D=\frac{3}{4E_1}(\frac{4}{E_2}-\frac{1}{E_1})$:

$$\begin{array}{rlrlrl}{S}_{11}^{\lg }& =\frac{1/(E_1{E}_2)}{D},& S_{12}^{\lg }& =\frac{-1/(2E_1^2)}{D},& S_{13}^{\lg }& =S_{12}^{\lg }\\ S_{22}^{\lg }& =\frac{3/(4E_1^2)}{D},& S_{23}^{\lg }& =\frac{1/(4E_1^2)-1/(E_1{E}_2)}{D},& S_{33}^{\lg }& =S_{22}^{\lg }\\ S_{44}^{\lg }& =\mu ,& S_{55}^{\lg }& =\mu ,& S_{66}^{\lg }& =\frac{1}{2E_{1}D}\end{array}$$

Note

Model 7 uses a pure Lagrange multiplier for incompressibility: ${\mathbf{K}}_{pp}=\mathbf{0}$. The bulk modulus $\kappa$ is not used in the pressure-pressure block.

Bond Rotation

For fiber-aligned models (5, 6, 7), the local stiffness ${\mathbf{S}}_{\lg }$ must be rotated to the global frame using the Bond transformation.

Rotation Matrix

Given fiber direction $\mathbf{n}=(n_1,n_2,n_3)$, construct a rotation $\mathbf{T}\in {\mathbb{R}}^{3\times 3}$ whose first column is the fiber direction. The remaining columns are chosen to form a right-handed orthonormal frame:

$$\mathbf{T}=\left[\begin{array}{ccc}\mathbf{n}& {\mathbf{t}}_1& {\mathbf{t}}_2\end{array}\right]$$

where ${\mathbf{t}}_1$ and ${\mathbf{t}}_2$ are computed from cross products with the global axes, with special handling when $\mathbf{n}$ is nearly aligned with $\hat{\mathbf{z}}$.

Bond Matrix

The $6\times 6$ Bond matrix $\mathbf{B}$ transforms Voigt-notation stress/strain vectors under rotation $\mathbf{T}$:

$$\mathbf{B}=\left[\begin{array}{cccccc}{T}_{11}^2& T_{12}^2& T_{13}^2& T_{12}{T}_{13}& T_{11}{T}_{13}& T_{11}{T}_{12}\\ T_{21}^2& T_{22}^2& T_{23}^2& T_{22}{T}_{23}& T_{21}{T}_{23}& T_{21}{T}_{22}\\ T_{31}^2& T_{32}^2& T_{33}^2& T_{32}{T}_{33}& T_{31}{T}_{33}& T_{31}{T}_{32}\\ 2T_{21}{T}_{31}& 2T_{22}{T}_{32}& 2T_{23}{T}_{33}& \cdots & \cdots & \cdots \\ 2T_{11}{T}_{31}& 2T_{12}{T}_{32}& 2T_{13}{T}_{33}& \cdots & \cdots & \cdots \\ 2T_{11}{T}_{21}& 2T_{12}{T}_{22}& 2T_{13}{T}_{23}& \cdots & \cdots & \cdots \end{array}\right]$$

The global stiffness is then:

$${\mathbf{S}}_{\mathrm{global}}={\mathbf{B}}^T{\mathbf{S}}_{\mathrm{local}}\mathbf{B}$$

Inverse Problem

Objective Function

The inverse problem minimizes the misfit between calculated and measured displacements:

$$\mathrm{\Phi }(\mathit{\theta })=\frac{1}{2}\parallel {\mathbf{u}}_{\mathrm{calc}}(\mathit{\theta })-{\mathbf{u}}_{\mathrm{meas}}{\parallel }^2+R(\mathit{\theta })$$

where $\mathit{\theta }$ is the vector of material parameters and $R(\mathit{\theta })$ collects regularization terms.

Adjoint Method

The gradient is computed efficiently via the adjoint method. The adjoint displacement $\mathit{\lambda }$ solves:

$$\bar{\mathbf{K}}\mathit{\lambda }=-({\mathbf{u}}_{\mathrm{calc}}-{\mathbf{u}}_{\mathrm{meas}})$$

where $\bar{\mathbf{K}}$ is the complex conjugate of the stiffness matrix (since $\mathbf{K}$ is complex symmetric, $\bar{\mathbf{K}}={\mathbf{K}}^H$).

Gradient Computation

The gradient of the data-misfit term with respect to material parameter ${\theta }_m$ is:

$$\frac{\partial \mathrm{\Phi }}{\partial {\theta }_m}=-\mathrm{Re}\left[\sum _e\sum _{\mathrm{gp}}{w}_{\mathrm{gp}}{\bar{\mathit{\lambda }}}_e^T\frac{\partial {\mathbf{K}}_e}{\partial {\theta }_m}{\mathbf{u}}_e\right]$$

where the sum runs over elements $e$ and Gauss points, $w_{\mathrm{gp}}$ is the quadrature weight times the Jacobian determinant, and $\partial {\mathbf{K}}_e/\partial {\theta }_m$ is the element stiffness derivative with respect to parameter ${\theta }_m$.

The parameter ${\theta }_m$ typically maps to a material property at a specific material mesh element via the GP2MTR interpolation. See the I/O page for details on the Gauss-point to material-point mapping.

Regularization Functionals

Tikhonov

Standard L2 penalty on deviation from reference:

$$R(\mathit{\theta })=\frac{1}{2}\alpha \parallel \mathit{\theta }-{\mathit{\theta }}_{\mathrm{ref}}{\parallel }^2$$

Gradient: $\mathrm{\nabla }R=\alpha (\mathit{\theta }-{\mathit{\theta }}_{\mathrm{ref}})$

Hessian diagonal: $H_{kk}=\alpha$

Total Variation

Smoothed TV functional on the material mesh:

$$R(\mathit{\theta })=\alpha \sum _e\sum _{\mathrm{gp}}{w}_{\mathrm{gp}}|J_e|\sqrt{|\mathrm{\nabla }\theta {|}^2+{\epsilon }^2}$$

where $\epsilon >0$ is a smoothing parameter and the gradient $\mathrm{\nabla }\theta$ is evaluated on the structured hex8 material mesh using trilinear basis functions.

Gradient: $\frac{\partial R}{\partial {\theta }_k}=\alpha \sum _e\sum _{\mathrm{gp}}{w}_{\mathrm{gp}}|J_e|\frac{\mathrm{\nabla }\theta \cdot \mathrm{\nabla }{N}_k}{\sqrt{|\mathrm{\nabla }\theta {|}^2+{\epsilon }^2}}$

Soft Prior

Region-based discrete Laplacian penalty:

$$R(\mathit{\theta })=\alpha \sum _j(L\theta {)}_j^2$$

where $L$ is a discrete Laplacian operator computed per material region.

Gradient: $\mathrm{\nabla }R=2\alpha L^{T}L\mathit{\theta }$

Exterior Constraint

Barrier penalty for out-of-domain parameters:

$$R(\mathit{\theta })=\alpha \sum _k[max(0,{\theta }_k-{\theta }_{max}{)}^2+max(0,{\theta }_{min}-{\theta }_k{)}^2]$$

Gradient: Only nonzero for parameters violating bounds:

$$\frac{\partial R}{\partial {\theta }_k}=\{\begin{array}{ll}2\alpha ({\theta }_k-{\theta }_{max}),& {\theta }_k>{\theta }_{max}\\ -2\alpha ({\theta }_{min}-{\theta }_k),& {\theta }_k<{\theta }_{min}\\ 0,& \text{otherwise}\end{array}$$

Hessian Modifiers

Joachimowicz: Scales diagonal Hessian by average value, error, and weight:

$$H_{kk}\leftarrow \overline{H}\cdot ϵ\cdot w$$

where $\overline{H}$ is the average diagonal Hessian, $ϵ$ is the current displacement error, and $w$ is a user-specified weight.

Marquardt: Normalized Levenberg-Marquardt damping. First normalize the Hessian:

$${\hat{H}}_{kk}=\frac{H_{kk}}{D_k^2},D_k=\sqrt{H_{kk}}$$

Then add damping ${\hat{H}}_{kk}\leftarrow {\hat{H}}_{kk}+\lambda$, solve for $\widehat{\mathrm{\Delta }\theta }$, and rescale: $\mathrm{\Delta }{\theta }_k={\widehat{\mathrm{\Delta }\theta }}_k/D_k$.

Post-Processing Regularizations

Van Houten: Clips power-law exponents to $[{\alpha }_{min},{\alpha }_{max}]$ after each iteration.

Spatial Filter: Applies Gaussian smoothing with kernel widths $({\sigma }_x,{\sigma }_y,{\sigma }_z)$ on the structured material mesh.

Bounds (McGarry): Enforces Poisson ratio bounds $\nu \in [{\nu }_{min},{\nu }_{max}]$ by clipping derived ratios.

Optimizer Update Rules

Conjugate Gradient (Polak-Ribiere+)

Search direction:

$${\mathbf{d}}_k=-{\mathbf{g}}_k+{\beta }_k{\mathbf{d}}_{k-1}$$

with the Polak-Ribiere+ update:

$${\beta }_k=max(0,\frac{{\mathbf{g}}_k^T({\mathbf{g}}_k-{\mathbf{g}}_{k-1})}{{\mathbf{g}}_{k-1}^T{\mathbf{g}}_{k-1}})$$

The $max(0,\cdot )$ ensures automatic reset to steepest descent when conjugacy is lost. Step size ${\alpha }_k$ is determined by Armijo backtracking line search (see below).

Gauss-Newton

Uses a diagonal Hessian approximation:

$$\mathrm{\Delta }{\theta }_k=-\frac{g_k}{H_{kk}}$$

where $H_{kk}$ is the diagonal of the Gauss-Newton Hessian (assembled from element-level contributions). A floor value prevents division by near-zero entries.

L-BFGS

Limited-memory BFGS quasi-Newton method using the two-loop recursion. Maintains $m$ curvature pairs $\{({\mathbf{s}}_i,{\mathbf{y}}_i)\}$ where:

$${\mathbf{s}}_i={\mathit{\theta }}_{i+1}-{\mathit{\theta }}_i,{\mathbf{y}}_i={\mathbf{g}}_{i+1}-{\mathbf{g}}_i$$

Initial Hessian scaling:

$$\gamma =\frac{{\mathbf{s}}_{k-1}^T{\mathbf{y}}_{k-1}}{{\mathbf{y}}_{k-1}^T{\mathbf{y}}_{k-1}}$$

The search direction is computed via the standard two-loop recursion with $\gamma \mathbf{I}$ as the initial inverse Hessian approximation.

Line Search

All optimizers use Armijo backtracking line search. Given sufficient decrease parameter $c_1$ and backtracking factor $\rho$:

$${\alpha }_{k+1}=\rho {\alpha }_k\text{while}\mathrm{\Phi }(\mathit{\theta }+{\alpha }_k{\mathbf{d}}_k)>\mathrm{\Phi }(\mathit{\theta })+c_1{\alpha }_k{\mathbf{g}}^T{\mathbf{d}}_k$$