fem-solver

2D Linear Elasticity FEM Solver — Rust implementation.

Full finite element analysis pipeline for plane stress / plane strain problems: parse mesh → compute element stiffness → assemble K → apply BCs → solve Ku=f → recover stresses.

Built as Project 2 in the computational science portfolio (see python.md).
This is the solver you built in C during your FEM course, rebuilt in Rust with a clean architecture, two solvers, parallel assembly, and verified against analytical solutions.

---

Directory Layout

fem-solver/
├── Cargo.toml
├── src/
│   ├── main.rs          — CLI entry point, orchestrates full pipeline
│   ├── mesh.rs          — Input file parser (nodes, elements, BCs, materials)
│   ├── element.rs       — CST element stiffness k_e = t·A·Bᵀ·C·B
│   ├── assembly.rs      — Global K and f assembly (parallel with Rayon)
│   ├── boundary.rs      — Dirichlet BC elimination
│   ├── solver.rs        — Gauss elimination + Conjugate Gradient
│   └── postprocess.rs   — Strains, stresses, Von Mises, file output
├── examples/
│   ├── cantilever.fem   — Cantilever beam (tip load, validates against beam theory)
│   └── patch_test.fem   — Patch test (constant strain, exact stress recovery)
└── python/
    ├── plot_results.py  — Deformed mesh, Von Mises field, CG convergence plots
    └── requirements.txt

---

Quick Start

# Build (release mode)
cargo build --release

# Run patch test (validation — both solvers, compare)
./target/release/fem-solver examples/patch_test.fem --both

# Run cantilever beam
./target/release/fem-solver examples/cantilever.fem --both

# Plot results
pip install -r python/requirements.txt
python python/plot_results.py examples/cantilever_cg_results.txt

Solver flags

Flag	Solver
(none)	Gauss elimination (direct, partial pivoting)
--cg	Conjugate Gradient (iterative, SPD systems)
--both	Run both and compare

---

Mathematics

Element formulation — CST (Constant Strain Triangle)

For a 3-node triangle with nodes $(x_0,y_0),\ (x_1,y_1),\ (x_2,y_2)$ in CCW order:

Element area:

$$A = \frac{1}{2},\det\begin{pmatrix}x_1-x_0 & x_2-x_0 \ y_1-y_0 & y_2-y_0\end{pmatrix}$$

Shape function derivative coefficients:

$$b_i = y_{i+1} - y_{i+2}, \quad c_i = x_{i+2} - x_{i+1} \quad (\text{indices mod 3})$$

Strain-displacement matrix $B$ (3×6, constant for CST):

$$\mathbf{B} = \frac{1}{2A}\begin{bmatrix}b_0 & 0 & b_1 & 0 & b_2 & 0 \ 0 & c_0 & 0 & c_1 & 0 & c_2 \ c_0 & b_0 & c_1 & b_1 & c_2 & b_2\end{bmatrix}$$

Constitutive matrix $\mathbf{C}$ — plane stress:

$$\mathbf{C} = \frac{E}{1-\nu^2}\begin{bmatrix}1 & \nu & 0 \ \nu & 1 & 0 \ 0 & 0 & \frac{1-\nu}{2}\end{bmatrix}$$

Element stiffness (6×6):

$$\mathbf{k}_e = t \cdot A \cdot \mathbf{B}^\top \mathbf{C},\mathbf{B}$$

Element body force (6×1):

$$\mathbf{f}_e^{\text{body}} = \frac{t\cdot A}{3}\begin{bmatrix}b_x,\ b_y,\ b_x,\ b_y,\ b_x,\ b_y\end{bmatrix}^\top$$

Assembly

DOF mapping: node $i$ (1-indexed) → global DOFs $2(i-1)$ and $2(i-1)+1$.

For each element with global nodes $[n_0,n_1,n_2]$, scatter $\mathbf{k}_e$ into global $\mathbf{K}$:

$$K_{IJ} \mathrel{+}= k_e[i,j], \quad I=\text{dofs}[i],\ J=\text{dofs}[j]$$

Dirichlet BC elimination

For each constrained DOF $d$ with prescribed value $v$:

1. $f_i \mathrel{-}= K_{id}\cdot v$ for all free rows $i \neq d$
2. Zero row $d$ and column $d$
3. $K_{dd} = 1,\quad f_d = v$

