Poisson Equation with MFEM + Ginkgo

Solves the Laplace equation $-\Delta u = 1$ with homogeneous Dirichlet boundary conditions using MFEM for finite element assembly and Ginkgo for the linear solve.

Based on MFEM Example 1.

Integration Overview

MFEM’s SparseMatrix stores CSR data directly. After calling Finalize(), the raw arrays are accessible via GetI(), GetJ(), and GetData(). These are wrapped into a gko::matrix::Csr using gko::array::view for zero-copy interoperability.

Step-by-Step

1. Set Up Mesh and Finite Element Space

    // Load mesh from file, or generate a default square mesh.
    Mesh* mesh;
    if (std::string(mesh_file).empty())
        mesh = new Mesh(Mesh::MakeCartesian2D(16, 16, Element::QUADRILATERAL));
    else
        mesh = new Mesh(mesh_file, 1, 1);

    int dim = mesh->Dimension();
    for (int l = 0; l < ref_levels; l++) mesh->UniformRefinement();

    std::cout << "Number of elements: " << mesh->GetNE() << std::endl;

    H1_FECollection fec(order, dim);
    FiniteElementSpace fespace(mesh, &fec);
    std::cout << "Number of unknowns: " << fespace.GetTrueVSize() << std::endl;

2. Assemble the Linear System

    // Determine essential (Dirichlet) boundary DOFs.
    Array<int> ess_tdof_list;
    if (mesh->bdr_attributes.Size()) {
        Array<int> ess_bdr(mesh->bdr_attributes.Max());
        ess_bdr = 1;
        fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }

    // Right-hand side: f = 1.
    LinearForm b(&fespace);
    ConstantCoefficient one(1.0);
    b.AddDomainIntegrator(new DomainLFIntegrator(one));
    b.Assemble();

    // Solution grid function, initialized to zero.
    GridFunction x(&fespace);
    x = 0.0;

    // Stiffness matrix.
    BilinearForm a(&fespace);
    a.AddDomainIntegrator(new DiffusionIntegrator(one));
    a.Assemble();

    // Apply boundary conditions and form the final linear system.
    SparseMatrix A;
    Vector B, X;
    a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
    std::cout << "System size: " << A.Height() << " x " << A.Width()
              << ", nnz: " << A.NumNonZeroElems() << std::endl;

3. Convert to Ginkgo Objects

MFEM exposes CSR pointers directly — GetI() returns row pointers, GetJ() returns column indices, GetData() returns values.

    // MFEM's SparseMatrix stores CSR data directly accessible via
    // GetI() (row pointers), GetJ() (column indices), GetData() (values).
    A.Finalize();
    const int n = A.Height();
    const int nnz = A.NumNonZeroElems();

    auto exec = gko::ReferenceExecutor::create();

    auto gko_A = gko::share(gko::matrix::Csr<double, int>::create(
        exec,
        gko::dim<2>(static_cast<gko::size_type>(n),
                    static_cast<gko::size_type>(n)),
        gko::array<double>::view(exec, nnz, A.GetData()),
        gko::array<int>::view(exec, nnz, A.GetJ()),
        gko::array<int>::view(exec, n + 1, A.GetI())));

    auto gko_b = gko::matrix::Dense<double>::create(
        exec, gko::dim<2>(static_cast<gko::size_type>(n), 1),
        gko::array<double>::view(exec, n, B.GetData()), 1);
    auto gko_x = gko::matrix::Dense<double>::create(
        exec, gko::dim<2>(static_cast<gko::size_type>(n), 1),
        gko::array<double>::view(exec, n, X.GetData()), 1);

4. Solve with Ginkgo

    auto solver =
        gko::solver::Cg<double>::build()
            .with_preconditioner(gko::preconditioner::Jacobi<double>::build())
            .with_criteria(
                gko::stop::Iteration::build().with_max_iters(500u),
                gko::stop::ResidualNorm<double>::build().with_reduction_factor(
                    1e-12))
            .on(exec)
            ->generate(gko_A);

    solver->apply(gko_b, gko_x);
    std::cout << "Ginkgo solver converged." << std::endl;

5. Recover Solution

    // Map algebraic solution back to the FE grid function.
    a.RecoverFEMSolution(X, b, x);

    std::ofstream sol_ofs("solution.gf");
    x.Save(sol_ofs);

    std::ofstream mesh_ofs("mesh.mesh");
    mesh->Print(mesh_ofs);
    std::cout << "Output written to solution.gf and mesh.mesh" << std::endl;

Building

cmake -S . -B build --preset default
cmake --build build
./build/poisson

A mesh file can be passed with -m:

./build/poisson -m /path/to/mfem/data/star.mesh -o 2

Expected Output

Options used:
   --mesh
   --order 1
   --refine 2
Number of elements: 4096
Number of unknowns: 4225
System size: 4225 x 4225, nnz: 37249
Ginkgo solver converged.
Output written to solution.gf and mesh.mesh