cleaned up the code, added more documentation and restructured sources

dune.module

@@ -7,6 +7,6 @@ Module: dune-curvedgeometry
 Version: 2.7
 Maintainer: simon.praetorius@tu-dresden.de
 # Required build dependencies
-Depends: dune-geometry dune-localfunctions
+Depends: dune-geometry
 # Optional build dependencies
-Suggests: dune-grid dune-vtk
+Suggests: dune-localfunctions dune-grid dune-vtk

dune/curvedgeometry/CMakeLists.txt

 #install headers
-install(FILES curvedgeometry.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/curvedgeometry)
+install(FILES
+  curvedgeometry.hh
+  curvedgeometrydatacollector.hh
+  DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/curvedgeometry)
+
+add_subdirectory(test)
\ No newline at end of file

dune/curvedgeometry/curvedgeometry.hh

@@ -13,11 +13,14 @@
 #include <dune/common/std/type_traits.hh>
 
 #include <dune/geometry/affinegeometry.hh>
+#include <dune/geometry/multilineargeometry.hh>
 #include <dune/geometry/quadraturerules.hh>
 #include <dune/geometry/referenceelements.hh>
 #include <dune/geometry/type.hh>
 
+#if HAVE_DUNE_LOCALFUNCTIONS
 #include <dune/localfunctions/lagrange/lagrangelfecache.hh>
+#endif
 
 namespace Dune
 {
@@ -59,7 +62,7 @@ namespace Dune
    *  This structure provides the default values.
    *
    *  \tparam  ct  coordinate type
-   *  \tparam  LFECache  A LocalFiniteElementVariantCache implementation, 
+   *  \tparam  LFECache  A LocalFiniteElementVariantCache implementation,
    *                     e.g. LagrangeLocalFiniteElementCache
    */
   template <class ct, class LFECache>
@@ -82,13 +85,16 @@ namespace Dune
   // CurvedGeometry
   // -------------------
 
-  /// \brief Curved geometry implementation based on lagrange basis function parametrization
+  /// \brief Curved geometry implementation based on local basis function parametrization
   /**
+   *  Parametrization of the geometry by any localfunction given by a finite-element cache
+   *  defined in the Traits class. See \ref LagrangeCurvedGeometry for an example instantiation
+   *  of this class with local lagrange basis functions.
    *
    *  \tparam  ct      coordinate type
    *  \tparam  mydim   geometry dimension
    *  \tparam  cdim    coordinate dimension
-   *  \tparam  TraitsType  Parametrization of the geometry, see \ref CurvedGeometryTraits
+   *  \tparam  TraitsType  Parameters of the geometry, see \ref CurvedGeometryTraits
    *
    *  The requirements on the traits are documented along with their default,
    *  CurvedGeometryTraits.
@@ -102,13 +108,16 @@ namespace Dune
 
     /// geometry dimension
     static const int mydimension = mydim;
+
     /// coordinate dimension
     static const int coorddimension = cdim;
 
     /// type of local coordinates
     using LocalCoordinate = FieldVector<ctype, mydimension>;
+
     /// type of global coordinates
     using GlobalCoordinate = FieldVector<ctype, coorddimension>;
+
     /// type of volume
     using Volume = ctype;
 
@@ -118,11 +127,9 @@ namespace Dune
     /// type of jacobian inverse transposed
     using JacobianInverseTransposed = FieldMatrix<ctype, coorddimension, mydimension>;
 
-  protected:
-    using ReferenceElements = Dune::ReferenceElements<ctype, mydimension>;
-
   public:
     /// type of reference element
+    using ReferenceElements = Dune::ReferenceElements<ctype, mydimension>;
     using ReferenceElement = typename ReferenceElements::ReferenceElement;
 
     /// Parametrization of the geometry
@@ -143,16 +150,12 @@ namespace Dune
       , localFE_(localFECache_.get(refElement.type()))
     {}
 
-
   public:
-    /// \brief constructor from a vector of coefficients of the LocalBasis parametrizing the geometry
+    /// \brief constructor from a vector of coefficients of the LocalBasis parametrizing
+    /// the geometry
     /**
      *  \param[in]  refElement  reference element for the geometry
      *  \param[in]  vertices    vertices to store internally
-     *
-     *  \note The type of vertices is actually a template argument.
-     *        It is only required that the internal vertex storage can be
-     *        constructed from this object.
      */
     CurvedGeometry (const ReferenceElement& refElement, std::vector<GlobalCoordinate> vertices)
       : CurvedGeometry(refElement)
@@ -161,7 +164,13 @@ namespace Dune
       assert(localFE_.size() == vertices_.size());
     }
 
