MFEM 开源有限元库使用指南

约 1725 字大约 6 分钟

MFEM 有限元 C++开源软件

2026-03-17

MFEM（Modular Finite Element Methods）是一个轻量级、通用、可扩展的 C++ 有限元库，用于求解偏微分方程。本文将详细介绍 MFEM 的特性、安装和使用方法。

MFEM 简介

MFEM 由劳伦斯利弗莫尔国家实验室（LLNL）开发，是一个现代化的有限元库，具有以下特点：

模块化设计：灵活的架构，易于扩展
高性能：支持 MPI 并行和 GPU 加速
丰富的有限元：支持 H1、H(curl)、H(div)、L2 等多种元
自适应网格细化：AMR（Adaptive Mesh Refinement）
多种求解器：集成 HYPRE、PETSc、SuperLU 等
开源免费：LGPL 2.1 许可证

官网：https://mfem.org/

安装 MFEM

从源码编译

# 下载源码
git clone https://github.com/mfem/mfem.git
cd mfem

# 串行版本
make serial -j4

# 并行版本（需要 MPI）
make parallel -j4

# 安装
sudo make install

使用包管理器

# Spack
spack install mfem

# Conda
conda install -c conda-forge mfem

CMake 构建

mkdir build && cd build
cmake .. \
  -DMFEM_USE_MPI=ON \
  -DMFEM_USE_METIS=ON \
  -DMFEM_USE_HYPRE=ON
make -j4
sudo make install

MFEM 核心概念

1. 网格（Mesh）

MFEM 支持多种网格格式：

#include "mfem.hpp"
using namespace mfem;

// 从文件读取网格
Mesh mesh("mesh.msh");

// 生成简单网格
Mesh mesh = Mesh::MakeCartesian2D(
    10, 10,           // 网格数量
    Element::TRIANGLE // 单元类型
);

// 细化网格
mesh.UniformRefinement();

2. 有限元空间（FiniteElementSpace）

定义离散空间：

// H1 空间（Lagrange 元）
H1_FECollection fec(2, mesh.Dimension());  // P2 元
FiniteElementSpace fespace(&mesh, &fec);

// H(div) 空间（RT 元）
RT_FECollection rt_fec(0, mesh.Dimension());  // RT0
FiniteElementSpace rt_fespace(&mesh, &rt_fec);

// L2 空间（DG 元）
L2_FECollection l2_fec(1, mesh.Dimension());  // P1 DG
FiniteElementSpace l2_fespace(&mesh, &l2_fec);

// H(curl) 空间（Nedelec 元）
ND_FECollection nd_fec(1, mesh.Dimension());  // Nedelec 1st kind
FiniteElementSpace nd_fespace(&mesh, &nd_fec);

3. 网格函数（GridFunction）

表示有限元解：

GridFunction u(&fespace);

// 初始化
u = 0.0;

// 从函数投影
FunctionCoefficient u0(initial_condition);
u.ProjectCoefficient(u0);

// 保存到文件
ofstream sol_ofs("solution.gf");
sol_ofs.precision(8);
u.Save(sol_ofs);

4. 双线性形式（BilinearForm）

定义刚度矩阵：

BilinearForm a(&fespace);

// 添加积分子
a.AddDomainIntegrator(new DiffusionIntegrator);  // -∇·(κ∇u)
a.AddDomainIntegrator(new MassIntegrator);       // u
a.AddDomainIntegrator(new ConvectionIntegrator); // b·∇u

// 组装矩阵
a.Assemble();

5. 线性形式（LinearForm）

定义右端项：

LinearForm b(&fespace);

// 添加积分子
FunctionCoefficient f(source_function);
b.AddDomainIntegrator(new DomainLFIntegrator(f));

// 组装向量
b.Assemble();

求解 Poisson 方程

完整示例：求解 $-\Delta u = f$ 在 $\Omega = [0,1]^2$ 上， $u = 0$ 在 $\partial\Omega$ 。

#include "mfem.hpp"
#include <fstream>
#include <iostream>

using namespace std;
using namespace mfem;

// 右端项
double rhs_func(const Vector &x) {
    return 2.0 * M_PI * M_PI * sin(M_PI * x(0)) * sin(M_PI * x(1));
}

// 精确解
double exact_sol(const Vector &x) {
    return sin(M_PI * x(0)) * sin(M_PI * x(1));
}

