Course: Materials Simulation Practical | FAU Erlangen-Nürnberg
Tools: Python · NumPy · Matplotlib
From-scratch implementation of the Finite Element Method for 1D linear elasticity in Python/NumPy. The project covers the full pipeline: constitutive tensor symmetry reduction, weak-form derivation, linear shape functions and B-matrix, Newton–Cotes numerical quadrature, global stiffness assembly, Dirichlet/Neumann boundary condition enforcement via permutation-matrix reduction, and element-wise stress recovery across four loading cases.
Starting from the static balance of linear momentum div(σ) = −ρb, Hooke's law σ = C:ε is applied. The full symmetry analysis of the stiffness tensor C_ijkl reduces 81 independent components to 21 via minor symmetries (stress/strain symmetry) and major symmetry (strain energy density), cast in Voigt notation. The 1D weak form is derived by multiplying by a test function and integrating by parts:
Linear Lagrange shape functions (N_1) and (N_2) are defined over each element. Their spatial derivatives form the strain-displacement matrix (B), which is constant for linear 1D elements. As a result, the strain and stress recovered from the displacement field are elementwise constant.
 |
| Shape functions N₁ and N₂ |
Interpolation accuracy is tested for three functions: a linear function (f(x)=2x+3), a sinusoidal function (f(x)=\sin(x)), and a quadratic function (f(x)=x^2).
 |
 |
 |
| Linear function | Sinusoidal function | Quadratic function |
The linear function is reproduced exactly, while the sinusoidal and quadratic functions are only approximated, producing a bounded interpolation error.
The same three functions are then examined at the derivative level by comparing the exact gradient with the gradient reconstructed from the shape-function derivatives (B = dN/dx).
 |
 |
 |
| Derivative of linear function | Derivative of sinusoidal function | Derivative of quadratic function |
Because linear shape functions have constant derivatives within each element, the gradient approximation is piecewise constant. This is why the stress field recovered from linear bar elements is also piecewise constant, with jumps possible across element boundaries.
Three Newton–Cotes rules are implemented and validated:
| Rule | Points | Exact for polynomial degree |
|---|---|---|
| Trapezoidal | 2 | ≤ 1 |
| Simpson's | 3 | ≤ 2 |
| 3/8 Rule | 4 | ≤ 3 |
The element stiffness matrix for a 1D bar (E = 210,000 N/mm², A = 25 mm², L = 50 mm) is computed via trapezoidal quadrature and verified analytically.
An assembly routine maps local element stiffness matrices into the global system via node connectivity, producing the symmetric tridiagonal global stiffness matrix K. Dirichlet boundary conditions are enforced using a permutation-matrix reduction: a matrix P of shape (n−m)×n maps the full displacement vector to the free-DOF subspace, solves the reduced system K_r·d_r = f_r, and reconstructs the full solution via back-substitution.
Element-wise stress is recovered as σ = E·B·û on a 6-element, 7-node bar with total length L = 50 mm (E = 210,000 N/mm², A = 25 mm²). Nodes are indexed 0–6.
 |
 |
| Displacement | Stress |
The displacement varies linearly and the stress remains uniform at σ = 0.2 N/mm² across all elements, consistent with the analytical solution for an axially loaded bar with one fixed end and an applied end force.
 |
 |
| Displacement | Stress |
Displacement ramps linearly to the prescribed value at x = 25 mm and plateaus thereafter. Stress is elevated in the constrained region and drops sharply in the free segment.
 |
 |
| Displacement | Stress |
A concentrated internal load at node 3 produces a stress jump at that location. Elements to the left carry the combined load while elements to the right carry only the end load, resulting in a piecewise-constant stress profile with a clear discontinuity.
 |
 |
| Displacement | Stress |
With both ends constrained, the displacement remains continuous but changes slope across the internal load location. The stress field is piecewise constant and exhibits a clear jump at the loaded node, increasing from approximately 4.7 N/mm² on the left to 7.9 N/mm² on the right.
pip install numpy matplotlibOpen fem_1d_elastic_bar.ipynb in Jupyter and run cells sequentially. All functions — shape functions, B-matrix, quadrature, assembly, BC enforcement, and stress recovery — are defined incrementally and reused across tasks.