Zone Decomposition

The zone decomposition system splits the global inverse problem into overlapping subdomains (zones) that can be solved independently. This is the primary approach for full-scale inverse reconstruction in Sentinel.

Overview

The algorithm proceeds as follows:

Generate a regular grid of zone centers (random seed placement via Park-Miller RNG)
Assign elements to zones based on centroid-in-box test with configurable overlap
Grow zones until each has sufficient interior DOF nodes (numsig target)
Extract zone meshes, BCs, and material from the global problem
Solve each zone with CG/GN/L-BFGS optimization
Consolidate zone results into global material (GP-weighted) and displacement (count-averaged)
Apply post-processing regularization (spatial filter, bounds)
Repeat for multiple global iterations

Configuration

Sentinel.ZoneConfig Type

julia

ZoneConfig

Configuration for zone-based domain decomposition.

Fields

zoneedge::SVector{3, Float64}: zone box half-extents in each direction [m]
overlap::Float64: overlap fraction (0.0–1.0); overlap region = overlap × zoneedge
optimizer::Symbol: zone-level optimizer (:cg, :gn, :lbfgs)
max_zone_iter::Int: maximum iterations per zone solve
zone_tol::Float64: convergence tolerance for zone solve
max_global_iter::Int: maximum global (outer) iterations
global_tol::Float64: convergence tolerance for global iteration
seed_element::Int: element index to seed zone grid generation (0 = auto-center)
randomize_zones::Bool: if true, regenerate zone grid each global iteration with random seed

source

Zone Grid Generation

Sentinel.generate_zone_grid Function

julia

generate_zone_grid(grid, config::ZoneConfig[, centroids];
                   seed_element_override=-1) -> ZoneGrid

Generate a regular grid of zone centers covering the FEM mesh.

Algorithm (ports Fortran zonegen):

Compute element centroids (or use precomputed)
Determine seed center (from seed_element, override, or mesh center)
Generate regular grid radiating outward with spacing = zoneedge
Prune zones with no elements
Precompute element lists for each zone

Arguments

grid: Ferrite Grid
config: ZoneConfig with zoneedge, overlap, seed_element
centroids: (optional) precomputed element centroids to avoid recomputation
seed_element_override: if >= 1, use this element as seed instead of config.seed_element

source

Sentinel.ZoneGrid Type

julia

ZoneGrid

Collection of zone centers with precomputed element centroids for fast element-in-zone queries.

Fields

centers::Vector{SVector{3, Float64}}: zone center coordinates
config::ZoneConfig: zone configuration
element_centroids::Vector{SVector{3, Float64}}: precomputed element centroids
zone_elements::Vector{Vector{Int}}: precomputed element lists per zone

source

Sentinel.ParkMillerRNG Type

julia

ParkMillerRNG

Port of the Fortran random(idum) function from zonefunctions.f90. Park-Miller "minimal" RNG with Bays-Durham shuffle (Numerical Recipes RAN1).

Initialized with idum=-1 and one throwaway call, matching the Fortran idum=-1; rr=random(idum) initialization sequence.

source

Sentinel.park_miller_random! Function

julia

park_miller_random!(rng::ParkMillerRNG) -> Float64

Generate the next uniform random deviate in (0, 1), matching the Fortran random(idum) function exactly.

source

Sentinel.compute_element_centroids Function

julia

compute_element_centroids(grid) -> Vector{SVector{3, Float64}}

Compute the centroid of each element by averaging its corner node coordinates. Works for any Ferrite Grid with Hexahedron cells (8 corner nodes).

source

Zone Extraction

Sentinel.build_zone Function

julia

build_zone(zone_id::Int, center::SVector{3, Float64},
           zone_elements::Vector{Int},
           global_grid, global_dh, global_dof_coords,
           global_material::Material, global_meshes::Vector{MaterialMesh},
           global_gp2mtrs::Vector{GP2Mtr},
           global_meas::Displacement, global_disp::Displacement,
           model::AbstractMaterialModel, config::ZoneConfig;
           fiber_field=nothing) -> Zone

