home.md

+# Examples
+
+## Poisson equation
+
+We want to solve the Poisson equation
+
+<img src="/uploads/a38ecd40a8e2300ca4d0fe576da47989/poisson.png" height="70" title="-\Delta u &= f,\quad \text{in }\Omega=\mathbb{S}^2\\u(x_0) &= 0,\quad f(x)\equiv 1" />
+
+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:
+
+```cpp
+#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:
+
+```cpp
+  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.
+
+```cpp
+  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:
+
+```cpp
+    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:
+
+```cpp
+  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:
+
+```cpp
+  Eigen::UmfPackLU<typename Matrix::Matrix> solver;
+
+  DOFVector<double> u(gv, 0);
+  u = solver.compute(A).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`:
+
+```cpp
+  Dune::VTKWriter<typename GridBase::LeafGridView> writer(gridBase->leafGridView());
+  writer.addVertexData(u, "u");
+  writer.write("poisson");
+```
+
+That's it!
\ No newline at end of file