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();
}

mentioned in issue #22

added 1 commit

• d60a97a0 - Some tweaking of the MappingGeometry, removing unused fields and types and added some documentation

Compare with previous version

marked this merge request as ready

changed the description

changed the description

changed the description

added 1 commit

• f102d2d1 - Improve the documentation

Compare with previous version

added feature label

mentioned in issue dune-common#224 (closed)

changed the description

mentioned in issue staging/dune-functions#69 (closed)

I don't understand the 'context' concept and why it is required. We have something like this in dune-fem and there we only need to bind to an entity - that should provide all the information. So what does 'lf.localContext()' return and can the 'context' simply be an entity?

Btw; I'm doing some higher order method and would

• like some higher order derivative information in the geometry.
• like an order method which I can use to increase the quadrature order appropriately.

Regarding higher-order derivative information in the geometry. Since a geometry is "just" a fancy interface for a local-function, this should be possible to incorporate, but the current interface definition has no notion of any higher-order derivatives of the geometry parametrization. In the mapped geometry introduced here, it would require to have higher-order differentiability of the underlying base geometry and of the mapping. It is a nice idea to think about higher-order derivatives especially in the context of higher-order geometries.

changed the description

changed the description

The algorithm should be properly documented. It seems that this is essentially a Gauss-Newton method without regularization and damping.

Yes. And maybe one could improve this also to make it a bit more stable. This was not necessary for me so far, since the initial guess from the flat geometry was already very close and thus the newton method converged very rapidly, but in other cases this might be more delicate.

A simple dumping strategy would be:

    LocalCoordinate xOld = flatGeometry().local(globalCoord), x;
    LocalCoordinate dx;
    GlobalCoordinate dglobal = global(xOld) - globalCoord;
    ctype alpha = 1.0, resOld = dglobal.two_norm(), res = 0.0;

    for (int i = 0; i < Traits::maxIteration(); ++i)
    {
      // Newton's method: DF^n dx^n = F^n, x^{n+1} -= dx^n
      const bool invertible = MatrixHelper::xTRightInvA(jacobianTransposed(x), dglobal, dx);

      // break if jacobian is not invertible
      if (!invertible)
        return std::nullopt;

      // update x with correction
      for (int j = 0; j < Traits::maxIteration(); ++j) {
        x = xOld - alpha * dx;
        dglobal = global(x) - globalCoord;
        res = dglobal.two_norm();
        alpha /= 2.0;

        if (res < resOld)
          break;
      }

      // break if tolerance is reached.
      if (dx.two_norm2() < Traits::tolerance())
        return x;

      alpha *= 2.0;
    }

    if (dx.two_norm2() > Traits::tolerance())
      return std::nullopt;

    return x;

But, maybe it would be required to construct a concrete example where this really matters. For curved geometry I would try to not use the local() method, if possible.

This was not necessary for me so far, since the initial guess from the flat geometry was already very close [...]

In cases similar to those of goemtrygrid, where you really have arbitrary lifts the initial guess will not help much.

Sure. I meant: while in my experiments I didn't face such a situation, there are other cases where a simple Newton method will not converge.

The question for me is: How much effort do we want to put into a stable local() method in case of "arbitrary" nonlinear geometries? A simple or slightly more complex dumping strategy would probably be useful, but this does not guarantee convergence. I have replaced the original local() method with a method called checkedLocal() (name is up to discussion), that tries to find a local coordinate but return a std::nullopt in case the method fails for some reason. So, the user can decide what to do with this information. The alternative is an exception that is thrown in the default local() method in case of failure.

added 33 commits

• f102d2d1...afa1aea7 - 29 commits from branch master
• 7d1baa2d - Add geometry mapped by an element-function
• 5d0c29e4 - Some tweaking of the MappingGeometry, removing unused fields and types and added some documentation
• c402a35e - Improve the documentation
• a52b8502 - Remove necessity to bind the derivative to an element

Compare with previous version

As far as I can see this essentially implements the concatenation f o g is terms of the Geometry interface where f and g are provided as differentiable function (in the dune-functions interface) and Geometry, respectively, with g = Id as a special case. The documentation could probably be more clear about this.

BTW: Conceptually, a Geometry is essentially a differentiable chart and not just a differentiable map. That's why you have to pass the reference element explicitly in the special case g = Id.

changed the description

Maybe we could introduce the concept of an invertible function. (similar to differentiable function) with a function inverse(f) that might be again a differentiable function. Then we call the full concept a diffeomorphic function.

Then we could also provide some standard implementations for the inverse based on numerical algorithms, like Newton methods. This way the user is responsible for providing function with all necessary information to implement the geometry interface and the MappedGeometry is "just" an interface-wrapper.

In !192 (merged) I made an attempt to implement the mapped geometry using a (local) differentiable function only, and no relation to any element. In order to make this usable in the same way as in this MR, two special basis geometries are provided: ReferenceElementGeometry (that is comparable with ReferenceElement::Codim<0>::Geometry, but with a slightly more efficient jacobian[Inverse]Transposed() method that just returns an identity matrix) and a LocalDerivativeGeometry that is like the identity mapping, but returns an actual (transposed) jacobian from a geometry class. With this, a local-functions derivative (that is computed w.r.t. global coordinates) can be transformed into a derivative w.r.t. local coordinates while preserving the rest of the functions.

marked this merge request as draft

Closed, since superseded !192 (merged)

closed