int main(int argc, char *argv[]) {
    // 创建网格
    int n = 32;
    Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::TRIANGLE);
    int dim = mesh.Dimension();

    // 定义有限元空间（P2 元）
    int order = 2;
    H1_FECollection fec(order, dim);
    FiniteElementSpace fespace(&mesh, &fec);
    cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;

    // 定义边界条件
    Array<int> ess_tdof_list;
    Array<int> ess_bdr(mesh.bdr_attributes.Max());
    ess_bdr = 1;  // 所有边界
    fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);

    // 定义解和右端项
    GridFunction u(&fespace);
    u = 0.0;

    LinearForm b(&fespace);
    FunctionCoefficient rhs(rhs_func);
    b.AddDomainIntegrator(new DomainLFIntegrator(rhs));
    b.Assemble();

    // 定义双线性形式
    BilinearForm a(&fespace);
    a.AddDomainIntegrator(new DiffusionIntegrator);
    a.Assemble();

    // 形成线性系统
    SparseMatrix A;
    Vector B, X;
    a.FormLinearSystem(ess_tdof_list, u, b, A, X, B);

    // 求解线性系统
    GSSmoother M(A);
    PCG(A, M, B, X, 1, 500, 1e-12, 0.0);

    // 恢复解
    a.RecoverFEMSolution(X, b, u);

    // 计算误差
    FunctionCoefficient u_exact(exact_sol);
    double l2_error = u.ComputeL2Error(u_exact);
    double h1_error = u.ComputeH1Error(&u_exact, NULL, 1, 1);
    cout << "L2 error: " << l2_error << endl;
    cout << "H1 error: " << h1_error << endl;

    // 保存解
    ofstream mesh_ofs("mesh.mesh");
    mesh_ofs.precision(8);
    mesh.Print(mesh_ofs);

    ofstream sol_ofs("sol.gf");
    sol_ofs.precision(8);
    u.Save(sol_ofs);

    return 0;
}

编译运行：

g++ -O3 poisson.cpp -o poisson -I/path/to/mfem -L/path/to/mfem -lmfem
./poisson

时间相关问题

求解热传导方程： $\frac{\partial u}{\partial t} - \Delta u = f$

#include "mfem.hpp"

using namespace mfem;

int main() {
    // 网格和有限元空间
    Mesh mesh = Mesh::MakeCartesian2D(32, 32, Element::QUADRILATERAL);
    H1_FECollection fec(2, mesh.Dimension());
    FiniteElementSpace fespace(&mesh, &fec);

    // 质量矩阵和刚度矩阵
    BilinearForm m(&fespace);
    m.AddDomainIntegrator(new MassIntegrator);
    m.Assemble();

    BilinearForm k(&fespace);
    k.AddDomainIntegrator(new DiffusionIntegrator);
    k.Assemble();

    // 初始条件
    GridFunction u(&fespace);
    FunctionCoefficient u0([](const Vector &x) {
        return exp(-10.0 * (pow(x(0) - 0.5, 2) + pow(x(1) - 0.5, 2)));
    });
    u.ProjectCoefficient(u0);

    // 时间积分（向后 Euler）
    double dt = 0.01;
    double t = 0.0;
    double t_final = 1.0;

    SparseMatrix M, K;
    m.FormSystemMatrix(Array<int>(), M);
    k.FormSystemMatrix(Array<int>(), K);

    // M + dt * K
    SparseMatrix A(M);
    A.Add(dt, K);

    Vector U, RHS;
    u.GetTrueDofs(U);

    while (t < t_final) {
        // RHS = M * U
        M.Mult(U, RHS);

        // 求解 (M + dt*K) * U_new = M * U
        GSSmoother prec(A);
        PCG(A, prec, RHS, U, 0, 200, 1e-12, 0.0);

        t += dt;
        u.SetFromTrueDofs(U);

        cout << "t = " << t << endl;
    }

    // 保存最终解
    ofstream sol_ofs("heat_solution.gf");
    u.Save(sol_ofs);

    return 0;
}

并行计算

使用 MPI 并行：

#include "mfem.hpp"

using namespace mfem;

