Linear Elasticity

Solve structural mechanics problems with finite elements

Run AllReset

This tutorial follows MFEM's Example 2 (ex2.cpp), which demonstrates solving a linear elasticity problem for a multi-material cantilever beam. We cover vector finite elements, the elasticity bilinear form, piecewise material properties, and both Dirichlet and Neumann boundary conditions.

Problem: Solve the equilibrium equation $-\nabla \cdot \sigma(\mathbf{u}) = 0$ for a cantilever beam with fixed left edge and downward traction on the right edge.

Governing Equations

The equilibrium equation for a solid body under load is:

$-\nabla \cdot \sigma(\mathbf{u}) = \mathbf{f}$

where $\mathbf{u}$ is the displacement field, $\sigma$ is the Cauchy stress tensor, and $\mathbf{f}$ is the body force per unit volume (zero in this example).

Constitutive Law (Hooke's Law)

For linear elastic materials, stress is related to strain through Hooke's law using the Lamé parameters $\lambda$ and $\mu$ :

$\sigma(\mathbf{u}) = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + \mu (\nabla \mathbf{u} + \nabla \mathbf{u}^T)$

The strain tensor $\varepsilon$ is the symmetric part of the displacement gradient:

$\varepsilon(\mathbf{u}) = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$

The Lamé parameters relate to Young's modulus $E$ and Poisson's ratio $\nu$ by:

$\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)} \qquad \mu = \frac{E}{2(1+\nu)}$

Cantilever Beam Problem

The classic benchmark problem is a 3D two-material cantilever beam: fixed on the left end (Dirichlet BC) with a downward force on the right end (Neumann BC):

                     +----------+----------+
        boundary --->| material | material |<--- boundary
        attribute 1  |    1     |    2     |     attribute 2
        (fixed)      +----------+----------+     (pull down)

Boundary conditions:

Left face (attribute 1): Fixed support, $\mathbf{u} = \mathbf{0}$
Right face (attribute 2): Applied traction, $\sigma \cdot \mathbf{n} = [0, 0, -f]^T$
Other faces (attribute 3): Traction-free, $\sigma \cdot \mathbf{n} = \mathbf{0}$

Material properties:

Material 1 (left half): Stiff, $\lambda = 50$ , $\mu = 50$
Material 2 (right half): Soft, $\lambda = 1$ , $\mu = 1$

Weak Formulation

Multiplying by a test function $\mathbf{v}$ and integrating by parts yields:

$\int_\Omega \sigma(\mathbf{u}) : \varepsilon(\mathbf{v}) \, d\mathbf{x} = \int_{\Gamma_N} \mathbf{t} \cdot \mathbf{v} \, ds$

This leads to the bilinear and linear forms:

Form	Expression
Bilinear	$a(\mathbf{u}, \mathbf{v}) = \int_\Omega \sigma(\mathbf{u}) : \varepsilon(\mathbf{v}) \, d\mathbf{x}$
Linear	$L(\mathbf{v}) = \int_{\Gamma_N} \mathbf{t} \cdot \mathbf{v} \, ds$

Step 1-3: Load the Mesh

We load a 3D beam mesh with 8 hexahedral elements. The mesh has two material regions (element attributes 1 and 2) and three boundary attributes (1=left fixed, 2=right load, 3=other faces free).

// C++ equivalent:
Mesh mesh(mesh_file, 1, 1);
int dim = mesh.Dimension();

Code [1]Run

# MFEM Example 2: Load beam mesh
# In C++: Mesh mesh(mesh_file, 1, 1);

# Load the beam mesh from MFEM mesh format
mesh = mfem.mesh.Mesh.fromString(`${BEAM_MESH}`);

# Get mesh information
dim = mesh.dimension;
num_elements = mesh.numElements;
num_vertices = mesh.numVertices;
num_boundary = mesh.numBoundaryElements;

disp("=== Beam Mesh Loaded ===")
disp(strcat("Dimension: ", num2str(dim)))
disp(strcat("Number of elements: ", num2str(num_elements)))
disp(strcat("Number of vertices: ", num2str(num_vertices)))
disp(strcat("Number of boundary elements: ", num2str(num_boundary)))

Step 4: Mesh Refinement

Refine the mesh to increase resolution. We choose the number of refinement levels to get a final mesh with no more than 5,000 elements.

