libfem

Finite element method library: a framework for structural mechanics, modal analysis, and nonlinear problems

libfem is a flexible, extensible FEM framework that provides simulation capabilities for structural mechanics, heat transfer, modal analysis, and nonlinear problems.

Features

12 Element Types - Full library of 1D, 2D, and 3D elements (linear and quadratic)
Multiple Physics - Diffusion, elasticity, and mass matrix integrators
Solver Wrappers - Direct, iterative, and eigenvalue solvers
Modal Analysis - Natural frequency extraction and mass participation factors
Dynamic Analysis - Newmark-β, HHT-α, and explicit time integration
Rayleigh Damping - Proportional damping with coefficient computation
Nonlinear Analysis - Neo-Hookean hyperelasticity, J2 plasticity, Newton-Raphson solver
Post-Processing - Stress recovery, failure criteria (von Mises, Tresca), safety factors
Results Export - VTK for visualization in a VTK viewer, CSV/JSON for reporting
Sparse Assembly - Efficient triplet-based sparse matrix construction
Quadrature Rules - Gauss-Legendre and simplex quadrature for accurate integration
Header-Only Dependencies - Eigen 3.4.0 and Spectra 1.0.1 bundled

Element Library

1D Elements

Element	Nodes	Order	Use Case
Bar2	2	Linear	Truss members, axial loading
Bar3	3	Quadratic	Higher accuracy 1D problems

2D Elements

Element	Nodes	Order	Use Case
Tri3	3	Linear	Simple triangular meshes
Tri6	6	Quadratic	Curved boundaries, higher accuracy
Quad4	4	Bilinear	Regular rectangular regions
Quad8	8	Serendipity	Curved edges without center node
Quad9	9	Biquadratic	Full quadratic accuracy

3D Elements

Element	Nodes	Order	Use Case
Tet4	4	Linear	Simple tetrahedral meshes
Tet10	10	Quadratic	Complex geometry, higher accuracy
Hex8	8	Trilinear	Regular hexahedral regions
Hex20	20	Serendipity	Curved surfaces
Hex27	27	Triquadratic	Full quadratic accuracy in 3D

Quick Start

Basic Mesh Setup

#include "fem/core/Mesh.hpp"
#include "fem/core/DOFHandler.hpp"
#include "fem/elements/Element.hpp"

using namespace fem;

// Create a 2D mesh
Mesh mesh(2);

// Add nodes
mesh.addNode(0.0, 0.0);  // Node 0
mesh.addNode(1.0, 0.0);  // Node 1
mesh.addNode(0.5, 1.0);  // Node 2

// Add a triangular element
mesh.addElement(ElementType::Tri3, {0, 1, 2});

Applying Boundary Conditions

// Dirichlet BC: fix temperature at node 0 to 0.0
mesh.addDirichletBC(0, 0, 0.0);

// Dirichlet BC: fix temperature at node 1 to 100.0
mesh.addDirichletBC(1, 0, 100.0);

// Neumann BC: apply flux at node 2
mesh.addNeumannBC(2, 0, 10.0);

Assembling a Stiffness Matrix

#include "fem/integrators/DiffusionIntegrator.hpp"

// Setup DOF handler (1 DOF per node for scalar problems)
DOFHandler dof(mesh, 1);
dof.applyDirichletBCs();

// Create assembler
Assembler assembler(dof.totalDOFs());

// Create integrator with conductivity k = 1.0
DiffusionIntegrator diffusion(1.0);

// Get element and compute element matrix
auto elem = createElement(ElementType::Tri3);

// Extract node coordinates for element 0
MatrixX nodeCoords(3, 3);
for (int i = 0; i < 3; ++i) {
    Vector3 pos = mesh.node(mesh.element(0).nodeIndices[i]);
    nodeCoords.col(i) = pos;
}

// Compute and assemble element matrix
MatrixX Ke = diffusion.computeElementMatrix(*elem, nodeCoords);
auto elemDOFs = dof.elementDOFs(0);
assembler.addElementMatrix(elemDOFs, Ke);