-    /// \brief constructor from a local parametrization function, mapping local to (curved) global coordinates
+    /// \brief constructor from a local parametrization function, mapping local to (curved)
+    /// global coordinates
+    /**
+     *  \param[in]  refElement  reference element for the geometry
+     *  \param[in]  param       parametrization function with signature
+     *                          `GlobalCoordiante(LocalCoordinate)`
+     */
     template <class Parametrization,
       std::enable_if_t<Std::is_callable<Parametrization(LocalCoordinate), GlobalCoordinate>::value, bool> = true>
     CurvedGeometry (const ReferenceElement& refElement, Parametrization&& param)
@@ -174,7 +183,8 @@ namespace Dune
     /// \brief constructor, forwarding to the other constructor that take a reference-element
     /**
      *  \param[in]  gt     geometry type
-     *  \param[in]  param  either a vector of vertices, or a functor that can be used to construct the vertices
+     *  \param[in]  param  either a vector of vertices, or a functor that can be used to
+     *                     construct the vertices
      */
     template <class Parametrization>
     CurvedGeometry (Dune::GeometryType gt, Parametrization&& param)
@@ -216,7 +226,9 @@ namespace Dune
     /// \brief is this mapping affine?
     bool affine () const
     {
-      return refElement_.type().isSimplex() && order() == 1;
+      if (!affine_)
+        affine_ = (order() == 1 && ( type().isSimplex() || flatGeometry().affine() ));
+      return *affine_;
     }
 
     /// \brief obtain the name of the reference element
@@ -239,7 +251,7 @@ namespace Dune
     GlobalCoordinate corner (int i) const
     {
       assert( (i >= 0) && (i < corners()) );
-      return vertices_[i];
+      return global(refElement_.position(i, mydimension));
     }
 
     /// \brief obtain the centroid of the mapping's image
@@ -287,7 +299,8 @@ namespace Dune
     LocalCoordinate local (const GlobalCoordinate& globalCoord) const
     {
       const ctype tolerance = Traits::tolerance();
-      LocalCoordinate x = refElement().position( 0, 0 );
+      LocalCoordinate x = flatGeometry().local(globalCoord);
+
       LocalCoordinate dx;
       const bool affineMapping = affine();
 
@@ -315,16 +328,19 @@ namespace Dune
       return x;
     }
 
-    /// \brief Construct a global coordinate normal of the curvilinear element evaluated at a given local coordinate
+    /// \brief Construct a global coordinate normal of the curvilinear 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
     {
-      if ((mydim == 1) && (cdim == 2)) { return normalEdge(local, jacobianTransposed(local)); }
-      if ((mydim == 2) && (cdim == 3)) { return normalTriangle(local, jacobianTransposed(local)); }
+      if ((mydim == 1) && (cdim == 2)) { return normal1D(local, jacobianTransposed(local)); }
+      if ((mydim == 2) && (cdim == 3)) { return normal2D(local, jacobianTransposed(local)); }
+
+      DUNE_THROW(Dune::NotImplemented,
+        "ERROR: normal() method only defined for edges in 2D and faces in 3D");
 
-      DUNE_THROW(Dune::NotImplemented, "ERROR: normal() method only defined for edges in 2D and faces in 3D");
       return GlobalCoordinate(0);
     }
 
@@ -337,10 +353,6 @@ namespace Dune
      *  \param[in]  local  local coordinate to evaluate the integration element in
      *
      *  \returns the integration element \f$\mu(x)\f$.
-     *
-     *  \note For affine mappings, it is more efficient to call
-     *        jacobianInverseTransposed before integrationElement, if both
-     *        are required.
      */
     ctype integrationElement (const LocalCoordinate& local) const
     {
@@ -367,9 +379,6 @@ namespace Dune
      *  \param[in]  local  local coordinate to evaluate Jacobian in
      *
      *  \returns a reference to the transposed of the Jacobian
-     *
-     *  \note The returned reference is reused on the next call to
-     *        JacobianTransposed, destroying the previous value.
      */
     JacobianTransposed jacobianTransposed (const LocalCoordinate& local) const
     {
@@ -415,12 +424,12 @@ namespace Dune
       return refElement_;
     }
 
-    const LocalBasis& localBasis() const
+    const LocalBasis& localBasis () const
     {
       return localFE_.localBasis();
     }
 
-    GlobalCoordinate normalEdge (const LocalCoordinate& local, const JacobianTransposed& J) const
+    GlobalCoordinate normal1D (const LocalCoordinate& local, const JacobianTransposed& J) const
     {
       GlobalCoordinate res{
          J[0][1],
@@ -429,7 +438,7 @@ namespace Dune
       return res /= res.two_norm();
     }
 
-    GlobalCoordinate normalTriangle (const LocalCoordinate& local, const JacobianTransposed& J) const
+    GlobalCoordinate normal2D (const LocalCoordinate& local, const JacobianTransposed& J) const
     {
       GlobalCoordinate res{
         J[0][1] * J[1][2] - J[0][2] * J[1][1],
@@ -439,20 +448,50 @@ namespace Dune
       return res /= res.two_norm();
     }
 
+    using FlatGeometry = MultiLinearGeometry<ctype, mydim, cdim>;
+    FlatGeometry const& flatGeometry () const
+    {
+      if (!flatGeometry_) {
+        std::vector<GlobalCoordinate> corners;
+        corners.reserve(refElement_.size(mydimension));
+        for (int i = 0; i < refElement_.size(mydimension); ++i)
+          corners.push_back(global(refElement_.position(i, mydimension)));
+
+        flatGeometry_ = FlatGeometry{refElement_, corners};
+      }
+
+      return *flatGeometry_;
+    }
+
   private:
     ReferenceElement refElement_;
     LocalFECache localFECache_;
     LocalFiniteElement const& localFE_;
 
     std::vector<GlobalCoordinate> vertices_;
+
+    // some data optionally provided
+    Std::optional<bool> affine_;
+    Std::optional<FlatGeometry> flatGeometry_;
   };
 
+
+#if HAVE_DUNE_LOCALFUNCTIONS
+
   /// \brief curved geometry parametrized with lagrange basis functions
+  /**
+   *  \tparam  ct      coordinate type
+   *  \tparam  mydim   geometry dimension
+   *  \tparam  cdim    coordinate dimension
+   *  \tparam  order   polynomial order of the lagrange functions for the geometry
+   *                   parametrization
+   **/
   template <class ctype, int mydim, int cdim, int order>
   using LagrangeCurvedGeometry = CurvedGeometry<ctype,mydim,cdim,
     CurvedGeometryTraits<ctype, LagrangeLocalFiniteElementCache<ctype, ctype, mydim, order>> >;
 
-} // namespace Dune
+#endif
 
+} // namespace Dune
 
 #endif // DUNE_CURVEDGEOMETRY_HH

