Cauchy: 2D Finite Element Structural Solver
How does a bridge bear its load, a wing resist lift, or a building withstand wind? To predict this without building expensive physical prototypes, engineers use Finite Element Analysis to discretize structures into thousands of small elements and solve for deformation and stress.
Commercial solvers like ANSYS and Abaqus are powerful tools, but their black-box nature hides the underlying mechanics. To better understand FEA at the implementation level, I set out to build a solver from scratch.
Cauchy is a custom 2D structural solver written in C++20. It uses bilinear quad (Q4) and serendipity (Q8) elements with Gauss quadrature, sparse matrix assembly, direct/iterative solvers, and adaptive mesh refinement driven by the Zienkiewicz-Zhu error estimator. The solver validates against 6 classical benchmark cases with analytical references.
Simulation Overview
All 7 validation cases solved with the same C++20 solver core. Each case targets a specific aspect of structural mechanics and validates against analytical or reference solutions. Mesh convergence studies are run at 5 resolutions (8x8 through 128x128).
| Case | Elements | Mesh | DOFs | Material | Plane | Max Disp. | Max Stress | Validates |
|---|---|---|---|---|---|---|---|---|
| Cantilever Beam | Q8 256 | 32x8 | 1,698 | Steel (E=200 GPa) | Stress | 0.139 mm | 8.20 MPa | PL^3/(3EI) deflection, My/I stress |
| Cook's Membrane | Q8 1,024 | 32x32 | 6,402 | Aluminum (E=1 GPa) | Stress | 9.50 mm | 0.21 MPa | Shear-bending coupling, ~13.68 mm reference |
| Patch Test | Q4 16 | 4x4 | 50 | Steel (E=200 GPa) | Stress | 1.00 mm | 208 MPa | Constant stress recovery (element pass/fail) |
| Plate with Hole | Q8 15,551 | 128x128 | 94,332 | Steel (E=200 GPa) | Stress | 0.028 mm | 7.66 MPa | Kirsch SCF~3.0 at circular hole edge |
| L-Bracket | Q4 12,159 | 128x128 | 24,832 | Steel (E=200 GPa) | Stress | 0.001 mm | 0.45 MPa | Re-entrant corner singularity, mesh refinement |
| Michell Truss | 3 | 3 bars | 6 | Steel (E=200 GPa) | Stress | 0.044 mm | -- | Bar element assembly, hand calculation check |
| Thick Cylinder | Q4 13,108 | 128x128 | 26,730 | Steel (E=200 GPa) | Strain | 0.007 mm | 1.80 MPa | Lame solution, internal pressure + thermal |
Element types: 2-node truss, Q4 4-node bilinear quad (2x2 Gauss), Q8 8-node serendipity (3x3 Gauss). All cases use steel (E=200 GPa, nu=0.3) unless noted. Convergence studies run at 8, 16, 32, 64, 128 resolution.
Validation Cases
Cantilever Beam
Clamped beam with tip load. Validates element formulation against PL^3/(3EI).
View case →
Cook's Membrane
Trapezoidal panel under distributed shear. Tests Q4 bending and shear locking.
View case →
Patch Test
Constant stress recovery. Mandatory element verification for any new formulation.
View case →
Plate with Hole
Kirsch solution: stress concentration factor of 3 at a circular hole edge.
View case →
L-Bracket
Re-entrant corner stress concentration. Demonstrates mesh refinement needs.
View case →
Michell Truss
3-bar truss structure. Validates bar element assembly and hand calculations.
View case →
Thick Cylinder (Lame)
Annular cylinder under internal pressure and thermal gradient. Validates plane strain formulation against Lame solution.
View case →
3D Cantilever Beam
3D hexahedral (H8) element validation against Euler-Bernoulli beam theory. Demonstrates 3D solid element formulation.
View case →
3D Plate with Hole
3D stress concentration around a circular hole in a thick plate. SCF varies through thickness.
View case →
3D Thick Cylinder (Lame)
3D Lame solution validation for thick-walled cylinder under internal pressure. All 6 stress components.
View case →
PINN Surrogate Model
Physics-Informed Neural Network for real-time 2D elasticity prediction. Browser deployment via ONNX Runtime.
View PINN →Key Features
Bar + Q4 + Q8 Elements
Two-node truss, 4-node bilinear quad, and 8-node serendipity quad with full Gauss integration (2x2 for Q4, 3x3 for Q8).
Dual Solvers
Cholesky direct solver for small/medium systems. Conjugate Gradient with auto-selected preconditioner (IC(0), SSOR, Block Jacobi) for large sparse systems.
Sparse Assembly
COO format for natural element assembly, converted to CSR for fast matrix-vector products.
Stress Recovery
Element-centered stress, node-averaged smoothing, Von Mises and principal stresses. Q8 uses 2x2 Gauss averaging.
ZZ Error Estimator
Zienkiewicz-Zhu superconvergent patch recovery (SPR) for a posteriori error estimation. Drives adaptive mesh refinement.
Adaptive h-Refinement
Red-green refinement with element marking. Concentrates DOFs where error is highest, achieving 10-14x node reduction.
Mesh Convergence
h-refinement studies with GCI (Grid Convergence Index) per ASME V&V 20-2009.
OpenMP Parallel
Element assembly and stress recovery loops parallelized with OpenMP.
Shear Locking Fix
Selective Reduced Integration (SRI) and B-Bar methods eliminate Q4 shear locking in bending-dominated problems.
Contact Mechanics
Node-to-surface frictionless contact with penalty constraint enforcement. Master-slave formulation with gap function detection.
Geometric Nonlinearity
Total Lagrangian Newton-Raphson with Green-Lagrange strain, consistent tangent stiffness, and load stepping.
Dynamic Analysis
Newmark-beta time integration, consistent mass matrices, Rayleigh damping, and subspace iteration modal analysis.
Architecture
| Module | Purpose | Key Types |
|---|---|---|
fea_types.hpp | Core types, enums, globals | Material, Node, Mesh, CaseType |
elements.hpp | Element stiffness matrices | BarElement, Q4Element, Q8Element, T3Element |
sparse.hpp | Sparse matrix formats | COOMatrix, CSRMatrix, DenseMatrix |
solver.hpp | Linear solvers | CholeskySolver, CGSolver |
mesh.hpp | Mesh generation | generate_structured_quad, generate_structured_quad8, generate_lbracket |
postprocess.hpp | Stress recovery, JSON output | ElementStress, write_meta_json |
adaptivity.hpp | ZZ error estimator, adaptive refinement | SPRResult, ElementError, refine_mesh, adaptive_loop |
convergence.hpp | GCI computation, h-refinement | GCIResult, run_study |
fea.hpp | Assembly + solve pipeline | assemble, solve, compute_strain_energy |
Quick Start
cmake -B build && cmake --build build -j$(sysctl -n hw.ncpu)
./build/FEA_Cantilever 32 — cantilever beam (32x8 mesh)
./build/FEA_Cook 32 — Cook's membrane (32x32 mesh)
./build/FEA_PlateHole 16 --q8 — plate with hole (Q8 elements)
./build/FEA_AdaptHole 8 --iters 4 — adaptive refinement
./build/FEA_Tests — run 57 Google Test cases