Master function that builds a complete Zone from global problem data.

Steps:

Extract zone grid (Ferrite Grid subset)
Create zone DofHandler
Compute DOF coordinate maps
For each material mesh: extract subset following GP2Mtr links
Build zone-local GP2Mtr
Extract zone material properties
Build zone boundary conditions
Extract measured displacement
Allocate stiffness matrix

source

julia

build_zone(zone_id, center, zone_elements, global_grid, global_dh,
           global_dof_coords, global_material, global_meshes, global_gp2mtrs,
           global_meas, global_disp, model, config,
           global_dof_elcon, global_bnod; fiber_field=nothing) -> Zone

Build a zone using topological boundary classification (Fortran-matching).

Uses element connectivity to determine boundary nodes: a DOF node is boundary if any element connected to it (in the global mesh) is NOT present in the zone, or if it's a global mesh boundary node.

Additional Arguments

global_dof_elcon: precomputed DOF-node-to-element connectivity from build_dof_node_element_connectivity
global_bnod: global boundary node flags from build_global_bnod

source

julia

build_zone(zone_id, center, zone_elements, state::GlobalProblemState,
           material, disp, config; fiber_field=nothing)

Convenience overload that unpacks a GlobalProblemState and delegates to the full 14-argument build_zone. Reduces argument count from 14 to 7.

source

Sentinel.Zone Type

julia

Zone

Complete zone-local problem data, ready for inverse solve.

Fields

zone_id::Int: zone identifier
center::SVector{3, Float64}: zone center coordinates
grid: zone-local Ferrite Grid
dh: zone-local DofHandler
cv_disp: CellValues for displacement
model::AbstractMaterialModel: material model
mapping::ZoneIndexMapping: index maps to global
material::Material: zone-local material
meshes::Vector{MaterialMesh}: zone-local material meshes
gp2mtrs::Vector{GP2Mtr}: zone-local GP-to-material mapping
bcs::BoundaryConditions: zone boundary conditions
meas::Displacement: zone-local measured displacement
K::SparseMatrixCSC{ComplexF64}: preallocated stiffness matrix
fiber_field: zone-local fiber directions (or nothing)

source

Sentinel.ZoneIndexMapping Type

julia

ZoneIndexMapping

Bidirectional index maps between global and zone-local meshes.

Fields

zne2glb_elem::Vector{Int}: zone element index → global element index
glb2zne_elem::Dict{Int, Int}: global element index → zone element index
zne2glb_grid_node::Vector{Int}: zone grid node → global grid node
glb2zne_grid_node::Dict{Int, Int}: global grid node → zone grid node
zne2glb_dof_node::Vector{Int}: zone DOF node → global DOF node
glb2zne_dof_node::Dict{Int, Int}: global DOF node → zone DOF node
zne2glb_mtrl_node::Vector{Vector{Int}}: per-mesh zone material node → global
glb2zne_mtrl_node::Vector{Dict{Int, Int}}: per-mesh global material node → zone
zne2glb_mtrl_elem::Vector{Vector{Int}}: per-mesh zone material elem → global
glb2zne_mtrl_elem::Vector{Dict{Int, Int}}: per-mesh global material elem → zone
boundary_dof_nodes::Vector{Int}: zone DOF nodes on boundary (prescribed)
interior_dof_nodes::Vector{Int}: zone DOF nodes in interior (free)

source

Sentinel.extract_zone_dof_mapping Function

julia

extract_zone_dof_mapping(global_dh, zone_dh, global_dof_coords, zone_dof_coords;
                          tol=1e-10, global_coord_lookup=nothing)
    -> (zne2glb_dof_node, glb2zne_dof_node)

Build bidirectional DOF node maps by matching physical coordinates.

