Notebook

MFEM Example 1¶

Adapted from PyMFEM/ex1.py. Compare with the original Example 1 in MFEM.

This example code demonstrates the use of MFEM to define a simple finite element discretization of the Laplace problem

$\begin{equation*} -\Delta x = 1 \label{laplace}\tag{1} \end{equation*}$

in a domain $\Omega$ with homogeneous Dirichlet boundary conditions

$\begin{equation*} x = 0 \label{laplace2}\tag{2} \end{equation*}$

on the boundary $\partial \Omega$ .

The problme is discretized on a computational mesh in either 2D or 3D using a finite elements space of the specified order (2 by default) resulting in the global sparse linear system

$\begin{equation*} A X = B. \label{laplace3}\tag{3} \end{equation*}$

The example highlights the use of mesh refinement, finite element grid functions, as well as linear and bilinear forms corresponding to the left-hand side and right-hand side of the discrete linear system. We also cover the explicit elimination of essential boundary conditions and using the GLVis tool for visualization.

In [ ]:

# Requires PyMFEM, see https://github.com/mfem/PyMFEM
import mfem.ser as mfem
from glvis import glvis, to_stream

In [ ]:

# Load the mesh from a local file
# meshfile = '../../mfem/data/star.mesh'
# mesh = mfem.Mesh(meshfile)

# Alternatively, create a simple square mesh and refine it
mesh = mfem.Mesh(5, 5, "TRIANGLE")
mesh.UniformRefinement()

# Create H1 finite element function space
fec = mfem.H1_FECollection(2, mesh.Dimension()) # order=2
fespace = mfem.FiniteElementSpace(mesh, fec)      

# Determine essential degrees of freedom (the whole boundary here)
ess_tdof_list = mfem.intArray()
ess_bdr = mfem.intArray([1]*mesh.bdr_attributes.Size())
fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list)

# Define Bilinear and Linear forms for the Laplace problem -Δu=1
one = mfem.ConstantCoefficient(1.0)
a = mfem.BilinearForm(fespace)
a.AddDomainIntegrator(mfem.DiffusionIntegrator(one))
a.Assemble()
b = mfem.LinearForm(fespace)
b.AddDomainIntegrator(mfem.DomainLFIntegrator(one))
b.Assemble()

# Create a grid function for the solution and initialize with 0
x = mfem.GridFunction(fespace);
x.Assign(0.0)

# Form the linear system, AX=B, for the FEM discretization
A = mfem.OperatorPtr()
B = mfem.Vector()
X = mfem.Vector()
a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
print("Size of the linear system: " + str(A.Height()))

# Solve the system using PCG solver and get the solution in x
Asm = mfem.OperatorHandle2SparseMatrix(A)
Msm = mfem.GSSmoother(Asm)
mfem.PCG(Asm, Msm, B, X, 1, 200, 1e-12, 0.0)
a.RecoverFEMSolution(X, b, x)

Plot the Solution with GLVis¶

In [ ]:

# Plot the mesh + solution (all GLVis keys and mouse commands work)
glvis((mesh, x), 400, 400)

In [ ]:

# Plot the mesh only
glvis(mesh)

In [ ]:

# Visualization with additional GLVis keys
g = glvis(to_stream(mesh,x) + 'keys ARjlmcbp*******')
g.set_size(600, 400)
g