// Build global stiffness matrix
SparseMatrix<double> K = assembler.buildMatrix();

Using Different Element Types

// Create elements via factory
auto bar2 = createElement(ElementType::Bar2);
auto tri6 = createElement(ElementType::Tri6);
auto tet10 = createElement(ElementType::Tet10);
auto hex27 = createElement(ElementType::Hex27);

// Or instantiate directly
Tri3Element tri3;
Quad8Element quad8;
Hex20Element hex20;

Physics Integrators

DiffusionIntegrator

For heat conduction and scalar diffusion problems:

DiffusionIntegrator diffusion(conductivity);
MatrixX Ke = diffusion.computeElementMatrix(element, nodeCoords);

ElasticityIntegrator

For structural mechanics:

Material mat(E, nu);  // Young's modulus, Poisson's ratio

// Choose analysis type
ElasticityIntegrator elasticity(mat, ElasticityType::PlaneStress);
ElasticityIntegrator elasticity(mat, ElasticityType::PlaneStrain);
ElasticityIntegrator elasticity(mat, ElasticityType::ThreeD);

MatrixX Ke = elasticity.computeElementMatrix(element, nodeCoords);

MassIntegrator

For dynamic analysis:

MassIntegrator mass(density);
MatrixX Me = mass.computeElementMatrix(element, nodeCoords);

Solvers

LinearSolver (Direct)

For symmetric positive-definite systems (SimplicialLDLT):

#include "fem/solvers/LinearSolver.hpp"

LinearSolver solver(K);  // Factorize
if (solver.success()) {
    VectorX u = solver.solve(F);  // Solve Ku = F
}

LUSolver (General)

For non-symmetric or indefinite systems:

LUSolver solver(A);
VectorX x = solver.solve(b);

IterativeSolver

Conjugate gradient with preconditioner support for large systems:

#include "fem/solvers/IterativeSolver.hpp"

IterativeSolver solver(Preconditioner::Diagonal);  // Or IncompleteCholesky
solver.setTolerance(1e-10);
solver.setMaxIterations(1000);
solver.compute(K);

VectorX u = solver.solve(F);

if (solver.converged()) {
    std::cout << "Converged in " << solver.iterations() << " iterations\n";
}

EigenSolver

Generalized eigenvalue problem Kx = λMx using Spectra:

#include "fem/solvers/EigenSolver.hpp"

EigenSolver solver(K, M);
solver.setNumModes(10);
solver.compute();

VectorX lambda = solver.eigenvalues();
VectorX omega = solver.angularFrequencies();  // rad/s
VectorX freq = solver.frequencies();          // Hz
MatrixX modes = solver.modeShapes();          // Mass-normalized

ModalAnalysis

High-level modal analysis with mass participation factors:

#include "fem/solvers/ModalAnalysis.hpp"

ModalAnalysis modal(K, M);
modal.setNumModes(10);
modal.compute();

// Get results
for (const auto& result : modal.results()) {
    std::cout << "Mode " << result.modeNumber
              << ": " << result.frequency << " Hz"
              << " (T = " << result.period << " s)\n";
}

// Mass participation in X direction (for 3-DOF per node structures)
MassParticipation mpf = modal.massParticipationX(3);
std::cout << "Cumulative mass participation: "
          << mpf.cumulative(mpf.cumulative.size()-1) * 100 << "%\n";

Dynamic Analysis

RayleighDamping

Proportional damping (C = αM + βK):

#include "fem/solvers/RayleighDamping.hpp"

// Method 1: Specify coefficients directly
RayleighDamping damping(alpha, beta);

// Method 2: Compute from target damping ratios at two frequencies
RayleighDamping damping = RayleighDamping::fromDampingRatios(
    omega1, xi1,   // 5% damping at 10 rad/s
    omega2, xi2);  // 3% damping at 100 rad/s