Hex27 DOF nodes include edge midpoints, face centers, and body centers that are NOT stored as Grid nodes. We match by physical position with tolerance.

Pass a pre-built global_coord_lookup from build_dof_coord_lookup to enable O(n_zone) lookup instead of O(n_zone × n_global) linear search.

source

Sentinel.classify_zone_dof_nodes Function

julia

classify_zone_dof_nodes(zne2glb_dof_node::Vector{Int},
                         zone_elements_set::Set{Int},
                         global_dof_elcon::Vector{Vector{Int}},
                         global_bnod::Union{BitVector, Nothing}=nothing)
                         -> (boundary, interior)

Classify each zone DOF node as boundary (prescribed) or interior (free) using element connectivity (topological approach), matching Fortran.

A node is boundary if:

Any element connected to it (in the global mesh) is NOT in the zone, OR
It is a global boundary node (global_bnod[gi] == true)

Arguments

zne2glb_dof_node: maps zone DOF node indices to global DOF node indices
zone_elements_set: Set of global element indices in the zone (for O(1) lookup)
global_dof_elcon: global DOF-node-to-element connectivity from build_dof_node_element_connectivity
global_bnod: global boundary node flags (from BC file). If nothing, uses bounding box detection on global coordinates.

Returns vectors of zone-local DOF node indices.

source

Sentinel.build_zone_bcs Function

julia

build_zone_bcs(boundary_dof_nodes::Vector{Int},
                global_disp::Displacement,
                zne2glb_dof_node::Vector{Int},
                zone_dh,
                zone_ne::Int,
                model::AbstractMaterialModel) -> BoundaryConditions

Build zone-level BoundaryConditions by prescribing the current global displacement estimate at all boundary DOF nodes.

For each boundary DOF node zi:

Map to global DOF node gi = zne2glb_dof_node[zi]
Get displacement (u, v, w) from global_disp.disp[gi, ds]
Set Dirichlet BCs on DOFs 3(zi-1)+1, 3(zi-1)+2, 3(zi-1)+3

source

Zone Solving

Sentinel.solve_zone! Function

julia

solve_zone!(zone::Zone, config::ZoneConfig, ω::Float64, ρ::Float64;
            solver::AbstractLinearSolver=DirectSolver(),
            regularizations=nothing,
            iteration_offset::Int=0) -> ZoneResult

Solve the inverse problem on a single zone.

Build ForwardProblemContext from zone data
Dispatch to configured optimizer (:cg, :gn, or :lbfgs)
Return ZoneResult with optimized material and convergence history

source

Sentinel.ZoneResult Type

julia

ZoneResult

Result from solving a single zone's inverse problem.

Fields

zone_id::Int: zone identifier
material::Material: zone-local optimized material
calc::Displacement: zone-local calculated displacement
history::ConvergenceHistory: convergence history
mapping::ZoneIndexMapping: index maps for scattering back to global

source

Consolidation

Sentinel.prepare_accumulation! Function

julia

prepare_accumulation!(material::Material)

Zero all property values and counts in preparation for zone accumulation.

source

Sentinel.prepare_displacement_accumulation! Function

julia

prepare_displacement_accumulation!(disp::Displacement)

Zero displacement values and counts for accumulation.

source

Sentinel.scatter_zone_material! Function

julia

scatter_zone_material!(global_material::Material, result::ZoneResult)

Accumulate zone-local property values into the global material arrays, using Gauss-point weighted accumulation (matching Fortran MRE-Zone consolidation).

Each zone's contribution is weighted by the number of Gauss points that map to each material node in that zone. The accumulated GP counts (rgp/igp) are used by consolidate_material! to compute the weighted average.

source

Sentinel.scatter_zone_displacement! Function

julia

scatter_zone_displacement!(global_disp::Displacement,
                            zone_calc::Displacement,
                            zne2glb_dof_node::Vector{Int})

Accumulate zone-local calculated displacement into the global displacement.

source

