# Add geometry based on a local-finite-element parametrization

requested to merge feature/curved-geometries into master

### Summary

This MR provides a parametrized geometry representing a local-finite-element function in terms of the local-finite-element basis functions. Note, the local finite-element is just a template parameter and thus not a dependency to dune-localfunctions.

#### Example 1

Parametrize a geometry by 2nd order Lagrange elements:

``````// construct a base grid e.g. with piecewise flat elements
YaspGrid<2> grid({1.0,1.0}, {10,10});

// geometry polynomial order
const int order = 2;
// define a local finite-element (from dune-localfunction)
using LocalFiniteElement = LagrangeCubeLocalFiniteElement<double,double,2,order>;
LocalFiniteElement localFE{};

for (const auto& e : elements(grid.leafGridView()))
{
// projection from local coordinates in the reference element
// to 2d coordinates in the base grid to 3d global coordinates
auto X = [geo=e.geometry()](const auto& local) -> FieldVector<double,3>
{
auto x = geo.global(local);
return {x[0], x[1], std::sin(x[0] * x[1])};
};

// construct the parametrized geometry
LocalFiniteElementGeometry geometry{referenceElement(e), localFE, X};
}``````

#### Example 2

First construct a global grid parametrization in form of a vector holding the coordinates of all Lagrange points, then locally build the corresponding geometry. Note that the example uses a global basis from dune-functions, but it can be implemented also low-level

``````// construct a base grid e.g. with piecewise flat elements
YaspGrid<2> grid({1.0,1.0}, {10,10});

using namespace Dune::Functions::BasisFactory;
auto parametrizationBasis
= makeBasis(grid.leafGridView(), power<3>(lagrange<2>(), blockedInterleaved()));

// create avector that holds all the Lagrange points on the curved grid
std::vector<FieldVector<double,3>> points(parametrizationBasis.size());
Functions::interpolate(parametrizationBasis, points,
[](auto const& x) -> FieldVector<double,3> {
return {x[0], x[1], std::sin(x[0] * x[1])};
});

auto localView = parametrizationBasis.localView();
for (const auto& e : elements(grid.leafGridView()))
{
localView.bind(e);
auto const& node = localView.tree().child(0);

// extract the points on the current elements
std::vector<FieldVector<double,3>> localPoints(node.size());
for (std::size_t i = 0; i < node.size(); ++i) {
std::size_t idx = localView.index(node.localIndex(i))[0];
localPoints[i] = points[idx];
}

// construct the parametrized geometry
LocalFiniteElementGeometry geometry{referenceElement(e), node.finiteElement(), localPoints};
}``````

### Benchmark

Using a multi-linear Lagrange basis in the `LocalFiniteElementGeometry` is equivalent to a `MultiLinearGeometry`. If we define a callable based on the linear-combination of local-basis functions, we could also use the `MappedGeometry` and the result is the same. This motivates a simple benchmark, implemented in `dune/geometry/test/benchmark-geometries.cc`.

GeometryType dim dow time(LocalFiniteElementGeometry) time(MultiLinearGeometry) time(MappedGeometry)
simplex 2 2 0.000100264 sec 2.4485e-05 sec 0.000120727 sec
simplex 2 3 0.000128236 sec 0.000144824 sec 0.000176881 sec
cube 2 2 0.000256232 sec 4.4235e-05 sec 0.000590434 sec
cube 2 3 0.000489962 sec 0.000196707 sec 0.00087769 sec
simplex 3 3 0.00184553 sec 0.00201224 sec 0.00332258 sec
cube 3 3 0.00663592 sec 0.00484587 sec 0.010294 sec

So, in most cases the optimized `MultiLinearGeometry` is the winner. In a few cases the `LocalFiniteElementGeometry` might be faster. This depends on the implementation of the local-basis functions and probably also on the benchmark system.

The `MappedGeometry` is typically the slowest, since also the most general geometry. Often it comes close the `LocalFiniteElementGeometry` when using the same local-basis functions.

### ToDo

• Improve the `affine()` function. It is currently too restrictive, since also a multi-linear geometry can be affine even on non-simplex geometries.
Edited by Simon Praetorius