// Method 3: Single-frequency approximations
RayleighDamping massD = RayleighDamping::massProportional(omega, 0.05);
RayleighDamping stiffD = RayleighDamping::stiffnessProportional(omega, 0.05);

// Build damping matrix
SparseMatrix C = damping.computeDampingMatrix(M, K);

// Query damping ratio at any frequency
Scalar xi = damping.dampingRatio(omega);

NewmarkSolver (Implicit)

Unconditionally stable implicit time integration:

#include "fem/solvers/NewmarkSolver.hpp"

NewmarkSolver solver(M, C, K);
solver.setMethod(NewmarkMethod::AverageAcceleration);  // β=1/4, γ=1/2
// Or: LinearAcceleration (β=1/6, γ=1/2)
// Or: Custom parameters
solver.setTimeStep(0.001);

VectorX u0 = VectorX::Zero(n);
VectorX v0 = VectorX::Zero(n);
solver.setInitialConditions(u0, v0);
solver.initialize();

for (double t = 0; t < tEnd; t += dt) {
    VectorX F = computeLoad(t + dt);
    solver.step(F);

    VectorX u = solver.displacement();
    VectorX v = solver.velocity();
    VectorX a = solver.acceleration();
}

HHTSolver (Numerical Damping)

HHT-α method with controllable numerical damping for high frequencies:

#include "fem/solvers/HHTSolver.hpp"

HHTSolver solver(M, C, K);
solver.setAlpha(-0.1);  // α ∈ [-1/3, 0], more negative = more damping
solver.setTimeStep(0.001);
solver.setInitialConditions(u0, v0);
solver.initialize();

// Query numerical properties
Scalar rho = solver.spectralRadius();  // ρ∞

for (double t = 0; t < tEnd; t += dt) {
    solver.step(F);
}

ExplicitSolver (Conditionally Stable)

Central difference method for fast explicit integration:

#include "fem/solvers/ExplicitSolver.hpp"

// Lump mass matrix for efficiency
SparseMatrix Mlumped = ExplicitSolver::lumpMass(M);

ExplicitSolver solver(Mlumped, C, K);

// Check stability requirement
Scalar dtCrit = ExplicitSolver::criticalTimeStep(omegaMax);
// Or estimate from element size
Scalar dtCrit = ExplicitSolver::criticalTimeStepFromElement(minSize, waveSpeed);

Scalar dt = 0.8 * dtCrit;  // Safety factor
solver.setTimeStep(dt);

if (!solver.isStable(omegaMax)) {
    // Time step too large!
}

solver.setInitialConditions(u0, v0);
solver.initialize();

// Explicit methods are very fast per step
for (int step = 0; step < nSteps; ++step) {
    solver.step(F);
}

Shape Functions

Each element provides:

shapeFunctions(xi) - Evaluate N_i at parametric coordinates
shapeFunctionDerivatives(xi) - Evaluate dN_i/dxi
quadraturePoints(order) - Get integration points

Quad9Element quad9;
VectorX xi(2);
xi << 0.5, 0.25;

VectorX N = quad9.shapeFunctions(xi);      // 9 values
MatrixX dN = quad9.shapeFunctionDerivatives(xi);  // 9x2 matrix
auto qpts = quad9.quadraturePoints(3);    // 3x3 = 9 points

Reference Domains

Element Type	Domain	Description
Bar	ξ ∈ [-1, 1]	Standard 1D interval
Quad	ξ,η ∈ [-1, 1]	Standard biunit square
Hex	ξ,η,ζ ∈ [-1, 1]	Standard biunit cube
Tri	ξ,η ≥ 0, ξ+η ≤ 1	Standard simplex
Tet	ξ,η,ζ ≥ 0, ξ+η+ζ ≤ 1	Standard simplex

Build & Test

# Build the test suite and run it to its contract
kiln test test/testFEM

# Run with verbose output
./debug/testFEM/testFEM -v2