Sentinel.consolidate_material! Function

julia

consolidate_material!(material::Material, prev_material::Material)

Divide accumulated property values by accumulated Gauss-point counts to complete GP-weighted averaging (matching Fortran MRE-Zone consolidation).

For nodes with no GP contributions (rgp=0), restore the previous iteration's value to avoid zeroing out uncovered regions.

source

Sentinel.consolidate_displacement! Function

julia

consolidate_displacement!(disp::Displacement, prev_disp::Displacement)

Divide accumulated displacement values by their counts to complete averaging. For nodes with count == 0 (uncovered by any zone), retain the previous iteration's displacement (matching Fortran where globcalc is NOT zeroed between iterations, so uncovered nodes keep their previous value).

source

Sentinel.consolidate_global! Function

julia

consolidate_global!(material::Material, disp::Displacement,
                     prev_material::Material, prev_disp::Displacement)

Complete the consolidation by averaging both material and displacement. Uncovered material nodes are restored from prev_material. Uncovered displacement nodes are restored from prev_disp (previous iteration's displacement, matching Fortran where globcalc starts zeroed and is never reset).

source

High-Level Driver

Sentinel.zone_decomposition_solve! Function

julia

zone_decomposition_solve!(material::Material, grid, dh,
                           gp2mtrs::Vector{GP2Mtr},
                           meshes::Vector{MaterialMesh},
                           meas::Displacement,
                           bcs::BoundaryConditions,
                           model::AbstractMaterialModel,
                           ω::Float64, ρ::Float64,
                           config::ZoneConfig;
                           solver::AbstractLinearSolver=DirectSolver(),
                           regularizations=nothing,
                           fiber_field=nothing,
                           verbose::Bool=true)
    -> (Material, Displacement, GlobalIterationState)

Run the zone decomposition inverse solver (serial version).

Algorithm

Generate zone grid covering the mesh
Initialize displacement estimate from measured data
For each global iteration: a. Build zone-local problems from current global state b. Solve each zone independently c. Scatter/consolidate results back to global arrays d. Check convergence (material_epsilon < global_tol)

Returns

Updated material, calculated displacement, and iteration state

source

Sentinel.zone_decomposition_solve_parallel! Function

julia

zone_decomposition_solve_parallel!(material::Material, grid, dh,
                                    gp2mtrs::Vector{GP2Mtr},
                                    meshes::Vector{MaterialMesh},
                                    meas::Displacement,
                                    bcs::BoundaryConditions,
                                    model::AbstractMaterialModel,
                                    ω::Float64, ρ::Float64,
                                    config::ZoneConfig;
                                    solver::AbstractLinearSolver=DirectSolver(),
                                    regularizations=nothing,
                                    fiber_field=nothing,
                                    verbose::Bool=true,
                                    pmap_func::Function=map)
    -> (Material, Displacement, GlobalIterationState)

Run zone decomposition with a user-provided parallel map function.

To use Distributed.jl, pass pmap_func=Distributed.pmap after loading the Distributed stdlib. With the default map, this runs serially.

Arguments

Same as zone_decomposition_solve! plus:
pmap_func: parallel map function (default: map for serial)

source

Sentinel.GlobalIterationState Type

julia

GlobalIterationState

Tracks the state of the global zone decomposition iteration loop.

Fields

iteration::Int: current global iteration number
material::Material: current global material estimate
global_obj::Float64: global objective value
zone_objectives::Vector{Float64}: per-zone objective values
mat_epsilon::Float64: relative material change from previous iteration
converged::Bool: whether global convergence criterion is met

source

Inverse Problem Setup

These functions handle the full inverse problem setup from a runfile:

Sentinel.InverseRunfileConfig Type

julia

InverseRunfileConfig

Complete parsed inverse problem runfile.

source

Sentinel.parse_inverse_runfile Function

julia

parse_inverse_runfile(filename) -> InverseRunfileConfig

Parse an inverse problem .dat runfile and return an InverseRunfileConfig.

The inverse runfile alternates between descriptive label lines and data lines. This parser reads all lines and processes them sequentially, skipping label lines explicitly to match the Fortran readinverse read order.

source

Sentinel.setup_inverse_problem Function

julia

setup_inverse_problem(config::InverseRunfileConfig,
                       basedir::AbstractString=".") -> InverseProblemSetup

Load all files referenced by an InverseRunfileConfig and return an InverseProblemSetup with the complete inverse problem setup.

Returns

InverseProblemSetup with fields:

grid: Ferrite Grid (hex27)
dh: DofHandler
cv_disp: CellValues for scalar hex27
model: AbstractMaterialModel
material: Material with initialized properties
meshes: Vector
gp2mtrs: Vector
bcs: BoundaryConditions
meas: Displacement (measured)
omega: angular frequency
K: preallocated stiffness matrix
zone_config: ZoneConfig
regularizations: Vector
config: the original InverseRunfileConfig
solver: DirectSolver
full_connectivity: ne×27 connectivity
dof_to_mesh: DOF-to-mesh-node mapping

source

Sentinel.run_inverse_solve! Function

julia

run_inverse_solve!(setup::InverseProblemSetup; verbose::Bool=true,
                    callback::Union{Nothing, Function}=nothing,
                    pmap_func::Function=map)
    -> (Material, Displacement, GlobalIterationState)

Run the zone decomposition inverse solver with escalating iteration structures.

This is a feature-complete version of zone_decomposition_solve! that handles:

Escalating zone iteration structures (CG count changes with global iteration)
Reconstruction indicators (non-reconstructed properties are preserved)
Post-processing regularizations (VanHouten, Bounds, SpatialFilter)

The optional callback(giter, material, state) is called after each global iteration, enabling per-iteration recording of material state, timing, or file output.

Pass pmap_func (e.g. Distributed.pmap) to parallelize zone solves across workers.

source

Sentinel.build_regularizations Function

julia

build_regularizations(config::InverseRunfileConfig,
                       material::Material,
                       meshes::Vector{MaterialMesh};
                       gp2mtrs::Union{Vector{GP2Mtr}, Nothing}=nothing,
                       region_info::Union{RegionInfo, Nothing}=nothing) -> Vector{AbstractRegularization}

Construct regularization objects based on the parsed regularization config.

When gp2mtrs is provided, computes a sensitivity mask (matching Fortran nodsense) for each spatial filter property: only material nodes covered by at least one Gauss point participate in Gaussian filtering.

source

Sentinel.InverseProblemSetup Type

julia

InverseProblemSetup

Complete inverse problem setup with compile-time field name checking.

Replaces the NamedTuple returned by setup_inverse_problem. All field names match the original NamedTuple, so existing setup.field access works without changes.

Fields

grid: Ferrite Grid (hex27)
dh: DofHandler
cv_disp: CellValues for displacement interpolation
model: material model singleton
omega: angular frequency [rad/s]
material: Material with initialized properties
meshes: vector of MaterialMesh
gp2mtrs: vector of GP2Mtr
bcs: BoundaryConditions
meas: measured Displacement data
K: preallocated stiffness matrix
zone_config: ZoneConfig for zone decomposition
regularizations: vector of AbstractRegularization
config: original InverseRunfileConfig
solver: linear solver backend
full_connectivity: full 27-node connectivity from mesh file
dof_to_mesh: mapping from Ferrite DOF nodes to mesh file nodes

source

Zone Decomposition ​

Overview ​

Configuration ​

Zone Grid Generation ​

Zone Extraction ​

Zone Solving ​

Consolidation ​

High-Level Driver ​

Inverse Problem Setup ​

Zone Decomposition

Overview

Configuration

Zone Grid Generation

Zone Extraction

Zone Solving

Consolidation

High-Level Driver

Inverse Problem Setup