MFEM Tutorial

Dohyun Kim

Appl. Math., Brown Univ.

APMA 2560
Feb 25, 2025

Table of Contents

  • Weak Formulation
  • Finite Element Formulation
  • Basic Structures for Linear Problem:
  • Technical Details
    • Mesh, Basis, Assembly, Quadrature rules,

Weak Formulation

Consider the Poisson equation with Dirichlet BC

The corresponding weak formulation is: Find such that

Here, .

Discrete Space

Let be a given mesh, and

We impose the boundary condition by setting

where is a projection onto and globally continuous if necessary.

w:640, center

Discrete Formulation

The corresponding discrete problem is: Find such that

Here, is a bilinear form, and is a linear form so that

  • for all and ,

  • for all and ,

  • for all and .

  • Note that for the most of conforming FEMs.

  • But, when you consider, e.g., Discontinuous Galerkin method.

Basis of

  • Let be a basis of . Then for given , there exists such that

  • Find or find such that

  • Since is finite dimensional and is a bilinear form, it is enough to test the equation against basis functions.

Linear System

  • Now we obtained a set of equations

  • This can be rewritten as

    where

  • Boundary Condition: are already known for where is the set of indices related to essential boundary.

  • Letting be the set of indices related to interior degrees of freedom, we obtain

    where and are submatrices of related to and , respectively.

Summary

  • For given finite element space with , we find solution by

  • Using the basis representation, we obtain the corresponding linear system

    where

Finite Element Space

Discrete Space

FiniteElementCollection

Mesh

FiniteElementSpace

GridFunction / Coefficient

Discrete Space

  • Recall that a finite space is defined by

  • We need three components to define
    • Local space:
      • Polynomial space , Raviart-Thomas space , etc.
    • Mesh:
      • Triangular mesh, Quadrilateral mesh, Tetahedral mesh, etc.
    • Global continuity or regularity
      • , , , etc.
  • The local space is FiniteElementCollection
  • The mesh is (Par)Mesh
  • Global continuity is encoded in FiniteElementCollection

FiniteElementCollection

  • We use FiniteElementCollection to represent a local finite element space.

  • Examples are: (: Complete polynomial, : Tensor product polynomial)

    • H1_FECollection - Continuous space
    • DG_FECollection/L2_FECollection - Discontinuous space
    • RT_FECollection - -conforming Raviart-Thomas space
    • ND_FECollection - -conforming Nedelec space
    • see, FiniteElementCollection for exhaustive list.

  • Code Example:

    FiniteElementCollection *fec = new DG_FECollection(order, dim);
FiniteElementCollection *fec = new H1_FECollection(order, dim);

Mesh

  • We use Mesh class to represent meshes.
  • A Mesh object can be created either from a file
    Mesh mesh("path/to/meshfile.mesh");
  • or
    Mesh mesh = Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL); // Uniform 2D rectangular mesh
  • You can also perform refinements or save the mesh
    Mesh mesh = Mesh::MakeCartesian2D(nx, ny, Element::TRIANGLE);
mesh.UniformRefinement();     // Uniform refine
mesh.RefineByError(err, tol); // Refine where err > tol
mesh.Save("path/to/meshfile.mesh");
  • It supports nonconforming and/or curved meshes.
  • The mesh file can be created 1) manually, 2) Gmsh, and so on. See, this guide.

Finite Element Space

  • With FiniteElementCollection and Mesh, we can define a finite element space .
  • We use a FiniteElementSpace class to store the information
  • For example, finite element space associated with a uniform rectangular mesh
    const int dim = 2;
const int order = 1;
const int nx = ny = 10;
FiniteElementCollection *fec = new H1_FECollection(order, dim);
Mesh mesh = Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL);
FiniteElementSpace fes(&mesh, fec); // P1-C0 Finite element space
  • A vector finite element space can be constructed either
    • with vector FiniteElementCollection, e.g., RT_FECollection
    • with tensor product space with a scalar FiniteElementCollection
    const int vdim = 2;
FiniteElementSpace vfes(&mesh, fec, vdim); // P1-C0 vector FE space

GridFunction

  • A discrete function can be created using
    GridFunction u(&fes);
  • For a discrete function , there exists a unique such that

  • A GridFunction stores the primal vector, [c_1, c_2, ..., c_N] associated with and (fes).
  • If we use nodal basis functions , then for each node .
  • It is important to distinguish "point value" and "coefficients" of a discrete function.
  • is well-defined within each element even if is not a node.

