From 7d1baa2d0ca6279bb7f39aec98d3c50757e09383 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sat, 19 Jun 2021 13:15:44 +0200 Subject: [PATCH 1/4] Add geometry mapped by an element-function --- dune/geometry/CMakeLists.txt | 1 + dune/geometry/mappedgeometry.hh | 568 ++++++++++++++++++++++ dune/geometry/test/CMakeLists.txt | 3 + dune/geometry/test/test-mappedgeometry.cc | 154 ++++++ 4 files changed, 726 insertions(+) create mode 100644 dune/geometry/mappedgeometry.hh create mode 100644 dune/geometry/test/test-mappedgeometry.cc diff --git a/dune/geometry/CMakeLists.txt b/dune/geometry/CMakeLists.txt index afd8bf2..fe789f4 100644 --- a/dune/geometry/CMakeLists.txt +++ b/dune/geometry/CMakeLists.txt @@ -8,6 +8,7 @@ install(FILES axisalignedcubegeometry.hh dimension.hh generalvertexorder.hh + mappedgeometry.hh multilineargeometry.hh quadraturerules.hh referenceelement.hh diff --git a/dune/geometry/mappedgeometry.hh b/dune/geometry/mappedgeometry.hh new file mode 100644 index 0000000..ac5ab4b --- /dev/null +++ b/dune/geometry/mappedgeometry.hh @@ -0,0 +1,568 @@ +#ifndef DUNE_GEOMETRY_MAPPEDGEOMETRY_HH +#define DUNE_GEOMETRY_MAPPEDGEOMETRY_HH + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Dune { + +/// \brief Default traits class for MappedGeometryTraits +/** + * The MappedGeometry allows tweaking + * some implementation details through a traits class. + * + * This structure provides the default values. + * + * \tparam ct coordinate type + **/ +template +struct MappedGeometryTraits +{ + /// \brief helper structure containing some matrix routines. See affinegeometry.hh + using MatrixHelper = Impl::FieldMatrixHelper; + + /// \brief tolerance to numerical algorithms + static ct tolerance () { return ct(16) * std::numeric_limits::epsilon(); } + + /// \brief maximal number of Newton iteration in `geometry.local(global)` + static int maxIteration () { return 100; } + + /// \brief Geometry is associated to just one GeometryType + template + struct hasSingleGeometryType + { + static const bool v = false; + static const unsigned int topologyId = ~0u; //< optionally, the topologyId of the single GeometryType + }; +}; + +// Placeholder type for a trivial identity matrix without any functionality +struct IdentityMatrix {}; + +/// \brief Efficient implementation of a trivial local geometry for elements +template +struct DefaultLocalGeometry +{ + using LocalCoordinate = FieldVector; + using GlobalCoordinate = FieldVector; + + template + decltype(auto) global (LocalCoordinate&& local) const + { + return std::forward(local); + } + + template + IdentityMatrix jacobianTransposed (LocalCoordinate&& local) const + { + return {}; + } +}; + + +/// \brief Geometry parametrized by a LocalFunction and a LocalGeometry +/** + * \tparam Map Mapping mapping coordinates, must be differentiable + * \tparam LG A local geometry type defining the local coordinates and the input coordinates + * for the local function. See \ref DefaultLocalGeometry + * \tparam TraitsType Parameters of the geometry, see \ref MappedGeometryTraits + * + * The requirements on the traits are documented along with their default, + * MappedGeometryTraits. + **/ +template > +class MappedGeometry +{ +public: + /// type of local coordinates + using LocalCoordinate = typename LG::LocalCoordinate; + + /// type of global coordinates + using GlobalCoordinate = std::result_of_t; + + /// coordinate type + using ctype = typename LocalCoordinate::value_type; + + /// geometry dimension + static const int mydimension = LocalCoordinate::size(); + + /// coordinate dimension + static const int coorddimension = GlobalCoordinate::size(); + + /// type of volume + using Volume = ctype; + + /// type of jacobian transposed + using JacobianTransposed = FieldMatrix; + + /// type of jacobian inverse transposed + using JacobianInverseTransposed = FieldMatrix; + + /// type of the extended Weingarten map + using NormalGradient = FieldMatrix; + +public: + /// type of reference element + using ReferenceElements = Dune::ReferenceElements; + using ReferenceElement = typename ReferenceElements::ReferenceElement; + + /// Parametrization of the geometry + using Traits = TraitsType; + +protected: + using MatrixHelper = typename Traits::MatrixHelper; + static const bool hasSingleGeometryType + = Traits::template hasSingleGeometryType::v; + static const bool isFlatAffine + = hasSingleGeometryType && ((Traits::template hasSingleGeometryType::topologyId) >> 1 == 0); + + /// type of coordinate transformation for subEntities to codim=0 entities + using LocalGeometry = LG; + + /// type of the mapping representation the geometry parametrization + using Mapping = Map; + + /// type of a mapping representing the derivative of `Map` + using DerivativeMapping = std::decay_t()))>; + +public: + /// \brief Constructor from mapping to parametrize the geometry + /** + * \param[in] refElement reference element for the geometry + * \param[in] mapping mapping for the parametrization of the geometry (stored by value) + * \param[in] localGeometry local geometry for local coordinate transformation + * + **/ + template + MappedGeometry (const ReferenceElement& refElement, Map_&& mapping, LG_&& localGeometry) + : refElement_(refElement) + , mapping_(std::forward(mapping)) + , localGeometry_(std::forward(localGeometry)) + {} + + /// \brief Constructor, forwarding to the other constructor that take a reference-element + /** + * \param[in] gt geometry type + * \param[in] mapping mapping for the parametrization of the geometry (stored by value) + * \param[in] localGeometry local geometry for local coordinate transformation + **/ + template + MappedGeometry (GeometryType gt, Map_&& mapping, LG_&& localGeometry) + : MappedGeometry(ReferenceElements::general(gt), + std::forward(mapping), + std::forward(localGeometry)) + {} + + /// \brief Constructor, forwarding to the other constructors, with LocalGeometry is + /// DefaultLocalGeometry. + /** + * \param[in] refGeo geometry type or reference element + * \param[in] mapping mapping for the parametrization of the geometry (stored by value) + **/ + template >, bool> = true> + MappedGeometry (RefGeo&& refGeo, Map_&& mapping) + : MappedGeometry(std::forward(refGeo), + std::forward(mapping), + LG_{}) + {} + + + /// \brief Is this mapping affine? No! Since we do not know anything about the mapping. + bool affine () const + { + return false; + } + + /// \brief Obtain the element type from the reference element + GeometryType type () const + { + return refElement_.type(); + } + + /// \brief Obtain number of corners of the corresponding reference element + int corners () const + { + return refElement_.size(mydimension); + } + + /// \brief Obtain coordinates of the i-th corner + GlobalCoordinate corner (int i) const + { + assert( (i >= 0) && (i < corners()) ); + return global(refElement_.position(i, mydimension)); + } + + /// \brief Obtain the centroid of the mapping's image + GlobalCoordinate center () const + { + return global(refElement_.position(0, 0)); + } + + /// \brief Evaluate the coordinate mapping + /** + * Evaluate the local function in local coordinates + * + * \param[in] local local coordinates + * \returns corresponding global coordinate + **/ + GlobalCoordinate global (const LocalCoordinate& local) const + { + return mapping()(localGeometry().global(local)); + } + + /// \brief Evaluate the inverse coordinate mapping + /** + * \param[in] globalCoord global coordinate to map + * \return corresponding local coordinate + * + * \throws in case of an error indicating that the local coordinate can not be obtained, + * an exception is thrown. See \ref checkedLocal for a variant that returns + * an optional instead. + * + * \note For given global coordinate `y` the returned local coordinate `x` that minimizes + * the following function over the local coordinate space spanned by the reference element. + * \code + * (global( x ) - y).two_norm() + * \endcode + **/ + LocalCoordinate local (const GlobalCoordinate& globalCoord) const + { + auto localCoord = checkedLocal(globalCoord); + if (!localCoord) + DUNE_THROW(Dune::Exception, + "local coordinate can not be recovered from global coordinate " << globalCoord); + + return *localCoord; + } + + /// \brief Evaluate the inverse coordinate mapping + /** + * \param[in] globalCoord global coordinate to map + * \return corresponding local coordinate or nothing, in case the local + * coordinate could not be calculated. The return value is wrapped + * in an optional. + **/ + std::optional checkedLocal (const GlobalCoordinate& globalCoord) const + { + LocalCoordinate x = flatGeometry().local(globalCoord); + LocalCoordinate dx; + + for (int i = 0; i < Traits::maxIteration(); ++i) + { + // Newton's method: DF^n dx^n = F^n, x^{n+1} -= dx^n + const GlobalCoordinate dglobal = global(x) - globalCoord; + const bool invertible = MatrixHelper::xTRightInvA(jacobianTransposed(x), dglobal, dx); + + // break if jacobian is not invertible + if (!invertible) + return std::nullopt; + + // update x with correction + x -= dx; + + // break if tolerance is reached. + if (dx.two_norm2() < Traits::tolerance()) + return x; + } + + if (dx.two_norm2() > Traits::tolerance()) + return std::nullopt; + + return x; + } + + /// \brief Construct a normal vector of the curved element evaluated at + /// a given local coordinate + /** + * \note Implemented for codim=1 entities only, i.e. edges in 2D and faces in 3D + **/ + GlobalCoordinate normal (const LocalCoordinate& local) const + { + GlobalCoordinate n = normalDirection(local); + return n / n.two_norm(); + } + + /// \brief Construct a normal direction (not normalized) of the curved element + /// evaluated at a given local coordinate + /** + * \note Implemented for codim=1 entities only, i.e. edges in 2D and faces in 3D + **/ + GlobalCoordinate normalDirection (const LocalCoordinate& local) const + { + assert(coorddimension == mydimension+1); + return [&]() { + if constexpr ((mydimension == 1) && (coorddimension == 2)) + return normalDirection1D(local); + else if constexpr ((mydimension == 2) && (coorddimension == 3)) + return normalDirection2D(local); + else + return GlobalCoordinate(0); + }(); + } + + /// \brief Obtain the integration element + /** + * If the Jacobian of the mapping is denoted by $J(x)$, the integration + * integration element \f$\mu(x)\f$ is given by + * \f[ \mu(x) = \sqrt{|\det (J^T(x) J(x))|}.\f] + * + * \param[in] local local coordinate to evaluate the integration element in + * + * \returns the integration element \f$\mu(x)\f$. + **/ + ctype integrationElement (const LocalCoordinate& local) const + { + return MatrixHelper::sqrtDetAAT(jacobianTransposed(local)); + } + + /// \brief Obtain the volume of the mapping's image + /** + * Calculates the volume of the entity by numerical integration. Since the polynomial order of the + * Volume element is not known, iteratively compute numerical integrals with increasing order + * of the quadrature rules, until tolerance is reached. + **/ + Volume volume () const + { + using std::abs; + Volume vol0 = volume(QuadratureRules::rule(type(), 1)); + for (int p = 2; p < 10; ++p) { + Volume vol1 = volume(QuadratureRules::rule(type(), p)); + if (abs(vol1 - vol0) < Traits::tolerance()) + return vol1; + + vol0 = vol1; + } + return vol0; + } + + /// \brief Obtain the volume of the mapping's image by given quadrature rules + Volume volume (const QuadratureRule& quadRule) const + { + Volume vol(0); + for (const auto& qp : quadRule) + vol += integrationElement(qp.position()) * qp.weight(); + return vol; + } + + /// \brief Obtain the transposed of the Jacobian + /** + * \param[in] local local coordinate to evaluate Jacobian in + * \returns the matrix corresponding to the transposed of the Jacobian + **/ + JacobianTransposed jacobianTransposed (const LocalCoordinate& local) const + { + if (!dMapping_) { + dMapping_.emplace(derivative(mapping())); + dMapping_->bind(mapping().localContext()); + } + + // coordinate in the localContext of the mapping + auto&& elementLocal = localGeometry().global(local); + + auto&& jLocal = localGeometry().jacobianTransposed(local); + auto&& jGlobal = hostGeometry().jacobianTransposed(elementLocal); + auto&& jacobian = (*dMapping_)(elementLocal); + return AB(jLocal, ABt(jGlobal,jacobian)); + } + + /// \brief obtain the transposed of 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] + **/ + JacobianInverseTransposed jacobianInverseTransposed (const LocalCoordinate& local) const + { + JacobianInverseTransposed out; + MatrixHelper::rightInvA(jacobianTransposed(local), out); + return out; + } + + /// \brief obtain the reference-element related to this geometry + friend ReferenceElement referenceElement (const MappedGeometry& geometry) + { + return geometry.refElement(); + } + +protected: + // the internal stored reference element + const ReferenceElement& refElement () const + { + return refElement_; + } + + // normal vector to an edge line-element + GlobalCoordinate normalDirection1D (const LocalCoordinate& local) const + { + auto J = jacobianTransposed(local); + return GlobalCoordinate{ + J[0][1], + -J[0][0]}; + } + + // normal vector to a triangle or quad face element + GlobalCoordinate normalDirection2D (const LocalCoordinate& local) const + { + auto J = jacobianTransposed(local); + return GlobalCoordinate{ + J[0][1] * J[1][2] - J[0][2] * J[1][1], + J[0][2] * J[1][0] - J[0][0] * J[1][2], + J[0][0] * J[1][1] - J[0][1] * J[1][0]}; + } + +private: + // internal reference to the stored mapping + const Mapping& mapping () const + { + return mapping_; + } + + // internal reference to the stored localgeometry + const LocalGeometry& localGeometry () const + { + return localGeometry_; + } + + // Construct a flat geometry from corner vertices + using FlatGeometry = std::conditional_t, + MultiLinearGeometry>; + const FlatGeometry& flatGeometry () const + { + if (!flatGeometry_) { + std::vector corners; + corners.reserve(refElement_.size(mydimension)); + for (int i = 0; i < refElement_.size(mydimension); ++i) + corners.push_back(global(refElement_.position(i, mydimension))); + + flatGeometry_.emplace(refElement_, std::move(corners)); + } + + return *flatGeometry_; + } + + // Return the geometry of the element the mapping is bound to + using HostGeometry = std::decay_t().localContext().geometry())>; + const HostGeometry& hostGeometry () const + { + if (!hostGeometry_) + hostGeometry_.emplace(mapping().localContext().geometry()); + + return *hostGeometry_; + } + +private: + // matrix-matrix multiplication A*B^t + template + static FieldMatrix ABt (const FieldMatrix& A, const FieldMatrix& B) + { + FieldMatrix ABt; + for (int i = 0; i < n; ++i) + for (int j = 0; j < m; ++j) { + ABt[i][j] = 0; + for (int k = 0; k < l; ++k) + ABt[i][j] += A[i][k] * B[j][k]; + } + + return ABt; + } + + // matrix-matrix multiplication A*B^t where A is a diagonal matrix + template + static FieldMatrix ABt (const DiagonalMatrix& A, const FieldMatrix& B) + { + FieldMatrix ABt; + for (int i = 0; i < n; ++i) + for (int j = 0; j < m; ++j) { + ABt[i][j] = A[i][i] * B[j][i]; + } + + return ABt; + } + + // matrix-matrix multiplication A*B + template + static FieldMatrix AB (const FieldMatrix& A, const FieldMatrix& B) + { + FieldMatrix AB; + for (int i = 0; i < n; ++i) + for (int j = 0; j < m; ++j) { + AB[i][j] = 0; + for (int k = 0; k < l; ++k) + AB[i][j] += A[i][k] * B[k][j]; + } + + return AB; + } + + // matrix-matrix multiplication A*B where A is a diagonal matrix + template + static FieldMatrix AB (const DiagonalMatrix& A, const FieldMatrix& B) + { + FieldMatrix AB; + for (int i = 0; i < n; ++i) + for (int j = 0; j < m; ++j) { + AB[i][j] = A[i][i] * B[i][j]; + } + + return AB; + } + + // multiplication of an identity-matrix with any real matrix + template + static const Mat& AB (const IdentityMatrix& A, const Mat& B) + { + return B; + } + +private: + /// Reference of the geometry + ReferenceElement refElement_; + + /// local parametrization of the element + Mapping mapping_; + + /// transformation of local coordinates to element-local coordinates + LocalGeometry localGeometry_; + + // some data optionally provided + mutable std::optional dMapping_; + mutable std::optional flatGeometry_; + mutable std::optional hostGeometry_; + mutable std::vector nCoefficients_; + mutable std::vector nGradients_; +}; + +// deduction guides +template +MappedGeometry (RefGeo, const Map&, const LG&) + -> MappedGeometry; + +template +MappedGeometry (Geo::ReferenceElement, const Map&) + -> MappedGeometry, + MappedGeometryTraits >; + +// typedef for geometries on grid elements with local geometry is identity +template +using ElementMappedGeometry = MappedGeometry, MappedGeometryTraits >; + +} // end namespace Dune + +#endif // DUNE_GEOMETRY_MAPPEDGEOMETRY_HH diff --git a/dune/geometry/test/CMakeLists.txt b/dune/geometry/test/CMakeLists.txt index 300fadc..c6d2b90 100644 --- a/dune/geometry/test/CMakeLists.txt +++ b/dune/geometry/test/CMakeLists.txt @@ -65,6 +65,9 @@ dune_add_test(SOURCES test-referenceelements.cc dune_add_test(SOURCES test-quadrature.cc LINK_LIBRARIES dunegeometry) +dune_add_test(SOURCES test-mappedgeometry.cc + LINK_LIBRARIES dunegeometry) + dune_add_test(SOURCES test-multilineargeometry.cc LINK_LIBRARIES dunegeometry) diff --git a/dune/geometry/test/test-mappedgeometry.cc b/dune/geometry/test/test-mappedgeometry.cc new file mode 100644 index 0000000..f5f0e90 --- /dev/null +++ b/dune/geometry/test/test-mappedgeometry.cc @@ -0,0 +1,154 @@ +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- +// vi: set et ts=4 sw=2 sts=2: +#include + +#include +#include + +#include +#include +#include +#include +#include + +// A mapping modelling a local element-function +template +struct AffineMapping +{ + // Minimal base element with a geometry + struct LocalContext + { + const Dune::MultiLinearGeometry& geometry_; + const Dune::MultiLinearGeometry& geometry () const { return geometry_; } + }; + + AffineMapping (const Dune::MultiLinearGeometry& geometry) + : localContext_{geometry} + { + for (int i = 0; i < cdim; ++i) + b[i] = i + 1; + for (int i = 0; i < mydim; ++i) + A[i][i] = i + 1; + } + + Dune::FieldVector operator() (const Dune::FieldVector& x) const + { + Dune::FieldVector result; + A.mv(x, result); + result += b; + return result; + } + + // the element, this mapping is bound to + const LocalContext& localContext () const { return localContext_; } + + // bind the mapping to a local element + void bind (const LocalContext&) {} + + + struct DerivativeAffineMapping + { + const LocalContext& localContext_; + const Dune::FieldMatrix& A; + + const auto& operator() (const Dune::FieldVector& x) const + { + return A; + } + + const LocalContext& localContext () const { return localContext_; } + void bind (const LocalContext&) {} + }; + + friend auto derivative (const AffineMapping& mapping) + { + return DerivativeAffineMapping{mapping.localContext_,mapping.A}; + } + +private: + LocalContext localContext_; + Dune::FieldMatrix A; + Dune::FieldVector b; +}; + + +template +static bool testMappedGeometry () +{ + bool pass = true; + + constexpr auto gt = Dune::GeometryType{id}; + auto refElem = Dune::referenceElement(gt); + + std::vector> corners(refElem.size(gt.dim())); + for (int i = 0; i < refElem.size(gt.dim()); ++i) + corners[i] = refElem.position(i,gt.dim()); + + auto localContext = Dune::MultiLinearGeometry{refElem, corners}; + auto mapping = AffineMapping{localContext}; + + // construct a geometry using given parametrization coefficients + auto geometry = Dune::MappedGeometry{refElem, mapping}; + pass &= checkGeometry(geometry); + + return pass; +} + + +template +static bool testMappedGeometry () +{ + bool pass = true; + + // pass &= testMappedGeometry(); + + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + + pass &= testMappedGeometry(); + pass &= testMappedGeometry(); + + return pass; +} + +int main ( int argc, char **argv ) +{ + bool pass = true; + + std::cout << ">>> Checking ctype = double" << std::endl; + pass &= testMappedGeometry< double >(); + //std::cout << ">>> Checking ctype = float" << std::endl; + //pass &= testMappedGeometry< float >(); + + return (pass ? 0 : 1); +} -- GitLab From 5d0c29e4ce87daaa70cc0ad9b092f03b8517263c Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sat, 19 Jun 2021 23:43:13 +0200 Subject: [PATCH 2/4] Some tweaking of the MappingGeometry, removing unused fields and types and added some documentation --- dune/geometry/mappedgeometry.hh | 43 +++++++++++++++++++++------------ 1 file changed, 28 insertions(+), 15 deletions(-) diff --git a/dune/geometry/mappedgeometry.hh b/dune/geometry/mappedgeometry.hh index ac5ab4b..21f0e7a 100644 --- a/dune/geometry/mappedgeometry.hh +++ b/dune/geometry/mappedgeometry.hh @@ -75,13 +75,32 @@ struct DefaultLocalGeometry /// \brief Geometry parametrized by a LocalFunction and a LocalGeometry /** - * \tparam Map Mapping mapping coordinates, must be differentiable - * \tparam LG A local geometry type defining the local coordinates and the input coordinates - * for the local function. See \ref DefaultLocalGeometry + * \tparam Map Mapping of element local coordinates to world coordinates, must be differentiable + * \tparam LG A local geometry type defining the local coordinates of the geometry and its + * transformation into local coordinates of the element the mapping the bound to. + * See \ref DefaultLocalGeometry * \tparam TraitsType Parameters of the geometry, see \ref MappedGeometryTraits * - * The requirements on the traits are documented along with their default, - * MappedGeometryTraits. + * This class represents a geometry that is parametrized by the chained mapping of a local geometry + * mapping and the given function. Thereby, the local geometry mapping related coordinates in + * subentities of an element to element coordinates and the function then maps element coordinates + * to the world coordinates representing the actual geometry. + * + * Typical choices for the local geometry are + * - An identity mapping, if the geometry represents an element geometry. + * - An intersection geometry-in-inside or geometry-in-outside mapping. + * - A reference-element sub-entity geometry mapping. + * + * The function (or geometry mapping) is a local function in terms of the dune-functions concept. It + * requires the following members: + * - `bind(LocalContext)`: Initialize the mapping on a given local context, i.e., an element or + * entity or intersection. + * - `localContext()`: Return the element (or entity or intersection) the mapping is bound to. + * - `operator()(LG::GlobalCoordinate)`: Evaluation in local-context coordinates. + * - `derivative(Map)`: A free function returning a local function that represents the first + * derivative of the mapping. + * + * The requirements on the traits are documented along with their default, `MappedGeometryTraits`. **/ template > @@ -112,9 +131,6 @@ public: /// type of jacobian inverse transposed using JacobianInverseTransposed = FieldMatrix; - /// type of the extended Weingarten map - using NormalGradient = FieldMatrix; - public: /// type of reference element using ReferenceElements = Dune::ReferenceElements; @@ -130,13 +146,13 @@ protected: static const bool isFlatAffine = hasSingleGeometryType && ((Traits::template hasSingleGeometryType::topologyId) >> 1 == 0); - /// type of coordinate transformation for subEntities to codim=0 entities + // type of coordinate transformation for subEntities to codim=0 entities using LocalGeometry = LG; - /// type of the mapping representation the geometry parametrization + // type of the mapping representation the geometry parametrization using Mapping = Map; - /// type of a mapping representing the derivative of `Map` + // type of a mapping representing the derivative of `Map` using DerivativeMapping = std::decay_t()))>; public: @@ -145,7 +161,6 @@ public: * \param[in] refElement reference element for the geometry * \param[in] mapping mapping for the parametrization of the geometry (stored by value) * \param[in] localGeometry local geometry for local coordinate transformation - * **/ template MappedGeometry (const ReferenceElement& refElement, Map_&& mapping, LG_&& localGeometry) @@ -182,7 +197,7 @@ public: {} - /// \brief Is this mapping affine? No! Since we do not know anything about the mapping. + /// \brief Is this mapping affine? Not in general, since we don't know anything about the mapping. bool affine () const { return false; @@ -544,8 +559,6 @@ private: mutable std::optional dMapping_; mutable std::optional flatGeometry_; mutable std::optional hostGeometry_; - mutable std::vector nCoefficients_; - mutable std::vector nGradients_; }; // deduction guides -- GitLab From c402a35e9f79c6308e366f552e8f17fa34b4dd01 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Sun, 18 Jul 2021 12:46:29 +0200 Subject: [PATCH 3/4] Improve the documentation --- dune/geometry/mappedgeometry.hh | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/dune/geometry/mappedgeometry.hh b/dune/geometry/mappedgeometry.hh index 21f0e7a..172321a 100644 --- a/dune/geometry/mappedgeometry.hh +++ b/dune/geometry/mappedgeometry.hh @@ -159,7 +159,8 @@ public: /// \brief Constructor from mapping to parametrize the geometry /** * \param[in] refElement reference element for the geometry - * \param[in] mapping mapping for the parametrization of the geometry (stored by value) + * \param[in] mapping a differentiable local function for the parametrization of the + * geometry (stored by value) * \param[in] localGeometry local geometry for local coordinate transformation **/ template @@ -172,7 +173,8 @@ public: /// \brief Constructor, forwarding to the other constructor that take a reference-element /** * \param[in] gt geometry type - * \param[in] mapping mapping for the parametrization of the geometry (stored by value) + * \param[in] mapping a differentiable local function for the parametrization of the + * geometry (stored by value) * \param[in] localGeometry local geometry for local coordinate transformation **/ template @@ -186,7 +188,8 @@ public: /// DefaultLocalGeometry. /** * \param[in] refGeo geometry type or reference element - * \param[in] mapping mapping for the parametrization of the geometry (stored by value) + * \param[in] mapping a differentiable local function for the parametrization of the + * geometry (stored by value) **/ template >, bool> = true> -- GitLab From a52b85028f2013fe9bb40989bb02a176dcebf7ab Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Wed, 23 Feb 2022 15:11:11 +0100 Subject: [PATCH 4/4] Remove necessity to bind the derivative to an element --- dune/geometry/mappedgeometry.hh | 4 +--- dune/geometry/test/test-mappedgeometry.cc | 17 +---------------- 2 files changed, 2 insertions(+), 19 deletions(-) diff --git a/dune/geometry/mappedgeometry.hh b/dune/geometry/mappedgeometry.hh index 172321a..9786a9d 100644 --- a/dune/geometry/mappedgeometry.hh +++ b/dune/geometry/mappedgeometry.hh @@ -384,10 +384,8 @@ public: **/ JacobianTransposed jacobianTransposed (const LocalCoordinate& local) const { - if (!dMapping_) { + if (!dMapping_) dMapping_.emplace(derivative(mapping())); - dMapping_->bind(mapping().localContext()); - } // coordinate in the localContext of the mapping auto&& elementLocal = localGeometry().global(local); diff --git a/dune/geometry/test/test-mappedgeometry.cc b/dune/geometry/test/test-mappedgeometry.cc index f5f0e90..a17bdbf 100644 --- a/dune/geometry/test/test-mappedgeometry.cc +++ b/dune/geometry/test/test-mappedgeometry.cc @@ -45,24 +45,9 @@ struct AffineMapping // bind the mapping to a local element void bind (const LocalContext&) {} - - struct DerivativeAffineMapping - { - const LocalContext& localContext_; - const Dune::FieldMatrix& A; - - const auto& operator() (const Dune::FieldVector& x) const - { - return A; - } - - const LocalContext& localContext () const { return localContext_; } - void bind (const LocalContext&) {} - }; - friend auto derivative (const AffineMapping& mapping) { - return DerivativeAffineMapping{mapping.localContext_,mapping.A}; + return [&A=mapping.A](const Dune::FieldVector& x) { return A; }; } private: -- GitLab