Finite strain finite element prototype

This repository implements the v1 contract in docs/IMPLEMENTATION_SPEC.md:

• full-integration Hex8 and Hex20 serendipity elements;
• centroidal F-bar Hex8 with its complete projection tangent;
• compressible neo-Hookean material through the endpoint F_n/F_np1 material-point API;
• a stateful multiplicative finite-strain J2/Voce plug-in with an analytic algorithmic tangent;
• sparse COO/CSR assembly and a whole-system sparse-LU KKT solve for affine constraints;
• exact and modified Newton, optional residual-based backtracking, recoverable trial rollback, adaptive cutback, and mandatory events;
• Gmsh Physical Groups and translational periodic node maps;
• append-only accepted-state HDF5 output, HDF5 restart, and separate temporal XDMF postprocessing.

The production element tangents use the specification's current-configuration Truesdell push-forward, geometric stiffness, and F-bar projection correction. An independent total-reference F-bar kernel and centered directional finite differences cross-check that decomposition.

Project status and TODO

FE solver

• Implement full-integration Hex8 and Hex20 and Hex-8 F-bar elements.
• Add neo-Hookean and finite-strain J2 plasticity with consistent algorithmic tangent.
• Implement cutback and adaptive time stepping.
• Add support for periodic boundary conditions (via Gmsh-based periodic boundary conditions).
• Implement HDF5 results/restart and temporal XDMF postprocessing for ParaView.
• Implement process-based parallel element evaluation.
• Validate J2 plasticity model against benchmarks from literature.
• Add timing information around time consuming parts of the solver for debugging and observability.
• Add selectable SciPy/SuperLU_MT sparse-direct solvers for the global KKT system.
• Integrate the external multiphase TRIP steel model through its Python bridge and solve periodic multi-grain RVEs.

Parallel global linear solver

• Capture FE KKT systems for isolated linear-solver testing.
• Evaluate SuperLU_MT against the existing SciPy/SuperLU reference solver.
• Develop a persistent native C and Python bridges for SuperLU_MT factorization and solve.
• Benchmark SuperLU_MT thread scaling on representative FE systems.
• Integrate SuperLU_MT as a lazy-loaded selectable backend in the production global solver.
• Add persistent numerical refactorization (refact=YES) so symbolic/sparsity information can be reused when the KKT structure is unchanged.
• Preserve factor reuse for modified Newton and refactor only when the tangent changes.

Follow-up development

• Review and document the exact Newton/modified-Newton conventions used by the implementation and keep formulation, specification, and code terminology synchronized.
• Revisit the symmetry/nonsymmetry assumptions for the consistent tangent and assembled stiffness.

Run

The checked-in acceptance meshes can be regenerated with:

python examples/generate_elongated_block_gmsh.py --element hex8 --out examples/elongated_hex8.msh
python examples/generate_elongated_block_gmsh.py --element hex20 --out examples/elongated_hex20.msh
python examples/generate_periodic_cube_gmsh.py --out examples/periodic_cube_8x8x8.msh
python examples/generate_heterogeneous_periodic_cube_gmsh.py --elements-per-edge 4 --out examples/heterogeneous_periodic_cube_4x4x4.msh
python examples/generate_heterogeneous_periodic_cube_gmsh.py --elements-per-edge 8 --out examples/heterogeneous_periodic_cube_8x8x8.msh
python examples/generate_heterogeneous_periodic_cube_gmsh.py --elements-per-edge 16 --out examples/heterogeneous_periodic_cube_16x16x16.msh
python examples/generate_necking_bar_gmsh.py --out examples/necking_bar_quarter_hex8.msh
python examples/generate_necking_bar_gmsh.py --core-divisions 8 --radial-divisions 6 --nz 48 \
  --out examples/necking_bar_quarter_refined_hex8.msh
python examples/generate_necking_prism_gmsh.py --element hex8 \
  --out examples/necking_prism_small_hex8.msh
