Examples

Poisson equation

We want to solve the Poisson equation

On the unit sphere, where the Laplacian is the Laplace-Beltrami operator. To close the equation, we need to fix at least one point x0.

For the grid, we use as input an ikosaeder mesh that is refined by division of each triangle into four parts, until the desired grid-width is reached. This grid is provided in the examples/grid directory, i.e. ikosaeder.3d in an ALBERTA grid format. For the quartering refinement of surface grids the Dune::FoamGrid backend is appropriate:

#include <dune/grid/albertagrid/albertareader.hh>
#include <dune/foamgrid/foamgrid.hh>

#include <dune/dec/DecGrid.hpp>

int main(int argc, char** argv)
{
  decpde::init(argc, argv);

  assert_msg( argc > 1, "usage: ", argv[0], " grid-filename [numRef]");

  using GridBase = Dune::FoamGrid<2, 3>;

  // read the grid from ALBERTA grid file
  Dune::GridFactory<GridBase> gridFactory;

  std::string filename = argv[1];
  Dune::AlbertaReader<GridBase>().readGrid(filename, gridFactory);
  std::unique_ptr<GridBase> gridBase( gridFactory.createGrid() );

  // create a grid adapter
  DecGrid<GridBase> grid(*gridBase);

Reading grid-files is done using a Dune::GridFactory. Thus, you can read GMsh files or ALBERTA files into any grid data-structure supporting your desired grid-entities. Dune-DEC requires an (oriented) well-centered simplicial grid. Starting from an ikosaeder grid a refinement by quartering guarantees the well-centeredness.

In order to refine a surface grid it is necessary to project the created vertices to the underlying surface. This needs a projection map. For some surfaces these are provided, e.g. for a sphere or a torus. Therefore, instead or reading into the grid directly, we have to pass the reader a projection adapter:

  using GridBase = Dune::FoamGrid<2, 3>;
  using Param = Parametrization<GridBase, SphereMapping<2, 3, UnitRadius>>;

  Dune::GridFactory<Param> gridFactory;
  Dune::AlbertaReader<Param>().readGrid(argv[1], gridFactory);
  std::unique_ptr<GridBase> gridBase( gridFactory.create() );

The class Parametrization is parametrized with the underlying grid-type and a mapping. This mapping is stored in the grid-elements to create projected points during refinement.

  DecGrid<GridBase> grid(*gridBase);
  grid.globalRefine(argc > 2 ? std::atoi(argv[2]) : 5);

Assembling the linear system

The Laplace-Beltrami operator is predefined in Dune-DEC and can be assembled directly. Therefore, we have to create a matrix first and an operator:

    DOFMatrix<double> A;
   
    // create a gridview, to specify for which element to assemble
    using GridView = typename Grid::LeafGridView;
    auto gv = grid.leafGridView();

    LaplaceBeltrami<GridView> laplace(gv);
    laplace.build(A);

    // create a rhs-vector with value 1 on all vertices
    DOFVector<double> f(gv, 0, 1.0);

The Laplace-Beltrami operator is defined on the vertices of the grid. This is a specialization of the more general Laplace-DeRham operator for 0-forms. The Matrix is initialized automatically during assembling of the linear system in the build() method. An operator provides additionally the methods init(), assemble(), and finish() to split up the build-process into three steps.

Boundary conditions

The enforcement of boundary-conditions, or in this case a point-wise condition, is a two-step procedure: 1. clear the row corresponding the DOF you want to specify the condition for. 2. Set the enforced value on the rhs-vector:

  auto boundary = {0};
  for (auto i : boundary) {
    A.clear_dirichlet_row(i);
    f[i] = 0.0;
  }

The resulting linear system A*u=f can now be solved using any linear-solver you prefer. Some solvers are provided in the linear_algebra submodule of Dune-DEC:

  Eigen::UmfPackLU<typename Matrix::Matrix> solver;

  DOFVector<double> u(gv, 0);
  solver.compute(A);
  u = solver.solve(f);

  msg("|A*u-f| = ", residuum(A, u, f));

Visualization

For the visualization, we write the DOFVector in a VTK-format, that can be plotted using ParaView:

  Dune::VTKWriter<GridView> writer(gv);
  writer.addVertexData(u, "u");
  writer.write("poisson");

That's it!