Skip to content

Add Geometry::jacobianInverse()

Carsten Gräser requested to merge feature/jacobianinverse into master

This adds a method Geometry::jacobianInverse(localCoordinate) and a corresponding typedef Geometry::JacobianInverse to the Geometry interface class. The method is default implemented using the TransposedMatrixWrapper from dune-common which internally captures (by value) the matrix returned by Geometry::jacobianInverseTransposed().

This is marked WIP, because a few things need to be discussed:

  • Currently the method is default-implemented in the interface class Dune::Geometry. Hence it will work for all grid implementations. But one may argue (I tend to do so) that it should be implemented directly in the geometry implementations. Otherwise a Geometry provided by a grid provides a different interface than e.g. a MultilinearGeometry. To avoid that all grid implementations need to be modified, one could conditionally activate the default-implementation depending on presence of an implemenation in the grid implementation.
  • Which interface do we want to provide? Currently you can do exactly one thing: Multiply the matrix from the right to a FieldMatrix using operator*. This is exactly what you need for the chain rule. This is in some sense the 'linear operator interface' proposed by @peter, but for the case where you want to chain the operator instead of directly applying it. Is this sufficient?
  • Maybe one also wants a cast to a suitable FieldMatrix?
  • Should Geometry::jacobian() also be implemented in a similar way? I currently don't have a use case for this, but having both, jacobian() and jacobianTransposed() seems to be reasonable. This would e.g. allow to write g.jacobianTransposed(x)*g.jacobian(x) to compute the metric tensor which is conceptually nice (but probably inefficient).
  • Tests are missing (in dune-grid) due to the open questions. But the TransposedMatrixWrapper is already tested in dune-common and I also tested with the dune-functions examples, where this allows to replacing the ridiculous transpose(geometry.jacobianInverseTransposed(x) by a plain geometry.jacobianInverse(x).

This is related to several other MRs that should be merged in appropriate order:

  • dune-common!1100 (merged) provides a utility needed for the subsequent extension of transpose().
  • dune-common!1101 (closed) extends the transpose() function to support capturing matrices by value.
  • dune-common!1138 (merged) extends support for transpose(matrix) and matrix.transposed()
  • !577 (merged) (present MR) extends the (grid-)Geometry interface class by the new interface. This is conditionally default-implemented: If the actual grid geometry provides the new interface it is used, otherwise a default implementation based on transpose() is used.
  • dune-geometry!193 (merged) extends the actual geometry check and the implementations provided in dune-geomtry by the new interface.
Edited by Carsten Gräser

Merge request reports