python examples/generate_necking_prism_gmsh.py --element hex20 \
  --out examples/necking_prism_small_hex20.msh
python examples/generate_necking_prism_gmsh.py --element hex8 \
  --nx 4 --ny 4 --nz 48 --out examples/necking_prism_refined_hex8.msh

Run a deck from the repository root:

python -m fe_solver examples/case_a_hex8.toml
python -m fe_solver examples/case_b_hex8_fbar.toml
python -m fe_solver examples/frame_objectivity_hex8.toml
python -m fe_solver examples/periodic_core_isochoric_hex8_fbar.toml
python -m fe_solver examples/periodic_core_shear_hex8_fbar.toml
python -m fe_solver examples/periodic_core_j2_matrix_isochoric_hex8_fbar.toml
python -m fe_solver examples/periodic_core_j2_matrix_shear_hex8_fbar.toml
python -m fe_solver examples/j2_necking_bar_hex8_fbar.toml
python -m fe_solver examples/j2_necking_prism_small_hex8.toml
python -m fe_solver examples/j2_necking_prism_small_hex8_fbar.toml

For a controlled partial run, add --stop-time 0.5. Relative mesh and output paths are resolved from the deck directory. The explicit equivalent command is python -m fe_solver solve DECK.

Decks use the serial scipy_splu reference backend by default. To build the native bridge and select SuperLU_MT for a run, use an explicit build parallelism and solver thread count:

cmake -S . -B build \
  -DCMAKE_BUILD_TYPE=Release \
  -Denable_internal_blaslib=ON \
  -DTHREAD_API=PTHREAD \
  -DPLAT=_PTHREAD

cmake --build build -j
ctest --test-dir build --output-on-failure

python -m fe_solver examples/case_a_hex8.toml \
  --solver-backend superlu_mt --num-threads N

The bridge is loaded from build/native/libsuperlu_mt_bridge.so only when superlu_mt is selected. SuperLU_MT has no implicit thread-count default; omitting --num-threads is an input error. The selected backend and its thread count are recorded in run_log.json.

Element evaluation defaults to the serial reference path. To reuse two local worker processes throughout the solve:

python -m fe_solver examples/j2_necking_prism_small_hex8_fbar.toml --num-processes 2

The unit of work is one complete element; the solver does not create a process per element or per quadrature point. Worker completion may be unordered, while global force and stiffness reduction remains in mesh order for reproducible floating-point behavior. The startup message and run_log.json report requested/effective process counts and available physical/logical CPUs. The solver never chooses all CPUs automatically, and each worker limits nested BLAS pools to one thread. For CPU-bound element kernels, physical-core count is the conservative starting point; benchmark before using additional hardware threads. Add --debug-timing when detailed element-phase and sparse-finalization timings are needed.

The startup summary also reports node and element counts, element-formulation totals, and the selected global-solution path. After the first real tangent assembly, it reports displacement/multiplier unknowns and the stored sparse entries in K and the complete KKT matrix, including the KKT fill percentage relative to dense storage. The same structured data is retained under analysis in run_log.json. Load and reaction fields remain separate from this model-size summary; full nodal constraint reactions are stored in run.h5.

Compare two completed HDF5 runs, including all numerical result fields and metadata, with:

python -m verification.compare_run_databases REFERENCE_RUN_H5 CANDIDATE_RUN_H5

The comparison normalizes only the output-directory entry in the embedded input deck and ignores the expected Git-revision provenance difference; both exceptions are listed in its JSON report.

Create one ParaView-readable temporal dataset after a run with:

python -m fe_solver postprocess examples/results/case_a_hex8/run.h5

This writes an .xdmf entry point and an .xdmf.h5 sidecar with connectivity and virtual-dataset views of time slices. Field values remain in the solver's run.h5; the views store selection metadata, not copies of the values. Keep all three files together (or retain their relative paths). No visualization files are produced during solution. Ordinary HDF references avoid XML HyperSlab compatibility problems in ParaView's XDMF3 reader.

