From 48cbe57bf48ab290d63910edc6612200b4d03417 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= Date: Wed, 9 Mar 2022 10:25:41 +0100 Subject: [PATCH 1/5] [test] Add test for jacobian(), jacobianInverse(), and corresponding type defs --- dune/geometry/test/checkgeometry.hh | 134 ++++++++++++++++++---------- 1 file changed, 89 insertions(+), 45 deletions(-) diff --git a/dune/geometry/test/checkgeometry.hh b/dune/geometry/test/checkgeometry.hh index 5b389b0..bc8324b 100644 --- a/dune/geometry/test/checkgeometry.hh +++ b/dune/geometry/test/checkgeometry.hh @@ -17,6 +17,39 @@ namespace Dune { + namespace Impl + { + + // Helper function to convert a matrix to a FieldMatrix + // using only the mv method. + template + Dune::FieldMatrix toFieldMatrix(const M& m){ + Dune::FieldMatrix result; + for (int j=0; j idColumn(0); + idColumn[j] = 1; + FieldVector column; + m.mv(idColumn,column); + for (int k=0; k + bool isIdentity(const Dune::FieldMatrix& m, double tolerance){ + if (rows != cols) + return false; + for(int i=0; i= tolerance) + return false; + return true; + } + + } + /** * \brief Static and dynamic checks for all features of a Geometry * @@ -60,6 +93,12 @@ namespace Dune // Matrix-like type for the return value of the jacobianInverseTransposed method typedef typename TestGeometry::JacobianInverseTransposed JacobianInverseTransposed; + // Matrix-like type for the return value of the jacobian method + typedef typename TestGeometry::Jacobian Jacobian; + + // Matrix-like type for the return value of the jacobianInverse method + typedef typename TestGeometry::JacobianInverse JacobianInverse; + const ctype tolerance = std::sqrt( std::numeric_limits< ctype >::epsilon() ); //////////////////////////////////////////////////////////////////////// @@ -161,60 +200,65 @@ namespace Dune pass = false; } - // Test whether the methods 'jacobianTransposed' and 'jacobianInverseTransposed' - // return matrices that are inverse to each other. - const JacobianTransposed &jt = geometry.jacobianTransposed( x ); - const JacobianInverseTransposed &jit = geometry.jacobianInverseTransposed( x ); + const JacobianTransposed &Jt = geometry.jacobianTransposed( x ); + const JacobianInverseTransposed &Jit = geometry.jacobianInverseTransposed( x ); + const Jacobian &J = geometry.jacobian( x ); + const JacobianInverse &Ji = geometry.jacobianInverse( x ); // Transform to FieldMatrix, so we can have coefficent access and other goodies - // We need some black magic for the transformation, because there is no - // official easy way yet. - // The following code does the transformation by multiplying jt and jit from - // the right by identity matrices. That way, only the mv method is used. - FieldMatrix< ctype, mydim, coorddim > jtAsFieldMatrix; - for (int j=0; j idColumn(0); - idColumn[j] = 1; + auto JtAsFieldMatrix = Impl::toFieldMatrix< ctype, mydim, coorddim >(Jt); + auto JitAsFieldMatrix = Impl::toFieldMatrix< ctype, coorddim, mydim >(Jit); + auto JAsFieldMatrix = Impl::toFieldMatrix< ctype, coorddim, mydim >(J); + auto JiAsFieldMatrix = Impl::toFieldMatrix< ctype, mydim, coorddim >(Ji); - FieldVector column; - jt.mv(idColumn,column); - - for (int k=0; k id = JtAsFieldMatrix * JitAsFieldMatrix; + if(not Impl::isIdentity(id, tolerance)) + { + std::cerr << "Error: jacobianTransposed and jacobianInverseTransposed are not inverse to each other." << std::endl; + std::cout << " id != [ "; + for( int j = 0; j < mydim; ++j ) + std::cout << (j > 0 ? " | " : "") << id[ j ]; + std::cout << " ]" << std::endl; + pass = false; + } } - FieldMatrix< ctype, coorddim, mydim > jitAsFieldMatrix; - for (int j=0; j idColumn(0); - idColumn[j] = 1; - - FieldVector column; - jit.mv(idColumn,column); - - for (int k=0; k id = JiAsFieldMatrix * JAsFieldMatrix; + if(not Impl::isIdentity(id, tolerance)) + { + std::cerr << "Error: jacobian and jacobianInverse are not inverse to each other." << std::endl; + std::cout << " id != [ "; + for( int j = 0; j < mydim; ++j ) + std::cout << (j > 0 ? " | " : "") << id[ j ]; + std::cout << " ]" << std::endl; + pass = false; + } } - - FieldMatrix< ctype, mydim, mydim > id; - FMatrixHelp::multMatrix( jtAsFieldMatrix, jitAsFieldMatrix, id ); - bool isId = true; - for( int j = 0; j < mydim; ++j ) - for( int k = 0; k < mydim; ++k ) - isId &= (std::abs( id[ j ][ k ] - (j == k ? 1 : 0) ) < tolerance); - if( !isId) + // Test whether the methods 'jacobianTransposed' and 'jacobianInverseTransposed' + // are the transposed of 'jacobian' and 'jacobianInverse', respectively. { - std::cerr << "Error: jacobianTransposed and jacobianInverseTransposed are not inverse to each other." << std::endl; - std::cout << " id != [ "; - for( int j = 0; j < mydim; ++j ) - std::cout << (j > 0 ? " | " : "") << id[ j ]; - std::cout << " ]" << std::endl; - pass = false; + if( (JtAsFieldMatrix - JAsFieldMatrix.transposed()).infinity_norm() != 0 ) + { + std::cerr << "Error: jacobian and jacobianTransposed are not transposed to each other." << std::endl; + pass = false; + } } + { + if( (JitAsFieldMatrix - JiAsFieldMatrix.transposed()).infinity_norm() != 0 ) + { + std::cerr << "Error: jacobianInverse and jacobianInverseTransposed are not transposed to each other." << std::endl; + pass = false; + } + } + // Test whether integrationElement returns something nonnegative if( geometry.integrationElement( x ) < 0 ) { @@ -226,7 +270,7 @@ namespace Dune for( int i = 0; i < mydim; ++i ) for( int j = 0; j < mydim; ++j ) for( int k = 0; k < coorddim; ++k ) - jtj[ i ][ j ] += jtAsFieldMatrix[ i ][ k ] * jtAsFieldMatrix[ j ][ k ]; + jtj[ i ][ j ] += JtAsFieldMatrix[ i ][ k ] * JtAsFieldMatrix[ j ][ k ]; if( abs( sqrt( jtj.determinant() ) - geometry.integrationElement( x ) ) > tolerance ) { std::cerr << "Error: integrationElement is not consistent with jacobianTransposed." << std::endl; -- GitLab From bf5e28be07d4bebb31dc7e14b12a8a32c9c68032 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= Date: Wed, 9 Mar 2022 10:26:50 +0100 Subject: [PATCH 2/5] Implement jacobian() and jacobianInverse() in AffineGeometry --- dune/geometry/affinegeometry.hh | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/dune/geometry/affinegeometry.hh b/dune/geometry/affinegeometry.hh index 6c22539..e381206 100644 --- a/dune/geometry/affinegeometry.hh +++ b/dune/geometry/affinegeometry.hh @@ -506,6 +506,12 @@ namespace Dune /** \brief Type for the transposed inverse Jacobian matrix */ typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed; + /** \brief Type for the Jacobian matrix */ + typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian; + + /** \brief Type for the inverse Jacobian matrix */ + typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianInverse; + private: //! type of reference element typedef Geo::ReferenceElement< Geo::ReferenceElementImplementation< ctype, mydimension > > ReferenceElement; @@ -640,6 +646,28 @@ namespace Dune return jacobianInverseTransposed_; } + /** \brief Obtain the Jacobian + * + * \param[in] local local coordinate to evaluate Jacobian in + * + * \returns a copy of the transposed of the Jacobian + */ + Jacobian jacobian ([[maybe_unused]] const LocalCoordinate &local) const + { + return jacobianTransposed_.transposed(); + } + + /** \brief Obtain the Jacobian's inverse + * + * The Jacobian's inverse is defined as a pseudo-inverse. If we denote + * the Jacobian by \f$J(x)\f$, the following condition holds: + * \f[J^{-1}(x) J(x) = I.\f] + */ + JacobianInverse jacobianInverse ([[maybe_unused]] const LocalCoordinate &local) const + { + return jacobianInverseTransposed_.transposed(); + } + friend ReferenceElement referenceElement ( const AffineGeometry &geometry ) { return geometry.refElement_; -- GitLab From 24faf399d4b63f0860488252087485c6a98a0026 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= Date: Wed, 9 Mar 2022 10:27:46 +0100 Subject: [PATCH 3/5] Implement jacobian() and jacobianInverse() in (Cached)MultiLinearGeometry --- dune/geometry/multilineargeometry.hh | 53 ++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/dune/geometry/multilineargeometry.hh b/dune/geometry/multilineargeometry.hh index ff21839..578cee5 100644 --- a/dune/geometry/multilineargeometry.hh +++ b/dune/geometry/multilineargeometry.hh @@ -201,6 +201,12 @@ namespace Dune //! type of jacobian inverse transposed class JacobianInverseTransposed; + /** \brief Type for the Jacobian matrix */ + typedef FieldMatrix< ctype, coorddimension, mydimension > Jacobian; + + /** \brief Type for the inverse Jacobian matrix */ + typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianInverse; + protected: typedef Dune::ReferenceElements< ctype, mydimension > ReferenceElements; @@ -389,6 +395,29 @@ namespace Dune return geometry.refElement(); } + + /** \brief Obtain the Jacobian + * + * \param[in] local local coordinate to evaluate Jacobian in + * + * \returns a copy of the transposed of the Jacobian + */ + Jacobian jacobian (const LocalCoordinate &local) const + { + return jacobianTransposed(local).transposed(); + } + + /** \brief Obtain the Jacobian's inverse + * + * The Jacobian's inverse is defined as a pseudo-inverse. If we denote + * the Jacobian by \f$J(x)\f$, the following condition holds: + * \f[J^{-1}(x) J(x) = I.\f] + */ + JacobianInverse jacobianInverse (const LocalCoordinate &local) const + { + return jacobianInverseTransposed(local).transposed(); + } + protected: ReferenceElement refElement () const @@ -513,6 +542,8 @@ namespace Dune typedef typename Base::JacobianTransposed JacobianTransposed; typedef typename Base::JacobianInverseTransposed JacobianInverseTransposed; + typedef typename Base::Jacobian Jacobian; + typedef typename Base::JacobianInverse JacobianInverse; template< class CornerStorage > CachedMultiLinearGeometry ( const ReferenceElement &referenceElement, const CornerStorage &cornerStorage ) @@ -660,6 +691,28 @@ namespace Dune return Base::jacobianInverseTransposed( local ); } + /** \brief Obtain the Jacobian + * + * \param[in] local local coordinate to evaluate Jacobian in + * + * \returns a copy of the transposed of the Jacobian + */ + Jacobian jacobian (const LocalCoordinate &local) const + { + return jacobianTransposed(local).transposed(); + } + + /** \brief Obtain the Jacobian's inverse + * + * The Jacobian's inverse is defined as a pseudo-inverse. If we denote + * the Jacobian by \f$J(x)\f$, the following condition holds: + * \f[J^{-1}(x) J(x) = I.\f] + */ + JacobianInverse jacobianInverse (const LocalCoordinate &local) const + { + return jacobianInverseTransposed(local).transposed(); + } + protected: using Base::refElement; -- GitLab From 9ef54000ca260ddaca1dd0586ed285767faa480c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= Date: Wed, 9 Mar 2022 10:28:43 +0100 Subject: [PATCH 4/5] Implement jacobian() and jacobianInverse() in AxisAlignedCubeGeometry --- dune/geometry/axisalignedcubegeometry.hh | 32 +++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/dune/geometry/axisalignedcubegeometry.hh b/dune/geometry/axisalignedcubegeometry.hh index 62ae06e..3967d65 100644 --- a/dune/geometry/axisalignedcubegeometry.hh +++ b/dune/geometry/axisalignedcubegeometry.hh @@ -88,6 +88,24 @@ namespace Dune { DiagonalMatrix, FieldMatrix >::type JacobianInverseTransposed; + /** + * \brief Return type of jacobian + * + * This is a fast DiagonalMatrix if dim==coorddim, and a FieldMatrix otherwise. + * The FieldMatrix will never contain more than one entry per row, + * hence it could be replaced by something more efficient. + */ + using Jacobian = std::conditional_t, FieldMatrix >; + + /** + * \brief Return type of jacobianInverse + * + * This is a fast DiagonalMatrix if dim==coorddim, and a FieldMatrix otherwise. + * The FieldMatrix will never contain more than one entry per row, + * hence it could be replaced by something more efficient. + */ + using JacobianInverse = std::conditional_t, FieldMatrix >; + /** \brief Constructor from a lower left and an upper right corner \note Only for dim==coorddim @@ -184,7 +202,7 @@ namespace Dune { return result; } - /** \brief Jacobian transposed of the transformation from local to global coordinates */ + /** \brief Inverse Jacobian transposed of the transformation from local to global coordinates */ JacobianInverseTransposed jacobianInverseTransposed([[maybe_unused]] const LocalCoordinate& local) const { JacobianInverseTransposed result; @@ -196,6 +214,18 @@ namespace Dune { return result; } + /** \brief Jacobian of the transformation from local to global coordinates */ + Jacobian jacobian([[maybe_unused]] const LocalCoordinate& local) const + { + return jacobianTransposed(local).transposed(); + } + + /** \brief Inverse Jacobian of the transformation from local to global coordinates */ + JacobianInverse jacobianInverse([[maybe_unused]] const LocalCoordinate& local) const + { + return jacobianInverseTransposed(local).transposed(); + } + /** \brief Return the integration element, i.e., the determinant term in the integral transformation formula */ -- GitLab From e9a9938fb078754c8769d6b3c0e19c49fbb93b28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= Date: Wed, 9 Mar 2022 11:09:23 +0100 Subject: [PATCH 5/5] [doc] Mention extended jacobian interface of geometry classes in the change log. --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index fee697e..fcd3c35 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,10 @@ # Master (will become release 2.9) +- The `Geometry` interface was extended by methods `jacobian(local)` and `jacobianInverse(local)` + and corresponding typedefs `Jacobian` and `JacobianInverse`. This is implemented by all geometry + implementations provided by dune-geometry. But external implementations need to be adjusted + to pass the interface check provided by `checkgeometry.hh`. + - The `Geometry::integrationElement` now needs to return the type `Volume` instead of `ctype`. Note that this may be different from `ctype` if the geometry supports typed dimensions. In such case, `ctype` is a length, and not -- GitLab