merged remote branch

dune/curvedsurfacegrid/CMakeLists.txt

-add_subdirectory(surfacedistance)
-add_subdirectory(test)
 
-set(HEADERS
+install(FILES 
   backuprestore.hh
   capabilities.hh
+  concepts.hh
   coordfunctionimplementations.hh
   coordprovider.hh
   curvedsurfacegrid.hh
@@ -20,7 +19,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/concepts.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_CONCEPTS_HH
+#define DUNE_CURVED_SURFACE_GRID_CONCEPTS_HH
+
+#include <dune/common/concept.hh>
+#include <dune/curvedsurfacegrid/gridfunctions/gridentityset.hh>
+#include <dune/functions/common/functionconcepts.hh>
+
+namespace Dune {
+namespace Concept {
+
+template< class LocalContext >
+struct LocalFunction
+{
+  using LocalCoordinate = typename LocalContext::Geometry::LocalCoordinate;
+
+  template< class LF >
+  auto require(LF&& lf) -> decltype(
+    lf.bind(std::declval<LocalContext>()),
+    lf.unbind(),
+    lf.localContext(),
+    requireConcept<Dune::Functions::Concept::Callable<LocalCoordinate>>(lf),
+    requireConvertible<LocalContext>(lf.localContext())
+  );
+};
+
+template< class LF, class LocalContext >
+constexpr bool isLocalFunction()
+{ return models<Concept::LocalFunction<LocalContext>, LF>(); }
+
+
+template< class HostGrid >
+struct GridFunction
+{
+  using EntitySet = GridEntitySet<HostGrid,0>;
+  using LocalContext = typename EntitySet::Element;
+  using GlobalCoordinate = typename EntitySet::GlobalCoordinate;
+
+  template< class GF >
+  auto require(GF&& gf) -> decltype(
+    localFunction(gf),
+    gf.entitySet(),
+    requireConcept<Dune::Functions::Concept::Callable<GlobalCoordinate>>(gf),
+    requireConcept<LocalFunction<LocalContext>>(localFunction(gf))
+  );
+};
+
+template< class GF, class HostGrid >
+constexpr bool isGridFunction()
+{ return models<Concept::GridFunction<HostGrid>, GF>(); }
+
+} // end namespace Concept
+} // end namespace Dune 
+
+#endif // DUNE_CURVED_SURFACE_GRID_CONCEPTS_HH
\ No newline at end of file

dune/curvedsurfacegrid/geometry.hh

@@ -576,7 +576,7 @@ namespace Dune
 
     protected:
       VertexGeometry (Construct, const ReferenceElement& refElement, const LocalFunction& localFunction,
-                const LocalGeometry& lg)
+                      const LocalGeometry& lg)
         : Super(refElement, std::vector<GlobalCoordinate>{localFunction(lg.global(refElement.position(0,0)))})
       {}
 

dune/curvedsurfacegrid/grid.hh

@@ -4,13 +4,17 @@
 #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/concepts.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/analyticgridfunction.hh>
+#include <dune/curvedsurfacegrid/gridfunctions/gridfunction.hh>
+#include <dune/functions/common/functionconcepts.hh>
 #include <dune/grid/geometrygrid/identity.hh>
 #include <dune/grid/geometrygrid/persistentcontainer.hh>
 #include <dune/grid/geometrygrid/grid.hh>
@@ -46,18 +50,42 @@ namespace Dune
    *  of a CurvedGeometry.
    *
    *  \tparam GF     GridViewFunction defined on a (flat) hostgrid
-   *  \tparam order  Polynomial order of local lagrange basis functions to use when approximating
+   *  \tparam ORDER  Polynomial order of local lagrange basis functions to use when approximating
    *                 the curved geometry. [optional]
    *
    *  \nosubgrouping
    */
-  template< class GF, int order = -1 >
+  template< class GF, int ORDER = -1 >
+  class CurvedSurfaceGrid;
+
+  //! Generator for CurvedSurfaceGrid from a grid-functions
+  template< class HostGrid, class GF, int ORDER = -1, 
+    std::enable_if_t<Dune::Concept::isGridFunction<GF, HostGrid>(), int> = 0>
+  auto curedSurfaceGrid (HostGrid& hostGrid, GF&& gridFunction) 
+  {
+    static_assert(std::is_same<HostGrid, GridOf_t<std::decay_t<GF>>>::value, "GridFunction must be defined on the HostGrid");
+    return CurvedSurfaceGrid<std::decay_t<GF>,ORDER>{hostGrid, std::forward<GF>(gridFunction)};
+  }
+
+  //! Generator for CurvedSurfaceGrid from a callable
+  template< class HostGrid, class F, int ORDER = -1, 
+    std::enable_if_t<not Dune::Concept::isGridFunction<F, HostGrid>(), int> = 0>
+  auto curedSurfaceGrid (HostGrid& hostGrid, F&& callable) 
+  {
+    using GlobalCoordinate = typename GridEntitySet<HostGrid,0>::GlobalCoordinate;
+    static_assert(Functions::Concept::isCallable<F, GlobalCoordinate>(), "Function must be callable");
+    auto gridFct = analyticGridFunction<HostGrid>(std::forward<F>(callable));
+    return CurvedSurfaceGrid<decltype(gridFct),ORDER>{hostGrid, std::move(gridFct)};
+  }
+
+
+  template< class GF, int ORDER >
   class CurvedSurfaceGrid
-      : public CurvedSurfaceGridBase<GF,order>
-      , public CGeo::BackupRestoreFacilities<CurvedSurfaceGrid<GF,order> >
+      : public CurvedSurfaceGridBase<GF,ORDER>
+      , public CGeo::BackupRestoreFacilities<CurvedSurfaceGrid<GF,ORDER> >
   {
     using Self = CurvedSurfaceGrid;
-    using Super = CurvedSurfaceGridBase<GF,order>;
+    using Super = CurvedSurfaceGridBase<GF,ORDER>;
 
     // friend declarations
     template< class, class > friend class CGeo::IdSet;
@@ -65,7 +93,7 @@ namespace Dune
 
   public:
     using GridFunction = GF;
-    using GridFamily = CGeo::GridFamily<GF,order>;
+    using GridFamily = CGeo::GridFamily<GF,ORDER>;
 
     /** \name Traits
      *  \{ */
@@ -186,17 +214,26 @@ namespace Dune
      *  The references to host grid is stored in the grid. It must be valid until this grid
      *  wrapper is destroyed.
      *
-     *  \param[in]  hostGrid   reference to the grid to wrap
-     *  \param[in]  param      mapping from global coordinates in the host geometry
-     *                         to global coordinates in the curved geometry
-     *  \param[in]  allocator  storage allocator
+     *  \param[in]  hostGrid      reference to the grid to wrap
+     *  \param[in]  gridFunction  mapping from global coordinates in the host geometry
+     *                            to global coordinates in the curved geometry
+     * 
+     *  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_,
+      std::enable_if_t<std::is_same<GF,std::decay_t<GF_>>::value, int> = 0>
+    CurvedSurfaceGrid (HostGrid& hostGrid, GF_&& gridFunction, std::integral_constant<int,ORDER> = {})
       : hostGrid_(hostGrid)
-      , gridFunction_(gridFunction)
+      , gridFunction_(wrap_or_move(std::forward<GF_>(gridFunction)))
       , levelIndexSets_(hostGrid_.maxLevel()+1, nullptr)
     {}
 
+    template <class F_,
+      std::enable_if_t<std::is_same<GF,AnalyticGridFunction<HostGrid,std::decay_t<F_>>>::value, int> = 0>
+    CurvedSurfaceGrid (HostGrid& hostGrid, F_&& callable, std::integral_constant<int,ORDER> = {})
+      : CurvedSurfaceGrid(hostGrid, analyticGridFunction<HostGrid>(std::forward<F_>(callable)))
+    {}
 
     //! destructor
     ~CurvedSurfaceGrid ()
@@ -487,7 +524,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 +556,7 @@ namespace Dune
 
   private:
     HostGrid& hostGrid_;
-    GridFunction gridFunction_;
+    std::shared_ptr<GridFunction> gridFunction_;
 
     bool useGeometryCaching_ = false;
 
@@ -538,13 +583,18 @@ namespace Dune
     mutable typename GeometryCache<std::make_integer_sequence<int,Traits::dimension+1>>::type geometryCache_;
   };
 
+  //! Deduction guide for CurvedSurfaceGrid from a grid-functions or callables
+  template< class HostGrid, class GF, int ORDER = -1 >
+  CurvedSurfaceGrid (HostGrid& hostGrid, GF&& gridFunction, std::integral_constant<int,ORDER> = {}) 
+    -> CurvedSurfaceGrid<GridFunctionOf_t<HostGrid,std::decay_t<GF>>,ORDER>;
+
 
   // CurvedSurfaceGrid::Codim
   // ------------------------
 
-  template< class GridFunction, int order >
+  template< class GridFunction, int ORDER >
   template< int codim >
-  struct CurvedSurfaceGrid<GridFunction,order>::Codim
+  struct CurvedSurfaceGrid<GridFunction,ORDER>::Codim
     : public Super::template Codim< codim >
   {
     /** \name Entity types

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>
@@ -33,7 +33,7 @@ namespace Dune
     template< class GF >
     struct DimRange
     {
-      using SigTraits = Functions::SignatureTraits<typename GF::Signature>;
+      using SigTraits = Functions::SignatureTraits<GF>;
       using Range = typename SigTraits::RawRange;
       static const int value = Range::size();
     };
@@ -44,7 +44,7 @@ namespace Dune
       using EntitySet = typename GF::EntitySet;
       using LocalContext = typename EntitySet::Element;
 
-      using Range = typename Functions::SignatureTraits<typename GF::Signature>::Range;
+      using Range = typename Functions::SignatureTraits<GF>::Range;
       using LocalSignature = Range(typename EntitySet::LocalCoordinate);
 
       static const bool value = Functions::Concept::isDifferentiableLocalFunction<LF,LocalSignature,LocalContext,Functions::LocalDerivativeTraits<EntitySet>::template Traits>();
@@ -62,9 +62,10 @@ namespace Dune
         using HostGrid = GridOf_t<GF>;
 
         using GridFunction = GF;
-        using LocalFunction = std::decay_t<decltype(localFunction(std::declval<GF>()))>;
+        using LocalFunction = std::decay_t<decltype(localFunction(std::declval<GF const&>()))>;
 
         static const bool differentiableLocalFunction = DifferentiableLocalFunction<GF, LocalFunction>::value;
+        static_assert((differentiableLocalFunction || order>0), "Either provide a differentiable GridFunction or set ORDER > 0");
 
         using ctype = typename HostGrid::ctype;
 

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;