When ParaView asks which implementation should open the file, select XDMF Reader T. XDMF Reader S flattens the spatial structure into one block, and the generic XDMF Reader is not compatible with these generated temporal files. Reader T presents one persistent multiblock branch per region/material, so Multiblock Inspector or Extract Block selections remain active as the time slider advances. VTK currently labels those temporal collection branches Block0, Block1, and so on rather than propagating their XDMF names; their deterministic order follows the HDF5 block keys shown under /mesh/blocks.

Stress and both strain measures are stored as six components in the order 11,22,33,12,23,13. Shear strains are tensor components, not doubled engineering shears. HDF5 datasets record this convention in their attributes. The XDMF exposes only these six-component arrays; select their individual components in ParaView's component selector. It does not add duplicate scalar aliases such as cauchy_stress_11. Six-component array magnitudes are not tensor Frobenius norms; contractions require double weighting of shear products and compatible stress/strain measures. Internal FE assembly conventions have not changed.

The results schema is now version 3; old nine-component databases must be regenerated by rerunning the solver. Regenerate the XDMF too, then reopen it in ParaView rather than reusing an old reader with cached field shapes.

Restart files now use schema version 3 and include accepted Gauss-point Cauchy stress for reporting at the restored time. Cold-start/restart reporting does not invoke constitutive evolution; the next positive-time increment obtains fresh stress and tangent normally. No material tangents or global stiffness are saved. Older restart files lack this stress cache and are rejected; new runs generate the required format without deck changes.

Periodic elastic-core/J2-matrix cube

The two periodic_core_j2_matrix_* decks use a 16 x 16 x 16 periodic cube with a centered 8 x 8 x 8 neo-Hookean core and a finite-strain J2/Voce matrix. Their reference elastic tangents are matched: both use mu = 80193.8 MPa and bulk modulus 164210 MPa; the neo-Hookean input therefore uses kappa = K - 2 mu/3 = 110747.4666666667 MPa. The core becomes the harder phase after the matrix yields without combining plastic contrast with an arbitrary initial elastic contrast. The J2 matrix uses centroidal Hex8-Fbar elements, while the compressible elastic core uses standard full-integration Hex8 elements on the same conforming nodal mesh. This intentionally exercises mixed material and element blocks; the core's matched reference constants correspond to approximately nu = 0.29, so standard Hex8 does not face a near-incompressible elastic-locking regime here.

The isochoric-uniaxial deck prescribes F11 = 1.2, F22 = F33 = 1/sqrt(1.2) at the endpoint. The shear deck prescribes F12 = 0.2. Both verify the reference-volume average of raw kinematic F against the imposed macro deformation. The full meshes are long-run examples rather than default acceptance jobs. From a synchronized repository root on a 16-core worker, run them one at a time with:

/usr/bin/time -v conda run --no-capture-output -n py3.14 \
  python -m fe_solver examples/periodic_core_j2_matrix_isochoric_hex8_fbar.toml \
  --num-processes 16 --debug-timing
/usr/bin/time -v conda run --no-capture-output -n py3.14 \
  python -m fe_solver examples/periodic_core_j2_matrix_shear_hex8_fbar.toml \
  --num-processes 16 --debug-timing

After either solve, create compact numerical histories and the temporal ParaView entry point with:

python -m verification.check_periodic_composite \
  examples/results/periodic_core_j2_matrix_isochoric_hex8_fbar/run.h5 \
  --output /tmp/periodic_isochoric.json --plot /tmp/periodic_isochoric.png
python -m fe_solver postprocess \
  examples/results/periodic_core_j2_matrix_isochoric_hex8_fbar/run.h5

The diagnostic plot contains macroscopic nominal stress derived from nodal constraint reactions and current-volume-weighted core, matrix, and whole-cell averages of stored centroidal Cauchy stress. These are response curves, not benchmark acceptance data. ParaView exposes J2 state arrays on the matrix block only; when both blocks are selected, coloring by such a field leaves the elastic core uncolored. Use Extract Block when inspecting the matrix state alone.

