Draft: Add geometry mapped by an element-function

Summary

A geometry that is defined as the mapping of a reference geometry

Details

A geometry can be seen as a class interface representing a reference domain R plus a mapping f, such that f(R) is of the same topology as R and f a bijective (diffeomorphic) mapping between R and f(R). The MappedGeometry provides this interface for a given ReferenceElement R and mapping f.

Additionally, the MappedGeometry allows to be defined on a different reference domain than associated to the mapping f. This could be, for example, the domain of a subEntity or element-intersection. In this case, we build a composition of two mappings f and g, where g represents a local-geometry transform between the two reference domains. Let I be the reference domain of the subEntity and g the associated geometry mapping (that is bijective between I and g(I)), then the MappedGeometry represents the composition f o g associated to the reference domain I.

We require for the local function f to fulfill some interface to be allowed in this wrapper class. Let lf be an object of local function, x be a local coordinate in the reference element of the geometry, then (formulated in c++-20 concepts)

template <class LF, class Range, class Domain, class LocalContext,
  class DLF = std::decay_t<decltype(derivative(std::declval<const LF&>())> >
concept LocalFunction = require (const LF& lf, const Domain& x)
{
  { lf(x) }             -> std::convertible_to<Range>;
  { lf.localContext() } -> std::convertible_to<LocalContext>;
  { derivative(lf) }    -> std::regular_invocable<Domain>;

  require (const DLF& dlf, const Domain& x)
  {
    { dlf(x) }          -> std::convertible_to<typename DerivativeTraits<Range(Domain)>::Range>;
  };
};

In this concept we make explicit notion of a "LocalContext" the local-function is associated to. This context refers to an entity or intersection the local-function was bound. This means that the local-function can be evaluated in the local-coordinates of that LocalContext. We thus require, that the local-context provides also a geometry.

Example 1

This example uses dune-functions, just to make the example short. One can easily define a local function following the test-mappedgeometry.cc

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

using Domain = FieldVector<double,2>;
using Range = FieldVector<double,3>;

auto sig = Functions::SignatureTag<Range(Domain)>{};
auto f = Functions::makeDifferentiableFunctionFromCallables(sig, 
  [](const Domain& x) -> FieldVector<double,3> {
    return {x[0], x[1], std::sin(x[0] * x[1])};
  },
  [](const Domain& x) -> FieldMatrix<double,3,2> {
    return {
      {1.0, 0.0},
      {0.0, 1.0},
      {x[1] * std::cos(x[0] * x[1]), x[0] * std::cos(x[0] * x[1])}
    };
  });

// create a differentiable analytic grid function
auto fGridFunction = Functions::makeAnalyticGridViewFunction(f, grid.leafGridView());

// build the corresponding local function
auto fLocalFunction = localFunction(fGridFunction);

for (const auto& e : elements(grid.leafGridView()))
{
  // bind the local function to the grid element
  fLocalFunction.bind(e);

  // construct the mapped geometry
  MappedGeometry geometry{referenceElement(e), fLocalFunction};

  // unbind the local function from the element
  fLocalFunction.unbind();
}

Example 2

In this example also intersection geometries are transformed. Since the local function must be bound to a grid element, we need an additional geometry transformation from the intersection geometry to the element geometry. The local-geometry transform is passed as additional argument. It can also be used to define a geometry on sub-entities of the grid element, by using a reference sub-element geometry mapping.

Let grid, f, fGridFunction, fLocalFunction be defined as in Example 1.

for (const auto& e : elements(grid.leafGridView()))
{
  // bind the local function to the grid element
  fLocalFunction.bind(e);

  for (const auto& is : intersections(element, grid.leafGridView())
  {
    // construct the mapped geometry for the intersection
    MappedGeometry isGeometry{referenceElement(is), fLocalFunction, is.geometryInInside()};
  }

  auto refElem = referenceElement(e);
  for (int i = 0; i < refElem.size(1); ++i) 
  {
    auto edgeRefElem = referenceElement<double,grid.dimension-1>(refElem.type(i,1));

    // construct a mapped geometry for the i'th subentity of codim 1
    // NOTE: this required a well-oriented (twist-free) grid
    MappedGeometry edgeGeometry{edgeRefElem, fLocalFunction, refElem.geometry<1>(i)};
  }

  // unbind the local function from the element
  fLocalFunction.unbind();
}