Summary

const int order = 2; // order of polynomial
const int nx = ny = 4; // number of elements in each axis

FiniteElementCollection *fec = new H1_FECollection(order, dim); // C0, Quadratic Polynomial
Mesh mesh = Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL); // Uniform quadrilateral mesh
const int dim  = mesh.Dimension(); // spatial dimension, 2
// You can also define a vector FE space
const int vdim = 2; // vector dimension
FiniteElementSpace fespace(&mesh, fec, vdim); // [C0-Q2 space]^2
GridFunction u(&fespace);

q2-fecmesh-and-fespace

BilinearForm and LinearForm

BilinearForm

LinearForm

NonlinearForm

mesh-and-fespace

mesh-and-fespace

BilinearForm

  • We want to express our bilinear form .
  • This operator is .
  • We use BilinearForm to represent such operator
    BilinearForm a_h(&fespace);
  • If , then
    MixedBilinearForm a_h(&fespace1, &fespace2);
  • We can add local bilinear forms to populate our bilinear form a_h.
    a_h.AddDomainIntegrator(new DiffusionIntegrator()); // a_h += (∇u, ∇v)
a_h.AddDomainIntegrator(new MassIntegrator());      // a_h += (u, v)
  • DiffusionIntegrator and MassIntegrator are BilinearFormIntegrator
  • This corresponds to a local bilinear form . See, fem/bilininteg.hpp

LinearForm

  • Similarly, we have a linear form .
  • This operator is .
  • We use LinearForm to represent such operator
    LinearForm b_h(&fespace);
  • We can add local linear forms to populate our linear form b_h.
    b_h.AddDomainIntegrator(new DomainLFIntegrator(f)); // b_h += (f,v)
b_h.AddBoundaryIntegrator(new BoundaryLFIntegrator(g_N));// b_h += (g_N, v)_e
  • DomainLFIntegrator and BoundaryLFIntegrator are LinearFormIntegrator
  • This corresponds to a local linear form . See, fem/lininteg.hpp

Assembly

  • With basis of , can be realized

  • For a given , we define

  • If the support of basis is local, .
  • Also, we have

    Similar condition holds for .
  • This implies sparsity of , and provides an assembly procedure.

Assembly

assembly

  • R_T, P_T: Mapping between the global and local
    • R_T*A = A(I_T, I_T) = A_loc
  • B, Q: Mapping between basis and values at points.
    • B can be function value, derivative value, etc.
  • D: Operation on quadrature values (dot products, weight, ...)
    • D(U, Ux, Uy, V, Vx, Vy) = Ux*Vx + Uy*Vy

Assembly - Pseudocode

pseudocode

Assembly - Pseudocode

pseudocode

How to write FormIntegrator

void MassIntegrator::AssembleElementMatrix(
    const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
{
   // INPUT el: Local finite element: basis, dof, ...
   // INPUT Trans: Local mesh information. e.g., element location, Jacobian
   // OUTPUT elmat: resulting (ϕ_j, ϕ_i) where 0 ≤ i,j ≤ nd is local dof
   int nd = el.GetDof(); // the number of basis related to current element
   double w; // store weight

   Vector shape; // store basis function values at an integration point
   shape.SetSize(nd);
   elmat.SetSize(nd); // set output size

   // Determine integration rule in the location
   const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el, Trans);

   elmat = 0.0; // initialize with 0
   for (int i = 0; i < ir->GetNPoints(); i++) // for each integration point
   {
      // Get current integration point (x_i, w_i)
      const IntegrationPoint &ip = ir->IntPoint(i);
      Trans.SetIntPoint (&ip);
      // 
      el.CalcPhysShape(Trans, shape);
      // update weight based on the element jacobian and integration weight
      w = Trans.Weight() * ip.weight;
      // elmat = elmat + w * (shape * shape^T)
      // Why? DIY
      AddMult_a_VVt(w, shape, elmat);
   }
} // Also see, DiffusionIntegrator, MixedVectorDivergenceOperator, ...

Coefficient

  • Recall that we used f when we define local linear form
    LinearForm b_h(&fespace);
