Add Geometry::jacobianInverse()
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 classDune::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 aGeometry
provided by a grid provides a different interface than e.g. aMultilinearGeometry
. 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 aFieldMatrix
usingoperator*
. 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 suitableFieldMatrix
?ShouldGeometry::jacobian()
also be implemented in a similar way? I currently don't have a use case for this, but having both,jacobian()
andjacobianTransposed()
seems to be reasonable. This would e.g. allow to writeg.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 theTransposedMatrixWrapper
is already tested in dune-common and I also tested with the dune-functions examples, where this allows to replacing the ridiculoustranspose(geometry.jacobianInverseTransposed(x)
by a plaingeometry.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 oftranspose()
.dune-common!1101 (closed) extends thetranspose()
function to support capturing matrices by value.-
dune-common!1138 (merged) extends support for
transpose(matrix)
andmatrix.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 ontranspose()
is used. - dune-geometry!193 (merged) extends the actual geometry check and the implementations provided in dune-geomtry by the new interface.
Merge request reports
Activity
- 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 whereoperator*
is not enough one can use that to get a richer interface. I needed to do that with ajacobianInverseTranspose
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 withjacT,jacInv,jacInvT
being available.
- Resolved by Carsten Gräser
mentioned in merge request dune-geometry!193 (merged)
The tests are now extended in dune-geometry!193 (merged).
added 1 commit
- 7aeaf07d - [doc] Mention extended jacobian interface of geometries in change log
- Resolved by Carsten Gräser
- Resolved by Carsten Gräser
mentioned in merge request dune-common!1101 (closed)
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
Toggle commit list-
7aeaf07d...ce5efe9e - 43 commits from branch
mentioned in merge request !588 (merged)
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
Toggle commit list-
87abeae3...a28dfdf4 - 4 commits from branch
- 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)
, andjacobianInverse(local)
in theDune::Geometry
interface class of grid geometries. TheGeometry
implementations fromdune-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 classDune::Geometry
and providing implementations of those members for all grids in dune-grid.Edited by Carsten Gräser
mentioned in merge request !590 (merged)
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 indune-spgrid
? Or should the implementation oftranspose
be more general (i.e. not requirerow_type
) in your opinion? I don't really see where therow_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
mentioned in merge request extensions/dune-spgrid!27 (merged)
@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.