int main(int argc, char *argv[]) {
    // 初始化 MPI
    MPI_Init(&argc, &argv);
    int num_procs, myid;
    MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);

    // 读取网格
    Mesh *mesh = new Mesh("mesh.mesh");
    ParMesh pmesh(MPI_COMM_WORLD, *mesh);
    delete mesh;

    // 并行有限元空间
    H1_FECollection fec(2, pmesh.Dimension());
    ParFiniteElementSpace fespace(&pmesh, &fec);

    if (myid == 0) {
        cout << "Number of unknowns: " << fespace.GlobalTrueVSize() << endl;
    }

    // 并行网格函数
    ParGridFunction u(&fespace);
    u = 0.0;

    // 并行双线性形式
    ParBilinearForm a(&fespace);
    a.AddDomainIntegrator(new DiffusionIntegrator);
    a.Assemble();

    // 并行线性形式
    ParLinearForm b(&fespace);
    FunctionCoefficient rhs([](const Vector &x) { return 1.0; });
    b.AddDomainIntegrator(new DomainLFIntegrator(rhs));
    b.Assemble();

    // 形成并行线性系统
    HypreParMatrix A;
    Vector B, X;
    Array<int> ess_tdof_list;
    a.FormLinearSystem(ess_tdof_list, u, b, A, X, B);

    // 使用 HYPRE 求解器
    HypreBoomerAMG prec(A);
    prec.SetPrintLevel(0);

    HyprePCG solver(A);
    solver.SetTol(1e-12);
    solver.SetMaxIter(500);
    solver.SetPrintLevel(2);
    solver.SetPreconditioner(prec);
    solver.Mult(B, X);

    // 恢复解
    a.RecoverFEMSolution(X, b, u);

    // 保存并行解
    ostringstream sol_name;
    sol_name << "solution." << setfill('0') << setw(6) << myid;
    ofstream sol_ofs(sol_name.str().c_str());
    sol_ofs.precision(8);
    u.Save(sol_ofs);

    MPI_Finalize();
    return 0;
}

编译并行版本：

mpicxx -O3 parallel_poisson.cpp -o parallel_poisson \
    -I/path/to/mfem -L/path/to/mfem -lmfem -lHYPRE
mpirun -np 4 ./parallel_poisson

可视化

使用 GLVis

MFEM 配套的可视化工具：

# 启动 GLVis 服务器
glvis

# 在代码中发送数据
socketstream sol_sock("localhost", 19916);
sol_sock.precision(8);
sol_sock << "solution\n" << mesh << u << flush;

导出到 ParaView

ParaViewDataCollection paraview_dc("output", &mesh);
paraview_dc.SetPrefixPath("ParaView");
paraview_dc.SetLevelsOfDetail(order);
paraview_dc.SetDataFormat(VTKFormat::BINARY);
paraview_dc.SetHighOrderOutput(true);
paraview_dc.SetCycle(0);
paraview_dc.SetTime(0.0);
paraview_dc.RegisterField("solution", &u);
paraview_dc.Save();

高级功能

1. 自适应网格细化

ThresholdRefiner refiner(error_estimator);
refiner.SetTotalErrorFraction(0.7);
refiner.Apply(mesh);
fespace.Update();
u.Update();

2. 非线性问题

class NonlinearOperator : public Operator {
public:
    void Mult(const Vector &x, Vector &y) const override {
        // 实现非线性算子
    }
    Operator &GetGradient(const Vector &x) const override {
        // 返回 Jacobian
    }
};

NewtonSolver newton_solver;
newton_solver.SetOperator(nonlinear_op);
newton_solver.Mult(b, u);

3. 混合有限元

// Darcy 流：u + K∇p = 0, ∇·u = f
RT_FECollection u_fec(0, dim);
L2_FECollection p_fec(0, dim);

ParFiniteElementSpace u_fespace(&pmesh, &u_fec);
ParFiniteElementSpace p_fespace(&pmesh, &p_fec);

// 构造鞍点系统
BlockOperator A(block_offsets);
A.SetBlock(0, 0, &M);  // 质量矩阵
A.SetBlock(0, 1, &B);  // 散度算子
A.SetBlock(1, 0, &Bt); // 散度算子转置

学习资源

官方文档：https://mfem.org/
示例程序：mfem/examples/ 目录包含 30+ 个示例
教程：https://mfem.org/tutorial/
论文：MFEM 相关的学术论文
GitHub：https://github.com/mfem/mfem

总结

MFEM 是一个功能强大、设计优雅的有限元库，适合从简单的教学示例到大规模并行计算的各种应用。其主要优势包括：

易用性：清晰的 API 设计，丰富的示例
灵活性：支持多种有限元类型和网格
性能：高效的并行实现和 GPU 支持
可扩展性：模块化设计，易于添加新功能

无论是学习有限元方法，还是开发高性能科学计算应用，MFEM 都是一个优秀的选择。建议从官方示例开始，逐步深入学习其高级功能。

版权所有

版权归属：Guisong Wu

许可证：署名 4.0 国际 (CC-BY-4.0)