Poisson Equation with deal.II + Ginkgo#
Solves the Laplace equation $-\Delta u = 1$ on $[-1,1]^2$ with homogeneous Dirichlet boundary conditions using deal.II for mesh management and assembly, and Ginkgo for the linear solve.
Based on deal.II tutorial step-3.
Integration Overview#
deal.II assembles the system matrix in its own SparseMatrix format. We
extract the CSR structure by iterating over rows, then wrap the data into
a gko::matrix::Csr using non-owning array views. deal.II vectors are
similarly wrapped as gko::matrix::Dense with zero-copy.
Step-by-Step#
1. Create the Mesh#
// Create a uniform mesh on [-1,1]^2 with 2^5 = 32 cells per direction.
Triangulation<2> triangulation;
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(5);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
2. Set Up Finite Elements#
// Use bilinear Q1 finite elements.
const FE_Q<2> fe(1);
DoFHandler<2> dof_handler(triangulation);
dof_handler.distribute_dofs(fe);
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
3. Assemble the Linear System#
deal.II assembles the stiffness matrix and right-hand side using cell-wise quadrature, then applies boundary conditions.
// Build the sparsity pattern and assemble the Laplace matrix and RHS.
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp);
SparsityPattern sparsity_pattern;
sparsity_pattern.copy_from(dsp);
SparseMatrix<double> system_matrix;
system_matrix.reinit(sparsity_pattern);
Vector<double> system_rhs(dof_handler.n_dofs());
Vector<double> solution(dof_handler.n_dofs());
const QGauss<2> quadrature_formula(fe.degree + 1);
FEValues<2> fe_values(fe, quadrature_formula,
update_values | update_gradients | update_JxW_values);
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto& cell : dof_handler.active_cell_iterators()) {
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
for (const unsigned int q : fe_values.quadrature_point_indices()) {
for (const unsigned int i : fe_values.dof_indices()) {
for (const unsigned int j : fe_values.dof_indices())
cell_matrix(i, j) += fe_values.shape_grad(i, q) *
fe_values.shape_grad(j, q) *
fe_values.JxW(q);
cell_rhs(i) +=
fe_values.shape_value(i, q) * 1.0 * fe_values.JxW(q);
}
}
cell->get_dof_indices(local_dof_indices);
for (const unsigned int i : fe_values.dof_indices()) {
for (const unsigned int j : fe_values.dof_indices())
system_matrix.add(local_dof_indices[i], local_dof_indices[j],
cell_matrix(i, j));
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
// Apply homogeneous Dirichlet boundary conditions.
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(
dof_handler, 0, Functions::ZeroFunction<2>(), boundary_values);
MatrixTools::apply_boundary_values(boundary_values, system_matrix, solution,
system_rhs);
4. Extract CSR Data for Ginkgo#
deal.II does not expose raw CSR pointers, so we iterate over the matrix rows to build standard CSR arrays.
// Extract CSR data from deal.II's SparseMatrix by iterating over rows.
// deal.II does not expose raw CSR pointers directly, so we build arrays
// from the row iterators.
const auto n = system_matrix.m();
const auto nnz = system_matrix.n_nonzero_elements();
std::vector<int> row_ptrs(n + 1);
std::vector<int> col_idxs(nnz);
std::vector<double> values(nnz);
row_ptrs[0] = 0;
for (unsigned int row = 0; row < n; ++row) {
int idx = row_ptrs[row];
for (auto it = system_matrix.begin(row); it != system_matrix.end(row);
++it) {
col_idxs[idx] = static_cast<int>(it->column());
values[idx] = it->value();
++idx;
}
row_ptrs[row + 1] = idx;
}
5. Set Up the Ginkgo Executor#
// Create a Ginkgo executor. Switch to CudaExecutor or HipExecutor for GPU.
auto exec = gko::ReferenceExecutor::create();
6. Wrap as Ginkgo Objects#
The CSR arrays and deal.II vectors are wrapped as Ginkgo objects using non-owning views — no data is copied.
// Wrap the CSR arrays as a Ginkgo matrix using non-owning views.
auto mtx = gko::share(gko::matrix::Csr<double, int>::create(
exec, gko::dim<2>(n, n),
gko::array<double>::view(exec, nnz, values.data()),
gko::array<int>::view(exec, nnz, col_idxs.data()),
gko::array<int>::view(exec, n + 1, row_ptrs.data())));
// Wrap deal.II vectors as Ginkgo Dense vectors (zero-copy).
auto b = gko::matrix::Dense<double>::create(
exec, gko::dim<2>(n, 1),
gko::array<double>::view(exec, n, system_rhs.begin()), 1);
auto x = gko::matrix::Dense<double>::create(
exec, gko::dim<2>(n, 1),
gko::array<double>::view(exec, n, solution.begin()), 1);
7. Solve with Ginkgo#
// Build a CG solver with Jacobi preconditioner.
auto solver =
gko::solver::Cg<double>::build()
.with_preconditioner(gko::preconditioner::Jacobi<double>::build())
.with_criteria(
gko::stop::Iteration::build().with_max_iters(1000u),
gko::stop::ResidualNorm<double>::build().with_reduction_factor(
1e-12))
.on(exec)
->generate(mtx);
solver->apply(b, x);
std::cout << "Solver converged." << std::endl;
Building#
cmake -S . -B build --preset default
cmake --build build
./build/poisson
To specify library paths:
cmake -S . -B build --preset default \
-DGinkgo_DIR=/path/to/ginkgo/lib/cmake/Ginkgo \
-DDEAL_II_DIR=/path/to/dealii
Expected Output#
Number of active cells: 1024
Number of degrees of freedom: 1089
Solver converged.
Output written to solution.vtk