# Run specific element tests
./debug/testFEM/testFEM testHex27ShapeFunctionsAtNodes -v2

The program's path is that of a kiln build with no output.dir (The build convention).

Dependencies

Both dependencies are bundled (header-only, MPL2 license):

Library	Version	Purpose
Eigen	3.4.0	Linear algebra, sparse matrices
Spectra	1.0.1	Large-scale eigenvalue problems

Architecture

src/libs/libfem/fem/
├── core/           # Core data structures
│   ├── Mesh.hpp/cpp
│   ├── DOFHandler.hpp/cpp
│   ├── Quadrature.hpp/cpp
│   └── Types.hpp
├── elements/       # Element implementations
│   ├── Element.hpp/cpp
│   ├── Bar2.hpp/cpp, Bar3.hpp/cpp
│   ├── Tri3.hpp/cpp, Tri6.hpp/cpp
│   ├── Quad4.hpp/cpp, Quad8.hpp/cpp, Quad9.hpp/cpp
│   ├── Tet4.hpp/cpp, Tet10.hpp/cpp
│   └── Hex8.hpp/cpp, Hex20.hpp/cpp, Hex27.hpp/cpp
├── integrators/    # Physics modules
│   ├── Integrator.hpp
│   ├── DiffusionIntegrator.hpp/cpp
│   ├── MassIntegrator.hpp/cpp
│   └── ElasticityIntegrator.hpp/cpp
└── solvers/        # Solver wrappers
    ├── LinearSolver.hpp/cpp
    ├── IterativeSolver.hpp/cpp
    ├── EigenSolver.hpp/cpp
    ├── ModalAnalysis.hpp/cpp
    ├── RayleighDamping.hpp/cpp
    ├── NewmarkSolver.hpp/cpp
    ├── HHTSolver.hpp/cpp
    └── ExplicitSolver.hpp/cpp

Nonlinear Analysis

Strain Measures

Utilities for large deformation analysis:

#include "fem/nonlinear/StrainMeasures.hpp"

Matrix3 F = /* deformation gradient */;

// Strain tensors
Matrix3 E = StrainMeasures::greenLagrangeStrain(F);      // E = 0.5(F^T F - I)
Matrix3 C = StrainMeasures::rightCauchyGreen(F);          // C = F^T F
Matrix3 b = StrainMeasures::leftCauchyGreen(F);           // b = F F^T

// Invariants and principal values
Vector3 inv = StrainMeasures::invariants(C);              // I1, I2, I3
Vector3 stretches = StrainMeasures::principalStretches(F); // λ1, λ2, λ3

// Voigt notation conversions (for constitutive relations)
Vector6 Ev = StrainMeasures::strainToVoigt(E);           // Engineering shear
Matrix3 E2 = StrainMeasures::voigtToStrain(Ev);

// Stress push-forward/pull-back
Matrix3 sigma = StrainMeasures::pushForwardStress(S, F);  // PK2 → Cauchy
Matrix3 S2 = StrainMeasures::pullBackStress(sigma, F);    // Cauchy → PK2

Neo-Hookean Material

Hyperelastic material model for rubber-like behavior:

#include "fem/nonlinear/NeoHookean.hpp"

// Create from Lamé parameters
NeoHookean mat(lambda, mu);

// Or from engineering constants
NeoHookean mat = NeoHookean::fromENu(1000.0, 0.3);  // E, nu

// Compute stress from deformation gradient
MaterialState state;
Matrix3 F = /* deformation gradient */;
Matrix3 S = mat.computeStressTensor(F);               // 2nd Piola-Kirchhoff
Vector6 Sv = mat.computeStressFromF(F, state);        // Voigt form

// Compute tangent stiffness
Matrix6 D = mat.computeTangentFromF(F);               // Material tangent

// Strain energy density
Scalar W = mat.computeStrainEnergyFromF(F);

J2 Plasticity

von Mises plasticity with isotropic hardening:

#include "fem/nonlinear/J2Plasticity.hpp"

// Create with E, nu, initial yield stress
J2Plasticity mat(200e3, 0.3, 250.0);

// Set hardening law
mat.setHardening(HardeningType::Linear, 1000.0);       // Linear: σ_y + H*ε_p
mat.setHardening(HardeningType::Saturation, 100.0, 10.0); // Voce: σ_y + Q(1-exp(-b*ε_p))
mat.setHardening(HardeningType::Perfect);              // No hardening

// Compute stress (performs return mapping internally)
MaterialState state;  // Tracks plastic strain history
Vector6 strain = /* total strain in Voigt form */;
Vector6 stress = mat.computeStress(strain, state);

// Check yield state
if (state.yielded) {
    Scalar epEq = state.equivalentPlasticStrain;
}

// Static utilities
Scalar sigmaEq = J2Plasticity::vonMisesStress(stress);
Vector6 s = J2Plasticity::deviatoricStress(stress);

Newton-Raphson Solver

Iterative nonlinear solver with line search:

#include "fem/nonlinear/NewtonRaphson.hpp"

NewtonRaphson nr;
nr.setTolerance(1e-6);
nr.setMaxIterations(20);
nr.setConvergenceCriterion(ConvergenceCriterion::Force);
// Options: Force, Displacement, Energy, ForceAndDisp, ForceOrDisp

nr.setLineSearch(LineSearchMethod::Backtracking);
// Options: None, Backtracking, Bisection

// Define problem
VectorX u = VectorX::Zero(n);

auto residual = [&](const VectorX& x) -> VectorX {
    return Fext - computeInternalForce(x);  // R = F_ext - F_int
};

auto tangent = [&](const VectorX& x) -> SparseMatrix {
    return computeTangentStiffness(x);      // K = dF_int/du
};

// Solve
NewtonResult result = nr.solve(u, residual, tangent);

if (result.converged) {
    std::cout << "Converged in " << result.iterations << " iterations\n";
} else {
    std::cout << "Failed to converge, residual = " << result.residualNorm << "\n";
}

Arc-Length Method

For tracing equilibrium paths through limit points:

#include "fem/nonlinear/ArcLength.hpp"

ArcLength arc;
arc.setArcLength(0.1);
arc.setTolerance(1e-6);
arc.setMaxIterations(20);
arc.setType(ArcLengthType::Spherical);  // Or Cylindrical, Updated

VectorX u = VectorX::Zero(n);
Scalar lambda = 0.0;  // Load factor
VectorX Fref = /* reference load */;

auto internalForce = [&](const VectorX& x) -> VectorX {
    return computeInternalForce(x);
};

auto tangent = [&](const VectorX& x) -> SparseMatrix {
    return computeTangentStiffness(x);
};

// Step along equilibrium path
while (lambda < lambdaMax) {
    ArcLengthResult result = arc.step(u, lambda, Fref, internalForce, tangent);
    if (!result.converged) break;

    // lambda and u are updated automatically
    recordEquilibriumPoint(lambda, u);
}

High-Level FEModel Pipeline

For complete end-to-end analysis workflows, the FEModel class provides a simplified interface:

Basic Usage

#include "fem/model/FEModel.hpp"

using namespace fem;
using namespace fem::model;

// Create model and load geometry
FEModel model;
model.loadSTL("part.stl");  // Automatically tetrahedralizes

// Set material
model.setMaterial(EngineeringMaterial::aluminum());

// Apply boundary conditions
auto bottomNodes = model.getNodesOnPlane(Vector3(0, 0, 0), Vector3(0, 0, 1));
auto topNodes = model.getNodesOnPlane(Vector3(0, 0, 10), Vector3(0, 0, 1));

model.addLoad(Load::fixNodes(bottomNodes));
model.addLoad(Load::applyForce(topNodes, Vector3(0, 0, -1000)));