b_h.AddDomainIntegrator(new DomainLFIntegrator(f)); // b_h += (f,v)
  • Here, f is a Coefficient. More specifically, FunctionCoefficient.
  • Coefficient represents a function
    ConstantCoefficient one_cf(1.0);
FunctionCoefficient f([](Vector &x){return sin(pi*x[0])*sin(pi*x[1]);});
FunctionCoefficient g([](Vector &x, double t){return exp(-t)*sin(pi*x[0]);});
GridFunctionCoefficient u_cf(&u_h); // Use GridFunction as a function
u_h.ProjectCoefficient(f); // interpolate given function f
b_h.AddDomainIntegrator(new DomainLFIntegrator(f)); // b_h += (f,v)

Ax=b and Solve

Essential Boundary

FormLinearSystem

Solvers

RecoverFEMSolution

Visualization

and

  • Now we described our discrete problem in MFEM.
  • Note that the essential BC has not been addressed yet.
  • Recall that for -finite element space,
    we impose Dirichlet boundary condition by
  • Degrees of freedom related to are NOT unknowns.
  • Therefore, we need to remove those degrees of freedom from the final linear system.

Essential Boundary

  • In a mesh file, you can specify boundary attribute
    # Mesh faces/edges on the boundary, e.g. triangles (2)
boundary
8 # 8 boundaries
1 1 0 1 # Bottom (1), Segment, vertex 0 - vertex 1
1 1 1 2 # Borrom (1), Segment, vertex 1 - vertex 2
2 1 2 4 # Right (2), Segment, vertex 2 - vertex 4
2 1 4 8 # Right (2), Segment, vertex 4 - vertex 8
...
  • When is bottom,
    Mesh mesh(mesh_file);
// create an integer array of size, the number of boundary attributes
Array<int> ess_bdr(mesh.bdr_attributes.Max()); // in this case, 4
ess_bdr = 0; // initialize with 0
ess_bdr[0] = 1; // first boundary, bottom, is essential boundary
// stores all degrees of freedom on the Dirichlet Boundary
Array<int> ess_tdof_list;
fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list); // Obtain related dofs

FormLinearSystem

  • To obtain the linear system from BilinearForm a_h and LinarForm b_h, we can
    OperatorPtr A;
Vector B, X;
a_h.FormLinearSystem(ess_tdof_list, u, b_h, A, X, B);
  • Here,
    • u: GridFunction, the solution with specified value on essential boundary
    • A: OperatorPtr or SparseMatrix, the resulting matrix
    • B: Vector, the resulting load vector (also called dual vector)
    • X: Vector, the vector containing coefficient, (also called primal vector)

Solve Linear System

  • To solve a system, there are many solvers in MFEM.
  • Direct Sparse Solvers: UMFPackSolver (from SuiteSparse), and MUMPS (from PETSc)
  • Iterative Solvers:
    • CG: (Preconditioned) Conjugate gradient method for SPD system
      CGSolver cg();
cg.SetOperator(*A);
cg.Mult(B, X);
    • GMRES: (Preconditioned) generalized minimal residual method for general system
      GMRESSolver gmres();
gmres.SetOperator(*A);
gmres.Mult(B, X);

Recover Solution

  • Now, your solution is saved in Vector X.
  • You can convert it back to a GridFunction by using
    a_h.RecoverFEMSolution(X, b, u);
  • Then you can save and visualize using
    u.Save(solution.gf);
sout << "solution\n" << mesh << u; // plot solution
sout << "view 0 0\n"; // top view
sout << flush; // draw
  • For more about visualization, see ex1.cpp and glvis.
  • Note that you need to run glvis in the background!
  • For 3D visualization or high-quality results, use ParaView.

Parallel Implementation

  • make parallel - Requires additional dependencies, see instruction.
  • Usually, you can just put Par before types. See, ex0p.cpp.
    Mpi::Init(); // Initialize MPI
H1_FECollection fec(order, dim);
ParMesh mesh(mesh_file); // Elements are distributed in each process

ParFiniteElementSpace fes(&mesh, &fec);
ParBilinearForm a_h(&fes);
ParLinearForm b_h(&fes);
ParGridFunction u_h(&fes);

a_h.AddDomainIntegrator(new DiffusionIntegrator());
ConstantCoefficient one_cf(1.0);
b_h.AddDomainIntegrator(new DomainLFIntegrator(one_cf));
...
  • Run with mpirun -np <num_process> <file>