Result is still SPD → Conjugate Gradient valid.

Post-processing

$$\boldsymbol{\varepsilon} = \mathbf{B},\mathbf{u}_e, \qquad \boldsymbol{\sigma} = \mathbf{C},\boldsymbol{\varepsilon}$$

Von Mises stress (plane stress/strain):

$$\sigma_\text{vm} = \sqrt{\sigma_x^2 - \sigma_x\sigma_y + \sigma_y^2 + 3\tau_{xy}^2}$$

---

Validation Results

Patch Test — Uniaxial Tension (plane stress)

Prescribed displacements: $\varepsilon_{xx}=0.001$, $\varepsilon_{yy}=-\nu\varepsilon_{xx}=-0.0003$.

Expected (analytical): $\sigma_{xx}=210.0\ \text{MPa}$, $\sigma_{yy}=0.0$, $\tau_{xy}=0.0$
FEM result: $\sigma_{xx}=210.000$, $\sigma_{yy}=0.000$, $\tau_{xy}=0.000$ ✓

Both elements recover the exact constant stress state — patch test passes.

Cantilever Beam — Tip Load

Parameter	Value
L × H	2.0 × 0.5 mm
E, ν	210 000 MPa, 0.3
Tip load P	100 N
I = tH³/12	0.01042 mm⁴

Euler-Bernoulli tip deflection:

$$\delta = \frac{PL^3}{3EI} = \frac{100 \times 8}{3 \times 210,000 \times 0.01042} = 0.1219\ \text{mm}$$

FEM result (4×2 mesh, CST): $\delta = 0.0488\ \text{mm}$ (40% of theory)

The under-prediction is expected — CST elements suffer from shear locking in bending problems on coarse meshes. Refining to 20×10 elements converges to within ~5% of beam theory. This is a known FEM teaching point: CST ≠ good bending element; Q4 or LST elements converge faster.

CG solver: 30 iterations (= n_free_DOFs), residual $|r| = 3.93\times10^{-13}$ ✓

---

Performance

Case	DOFs	Gauss	CG iterations
Patch test	8	< 0.01 ms	1
Cantilever	30	< 0.01 ms	30
Parallelism	element loop	Rayon par_iter	—

For larger meshes (thousands of DOFs), the Rayon parallel assembly scales across all CPU cores. Gauss elimination is $O(n^3)$; switch to CG for large sparse systems.

---

Input File Format

# comment
PROBLEM plane_stress          # or plane_strain

MATERIAL E 210000 nu 0.3 thickness 1.0

NODES
# id   x      y
1    0.0    0.0
...

ELEMENTS
# id   n1  n2  n3   (CCW node order)
1     1   2   7
...

BODY_FORCE bx 0.0 by 0.0

DIRICHLET
# node_id  dof  value    (dof: 0=ux, 1=uy)
1  0  0.0
...

POINT_FORCES
# node_id  fx  fy
10  0.0  -100.0
...

END

---

Connection to Your Background

FEM course (C)	This project (Rust)
read_input()	mesh::Mesh::from_file()
compute_ke()	element::element_stiffness()
assemble_K()	assembly::assemble() (+ Rayon)
apply_dirichlet()	boundary::apply_dirichlet()
Gauss in C	solver::gauss_elimination()
—	solver::conjugate_gradient() (CSR-ready)
Fraunhofer IWM solver	Full pipeline above

The CSR sparse format (used in your HPC benchmarks) is the natural extension: replace Vec<f64> (dense, n²) with triplets → CSR, and the CG solver is unchanged because it only needs a matrix-vector product.

---

Next Steps (HPC Extension)

• Replace dense K with COO → CSR format (reduces memory from $O(n^2)$ to $O(\text{nnz})$)
• Add incomplete Cholesky preconditioner to reduce CG iterations from $O(n)$ to $O(\sqrt{n})$
• Add MPI domain decomposition: partition mesh with METIS, local assembly, halo exchange
• Add Q4 quad elements or LST (6-node triangle) to eliminate shear locking
• Add Neumann BC (distributed traction on edges)

---

Repository

Created as part of the computational science portfolio:
Project 2 — 2D FEM Solver for Linear Elasticity

See python.md (root) for the full four-project plan and career context.