dune/curvedgeometry/curvedgeometrydatacollector.hh

@@ -14,8 +14,15 @@
 
 namespace Dune
 {
-  // Extract point coordinates from parametrized geometry
-  // and write nonlinear vtk elements
+  /// \brief DataCollector for dune-vtk writers extracting the parametrized geometry from
+  /// a local-to-global coordinate mapping and writing quadratic VTK elements
+  /**
+   * To write curved elements, a \ref CurvedGeometry is constructed on-the-fly from the given
+   * mappings in the constructor.
+   *
+   * \tparam GridView  The grid-view to write
+   * \tparam Geometry  The CurvedGeometry to use for element parametrization
+   **/
   template <class GridView, class Geometry>
   class CurvedGeometryDataCollector
       : public UnstructuredDataCollectorInterface<GridView, CurvedGeometryDataCollector<GridView,Geometry>, Partitions::All>
@@ -34,21 +41,33 @@ namespace Dune
     using Super::partition;
 
   public:
-    // Construct the curved geometry from a coordinate mapping
+    /// \brief Construct the curved geometry from a coordinate mapping
+    /**
+     * \param gridView  An instance of the GridView type
+     * \param mapping   A coordinate mapping from global flat to global curved coordinates, i.e.
+     *                  with signature `typename Geometry::GlobalCoordinate(GlobalCoord)`,
+     *                  see \ref Mapping
+     **/
     CurvedGeometryDataCollector (GridView const& gridView, Mapping mapping)
       : Super(gridView)
       , dataCollector_(gridView)
       , mapping_(mapping)
     {}
 
-    // construct the curved geometry from a stored parametrization, i.e. a vector of coefficient vectors
+    /// \brief construct the curved geometry from a stored parametrization, i.e. a vector of
+    /// coefficient vectors
+    /**
+     * \param gridView         An instance of the GridView type
+     * \param parametrization  Vector (accessed by element index) of element parametrization
+     *                         coefficients, see \ref Parametrization
+     **/
     CurvedGeometryDataCollector (GridView const& gridView, Parametrization const& parametrization)
       : Super(gridView)
       , dataCollector_(gridView)
       , parametrization_(&parametrization)
     {}
 
-    // collect the points of the geometry. This uses the curved geometry instead of the element geometry
+    /// collect the points of the geometry. This uses the curved geometry instead of the element geometry
     template <class T>
     std::vector<T> pointsImpl () const
     {
@@ -82,7 +101,7 @@ namespace Dune
 
   private:
 
-    // construct the geometry object, either from the mapping or the coefficient vector
+    /// construct the geometry object, either from the mapping or the coefficient vector
     template <class IndexSet, class Element>
     Geometry makeGeometry(IndexSet const& indexSet, Element const& element) const
     {

dune/curvedgeometry/test/CMakeLists.txt

+dune_add_test(SOURCES curvedgeometry.cc)
+dune_add_test(SOURCES parametrization.cc)

src/CMakeLists.txt

-add_executable("curvedgeometry" curvedgeometry.cc)
-target_link_dune_default_libraries("curvedgeometry")
-
-add_executable("parametrization" parametrization.cc)
-target_link_dune_default_libraries("parametrization")
-
 add_executable("lagrangesimplex" lagrangesimplex.cc)
 target_link_dune_default_libraries("lagrangesimplex")