Skip to content
Snippets Groups Projects

Add Geometry::jacobianInverse()

Merged Carsten Gräser requested to merge feature/jacobianinverse into master
All threads resolved!

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

Pipeline #45726 passed

Pipeline passed for c24f6ef4 on feature/jacobianinverse

Approval is optional

Merged by Carsten GräserCarsten Gräser 2 years ago (Jun 10, 2022 2:24pm UTC)

Merge details

  • Changes merged into master with e20fce13.
  • Deleted the source branch.

Pipeline #45863 failed

Pipeline failed for e20fce13 on master

Activity

Filter activity
  • Approvals
  • Assignees & reviewers
  • Comments (from bots)
  • Comments (from users)
  • Commits & branches
  • Edits
  • Labels
  • Lock status
  • Mentions
  • Merge request status
  • Tracking
  • I canceled the CI pipeline, because is will fail as long as the corresponding MR in dune have not been merged.

    • Resolved by Carsten Gräser

      Maybe one also wants a cast to a suitable FieldMatrix?

      Isn't this what we except the return type of jacobianInverseTranspose to do as well? I think that is reasonable. In the cases where operator* is not enough one can use that to get a richer interface. I needed to do that with a jacobianInverseTranspose as well in some code I had if I remember correctly.

      A question concerning your implementation: since any grid geometry should use the wrapper class it in this form not possible for a given grid implementation to provide that method, is that right? If possible I would suggest to have this default and test if the implementation provides it. We can then decide to deprecate the default implementation at some point if we want to.

      I think providing jacobian makes sense - that seems expected with jacT,jacInv,jacInvT being available.

  • Carsten Gräser added 2 commits

    added 2 commits

    • 561257f8 - Add Geometry::jacobian() and Geometry::jacobianInverse()
    • f9f357fa - Add Geometry::jacobian() and Geometry::jacobianInverse() to static grid check

    Compare with previous version

  • I updated the proposal as follows:

    • The default-implementation is now only activated if the grid implementation does not provide one.
    • The method is added to the static grid check.
    • The method jacobian() and typedef Jacobian have been added analogously.
  • Carsten Gräser changed the description

    changed the description

  • Carsten Gräser changed the description

    changed the description

  • mentioned in merge request dune-geometry!193 (merged)

  • Carsten Gräser changed the description

    changed the description

  • Carsten Gräser resolved all threads

    resolved all threads

  • The tests are now extended in dune-geometry!193 (merged).

  • Carsten Gräser marked this merge request as ready

    marked this merge request as ready

  • added 1 commit

    • 7aeaf07d - [doc] Mention extended jacobian interface of geometries in change log

    Compare with previous version

  • IMO the feature is now ready.

  • Christian Engwer
  • Simon Praetorius
  • mentioned in merge request dune-common!1101 (closed)

  • Carsten Gräser added 46 commits

    added 46 commits

    • 7aeaf07d...ce5efe9e - 43 commits from branch master
    • 5ff53085 - Add Geometry::jacobian() and Geometry::jacobianInverse()
    • 596ea49e - Add Geometry::jacobian() and Geometry::jacobianInverse() to static grid check
    • 87abeae3 - [doc] Mention extended jacobian interface of geometries in change log

    Compare with previous version

  • Carsten Gräser resolved all threads

    resolved all threads

  • Carsten Gräser resolved all threads

    resolved all threads

  • Carsten Gräser changed the description

    changed the description

  • Carsten Gräser mentioned in merge request !588 (merged)

    mentioned in merge request !588 (merged)

  • Carsten Gräser added 7 commits

    added 7 commits

    • 87abeae3...a28dfdf4 - 4 commits from branch master
    • 30c840d7 - Add Geometry::jacobian() and Geometry::jacobianInverse()
    • 0499e3f7 - Add Geometry::jacobian() and Geometry::jacobianInverse() to static grid check
    • 3155cdf1 - [doc] Mention extended jacobian interface of geometries in change log

    Compare with previous version

    • Resolved by Carsten Gräser

      This implements what was decided on the last meeting. Any objections to merging this?

      Notice that there's a follow-up question: Currently there are default implementation of Jacobian, JacobianInverse, jacobian(local), and jacobianInverse(local) in the Dune::Geometry interface class of grid geometries. The Geometry implementations from dune-geometry should work independently and thus provide these on their own (cf. dune-geometry!193 (merged)). Should grid geometry implementations also provide these members by themselves in the future?

      I would strongly opt for yes because otherwise these no longer implement the Geometry interface from dune-geometry. If we agree on this I would open another MR deprecating the default implementation in the interface class Dune::Geometry and providing implementations of those members for all grids in dune-grid.

      Edited by Carsten Gräser
  • Carsten Gräser added 3 commits

    added 3 commits

    • a5308a3a - Add Geometry::jacobian() and Geometry::jacobianInverse()
    • e142d12f - Add Geometry::jacobian() and Geometry::jacobianInverse() to static grid check
    • c24f6ef4 - [doc] Mention extended jacobian interface of geometries in change log

    Compare with previous version

  • Carsten Gräser mentioned in merge request !590 (merged)

    mentioned in merge request !590 (merged)

  • Carsten Gräser resolved all threads

    resolved all threads

  • Carsten Gräser mentioned in commit e20fce13

    mentioned in commit e20fce13

    • Resolved by Carsten Gräser

      @carsten.graeser This breaks at least dune-spgrid (in code not using any of the new interfaces). Seems to have some hidden requirement that wasn't there before. Should this be fixed in dune-spgrid? Or should the implementation of transpose be more general (i.e. not require row_type) in your opinion? I don't really see where the row_type is required other than in the SFINAE construct.

      /dune/modules/dune-grid/dune/grid/common/geometry.hh:136:54: error: no matching function for call to transpose(Dune::SPJacobianInverseTransposed<double, 2, 2>)
        136 |     using JacobianInverseDefault = decltype(transpose(std::declval<JacobianInverseTransposed>()));
            |                                             ~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      In file included from /dune/modules/dune-grid/dune/grid/common/geometry.hh:14,
                       from /dune/modules/dune-grid/dune/grid/common/grid.hh:1095,
                       from /dune/modules/dune-grid/dune/grid/common/intersection.hh:7,
                       from /dune/modules/dune-grid/dune/grid/io/file/dgfparser/gridptr.hh:23,
                       from /dune/modules/dune-grid/dune/grid/io/file/dgfparser/dgfparser.hh:28,
                       from /builds/dumux-repositories/dumux/dumux/io/grid/gridmanager_base.hh:37,
                       from /builds/dumux-repositories/dumux/dumux/io/grid/gridmanager.hh:29,
                       from /builds/dumux-repositories/dumux/test/freeflow/navierstokes/periodic/main.cc:37:
      /dune/modules/dune-common/dune/common/transpose.hh:182:6: note: candidate: template<class Matrix, class, typename std::enable_if<Dune::Impl::HasMemberFunctionTransposed<M>::value, int>::type <anonymous> > auto Dune::transpose(const Matrix&)
        182 | auto transpose(const Matrix& matrix) {
            |      ^~~~~~~~~
      /dune/modules/dune-common/dune/common/transpose.hh:182:6: note:   template argument deduction/substitution failed:
      /dune/modules/dune-common/dune/common/transpose.hh:180:3: error: no type named row_type in class Dune::SPJacobianInverseTransposed<double, 2, 2>
        180 |   typename = std::void_t<typename Matrix::row_type>,
            |   ^~~~~~~~
      Edited by Timo Koch
  • Timo Koch resolved all threads

    resolved all threads

  • @carsten.graeser: I just saw that the jacobianInverseImpl and jacobianImpl methods are dprecated. I guess this is to notify of the change. But why not keeping these as a default methods? The code looks good and it mainly would be exactly what I would implement in a grid until it becomes a performance problem.

  • Please register or sign in to reply
    Loading