Uniaxial stretch followed by rotation

frame_objectivity_hex8.toml stretches the rod by 10% through t=0.1, then rotates it counterclockwise about z by 90 degrees through t=1. Only the two end faces are driven. They contract according to the exact neo-Hookean uniaxial-stress solution; the four lateral faces are traction-free. The interior and free-side displacements are solved, not prescribed. The material-specific neo_hook_uniaxial_rotation boundary path uses the named material's properties and remains exact between output events, including cutbacks.

python -m fe_solver examples/frame_objectivity_hex8.toml
python -m fe_solver postprocess examples/results/frame_objectivity_hex8/run.h5

Open examples/results/frame_objectivity_hex8/run.xdmf. Use Warp By Vector with displacement and scale 1. At t=0.1, component 11 of cauchy_stress is positive and component 22 is approximately zero. At t=1, those values exchange places. In the six-component cauchy_stress array, these are indices 0 and 1. Spatial Almansi strain rotates too, but transverse contraction means its lateral normal components are negative, not zero. Green-Lagrange strain remains fixed in the material frame during rotation.

Finite-strain J2 necking benchmark

j2_necking_bar_hex8_fbar.toml models one eighth of the classical imperfect circular tensile bar. Its quarter cross-section uses a regular central square and two transfinite outer sectors, avoiding collapsed axis elements and the skewed 45-degree surface cells of a one-block square-to-disk map. Half of the axial layers lie in the central third where necking begins. The prescribed 7 mm motion is applied to the end of the half-model. Material properties are written in the consistent mm--N--MPa system.

After solving, extract end reaction, middle radius, monitor-point radial displacement, and maximum equivalent plastic strain with:

python -m verification.check_j2_necking examples/results/j2_necking_bar_hex8_fbar/run.h5
python -m verification.check_j2_necking examples/results/j2_necking_bar_hex8_fbar/run.h5 \
  --plot /tmp/j2_necking_bar.png

The HDF5 database stores the inverse plastic metric and equivalent plastic strain at element centroids. J2 equivalent stress is intentionally not stored because it can be calculated from the six Cauchy-stress components in ParaView.

The 96-element j2_necking_prism_small_hex8_fbar.toml case is a qualitative nonlinear-solver diagnostic, not a replacement for the circular benchmark. j2_necking_prism_small_hex8.toml uses the identical mesh, material, loading, and solver controls as an unstabilized standard-Hex8 comparison. All J2 necking decks enable residual-based Newton backtracking. Invalid intermediate configurations reject only the current search length; pseudo-time cutback remains the fallback when no admissible decreasing search step exists.

Extract the square-prism response or compare the two formulations with:

python -m verification.check_j2_prism \
  examples/results/j2_necking_prism_small_hex8_fbar/run.h5 --summary-only
python -m verification.check_j2_prism \
  examples/results/j2_necking_prism_small_hex8_fbar/run.h5 \
  --compare examples/results/j2_necking_prism_small_hex8/run.h5 \
  --summary-only --output /tmp/j2_prism_comparison.json \
  --plot /tmp/j2_prism_comparison.png

The comparison reports reaction, transverse contraction, plastic localization, raw and material-seen Jacobian ranges, and cross-section mean-stress variation. Its optional six-panel plot places the physical response and volumetric diagnostics against prescribed end elongation. The standard Hex8 result is a locking control; convergence alone does not make it a reference solution. A plot from the first completed paired 96-element runs is retained under verification/j2_prism_results/ as a qualitative regression artifact.

The follow-on formulation study is expressed as orthogonal variants of that base deck rather than duplicated TOML files. List the cases with:

python -m verification.run_j2_formulation_study --list-cases

The matrix contains the completed 2x2x24 Hex8/Hex8-Fbar pair, a paired 4x4x48 mesh refinement, a coarse pair with the bulk modulus halved from 164210 to 82105 MPa, a 2x2x24 serendipity Hex20 case, and the 960/7680-element circular-bar mesh pair. The half-bulk case changes the elastic Poisson ratio from approximately 0.290 to 0.132 while retaining every shear/plastic parameter. Run one new case with, for example:

