User Guide

This guide covers practical workflows for using Sentinel.jl, from selecting a material model to running a full reconstruction and visualizing results.

Material Model Selection

Sentinel implements 7 material models. Choose based on your tissue type and available data:

Model	Type	When to Use	Parameters
1	IsotropicIncompressible	Most soft tissues (brain, liver). Default choice.	${\mu }^{\ast },{\kappa }^{\ast },\rho$
2	Orthotropic	Tissues with 3 distinct material axes	$E_1^{\ast },E_2^{\ast },E_3^{\ast },\rho$
3	IsotropicCompressible	When compressibility matters (lung)	${\mu }^{\ast },{\kappa }^{\ast },\rho$
4	Poroelastic	Fluid-saturated tissues (cartilage)	${\mu }_s,{\kappa }_s,{\kappa }_f,k,\varphi ,\eta ,{\rho }_s,{\rho }_f,{\rho }_a$
5	TransverseIsotropicV1	Fiber-reinforced tissue (muscle, white matter)	${\mu }^{\ast },{\mu }_s^{\ast },E_t^{\ast },{\kappa }^{\ast },\rho ,\mathbf{n}$
6	TransverseIsotropicV2	Same as V1, alternate parameterization	${\mu }^{\ast },\varphi ,\zeta ,{\kappa }^{\ast },\rho ,\mathbf{n}$
7	GeneralizedAnisotropic	General TI with Itzkov-Askel incompressibility	$E_1^{\ast },E_2^{\ast },{\mu }^{\ast },{\kappa }^{\ast },\rho ,\mathbf{n}$

Key points:

• Models 1-3 and 5-7 use Hex27 elements; Model 4 uses Tet4

• Models 1, 2, 5, 6, 7 have pressure DOFs (mixed formulation)

• Models 5, 6, 7 require fiber directions from DTI data

• All complex moduli use the convention ${\mu }^{\ast }={\mu }^{\prime }+i{\mu }^{″}$ where ${\mu }^{″}<0$ for physical damping

See the Mathematical Reference for complete constitutive equations.

Mesh Generation from MRI

Input Data

Sentinel reads MRI displacement data from three formats:

# Siemens DICOM (phase images)
mre_data = read_mre_siemens("dicom_dir/", 60.0, 1.0)

# WashU .mat format
mre_data = read_mre_washu("washu_data.mat")

# UIUC .mat format
mre_data = read_mre_uiuc("uiuc_data.mat")

All readers return an MREData struct containing complex displacement arrays, anatomical data, mask, and voxel spacing.

Mesh Configuration

The HexMeshConfig struct controls mesh generation:

config = HexMeshConfig(
    mesh_strategy = 2,     # 1=voxel, 2=wavelength, 3=target resolution
    mu_estimate = 3000.0,  # for wavelength calculation [Pa]
    rho_estimate = 1000.0, # density [kg/m^3]
    nodes_per_wavelength = 8.0,
    buffer_size = 2,       # boundary padding elements
    spline_bc_type = :natural
)

Mesh strategies:

1. One node per voxel: direct mapping, preserves MRI resolution

2. Wavelength-based: automatically sizes elements to resolve shear wavelength (recommended)

3. Target resolution: user specifies element size in mm

Generate and Export

result = generate_hex_mesh(mre_data, config)

# Write NLI-format files (.nod, .elm, .dsp, .bnd)
write_hex_mesh_files(result, "output/mesh")

# Write VTK for ParaView visualization
export_vtk("output/mesh_preview", result)

The HexMeshResult contains the Ferrite grid, full hex27 connectivity, interpolated displacement data, boundary nodes, and metadata.

DTI Fiber Processing

For transverse isotropic models (5, 6, 7), fiber directions are needed:

# Process DTI data
dti_data = process_dti("dti_data.mat", size(mre_data.Ur)[1:3])

# Interpolate fiber directions to FEM mesh elements
fiber_dirs = interpolate_fiber_directions(dti_data, hex_mesh_result)

The DTIData struct contains the principal eigenvector field (V1) and fractional anisotropy (FA). Fiber directions are interpolated from the DTI grid to each element centroid using trilinear interpolation.

