Add geometry based on a local-finite-element parametrization

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.