python -m verification.run_j2_formulation_study \
  --case soft_bulk_hex8_fbar --num-processes 2 --debug-timing --stop-time 0.5

New variant outputs are isolated under examples/results/j2_formulation_study/; the three baseline cases reuse their existing result directories. Use check_j2_prism to compare any completed pair directly from HDF5; this avoids ParaView interpolation, warp scaling, camera, and color-range ambiguity. The refined cases are intentionally expensive and should use the long-run workflow rather than the ordinary test suite.

Standalone material-point characterization and evolved-state tangent diagnostics are available through:

python -m verification.run_material_point_characterization
python -m verification.check_evolved_tangents DECK RESTART RUN_H5 --trial-time 0.405 --output report.json
python -m verification.diagnose_newton_conditioning DECK RESTART RUN_H5 --committed-time 0.43 --attempt-time 0.44 --line-search none --output report.json

The first command writes reviewable J2 uniaxial and simple-shear CSV, compressed NumPy, PNG, and JSON results under verification/material_point_results/. The conditioning utility deliberately uses a dense null-space/SVD and is restricted to small diagnostics; production assembly and KKT solution remain sparse.

Numba is a required project dependency. The shared material-to-spatial Truesdell tangent transform and J2 return map are always compiled in production. The synthetic material benchmark compares the J2 dispatcher with its uncompiled .py_func reference:

python -m verification.run_j2_numba_feasibility material

The earlier full-solver Python/Numba comparison databases and their compare-prisms reader remain as historical adoption evidence. Backend selection is no longer exposed by the solver or formulation-study commands.

Multiphase TRIP results

The accepted 8A/56F 8³ periodic simple-shear and isochoric-uniaxial runs, input decks, mesh/property data, and SVG/numerical reports are archived in verification/multiphase_trip_results/. Git LFS stores the simulation databases/logs and binary property/history files. Retrieve them with git lfs pull when using a clone. The external constitutive source and compiled library are not distributed in this repository; reading results and regenerating reports does not require them. The archive README records runtime solver overrides, provenance limitations, and reproduction commands.

Verify

The test suite uses the standard library runner, so no separate test dependency is required:

python -m unittest discover -s tests -v

It checks shape functions, material and element tangents, homogeneous F-bar equivalence, Gmsh periodic maps, spanning-tree constraints, nonsymmetric whole-KKT solves, event merging, restart compatibility, forced-failure cutback/rollback, and serial-versus-process assembly equivalence.

Run the complete regression and featured-example acceptance suite with:

python -m verification.run_acceptance

The substantially larger necking benchmark is selectable explicitly and is not part of the default quick suite:

python -m verification.run_acceptance --deck j2_necking_bar_hex8_fbar.toml

When ParaView is installed, also test its real readers and Warp By Vector against the HDF5 database:

python -m verification.check_paraview examples/results/frame_objectivity_hex8/run.h5 --pvpython /path/to/paraview/bin/pvpython

This checks both XDMF3 reader variants at every saved time, all field components, and warped node positions.

Package layout

• fe_solver/shape.py, quadrature.py: interpolation and integration
• fe_solver/materials.py, material_point.py: named material API, material models, and prescribed-F verification driver
• fe_solver/tangents.py: shared material-to-spatial tangent transformation
• fe_solver/elements.py: standard and F-bar element kernels
• fe_solver/mesh.py, constraints.py: Gmsh import and affine constraints
• fe_solver/preprocess.py, assembly.py, execution.py, solver.py: validated setup, deterministic serial/process element execution, sparse assembly, and nonlinear solution
• fe_solver/io.py, config.py: HDF5 result/restart I/O and TOML/event handling
• fe_solver/postprocess.py: temporal XDMF construction from a completed HDF5 run