// Solve and check results
if (model.solve()) {
    std::cout << model.summary().toStdString() << std::endl;
    std::cout << "Max displacement: " << model.results().maxDisplacement() << " m" << std::endl;
    std::cout << "Min FOS: " << model.results().safety.minFOS << std::endl;

    // Export for visualization
    model.exportVTK("results.vtk");
    model.exportCSV("results.csv");
}

Loading a mesh from libmesh

A mesh that already exists in memory, for example the Tet4 volume mesh that the Construct SDF export writes, goes in through loadMesh(). This skips the STL import and the tetrahedralization that the export already did. Both precisions are accepted. A float mesh is widened to double before conversion.

#include "fem/model/FEModel.hpp"
#include "mesh/io/IORegistry.hpp"

auto loaded = mesh::io::IORegistry<float>::instance().load("volume.vtk");
fem::FEModel model;
if (!model.loadMesh(loaded.value())) {
    qWarning() << model.errorString();
}

loadMesh() returns false and sets errorString() if the source mesh has no vertices or no elements. A surface mesh loads, but it will not solve as a volume. The end-to-end case is testMeshExport::testFEModelFromSDF.

Engineering Materials

Pre-defined material properties with yield strength for safety factor calculations:

// Built-in materials
EngineeringMaterial steel = EngineeringMaterial::steel();     // E=210 GPa, nu=0.3, σy=250 MPa
EngineeringMaterial alum = EngineeringMaterial::aluminum();   // E=70 GPa, nu=0.33, σy=270 MPa
EngineeringMaterial titan = EngineeringMaterial::titanium();  // E=116 GPa, nu=0.34, σy=880 MPa
EngineeringMaterial pla = EngineeringMaterial::pla();         // E=3.5 GPa, nu=0.36, σy=50 MPa
EngineeringMaterial abs = EngineeringMaterial::abs();         // E=2.3 GPa, nu=0.35, σy=40 MPa

// Custom material
EngineeringMaterial custom(E, nu, yieldStrength, density, "Custom Material");

Load Types

// Fix all DOFs at nodes
Load::fixNodes(nodeIndices);

// Apply force distributed across nodes
Load::applyForce(nodeIndices, Vector3(fx, fy, fz));

// Prescribed displacement
Load::applyDisplacement(nodeIndices, Vector3(dx, dy, dz));

Node Selection Helpers

// Get nodes in a bounding box
auto nodes = model.getNodesInBox(Vector3(xMin, yMin, zMin), Vector3(xMax, yMax, zMax));

// Get nodes on a plane (with tolerance)
auto planeNodes = model.getNodesOnPlane(pointOnPlane, planeNormal, tolerance);

STL Import

For loading surface meshes from STL files:

#include "fem/io/STLImporter.hpp"

STLImporter importer;

// Import as surface mesh (Tri3 shell elements)
auto surfaceMesh = importer.importSurface("part.stl");

// Import and tetrahedralize to volume mesh (Tet4 elements)
auto volumeMesh = importer.importVolume("part.stl");

// With tetrahedralization options
TetrahedralizeOptions opts;
opts.maxEdgeLength = 0.01;    // Maximum edge length
opts.qualityRatio = 2.0;      // Quality ratio
opts.preserveSurface = true;  // Keep original surface triangles
auto mesh = importer.importVolume("part.stl", opts);

// Check for errors
if (!mesh) {
    std::cerr << importer.errorString().toStdString() << std::endl;
}

Note: The built-in tetrahedralization uses a centroid-fan method suitable for convex shapes. For complex non-convex geometries, consider using an external tetrahedral mesher.

Post-Processing & Results

Stress Recovery

Compute stress and strain fields from displacement solutions:

#include "fem/postprocess/StressRecovery.hpp"

// Create stress recovery with material properties
StressRecovery sr(E, nu);  // Young's modulus, Poisson's ratio

// For 2D analysis
StressRecovery sr2D(E, nu, ElasticityIntegrator::ProblemType::PlaneStress);