Questions

structure_all

https://dohyun-cse.github.io/mfem-tutorial

Hands-on session: ex0.cpp

  • feel free to ask C++ related questions!

# FunctionCoefficient * On the other hand, `FunctionCoefficient` represents a function $\mathbb{R}^n\rightarrow\mathbb{R}$. ```cpp // 2D coordinates x=(x,y) ↦ sin(x)*sin(y) FunctionCoefficient f([](const Vector& x){return sin(x[0])*sin(x[1])}); ``` * For a vector-valued function, we use `VectorFunctionCoefficient` ```cpp // 2D coordinates x=(x,y) ↦ (sin(x), sin(y)) VectorFunctionCoefficient f(2, [](const Vector &x, Vector y){y[0] = sin(x[0]); y[1] = sin(x[1])}); ``` * When you want to compute its value at `x`, then use ```cpp Vector x(2); x = 0.5; // equivalent to x[0] = x[1] = 0.5; const double value = f.Eval(x); // if f:R^n → R Vector y(2); f.Eval(x, y); // if f:R^n → R^d ``` ---

# FormIntegrator * Now you have a bilinear form $a_h(\cdot,\cdot)$, `BilinearForm a` * But you did not specify its actual implementation $$ a_h(u, v) = \sum_{T\in\mathcal{T}_h}(∇ u, ∇ v)_T $$ * The local bilinear form $(∇ u, ∇ v)_T$ is called `BilinearFormIntegrator`. * This class should contains: * `AssembleElementMatrix`: if it is an element integrator, $(\cdot,\cdot)_T$ * `AssembleFaceMatrix`: if it is a face integrator, $(\cdot,\cdot)_F$ * ... see, <a href="https://docs.mfem.org/html/bilininteg_8hpp_source.html" target="_blank">`fem/bilininteg.hpp`</a> * Examples are: * `MassIntegrator`: $(u, v)_T$ * `DiffusionIntegrator`: $(∇ u, ∇ v)_T$ * `MixedVectorDivergenceIntegrator`: $(∇ ⋅ u, v)_T$ * `TransposeIntegrator`: convert $(f(u), g(v))$ to $(f(v), g(u))$ --- # FormIntegrator - continued * You can add `BilinearFormIntegrator` to `BilinearForm` ```cpp BilinearForm a_h(&fespace); a_h.AddDomainIntegrator(new DiffusionIntegrator); // a_h = ∑_T(∇ u, ∇ v)_T a_h.AddDomainIntegrator(new MassIntegrator); // a_h = ∑_T [ (∇ u, ∇ v)_T + (u, v)_T ] a_h.AddInteriorFaceIntegrator(new SomeFaceIntegrator); // a_h = ... + ∑_{F_i}(⋅⋅⋅, ⋅⋅⋅)_F a_h.AddBoundaryIntegrator(new SomeFaceIntegrator); // a_h = ... + ∑_{F_b}(⋅⋅⋅, ⋅⋅⋅)_F ``` * For a `LinearForm`, you can use `LinearFormIntegrator` * The most useful one is `DomainLFIntegrator` and `BoundaryLFIntegrator` ```cpp LinearForm b_h(&fespace); b_h.AddDomainIntegrator(new DomainLFIntegrator(f)); // b_h = ∑_T (f, v) b_h.AddBoundaryIntegrator(new BoundaryLFIntegrator(g, nbc_bdr)); // b_h = ... + ∑_{F_{b,N}}(g_N, v)_F ``` ---

---

# Assembly Level

```cpp

a_h.SetAssemblyLevel(AssemblyLevel::<LEVEL>)

```

* `FULL` (or `LEGACY`)

* Construct global sparse matrix

* `ELEMENT`

* Assemble only local matrices `A_loc`. `Ax` can be computed element-wisely,

$$

\texttt{A*x}=\sum_T P_T*(\texttt{A\_loc}*\texttt{x\_loc})

$$

* `PARTIAL`

* Assemble only on a "reference" element, and other computations are done on-the-fly. See, <a href="http://www.lcad.icmc.usp.br/~buscaglia/teaching/fem/2015/slides20150429.pdf" target="_blank">here</a>

* `NONE`

* Store nothing, everything is on-the-fly.