// C++ equivalent:
int ref_levels = (int)floor(log(5000./mesh.GetNE())/log(2.)/dim);
for (int l = 0; l < ref_levels; l++)
{
   mesh.UniformRefinement();
}

Code [2]Run

# MFEM ex2 Step 4: Refine mesh to increase resolution
# Target: no more than 5,000 elements

dim = mesh.dimension;
num_elements = mesh.numElements;

# Calculate refinement levels (from ex2.cpp)
max_elements = 5000;
ref_levels = floor(log(max_elements / num_elements) / log(2) / dim);

disp(strcat("Initial elements: ", num2str(num_elements)))
disp(strcat("Refinement levels needed: ", num2str(ref_levels)))

# Perform uniform refinement
for l = 1:ref_levels
    mesh.refineUniform()
    disp(strcat("After level ", num2str(l), ": ", num2str(mesh.numElements), " elements"))
end

disp("")
disp(strcat("Final mesh: ", num2str(mesh.numElements), " elements"))

Mesh Visualization

Visualize the refined beam mesh to understand its structure. The mesh shows the cantilever beam geometry with quadrilateral elements.

Code [3]Run

# Visualize the refined beam mesh using VTK

disp(strcat("Visualizing mesh: ", num2str(mesh.numElements), " elements, ", num2str(mesh.numVertices), " vertices"))

# Visualize with VTK
vtk.colormap("viridis")
vtk.edges("on")
vtk.mesh(mesh)
title("Cantilever Beam Mesh (MFEM ex2)")

Step 5: Define the Vector Finite Element Space

Unlike scalar PDEs (like Laplace), elasticity involves a vector unknown: $\mathbf{u} = [u_x, u_y, u_z]^T$ in 3D. We create an H1 space with vdim=3 (vector dimension = 3).

// C++ equivalent:
FiniteElementCollection *fec = new H1_FECollection(order, dim);
FiniteElementSpace fespace(&mesh, fec, dim);  // dim copies of scalar space

Each node has $d$ degrees of freedom (DOFs). For a 3D mesh with $N$ nodes, the global system has size $3N \times 3N$ .

Code [4]Run

# MFEM ex2 Step 5: Define Vector Finite Element Space
# Create H1 space with vdim=3 for 3D displacement field

order = 1;  # Polynomial order (linear elements)
vdim = 3;   # Vector dimension (displacement has 3 components in 3D)

# Create H1 finite element space for vector field
fespace = mfem.fespace.createH1(mesh, order, vdim);

# Get the number of degrees of freedom
ndofs = fespace.ndofs;

disp("=== Vector Finite Element Space ===")
disp(strcat("Polynomial order: ", num2str(order)))
disp(strcat("Vector dimension: ", num2str(vdim)))
disp(strcat("Number of DOFs: ", num2str(ndofs)))
disp(strcat("DOFs per node: ", num2str(vdim)))

Step 6: Set Up Material Coefficients

Define piecewise constant Lamé parameters for the two materials:

Material 1 (element attribute 1): $\lambda = 50$ , $\mu = 50$ (stiff)
Material 2 (element attribute 2): $\lambda = 1$ , $\mu = 1$ (soft)

// C++ equivalent:
Vector lambda(mesh.attributes.Max());
lambda = 1.0;
lambda(0) = lambda(1)*50;  // Material 1 is 50x stiffer
PWConstCoefficient lambda_func(lambda);

Code [5]Run

# MFEM ex2 Step 6: Set up piecewise constant material coefficients
# Two materials with different stiffness

# Material 1: stiff (left half of beam)
# Material 2: soft (right half of beam)
lambda_vals = [50.0, 1.0];  # Lame first parameter
mu_vals = [50.0, 1.0];      # Lame second parameter (shear modulus)

# Create piecewise constant coefficients
lambda = mfem.coefficients.PWConstCoefficient.create(lambda_vals);
mu = mfem.coefficients.PWConstCoefficient.create(mu_vals);

disp("=== Material Properties ===")
disp("Material 1 (stiff): lambda = 50, mu = 50")
disp("Material 2 (soft):  lambda = 1,  mu = 1")
disp("Stiffness ratio: 50:1")

Step 7: Determine Essential Boundary DOFs

