|
|
# 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 |