Regularization Selection

Sentinel provides 10 regularization methods in three categories:

Penalty Methods (add to objective function)

Method	Use Case	Key Parameter
TikhonovReg	General smoothing, good default	weight $\alpha$
TotalVariationReg	Edge-preserving reconstruction	weight $\alpha$, smoothing $\epsilon$
SoftPriorReg	Region-based prior information	weight $\alpha$, region assignments
ExteriorConstraintReg	Keep properties in physical range	weight $\alpha$, bounds

Hessian Modifiers (modify update direction)

Method	Use Case
JoachimowiczReg	Error-scaled Hessian damping
MarquardtReg	Levenberg-Marquardt damping with normalization

Post-Processing (applied after each iteration)

Method	Use Case
VanHoutenReg	Clip power-law exponents to physical range
SpatialFilterReg	Gaussian smoothing of property maps
BoundsReg	Enforce Poisson ratio bounds

Weight Scheduling

All penalty regularizations support scheduled weights via WeightSchedule:

ws = WeightSchedule(start_weight=1e-2, end_weight=1e-4, max_iter=20, delay=0)
tikhonov = TikhonovReg(1, ws, ws, material)

The weight interpolates linearly from start_weight to end_weight over max_iter iterations.

Optimizer Selection

Optimizer	Best For	Memory	Convergence
conjugate_gradient!	General use, large problems	O(n)	Linear
gauss_newton!	Well-conditioned problems	O(n)	Superlinear
quasi_newton!	Fast convergence, moderate problems	O(mn)	Superlinear

Recommendations:

• Start with CG — robust, low memory, good baseline

• Try Gauss-Newton with Marquardt regularization for faster convergence

• Use L-BFGS (quasi_newton!) when CG is too slow and memory permits (controlled by memory parameter, default 10)

All optimizers use Armijo backtracking line search. See line_search for parameters.

Forward Problem Workflow

Manual Setup

using Sentinel, Ferrite

# Load mesh
grid, full_conn = read_mesh_legacy("mesh.nod", "mesh.elm")
dh = setup_hex27_dofhandler(grid)

# Cell values
ip = Lagrange{RefHexahedron, 2}()
qr = QuadratureRule{RefHexahedron}(3)
cv_disp = CellValues(qr, ip)
cv_press = CellValues(qr, Lagrange{RefHexahedron, 1}())

# Material
model = IsotropicIncompressible()
# ... set up GaussPointMaterial or Material + GP2Mtr ...

# Assemble
K = allocate_stiffness(dh, model)
assemble_stiffness!(K, dh, cv_disp, cv_press, model, props, omega)

# Boundary conditions
bcs = read_bcs_file("mesh.bcs")
f = zeros(ComplexF64, size(K, 1))
apply_dirichlet!(K, f, bcs, 1)  # dispset 1

# Solve
disp = Displacement()
init_displacement!(disp, ndofs(dh) ÷ 3, getncells(grid), 4, 1)
forward_solve!(disp, K, f, model, 1; solver=DirectSolver())

Runfile Setup (Recommended)

config = parse_runfile("problem.dat")
setup = setup_forward_problem(config, basedir)
# setup contains: grid, dh, cv_disp, cv_press, model, material, bcs, K, ...

Postprocessing

VTK Export

# From HexMeshResult (preprocessing)
export_vtk("mesh_output", hex_mesh_result;
    write_disp=true, write_anatomical=true, write_region=true)

Convergence Analysis

# Print formatted table
print_convergence_table(history)

# Get statistics
stats = convergence_statistics(history)
# stats.initial_objective, stats.final_objective, stats.objective_reduction, ...

# Export to CSV
export_convergence_csv("convergence.csv", history)

Property Interpolation

Interpolate reconstructed properties from the material mesh to arbitrary points or a structured grid:

# To a structured grid
result = interpolate_to_grid(material, mesh, origin, resolution, dims)
# result.data["prop1_real"] is a 3D array

# To arbitrary points
values = interpolate_properties(material, mesh, target_coords;
    prop_indices=[1], val_idx=1)