Mark the left edge (boundary attribute 1) as essential (Dirichlet). This enforces $\mathbf{u} = \mathbf{0}$ on the fixed support.

// C++ equivalent:
Array<int> ess_tdof_list, ess_bdr(mesh.bdr_attributes.Max());
ess_bdr = 0;
ess_bdr[0] = 1;  // Mark boundary attribute 1
fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);

Code [6]Run

# MFEM ex2 Step 7: Determine essential boundary DOFs
# Mark left edge (boundary attribute 1) as fixed

# Boundary attribute 1 = left edge (fixed support)
ess_bdr_attrs = [1];

# Get the essential true DOFs
ess_tdof_list = mfem.boundaryConditions.getEssentialTrueDofs(fespace, ess_bdr_attrs);

disp("=== Essential Boundary Conditions ===")
disp("Fixed support on left edge (boundary attribute 1)")
disp(strcat("Number of essential DOFs: ", num2str(length(ess_tdof_list))))

Step 8: Set Up the Linear Form (RHS)

The linear form represents the applied traction on the right face (boundary attribute 2). We apply a downward pull force $\mathbf{t} = [0, 0, -0.01]^T$ .

// C++ equivalent:
VectorArrayCoefficient f(dim);
f.Set(dim-1, new PWConstCoefficient(pull_force));  // Pull down in z-direction

LinearForm *b = new LinearForm(fespace);
b->AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(f));
b->Assemble();

Code [7]Run

# MFEM ex2 Step 8: Set up the linear form (RHS)
# Traction force on right face (boundary attribute 2)

# Define traction vector: pull down in z-direction
traction = mfem.coefficients.vectorConstant([0.0, 0.0, -0.01]);

# Create linear form
b = mfem.linearForms.create(fespace);

# Add boundary integrator for traction on right face
# Boundary attribute 2 = right face (loaded)
b.addVectorBoundaryIntegrator(traction, [2]);

# Assemble the linear form
b.assemble();

disp("=== Linear Form (RHS) ===")
disp("Traction force: [0, 0, -0.01] on right face")
disp(strcat("Linear form size: ", num2str(b.size)))

Step 9: Define the Solution Vector

Create the solution grid function for the displacement field $\mathbf{u}$ . Initialize to zero, which satisfies the essential boundary conditions.

// C++ equivalent:
GridFunction x(&fespace);
x = 0.0;

Code [8]Run

# MFEM ex2 Step 9: Create solution grid function
# Initialize displacement field to zero

# Create grid function on vector finite element space
u = mfem.gridfunction.GridFunction.create(fespace);

# Initialize to zero (satisfies homogeneous Dirichlet BC)
u.projectConstant(0.0);

disp("=== Solution GridFunction ===")
disp("Displacement field u = [ux, uy, uz]")
disp("Initial guess: u = 0")
disp(strcat("GridFunction size: ", num2str(u.size)))

Step 10: Set Up the Bilinear Form (LHS)

The bilinear form uses the ElasticityIntegrator to compute $a(\mathbf{u}, \mathbf{v}) = \int_\Omega \sigma(\mathbf{u}) : \varepsilon(\mathbf{v}) \, d\mathbf{x}$ .

// C++ equivalent:
BilinearForm *a = new BilinearForm(fespace);
a->AddDomainIntegrator(new ElasticityIntegrator(lambda_func, mu_func));
a->Assemble();

Code [9]Run

# MFEM ex2 Step 10: Set up the bilinear form (LHS)
# Elasticity integrator: integral(sigma(u) : epsilon(v)) dx

# Create bilinear form on the vector finite element space
a = mfem.bilinearForms.create(fespace);

# Add elasticity integrator with piecewise Lame parameters
a.addElasticityIntegrator(lambda, mu);

# Assemble the bilinear form
a.assemble();

disp("=== Bilinear Form (LHS) ===")
disp("Elasticity integrator: a(u,v) = integral(sigma(u) : epsilon(v)) dx")
disp("Bilinear form assembled")

Step 11: Apply Essential Boundary Conditions

Eliminate essential boundary conditions directly on the bilinear form and linear form. This modifies the system matrix rows/columns for essential DOFs (setting diagonal to 1, off-diagonal to 0) and adjusts the RHS accordingly.

// C++ equivalent:
a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);

Note: We use applyEssentialBC which modifies the bilinear form and RHS in-place, then finalize to get the sparse matrix ready for solving.

