added two gridfunctions for normals and coords, changed the storage type for...

added two gridfunctions for normals and coords, changed the storage type for gridfunctions in the grid to shared_ptr

dune/curvedsurfacegrid/CMakeLists.txt

-add_subdirectory(surfacedistance)
-add_subdirectory(test)
 
-set(HEADERS
+install(FILES 
   backuprestore.hh
   capabilities.hh
   coordfunctionimplementations.hh
@@ -20,7 +18,9 @@ set(HEADERS
   intersection.hh
   intersectioniterator.hh
   iterator.hh
-  )
-
-install(FILES ${HEADERS}
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/curvedsurfacegrid)
+
+add_subdirectory(gridfunctions)
+add_subdirectory(surfacedistance)
+add_subdirectory(test)
+  
\ No newline at end of file

dune/curvedsurfacegrid/capabilities.hh

@@ -11,7 +11,7 @@
 #include <dune/grid/common/capabilities.hh>
 #include <dune/grid/geometrygrid/capabilities.hh>
 #include <dune/curvedsurfacegrid/declaration.hh>
-#include <dune/curvedsurfacegrid/gridfunction/gridfunction.hh>
+#include <dune/curvedsurfacegrid/gridfunctions/gridfunction.hh>
 
 namespace Dune
 {

dune/curvedsurfacegrid/grid.hh

@@ -4,13 +4,14 @@
 #define DUNE_CURVED_SURFACE_GRID_GRID_HH
 
 #include <dune/common/deprecated.hh>
+#include <dune/common/shared_ptr.hh>
 
 #include <dune/grid/common/grid.hh>
 
 #include <dune/curvedsurfacegrid/gridfamily.hh>
 #include <dune/curvedsurfacegrid/backuprestore.hh>
 #include <dune/curvedsurfacegrid/datahandle.hh>
-#include <dune/curvedsurfacegrid/gridfunction/gridfunction.hh>
+#include <dune/curvedsurfacegrid/gridfunctions/gridfunction.hh>
 #include <dune/grid/geometrygrid/identity.hh>
 #include <dune/grid/geometrygrid/persistentcontainer.hh>
 #include <dune/grid/geometrygrid/grid.hh>
@@ -187,13 +188,16 @@ namespace Dune
      *  wrapper is destroyed.
      *
      *  \param[in]  hostGrid      reference to the grid to wrap
-     *  \param[in]  param      mapping from global coordinates in the host geometry
+     *  \param[in]  gridFunction  mapping from global coordinates in the host geometry
      *                            to global coordinates in the curved geometry
-     *  \param[in]  allocator  storage allocator
+     * 
+     *  If the gridFunction is passed by (non-const) reference it is stored as a non-destroying shared_ptr.
+     *  Otherwise it is copied or moved into a new object (stored as shared_ptr as well).
      */
-    CurvedSurfaceGrid (HostGrid& hostGrid, const GridFunction& gridFunction)
+    template <class GF_>
+    CurvedSurfaceGrid (HostGrid& hostGrid, GF_&& gridFunction)
       : hostGrid_(hostGrid)
-      , gridFunction_(gridFunction)
+      , gridFunction_(wrap_or_move(std::forward<GF_>(gridFunction)))
       , levelIndexSets_(hostGrid_.maxLevel()+1, nullptr)
     {}
 
@@ -487,7 +491,15 @@ namespace Dune
     }
 
     //! obtain constant reference to the coordinate function
-    GridFunction const& gridFunction () const { return gridFunction_; }
+    GridFunction const& gridFunction () const 
+    { 
+      return *gridFunction_; 
+    }
+
+    GridFunction& gridFunction () 
+    { 
+      return *gridFunction_; 
+    }
 
     bool useGeometryCaching () const
     {
@@ -511,7 +523,7 @@ namespace Dune
 
   private:
     HostGrid& hostGrid_;
-    GridFunction gridFunction_;
+    std::shared_ptr<GridFunction> gridFunction_;
 
     bool useGeometryCaching_ = false;
 

dune/curvedsurfacegrid/gridfamily.hh

@@ -16,7 +16,7 @@
 #include <dune/curvedsurfacegrid/iterator.hh>
 #include <dune/curvedsurfacegrid/idset.hh>
 #include <dune/curvedsurfacegrid/indexsets.hh>
-#include <dune/curvedsurfacegrid/gridfunction/gridfunction.hh>
+#include <dune/curvedsurfacegrid/gridfunctions/gridfunction.hh>
 
 #include <dune/functions/common/functionconcepts.hh>
 #include <dune/functions/common/signature.hh>

dune/curvedsurfacegrid/gridfunctions/CMakeLists.txt

+install(FILES 
+  analyticgridfunction.hh
+  discretegridviewfunction.hh
+  gridentityset.hh
+  gridfunction.hh
+  spheregridfunction.hh
+  torusgridfunction.hh
+  DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/curvedsurfacegrid/gridfunctions)

dune/curvedsurfacegrid/gridfunction/analyticgridfunction.hh → dune/curvedsurfacegrid/gridfunctions/analyticgridfunction.hh

@@ -3,33 +3,16 @@
 #ifndef DUNE_CURVED_SURFACE_GRID_ANALYTIC_GRIDFUNCTION_HH
 #define DUNE_CURVED_SURFACE_GRID_ANALYTIC_GRIDFUNCTION_HH
 
-namespace Dune
-{
-  //! A set of entities of given `codim` of a `Grid`
-  /**
-   * \tparam G      The grid type
-   * \tparam codim  Codimension of the entities to define the set of.
-   *
-   * \note This entityset just defines types
-   **/
-  template< class G, int codim >
-  class GridEntitySet
-  {
-  public:
-    //! Type of the grid
-    using Grid = G;
+#include <type_traits>
+#include <utility>
 
-    //! Type of Elements contained in this EntitySet
-    using Element = typename Grid::template Codim<codim>::Entity;
-
-    //! Type of local coordinates with respect to the Element
-    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
-
-    //! Type of global coordinates with respect to the Element
-    using GlobalCoordinate = typename Element::Geometry::GlobalCoordinate;
-  };
+#include <dune/common/typeutilities.hh>
+#include <dune/common/std/optional.hh>
 
+#include "gridentityset.hh"
 
+namespace Dune
+{
   //! LocalFunction associated to the \ref AnalyticGridFunction
   /**
    * \tparam LocalContext  Context this localFunction can be bound to, e.g. a grid element

dune/curvedsurfacegrid/gridfunctions/discretegridviewfunction.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_CURVED_SURFACE_GRID_DISCRETE_GRIDVIEWFUNCTION_HH
+#define DUNE_CURVED_SURFACE_GRID_DISCRETE_GRIDVIEWFUNCTION_HH
+
+#include <array>
+#include <vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+#include <dune/functions/common/defaultderivativetraits.hh>
+#include <dune/functions/functionspacebases/basistags.hh>
+#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+#include <dune/functions/gridfunctions/localderivativetraits.hh>
+#include <dune/grid/utility/hierarchicsearch.hh>
+#include <dune/istl/bvector.hh>
+
+
+namespace Dune
+{
+  namespace Impl
+  {
+    template <class Sig, int degree, template <class> class Traits>
+    struct DerivativeRangeType;
+
+    template <class R, class D, int degree, template <class> class DerivativeTraits>
+    struct DerivativeRangeType<R(D), degree, DerivativeTraits>
+    {
+      using DerivativeRange = typename DerivativeTraits<R(D)>::Range;
+      using type = typename DerivativeRangeType<DerivativeRange(D), degree-1, DerivativeTraits>::type;
+    };
+
+    template <class R, class D, template <class> class Traits>
+    struct DerivativeRangeType<R(D), 0, Traits>
+    {
+      using type = R;
+    };
+  }
+
+  //! Grid-view function representing a coordinate map from a reference grid to world
+  /**
+   * \tparam GridView   The grid-view type this coordinate function is defined on
+   * \tparam dow        Dimension of the world this coordfunction mapps into
+   * \tparam ORDER      Polynomial order of the local parametrization
+   * \tparam T          Type of the basis range-type and coefficient type
+   **/
+  template< class GridView, int dow = GridView::dimensionworld, int ORDER = -1, class T = double >
+  class DiscreteGridViewFunction;
+
+  //! Generator for \ref DiscreteGridViewFunction
+  /**
+   * \param gridView   Function that can be evaluated at global coordinates of the \ref Grid
+   * \param order      Polynomial order of the local parametrization [ORDER]
+   * 
+   * \tparam ORDER  Polynomial order of the local parametrization [-1]
+   * \tparam T      Range-type of the basis and coefficient value-type [double]
+   **/
+  template< int dow, int ORDER = -1, class T = double, class GridView >
+  auto discreteGridViewFunction (const GridView& gridView, int order = ORDER)
+  {
+    return DiscreteGridViewFunction<GridView,dow,ORDER,T>{gridView, order};
+  }
+
+
+  template< class GridView, int dow, int ORDER, class T >
+  class DiscreteGridViewFunction
+  {
+    static auto makeBasis (const GridView& gridView, int order)
+    {
+      namespace BF = Functions::BasisFactory;
+      return BF::makeBasis(gridView, BF::power<dow>(BF::lagrange<T>(order), BF::blockedInterleaved()));
+    }
+
+    using Basis = decltype(makeBasis(std::declval<GridView>(), ORDER));
+
+  public:
+    using EntitySet = Functions::GridViewEntitySet<GridView,0>;
+
+    using Domain = typename EntitySet::GlobalCoordinate;
+    using Range = FieldVector<T,dow>;
+    using Signature = Range(Domain);
+
+  private:
+    using VectorType = BlockVector<Range>;
+
+    template< int derivativeOrder = 0 >
+    class LocalFunction
+    {
+      using LocalView = typename Basis::LocalView;
+      using LocalFiniteElement = typename LocalView::Tree::ChildType::FiniteElement;
+      using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
+
+      template <class Sig>
+      using DerivativeTraits = typename Functions::LocalDerivativeTraits<EntitySet, Functions::DefaultDerivativeTraits>::template Traits<Sig>;
+
+    public:
+      using LocalContext = typename LocalView::Element;
+      using Geometry = typename LocalContext::Geometry;
+
+      using RangeType = typename DiscreteGridViewFunction::Range;
+
+      using Domain = typename EntitySet::LocalCoordinate;
+
+      template <int degree>
+      using DerivativeRange = typename Impl::DerivativeRangeType<RangeType(Domain), degree, DerivativeTraits>::type;
+
+      using Range = DerivativeRange<derivativeOrder>;
+      using Signature = Range(Domain);
+
+    public:
+      template <class LV>
+      LocalFunction (LV&& localView, const VectorType& coords)
+        : localView_(std::forward<LV>(localView))
+        , coords_(coords)
+      {}
+
+      //! Collect the coords from all element DOFs into a local
+      //! vector that can be accessed in the operator() for interpolation
+      void bind (const LocalContext& element)
+      {
+        localView_.bind(element);
+
+        const auto& leafNode = localView_.tree().child(0);
+        localCoords_.resize(leafNode.size());
+
+        // collect local coordinate vectors
+        for (std::size_t i = 0; i < localCoords_.size(); ++i) {
+          auto idx = localView_.index(leafNode.localIndex(i));
+          localCoords_[i] = coords_[idx[0]];
+        }
+
+        if constexpr (derivativeOrder == 1)
+          geometry_.emplace(element.geometry());
+
+        bound_ = true;
+      }
+
+      void unbind ()
+      {
+        localView_.unbind();
+        bound_ = false;
+      }
+
+      LocalContext const& localContext () const
+      {
+        assert(bound_);
+        return localView_.element();
+      }
+
+      //! Evaluate coordinates in local coordinates
+      //! by interpolation of stored coords in \ref localCoords_.
+      Range operator() (const Domain& local) const 
+      {
+        static_assert(derivativeOrder < 2, "Higher-order derivatives not implemented");
+
+        if constexpr (derivativeOrder == 0)
+          return evaluateFunction(local);
+        else if constexpr (derivativeOrder == 1)
+          return evaluateJacobian(local);
+
+        return Range(0);
+      }
+
+      friend LocalFunction<derivativeOrder+1> derivative (LocalFunction const& lf)
+      {
+        return LocalFunction<derivativeOrder+1>{lf.localView_, lf.coords_};
+      }
+
+    private:
+      DerivativeRange<0> evaluateFunction (const Domain& local) const 
+      {
+        assert(bound_);
+
+        const auto& leafNode = localView_.tree().child(0);
+        const auto& lfe = leafNode.finiteElement();
+
+        // evaluate basis functions in local coordinate
+        lfe.localBasis().evaluateFunction(local, shapeValues_);
+        assert(localCoords_.size() == shapeValues_.size());
+
+        DerivativeRange<0> x(0);
+        for (std::size_t i = 0; i < localCoords_.size(); ++i)
+          x.axpy(shapeValues_[i], localCoords_[i]);
+
+        return x;
+      }
+
+      DerivativeRange<1> evaluateJacobian (const Domain& local) const 
+      {
+        assert(bound_);
+
+        const auto& leafNode = localView_.tree().child(0);
+        const auto& lfe = leafNode.finiteElement();
+
+        // evaluate basis functions in local coordinate
+        lfe.localBasis().evaluateJacobian(local, shapeGradients_);
+        assert(localCoords_.size() == shapeGradients_.size());
+
+        // transform gradients to global coordinates
+        auto jit = geometry_->jacobianInverseTransposed(local);
+        gradients_.resize(shapeGradients_.size());
+        for (std::size_t i = 0; i < shapeGradients_.size(); ++i)
+          jit.mv(shapeGradients_[i][0], gradients_[i]);
+
+        DerivativeRange<1> J(0);
+        for (std::size_t i = 0; i < localCoords_.size(); ++i)
+          for (int j = 0; j < J.N(); ++j)
+            J[j].axpy(localCoords_[i][j], gradients_[i]);
+
+        return J;
+      }
+
+    private:
+      LocalView localView_;
+      const VectorType& coords_;
+
+      std::vector<RangeType> localCoords_;
+      Std::optional<Geometry> geometry_;
+      mutable std::vector<typename LocalBasis::Traits::RangeType> shapeValues_;
+      mutable std::vector<typename LocalBasis::Traits::JacobianType> shapeGradients_;
+      mutable std::vector<DerivativeRange<0>> gradients_;
+
+      bool bound_ = false;
+    };
+
+  public:
+    //! Constructor.
+    DiscreteGridViewFunction (const GridView& gridView, int order = (ORDER > 0 ? ORDER : 1))
+      : entitySet_(gridView)
+      , basis_(makeBasis(gridView, order))
+    {
+      update(gridView);
+    }
+
+    void update (const GridView& gridView)
+    {
+      entitySet_ = EntitySet{gridView};
+      basis_.update(gridView);
+
+      coords_.resize(basis_.size());
+    }
+
+    //! evaluate in global coordinates
+    Range operator() (const Domain& x) const
+    {
+      using Grid = typename GridView::Grid;
+      using IS = typename GridView::IndexSet;
+
+      const auto& gv = entitySet_.gridView();
+      HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};
+
+      auto element = hsearch.findEntity(x);
+      auto geometry = element.geometry();
+      auto localFct = localFunction(*this);
+      localFct.bind(element);
+      return localFct(geometry.local(x));
+    }
+
+    //! Create a local function of this grifunction
+    friend LocalFunction<0> localFunction (const DiscreteGridViewFunction& gf)
+    {
+      return LocalFunction<0>{gf.basis_.localView(), gf.coords_};
+    }
+
+    //! obtain the stored \ref GridViewEntitySet
+    const EntitySet& entitySet () const
+    {
+      return entitySet_;
+    }
+
+    const Basis& basis () const 
+    {
+      return basis_;
+    }
+
+    const VectorType& coefficients () const 
+    {
+      return coords_;
+    }
+
+    VectorType& coefficients ()
+    {
+      return coords_;
+    }
+
+  private:
+    EntitySet entitySet_;
+    Basis basis_;
+    VectorType coords_;
+  };
+
+} // end namespace Dune
+
+#endif // DUNE_CURVED_SURFACE_GRID_DISCRETE_GRIDVIEWFUNCTION_HH

dune/curvedsurfacegrid/gridfunctions/gridentityset.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_CURVED_SURFACE_GRID_GRID_ENTITYSET_HH
+#define DUNE_CURVED_SURFACE_GRID_GRID_ENTITYSET_HH
+
+namespace Dune
+{
+  //! A set of entities of given `codim` of a `Grid`
+  /**
+   * \tparam G      The grid type
+   * \tparam codim  Codimension of the entities to define the set of.
+   *
+   * \note This entityset just defines types
+   **/
+  template< class G, int codim >
+  class GridEntitySet
+  {
+  public:
+    //! Type of the grid
+    using Grid = G;
+
+    //! Type of Elements contained in this EntitySet
+    using Element = typename Grid::template Codim<codim>::Entity;
+
+    //! Type of local coordinates with respect to the Element
+    using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
+
+    //! Type of global coordinates with respect to the Element
+    using GlobalCoordinate = typename Element::Geometry::GlobalCoordinate;
+  };
+
+} // end namespace Dune 
+
+#endif // DUNE_CURVED_SURFACE_GRID_GRID_ENTITYSET_HH

dune/curvedsurfacegrid/gridfunction/gridfunction.hh → dune/curvedsurfacegrid/gridfunctions/gridfunction.hh

@@ -9,25 +9,25 @@ namespace Dune
 {
   namespace Impl
   {
-    template< class GF, class = void >
+    template< class ES, class = void >
     struct GridOf;
 
-    template< class GF >
-    struct GridOf<GF, std::void_t<typename GF::GridView>>
+    template< class ES >
+    struct GridOf<ES, std::void_t<typename ES::GridView>>
     {
-      using type = typename GF::GridView::Grid;
+      using type = typename ES::GridView::Grid;
     };
 
-    template< class GF >
-    struct GridOf<GF, std::void_t<typename GF::Grid>>
+    template< class ES >
+    struct GridOf<ES, std::void_t<typename ES::Grid>>
     {
-      using type = typename GF::Grid;
+      using type = typename ES::Grid;
     };
 
   } // end namespace Impl
 
   template< class GF >
-  using GridOf_t = typename Impl::GridOf<GF>::type;
+  using GridOf_t = typename Impl::GridOf<typename GF::EntitySet>::type;
 
 } // end namespace Dune
 

dune/curvedsurfacegrid/gridfunctions/normalgridviewfunction.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_CURVED_SURFACE_GRID_NORMAL_GRIDVIEWFUNCTION_HH
+#define DUNE_CURVED_SURFACE_GRID_NORMAL_GRIDVIEWFUNCTION_HH
+
+#include <array>
+#include <vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+#include <dune/functions/functionspacebases/basistags.hh>
+#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+#include <dune/grid/utility/hierarchicsearch.hh>
+#include <dune/istl/bvector.hh>
+
+
+namespace Dune 
+{
+  //! Grid-view function representing averaged normal vector
+  /**
+   * \tparam GridView   The grid-view this grid-view-function is defined on
+   * \tparam ORDER      Polynomial order of the lagrange bases used for representing the normals
+   * \tparam T          Value type used for the basis and the coefficients
+   **/
+  template< class GridView, int ORDER = -1, class T = double >
+  class NormalGridViewFunction
+  {
+    static auto makeBasis (const GridView& gridView, int order)
+    {
+      namespace BF = BasisFactory;
+      return BF::makeBasis(gridView, BF::power<GridView::dimensionworld>(BF::lagrange<T>(order), BF::blockedInterleaved()));
+    }
+
+    using Basis = decltype(makeBasis(std::declval<GridView>(), ORDER));
+
+  public:
+    using EntitySet = GridViewEntitySet<GridView,0>;
+
+    using Domain = typename EntitySet::GlobalCoordinate;
+    using Range = FieldVector<T,GridView::dimensionworld>;
+
+    using VectorType = BlockVector<Range>;
+
+  private:
+    class LocalFunction
+    {
+      using LocalView = typename Basis::LocalView;
+      using LocalContext = typename LocalView::Element;
+
+      using Domain = typename EntitySet::LocalCoordinate;
+      using Range = typename NormalGridViewFunction::Range;
+
+    public:
+      LocalFunction (LocalView&& localView, const VectorType& normals)
+        : localView_(std::move(localView))
+        , normals_(normals)
+      {}
+
+      //! Collect the normal vector from all element DOFs into a local
+      //! vector that can be accessed in the operator() for interpolation
+      void bind (const LocalContext& element)
+      {
+        localView_.bind(element);
+