MFEM Example 1: Laplace Problem

Solving PDEs with the Finite Element Method

Run AllReset

This tutorial follows MFEM's Example 1 (ex1.cpp), which demonstrates the basic usage of MFEM to solve a simple Laplace problem. We cover finite element grid functions, linear and bilinear forms, elimination of essential boundary conditions, and visualization.

Problem: Solve the Laplace equation $-\Delta u = 1$ with homogeneous Dirichlet boundary conditions on a star-shaped domain.

Step 1-3: Load the Mesh

MFEM can handle triangular, quadrilateral, tetrahedral, hexahedral, surface and volume meshes. Here we load the star mesh from MFEM's data directory. The mesh file format specifies dimension, elements, boundary markers, and vertex coordinates.

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

Code [1]Run

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

# Load the star mesh from MFEM mesh format
mesh = mfem.mesh.Mesh.fromString(`${STAR_MESH}`)

# Get mesh information (using mfemwasm property names)
dim = mesh.dimension
num_elements = mesh.numElements
num_vertices = mesh.numVertices
num_boundary = mesh.numBoundaryElements

disp("=== Star 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 10,000 elements. Each uniform refinement in 2D splits each element into 4.

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

Note: The mesh variable persists from the previous cell, so we can continue working with it.

Code [2]Run

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

# Use the mesh from the previous cell
dim = mesh.dimension
num_elements = mesh.numElements

# Calculate refinement levels (from ex1.cpp)
# ref_levels = floor(log(10000/NE)/log(2)/dim)
max_elements = 10000
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(strcat("Final mesh: ", num2str(mesh.numElements), " elements"))

Mesh Visualization

Visualize the refined star mesh to understand its structure. The vtk.mesh function takes an MFEM mesh object directly and renders it using VTK.js.

Since the mesh variable persists from the previous cell, we can visualize it directly.

Code [3]Run

# Visualize the refined star mesh using VTK
# The mesh was already refined in the previous cell

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("Refined Star Mesh (MFEM ex1)")

Step 5: Define the Finite Element Space

Define a finite element space on the mesh. For the Laplace problem, we use H1 (continuous) finite elements of a given polynomial order.

// C++ equivalent:
FiniteElementCollection *fec = new H1_FECollection(order, dim);
FiniteElementSpace fespace(&mesh, fec);
cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;

Code [4]Run

# MFEM ex1 Step 5: Define the Finite Element Space
# Create H1 (continuous) finite element space of order 1

order = 1  # Polynomial order

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

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

disp("=== Finite Element Space ===")
disp(strcat("Polynomial order: ", num2str(order)))
disp(strcat("Number of DOFs: ", num2str(ndofs)))

Step 6: Determine Essential Boundary DOFs

Determine the list of essential (Dirichlet) boundary degrees of freedom. In this example, we mark all external boundary attributes as essential, which corresponds to homogeneous Dirichlet boundary conditions ( $u = 0$ on the boundary).

// C++ equivalent:
Array<int> ess_tdof_list;
if (mesh.bdr_attributes.Size())
{
   Array<int> ess_bdr(mesh.bdr_attributes.Max());
   ess_bdr = 0;
   // Apply boundary conditions on all external boundaries:
   mesh.MarkExternalBoundaries(ess_bdr);
   fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
}

Code [5]Run

# MFEM ex1 Step 6: Determine essential boundary DOFs
# Mark all boundary attributes as essential (Dirichlet)

# Get all boundary attributes (in our star mesh, all boundaries have attribute 1)
bdr_attrs = [1]  # Array of boundary attributes to mark as essential

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

disp("=== Essential Boundary Conditions ===")
disp(strcat("Boundary attributes marked: ", num2str(bdr_attrs)))
disp(strcat("Number of essential DOFs: ", num2str(length(ess_tdof_list))))

Step 7: Set Up the Linear Form (RHS)

Set up the linear form $b(v)$ which corresponds to the right-hand side of the FEM linear system. In this case, $b(v) = (1, \phi_i)$ where $\phi_i$ are the basis functions in the finite element space.

// C++ equivalent:
LinearForm b(&fespace);
ConstantCoefficient one(1.0);
b.AddDomainIntegrator(new DomainLFIntegrator(one));
b.Assemble();

Code [6]Run

# MFEM ex1 Step 7: Set up the linear form (RHS)
# b(v) = integral(f * v) where f = 1

# Create constant coefficient f = 1
one = mfem.coefficients.constant(1.0)

# Create linear form on the finite element space
b = mfem.linearForms.create(fespace)

# Add domain integrator: integral(f * v) dx
b.addDomainIntegrator(one)

# Assemble the linear form
b.assemble()

disp("=== Linear Form (RHS) ===")
disp("Source term f = 1 (constant)")
disp(strcat("Linear form size: ", num2str(b.size)))

Step 8: Define the Solution Vector

Define the solution vector $x$ as a finite element grid function corresponding to the finite element space. Initialize $x$ with an initial guess of zero, which satisfies the boundary conditions.

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

Code [7]Run

# MFEM ex1 Step 8: Define the solution grid function
# Initialize with zero (satisfies homogeneous Dirichlet BC)

# Create grid function on the finite element space
x = mfem.gridfunction.GridFunction.create(fespace)

# Initialize to zero
x.projectConstant(0.0)

disp("=== Solution GridFunction ===")
disp("Initial guess: x = 0")
disp(strcat("GridFunction size: ", num2str(x.size)))

Step 9: Set Up the Bilinear Form (LHS)

Set up the bilinear form $a(u,v)$ on the finite element space corresponding to the Laplacian operator $-\Delta$ , by adding the Diffusion domain integrator.

The weak form of the Laplace equation $-\Delta u = f$ is: $a(u,v) = \int_\Omega \nabla u \cdot \nabla v \, dx = \int_\Omega f v \, dx = b(v)$

// C++ equivalent:
BilinearForm a(&fespace);
a.AddDomainIntegrator(new DiffusionIntegrator(one));
a.Assemble();

Code [8]Run

# MFEM ex1 Step 9: Set up the bilinear form (LHS)
# a(u,v) = integral(grad(u) . grad(v)) dx

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

# Add diffusion integrator: integral(k * grad(u) . grad(v)) with k = 1
a.addDiffusionIntegrator(one)

# Assemble the bilinear form
a.assemble()

disp("=== Bilinear Form (LHS) ===")
disp("Diffusion operator: -Delta u")
disp("Bilinear form assembled")

Step 10: Form Linear System and Apply Boundary Conditions

Apply the essential boundary conditions to the system. This modifies both the matrix $A$ and the right-hand side vector $B$ to enforce the Dirichlet boundary conditions ( $u = 0$ on the boundary).

// C++ equivalent:
OperatorPtr A;
Vector B, X;
a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
cout << "Size of linear system: " << A->Height() << endl;

Code [9]Run

# MFEM ex1 Step 10: Apply essential BCs to the system
# This modifies the bilinear form and linear form for Dirichlet BCs

# Apply essential boundary conditions
# This eliminates BC DOFs from the system
mfem.boundaryConditions.applyEssentialBC(a, ess_tdof_list, x, b)

disp("=== Linear System Formed ===")
disp("Essential BCs applied to system")
disp(strcat("System size: ", num2str(ndofs), " DOFs"))
disp(strcat("Essential (fixed) DOFs: ", num2str(length(ess_tdof_list))))

Step 11: Solve the Linear System

Solve the linear system $Ax = B$ using a Conjugate Gradient (CG) solver. CG is appropriate here because the diffusion operator produces a symmetric positive definite matrix.

// C++ equivalent:
GSSmoother M((SparseMatrix&)(*A));
PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
a.RecoverFEMSolution(X, b, x);

Note: The solution is recovered back into the GridFunction x, which can then be visualized or post-processed.

Code [10]Run

# MFEM ex1 Step 11: Solve the linear system
# Use CG (Conjugate Gradient) solver for symmetric positive definite system

# Finalize the bilinear form (prepares matrix for solving)
a.finalize()

# Get the system matrix pointer for solving
A_ptr = a.getSystemMatrixPointer()

# Get RHS vector from linear form
b_vec = b.getVectorPointer()

# Get solution vector pointer from grid function
x_ptr = x.getPointer()

disp("=== Solving Ax = b ===")

# Create CG solver (from the linear solvers namespace)
solver = mfem.linear.CGSolver.create()

# Configure solver tolerances and iterations
solver.setRelTol(1e-10)
solver.setAbsTol(0.0)
solver.setMaxIter(200)
solver.setPrintLevel(1)

# Set the system matrix as the operator
solver.setOperator(A_ptr)

# Solve: find x such that A*x = b
disp("Running CG solver...")
solver.mult(b_vec, x_ptr)

# 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 12: Visualize the Solution

Visualize the computed solution using VTK. The solution shows the temperature/potential distribution on the star-shaped domain, with zero values on the boundary (Dirichlet BC) and maximum values in the interior where the source term $f=1$ pushes the solution up.

Code [11]Run

# MFEM ex1 Step 12: Visualize the solution with VTK

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

# Get solution data and compute min/max
sol_arr = reshape(x.getData(), x.size);
min_val = min(sol_arr);
max_val = max(sol_arr);
disp(strcat("Solution range: [", num2str(min_val), ", ", num2str(max_val), "]"));

# Visualize the solution using VTK
vtk.colormap("viridis")
vtk.gridfunction(mesh, x)
title("Laplace Solution on Star Domain")

Summary

We have completed MFEM Example 1, which solves the Laplace equation $-\Delta u = 1$ with homogeneous Dirichlet boundary conditions on a star-shaped domain.

Steps Covered:

Load Mesh - Read the star mesh from MFEM format
Refine Mesh - Uniform refinement for better resolution
Visualize Mesh - VTK rendering of the refined mesh
Create FE Space - H1 continuous finite elements (order 1)
Mark Boundary DOFs - Identify essential (Dirichlet) boundary
Set Up Linear Form - RHS: $b(v) = \int f v \, dx$ with $f=1$
Create Solution - Initialize GridFunction to zero
Set Up Bilinear Form - LHS: $a(u,v) = \int \nabla u \cdot \nabla v \, dx$
Apply BCs - Eliminate essential DOFs from system
Solve - CG solver for $Ax = b$
Visualize - Display solution with VTK

Workbench

Clear

No variables in workbench

Next Steps

VTK.js Visualization

Advanced 3D visualization for volumes and FEM meshes

Working with Matrices

Matrix operations essential for FEM assembly and solving

Equana