Code [10]Run

# MFEM ex2 Step 11: Apply essential boundary conditions
# Modifies bilinear form and RHS in-place

disp(strcat("Total DOFs: ", num2str(fespace.trueVSize)))
disp(strcat("Essential DOFs count: ", num2str(length(ess_tdof_list))))

# Apply essential BCs: sets essential rows/cols in A, adjusts b
mfem.boundaryConditions.applyEssentialBC(a, ess_tdof_list, u, b);

# Finalize the bilinear form (convert to sparse matrix)
a.finalize();

disp("=== Essential BCs Applied ===")
disp("System matrix and RHS modified for Dirichlet conditions")

Step 12: Solve the Linear System

Solve the system $A\mathbf{u} = \mathbf{b}$ using Preconditioned Conjugate Gradient (PCG) with a Gauss-Seidel smoother. The solution is computed directly into the grid function $\mathbf{u}$ .

// C++ equivalent:
GSSmoother M((SparseMatrix&)(*A));
PCG(*A, M, B, X, 1, 500, 1e-8, 0.0);

Code [11]Run

# MFEM ex2 Step 12: Solve the system
# Uses PCG with Gauss-Seidel preconditioner

disp("=== Solving Au = b with PCG ===")

# Create CG solver
solver = mfem.linear.CGSolver.create();

# Configure solver
solver.setRelTol(1e-8);
solver.setAbsTol(0.0);
solver.setMaxIter(500);
solver.setPrintLevel(1);

# Set the system matrix (from the bilinear form)
solver.setOperator(a.getSystemMatrixPointer());

# Solve: A * u = b
disp("Running PCG solver...")
solver.mult(b.getVectorPointer(), u.getPointer());

# Report solver results
num_iters = solver.getNumIterations();
converged = solver.getConverged();
final_norm = solver.getFinalNorm();
disp(strcat("Iterations: ", num2str(num_iters)))
disp(strcat("Converged: ", num2str(converged)))
disp(strcat("Final norm: ", num2str(final_norm)))

Step 13: Visualize the Displacement Field

Visualize the computed displacement field. The beam bends downward due to the applied traction, with the stiff material (left) deforming less than the soft material (right).

Code [12]Run

# MFEM ex2 Step 13: Visualize the displacement field

disp("=== Visualization ===")

# Get displacement data and statistics
disp_arr = u.getData();
disp(strcat("Number of DOFs: ", num2str(length(disp_arr))))
disp(strcat("Number of vertices: ", num2str(mesh.numVertices)))

min_val = min(disp_arr);
max_val = max(disp_arr);
disp(strcat("Displacement range: [", num2str(min_val), ", ", num2str(max_val), "]"))

# Visualize displacement field with VTK
# Use vtk.displacement with scale factor to amplify small displacements
vtk.colormap("coolwarm")
vtk.edges("on")
vtk.displacement(mesh, u, 100)  # Scale deformation 100x for visibility
title("Cantilever Beam Displacement (MFEM ex2)")

Summary

We have completed MFEM Example 2, solving a linear elasticity problem for a multi-material cantilever beam.

Steps Covered:

Load Mesh - 3D beam mesh with 8 hex elements, 2 materials
Refine Mesh - Uniform refinement for better resolution
Visualize Mesh - VTK rendering of the 3D beam mesh
Create Vector FE Space - H1 space with vdim=3 for 3D displacement
Set Up Materials - Piecewise constant Lamé parameters
Mark Essential BCs - Fixed left face (Dirichlet)
Create Linear Form - Traction force on right face (Neumann)
Create Solution - Initialize displacement to zero
Create Bilinear Form - Elasticity integrator
Apply Essential BCs - Eliminate fixed DOFs
Solve - CG solver for the linear system
Visualize - Display 3D displacement field

Key Concepts:

Vector finite elements - Each node has 3 DOFs in 3D
Elasticity integrator - Computes stress-strain relationship
Piecewise coefficients - Different material properties per region
Mixed BCs - Dirichlet (fixed) + Neumann (traction) boundaries

Workbench

Clear

No variables in workbench

Next Steps

Finite Element Method

Learn FEM fundamentals: mesh generation, basis functions, assembly

Heat Conduction

Time-dependent problems with ODE solvers

Equana