// Compute stress from displacement solution
StressResult stress = sr.compute(mesh, displacement);

// Access results
VectorX vonMises = stress.vonMises;           // Per-node von Mises stress
VectorX tresca = stress.tresca;               // Per-node Tresca stress
MatrixX sigma = stress.stress;                // Full stress tensor (nodes x 6, Voigt)
MatrixX eps = stress.strain;                  // Full strain tensor
MatrixX principals = stress.principalStress;  // Principal stresses (nodes x 3)

// Static utilities for stress analysis
Scalar vm = StressRecovery::vonMisesStress(stressVoigt);
Scalar tr = StressRecovery::trescaStress(stressVoigt);
Vector3 sigmaP = StressRecovery::principalStresses(stressVoigt);

// Get principal directions
Vector3 principals;
Matrix3 directions;
StressRecovery::principalStressesWithDirections(stressVoigt, principals, directions);

Failure Criteria

Check structural integrity using standard failure theories:

#include "fem/postprocess/FailureCriteria.hpp"

// Available criteria
FailureCriteria::Criterion::VonMises      // Distortion energy (ductile metals)
FailureCriteria::Criterion::Tresca        // Maximum shear (conservative)
FailureCriteria::Criterion::MaxPrincipal  // Maximum tensile stress (brittle)
FailureCriteria::Criterion::MinPrincipal  // Maximum compressive stress
FailureCriteria::Criterion::Rankine       // Maximum normal stress
FailureCriteria::Criterion::MaxStrain     // Maximum principal strain

// Check if failure occurred
bool failed = FailureCriteria::checkFailure(stressVoigt, criterion, yieldStrength);

// Compute equivalent stress
Scalar sigmaEq = FailureCriteria::equivalentStress(stressVoigt, criterion);

// Compute safety margin (positive = safe)
Scalar margin = FailureCriteria::computeMargin(stressVoigt, criterion, limit);

// Compute factor of safety
Scalar fos = FailureCriteria::computeFOS(stressVoigt, criterion, limit);

Safety Factor Analysis

Compute factor of safety across the entire structure:

#include "fem/postprocess/SafetyFactor.hpp"

// Create with yield and ultimate strengths
SafetyFactor sf(yieldStrength, ultimateStrength);

// Compute yield FOS using von Mises criterion
SafetyFactorResult fosResult = sf.computeYieldFOS(stress, mesh);

// Or use Tresca criterion
SafetyFactorResult fosResult = sf.computeYieldFOS(stress, mesh,
    FailureCriteria::Criterion::Tresca);

// Compute ultimate FOS
SafetyFactorResult ultResult = sf.computeUltimateFOS(stress, mesh);

// Access results
Scalar minFOS = fosResult.minFOS;              // Minimum FOS in structure
Index critNode = fosResult.criticalNode;       // Node with minimum FOS
Vector3 critLoc = fosResult.criticalLocation;  // Position of critical node
VectorX allFOS = fosResult.fosPerNode;         // FOS at every node
bool passed = fosResult.passed;                // True if minFOS >= required

// Set required FOS for pass/fail check
sf.setRequiredFOS(1.5);  // Require 50% safety margin

VTK Export

Export results for visualization in a VTK viewer:

#include "fem/io/VTKExporter.hpp"

VTKExporter vtk;

// Set deformation scale for visualization
vtk.setDeformationScale(100.0);  // Amplify deformation 100x

// Export mesh only
vtk.exportMesh("mesh.vtk", mesh);

// Export with displacement
vtk.exportDisplacement("displaced.vtk", mesh, displacement);

// Export with full stress results
vtk.exportResults("stress.vtk", mesh, displacement, stress);

// Export with safety factors
vtk.exportSafetyFactors("fos.vtk", mesh, displacement, fosResult);

// Export everything
vtk.exportAll("complete.vtk", mesh, displacement, stress, fosResult);

VTK output includes:

Deformed mesh geometry
Displacement vectors and magnitude
von Mises and Tresca stress
Principal stresses (σ₁, σ₂, σ₃)
Full stress tensor
Factor of safety per node

CSV/JSON Export

Export for reporting and data analysis:

#include "fem/io/ResultsExporter.hpp"

ResultsExporter exporter;
exporter.setPrecision(8);  // Decimal places

// CSV export (spreadsheet-compatible)
exporter.exportCSV("results.csv", mesh, displacement, stress, fosResult);

// JSON export (programmatic access)
exporter.exportJSON("results.json", mesh, displacement, stress, fosResult);

// Generate text summary
std::string summary = exporter.generateSummary(mesh, stress, fosResult);
std::cout << summary;

Example summary output:

=== FEM Analysis Summary ===

Mesh Statistics:
  Number of Nodes: 1247
  Number of Elements: 4823
  Spatial Dimension: 3
  Bounding Box:
    Min: (0, 0, 0)
    Max: (100, 50, 25)

Safety Factor Results:
  Minimum FOS: 1.42
  Critical Node: 847
  Critical Location: (45.2, 12.8, 24.9)
  Status: PASSED

Stress Results:
  Max von Mises Stress: 175.6 MPa
  Max Tresca Stress: 182.3 MPa
  Max Principal Stress (σ1): 168.4 MPa
  Min Principal Stress (σ3): -42.1 MPa

Architecture

src/libs/libfem/fem/
├── core/           # Core data structures
│   ├── Mesh.hpp/cpp
│   ├── DOFHandler.hpp/cpp
│   ├── Quadrature.hpp/cpp
│   └── Types.hpp
├── elements/       # Element implementations
│   ├── Element.hpp/cpp
│   ├── Bar2.hpp/cpp, Bar3.hpp/cpp
│   ├── Tri3.hpp/cpp, Tri6.hpp/cpp
│   ├── Quad4.hpp/cpp, Quad8.hpp/cpp, Quad9.hpp/cpp
│   ├── Tet4.hpp/cpp, Tet10.hpp/cpp
│   └── Hex8.hpp/cpp, Hex20.hpp/cpp, Hex27.hpp/cpp
├── integrators/    # Physics modules
│   ├── Integrator.hpp
│   ├── DiffusionIntegrator.hpp/cpp
│   ├── MassIntegrator.hpp/cpp
│   └── ElasticityIntegrator.hpp/cpp
├── nonlinear/      # Nonlinear analysis
│   ├── StrainMeasures.hpp/cpp
│   ├── NonlinearMaterial.hpp/cpp
│   ├── NeoHookean.hpp/cpp
│   ├── J2Plasticity.hpp/cpp
│   ├── NewtonRaphson.hpp/cpp
│   └── ArcLength.hpp/cpp
├── postprocess/    # Post-processing
│   ├── StressRecovery.hpp/cpp
│   ├── FailureCriteria.hpp/cpp
│   └── SafetyFactor.hpp/cpp
├── io/             # Import/Export
│   ├── STLImporter.hpp/cpp
│   ├── VTKExporter.hpp/cpp
│   └── ResultsExporter.hpp/cpp
├── model/          # High-level pipeline
│   └── FEModel.hpp/cpp
└── solvers/        # Solver wrappers
    ├── LinearSolver.hpp/cpp
    ├── IterativeSolver.hpp/cpp
    ├── EigenSolver.hpp/cpp
    ├── ModalAnalysis.hpp/cpp
    ├── RayleighDamping.hpp/cpp
    ├── NewmarkSolver.hpp/cpp
    ├── HHTSolver.hpp/cpp
    └── ExplicitSolver.hpp/cpp

References

Eigen Documentation
Spectra Documentation
Hughes, T.J.R. "The Finite Element Method: Linear Static and Dynamic Finite Element Analysis"
Belytschko, Liu, Moran. "Nonlinear Finite Elements for Continua and Structures"