added gridfunction for functors, and specializations for the sphere and torus

dune/curvedsurfacegrid/capabilities.hh

@@ -11,6 +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>
 
 namespace Dune
 {
@@ -27,7 +28,7 @@ namespace Dune
     template< class GridFunction, int order >
     struct hasSingleGeometryType< CurvedSurfaceGrid<GridFunction,order> >
     {
-      using HostGrid = typename GridFunction::GridView::Grid;
+      using HostGrid = GridOf_t<GridFunction>;
       static const bool v = hasSingleGeometryType< HostGrid > :: v;
       static const unsigned int topologyId = hasSingleGeometryType< HostGrid > :: topologyId;
     };
@@ -50,7 +51,7 @@ namespace Dune
     template< class GridFunction, int order, int codim >
     struct canCommunicate< CurvedSurfaceGrid<GridFunction,order>, codim >
     {
-      using HostGrid = typename GridFunction::GridView::Grid;
+      using HostGrid = GridOf_t<GridFunction>;
       static const bool v = canCommunicate< HostGrid, codim >::v && hasEntity< HostGrid, codim >::v;
     };
 
@@ -58,21 +59,21 @@ namespace Dune
     template< class GridFunction, int order >
     struct hasBackupRestoreFacilities< CurvedSurfaceGrid<GridFunction,order> >
     {
-      using HostGrid = typename GridFunction::GridView::Grid;
+      using HostGrid = GridOf_t<GridFunction>;
       static const bool v = hasBackupRestoreFacilities< HostGrid >::v;
     };
 
     template< class GridFunction, int order >
     struct isLevelwiseConforming< CurvedSurfaceGrid<GridFunction,order> >
     {
-      using HostGrid = typename GridFunction::GridView::Grid;
+      using HostGrid = GridOf_t<GridFunction>;
       static const bool v = isLevelwiseConforming< HostGrid >::v;
     };
 
     template< class GridFunction, int order >
     struct isLeafwiseConforming< CurvedSurfaceGrid<GridFunction,order> >
     {
-      using HostGrid = typename GridFunction::GridView::Grid;
+      using HostGrid = GridOf_t<GridFunction>;
       static const bool v = isLeafwiseConforming< HostGrid >::v;
     };
 
@@ -96,7 +97,7 @@ namespace Dune
     template< class GridFunction, int order, int codim >
     struct hasHostEntity< CurvedSurfaceGrid<GridFunction,order>, codim >
     {
-      using HostGrid = typename GridFunction::GridView::Grid;
+      using HostGrid = GridOf_t<GridFunction>;
       static const bool v = hasEntity< HostGrid, codim >::v;
     };
 

dune/curvedsurfacegrid/grid.hh

@@ -10,6 +10,7 @@
 #include <dune/curvedsurfacegrid/gridfamily.hh>
 #include <dune/curvedsurfacegrid/backuprestore.hh>
 #include <dune/curvedsurfacegrid/datahandle.hh>
+#include <dune/curvedsurfacegrid/gridfunction/gridfunction.hh>
 #include <dune/grid/geometrygrid/identity.hh>
 #include <dune/grid/geometrygrid/persistentcontainer.hh>
 #include <dune/grid/geometrygrid/grid.hh>
@@ -21,7 +22,7 @@ namespace Dune
     template< class GridFunction, int order >
     struct CurvedSurfaceGridBase
     {
-      using HG = typename GridFunction::GridView::Grid;
+      using HG = GridOf_t<GridFunction>;
 
       using type = GridDefaultImplementation<
         HG::dimension, CGeo::DimRange<GridFunction>::value, typename HG::ctype,
@@ -59,20 +60,12 @@ namespace Dune
     using Super = CurvedSurfaceGridBase<GF,order>;
 
     // friend declarations
-    // friend class CGeo::HierarchicIterator<const Self,GridFunction>;
-    // template< int, class, bool > friend class CGeo::EntityBase;
-    // template< class, class, class > friend class CGeo::GridView;
-    // template< class, class > friend class CGeo::Intersection;
-    // template< class, class > friend class CGeo::IntersectionIterator;
     template< class, class > friend class CGeo::IdSet;
     template< class, class > friend class CGeo::IndexSet;
-    // template< class > friend class HostGridAccess;
-    // template< class, class > friend class CGeo::CommDataHandle;
 
   public:
     using GridFunction = GF;
     using GridFamily = CGeo::GridFamily<GF,order>;
-    using HostGrid = typename GridFunction::GridView::Grid;
 
     /** \name Traits
      *  \{ */
@@ -80,6 +73,8 @@ namespace Dune
     //! type of the grid traits
     using Traits = typename GridFamily::Traits;
 
+    using HostGrid = typename Traits::HostGrid;
+
     //! traits structure containing types for a codimension
     /**
      *  \tparam codim  codimension

dune/curvedsurfacegrid/gridfamily.hh

@@ -3,6 +3,8 @@
 #ifndef DUNE_CURVED_SURFACE_GRID_GRIDFAMILY_HH
 #define DUNE_CURVED_SURFACE_GRID_GRIDFAMILY_HH
 
+#include <type_traits>
+
 #include <dune/grid/common/grid.hh>
 #include <dune/curvedsurfacegrid/capabilities.hh>
 #include <dune/curvedsurfacegrid/entity.hh>
@@ -14,10 +16,10 @@
 #include <dune/curvedsurfacegrid/iterator.hh>
 #include <dune/curvedsurfacegrid/idset.hh>
 #include <dune/curvedsurfacegrid/indexsets.hh>
+#include <dune/curvedsurfacegrid/gridfunction/gridfunction.hh>
 
 #include <dune/functions/common/functionconcepts.hh>
 #include <dune/functions/common/signature.hh>
-#include <dune/functions/gridfunctions/gridviewentityset.hh>
 #include <dune/functions/gridfunctions/localderivativetraits.hh>
 
 namespace Dune
@@ -39,8 +41,7 @@ namespace Dune
     template< class GF, class LF >
     struct DifferentiableLocalFunction
     {
-      using GridView = typename GF::GridView;
-      using EntitySet = Functions::GridViewEntitySet<GridView,0>;
+      using EntitySet = typename GF::EntitySet;
       using LocalContext = typename EntitySet::Element;
 
       using Range = typename Functions::SignatureTraits<typename GF::Signature>::Range;
@@ -58,13 +59,12 @@ namespace Dune
       struct Traits
       {
         using Grid = CurvedSurfaceGrid<GF,order>;
-        using HostGrid = typename GF::GridView::Grid;
+        using HostGrid = GridOf_t<GF>;
 
         using GridFunction = GF;
         using LocalFunction = std::decay_t<decltype(localFunction(std::declval<GF>()))>;
 
         static const bool differentiableLocalFunction = DifferentiableLocalFunction<GF, LocalFunction>::value;
-        static_assert(differentiableLocalFunction, "");
 
         using ctype = typename HostGrid::ctype;
 

dune/curvedsurfacegrid/gridfunction/analyticgridfunction.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_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;
+
+    //! 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;
+  };
+
+
+  //! LocalFunction associated to the \ref AnalyticGridFunction
+  /**
+   * \tparam LocalContext  Context this localFunction can be bound to, e.g. a grid element
+   * \tparam F             Type of a function that can be evaluated in global coordinates.
+   **/
+  template< class LocalContext, class F >
+  class LocalAnalyticGridFunction;
+
+  //! Generator for \ref LocalAnalyticGridFunction
+  /**
+   * \param  ff            Function that can be evaluated at global coordinates of the \ref LocalContext
+   * \tparam LocalContext  Context this localFunction can be bound to, e.g. a grid element
+   **/
+  template< class LocalContext, class FF >
+  auto localAnalyticGridFunction (FF&& ff)
+  {
+    using F = std::decay_t<FF>;
+    return LocalAnalyticGridFunction<LocalContext, F>{std::forward<FF>(ff)};
+  }
+
+
+  template< class LC, class Functor >
+  class LocalAnalyticGridFunction
+  {
+  public:
+    using LocalContext = LC;
+    using Geometry = typename LocalContext::Geometry;
+
+    using Domain = typename Geometry::GlobalCoordinate;
+    using LocalDomain = typename Geometry::LocalCoordinate;
+    using Range = std::result_of_t<Functor(Domain)>;
+    using Signature = Range(LocalDomain);
+
+  public:
+    //! Constructor. Stores the functor f by value
+    template< class FF,
+      disableCopyMove<LocalAnalyticGridFunction, FF> = 0>
+    LocalAnalyticGridFunction (FF&& f)
+      : f_(std::forward<FF>(f))
+    {}
+
+    //! bind the LocalFunction to the local context
+    /**
+     * Stores the localContext and its geometry in a cache variable
+     **/
+    void bind (const LocalContext& localContext)
+    {
+      localContext_.emplace(localContext);
+      geometry_.emplace(localContext.geometry());
+    }
+
+    //! unbind the localContext from the localFunction
+    /**
+     * Reset the geometry
+     **/
+    void unbind ()
+    {
+      geometry_.reset();
+      localContext_.reset();
+    }
+
+    //! evaluate the stored function in local coordinates
+    /**
+     * Transform the local coordinates to global coordinates first,
+     * then evalute the stored functor.
+     **/
+    Range operator() (const LocalDomain& x) const
+    {
+      assert(!!geometry_);
+      return f_(geometry_->global(x));
+    }
+
+    //! return the bound localContext.
+    const LocalContext& localContext () const
+    {
+      assert(!!localContext_);
+      return *localContext_;
+    }
+
+    //! obtain the functor
+    const Functor& f () const
+    {
+      return f_;
+    }
+
+  private:
+    Functor f_;
+
+    // some caches
+    Std::optional<LocalContext> localContext_;
+    Std::optional<Geometry> geometry_;
+  };
+
+  //! Derivative of a \ref LocalAnalyticGridFunction
+  /**
+   * Participates in overload resolution only if the functor `F` is differentiable
+   **/
+  template< class LC, class F >
+  auto derivative (const LocalAnalyticGridFunction<LC,F>& t)
+    -> decltype(localAnalyticGridFunction<LC>(derivative(t.f())))
+  {
+    return localAnalyticGridFunction<LC>(derivative(t.f()));
+  }
+
+
+  //! GridFunction associated to the mapping F
+  /**
+   * \tparam Grid       The grid type with elements the corresponding LocalFunction can be bound to
+   * \tparam F          Type of a function that can be evaluated in global coordinates.
+   **/
+  template< class Grid, class F >
+  class AnalyticGridFunction;
+
+  //! Generator for \ref AnalyticGridFunction
+  /**
+   * \param  ff    Function that can be evaluated at global coordinates of the \ref Grid
+   * \tparam Grid  The grid type with elements the corresponding LocalFunction can be bound to
+   **/
+  template< class Grid, class FF >
+  auto analyticGridFunction (FF&& ff)
+  {
+    using F = std::decay_t<FF>;
+    return AnalyticGridFunction<Grid, F>{std::forward<FF>(ff)};
+  }
+
+
+  template< class GridType, class Functor >
+  class AnalyticGridFunction
+  {
+  public:
+    using Grid = GridType;
+    using EntitySet = GridEntitySet<Grid,0>;
+
+    using Domain = typename EntitySet::GlobalCoordinate;
+    using Range = std::result_of_t<Functor(Domain)>;
+    using Signature = Range(Domain);
+
+  public:
+    //! Constructor. Stores the functor f by value
+    template< class FF,
+      disableCopyMove<AnalyticGridFunction, FF> = 0>
+    AnalyticGridFunction (FF&& f)
+      : f_(std::forward<FF>(f))
+    {}
+
+    //! evaluate the stored function in global coordinates
+    Range operator() (const Domain& x) const
+    {
+      return f_(x);
+    }
+
+    //! construct the \ref LocalAnalyticGridFunction
+    using LocalFunction = LocalAnalyticGridFunction<typename EntitySet::Element, Functor>;
+    friend LocalFunction localFunction (const AnalyticGridFunction& t)
+    {
+      return LocalFunction(t.f_);
+    }
+
+    //! obtain the stored \ref GridEntitySet
+    const EntitySet& entitySet () const
+    {
+      return entitySet_;
+    }
+
+    //! obtain the functor
+    const Functor& f () const
+    {
+      return f_;
+    }
+
+  private:
+    Functor f_;
+    EntitySet entitySet_;
+  };
+
+  //! Derivative of an \ref AnalyticGridFunction
+  /**
+   * Participates in overload resolution only if the functor `F` is differentiable
+   **/
+  template< class Grid, class F >
+  auto derivative (const AnalyticGridFunction<Grid,F>& t)
+    -> decltype(analyticGridFunction<Grid>(derivative(t.f())))
+  {
+    return analyticGridFunction<Grid>(derivative(t.f()));
+  }
+
+} // end namespace Dune
+
+#endif // DUNE_CURVED_SURFACE_GRID_ANALYTIC_GRIDFUNCTION_HH

dune/curvedsurfacegrid/gridfunction/gridfunction.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_GRIDFUNCTION_HH
+#define DUNE_CURVED_SURFACE_GRID_GRIDFUNCTION_HH
+
+#include <type_traits>
+
+namespace Dune
+{
+  namespace Impl
+  {
+    template< class GF, class = void >
+    struct GridOf;
+
+    template< class GF >
+    struct GridOf<GF, std::void_t<typename GF::GridView>>
+    {
+      using type = typename GF::GridView::Grid;
+    };
+
+    template< class GF >
+    struct GridOf<GF, std::void_t<typename GF::Grid>>
+    {
+      using type = typename GF::Grid;
+    };
+
+  } // end namespace Impl
+
+  template< class GF >
+  using GridOf_t = typename Impl::GridOf<GF>::type;
+
+} // end namespace Dune
+
+#endif // DUNE_CURVED_SURFACE_GRID_GRIDFUNCTION_HH

dune/curvedsurfacegrid/gridfunction/spheregridfunction.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_SPHERE_GRIDFUNCTION_HH
+#define DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH
+
+#include <type_traits>
+
+#include <dune/curvedsurfacegrid/gridfunction/analyticgridfunction.hh>
+#include <dune/functions/common/defaultderivativetraits.hh>
+
+namespace Dune
+{
+  // Sphere functor
+  template< class T >
+  struct SphereProjection
+  {
+    T radius_;
+
+    SphereProjection (T radius)
+      : radius_(radius)
+    {}
+
+    //! project the coordinate to the sphere at origin with `radius`
+    template< class Domain >
+    Domain operator() (const Domain& x) const
+    {
+      return x * (radius_ / x.two_norm());
+    }
+
+    //! derivative of the projection
+    friend auto derivative (const SphereProjection& sphere)
+    {
+      return [r=sphere.radius_](auto const& x)
+      {
+        using Domain = std::decay_t<decltype(x)>;
+        using DerivativeTraits = Functions::DefaultDerivativeTraits<Domain(Domain)>;
+        typename DerivativeTraits::Range out;
+
+        auto nrm = x.two_norm();
+        for (int i = 0; i < out.N(); ++i)
+          for (int j = 0; j < out.M(); ++j)
+            out[i][j] = r * ((i == j ? 1 : 0) - (x[i]/nrm) * (x[j]/nrm)) / nrm;
+        return out;
+      };
+    }
+  };
+
+  //! construct a grid function representing a sphere parametrization
+  template< class Grid, class T >
+  auto sphereGridFunction (T radius)
+  {
+    return analyticGridFunction<Grid>(SphereProjection<T>{radius});
+  }
+
+} // end namespace Dune
+
+#endif // DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH

dune/curvedsurfacegrid/gridfunction/torusgridfunction.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_SPHERE_GRIDFUNCTION_HH
+#define DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH
+
+#include <type_traits>
+
+#include <dune/common/math.hh>
+#include <dune/curvedsurfacegrid/gridfunction/analyticgridfunction.hh>
+#include <dune/functions/common/defaultderivativetraits.hh>
+
+namespace Dune
+{
+  // torus functor
+  template< class T >
+  struct TorusProjection
+  {
+    T radius1_;
+    T radius2_;
+
+    TorusProjection (T radius1, T radius2)
+      : radius1_(radius1)
+      , radius2_(radius2)
+    {}
+
+    //! project the coordinate to the sphere at origin with `radius`
+    template< class Domain >
+    Domain operator() (const Domain& x) const
+    {
+      using std::sqrt;
+      auto scale1 = radius1_ / sqrt(x[0]*x[0] + x[1]*x[1]);
+      Domain center{x[0] * scale1, x[1] * scale1, 0};
+      Domain out{x[0] - center[0], x[1] - center[1], x[2]};
+      out *= radius2_ / out.two_norm();
+
+      return out + center;
+    }
+
+    //! derivative of the projection
+    friend auto derivative (const TorusProjection& t)
+    {
+      return [r1=t.radius1_,r2=t.radius2_](auto const& x)
+      {
+        using std::sqrt;
+        using Domain = std::decay_t<decltype(x)>;
+        using DerivativeTraits = Functions::DefaultDerivativeTraits<Domain(Domain)>;
+        using Matrix = typename DerivativeTraits::Range;
+
+        T x0 = x[0]*x[0], x1 = x[1]*x[1], x2 = x0+x1, x3 = sqrt(x2), x4 = T(1)/x3, x5 = r1*x4;
+        T x6 = sqrt(x2*x2*x2), x7 = r1/x6, x8 = x0*x7, x9 = T(1)-x5, x10 = x8+x9;
+        T x11 = x[2]*x[2], x12 = x9*x9, x13 = x0*x12, x14 = x1*x12, x15 = x11+x13+x14;
+        T x16 = r2/sqrt(x15), x17 = x1*x7, x18 = sqrt(power(r2*(x10 + x17)/x15, 3));
+        T x19 = sqrt(power(x2,5)), x20 = r1*x19, x21 = x2*x2, x22 = r1*x3;
+        T x23 = x21-r1*x6; x24 = r1-x3, x25 = x24*x24, x26 = x0*x25+x1*x25+x11*x2;
+        T x27 = x26/x2, x28 = sqrt(x27), x29 = r2*x25*x28, x30 = x2*x2*x2*x29, x31 = x26*x26;
+        T x32 = sqrt(x27*x27*x27), x33 = -r1*r2*sqrt(power(x2,13))*x32 + r1*sqrt(power(x2,9))*x31;
+        T x34 = T(1)/x31, x35 = x[0]*x[1]*x34/power(x2,6), x36 = r2*x[2], x37 = x[0]*x36;
+        T x38 = (x21*(x0*x24 + x1*x24 - x3*(x11 + x25)) + x26*x6)/(x19*x26*x28), x39 = x36*x[1];
+        T x40 = x24*x4/x32;
+
+        return Matrix{
+          {x10*x16 - x13*x18 + x5 - x8, -x35*(x30*(x1*x20 + x21*(x0*x22 + x23)) + x33), x37*x38},
+          {-x35*(x30*(x0*x20 + x21*(x1*x22 + x23)) + x33), -x14*x18 + x16*(x17 + x9) - x17 + x5, x38*x39},
+          {x37*x40, x39*x40, x21*x29*x34}
+        };
+      };
+    }
+  };
+
+  //! construct a grid function representing a torus parametrization
+  template< class Grid, class T >
+  auto torusGridFunction (T radius1, T radius2)
+  {
+    return analyticGridFunction<Grid>(TorusProjection<T>{radius1,radius2});
+  }
+
+} // end namespace Dune
+
+#endif // DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH

dune/curvedsurfacegrid/test/convergence.cc

@@ -7,6 +7,7 @@
 #include <iostream>
 #include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
 #include <dune/curvedsurfacegrid/curvedsurfacegrid.hh>
+#include <dune/curvedsurfacegrid/gridfunction/analyticgridfunction.hh>
 #include <dune/geometry/quadraturerules.hh>
 #include <dune/grid/io/file/gmshreader.hh>
 #include <dune/localfunctions/lagrange/pk.hh>
@@ -14,7 +15,7 @@
 #include <dune/foamgrid/foamgrid.hh>
 
 #include <dune/functions/common/differentiablefunctionfromcallables.hh>
-#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
+//#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
 
 #define STR(s) STR_HELPER(s)
 #define STR_HELPER(s) #s
@@ -177,7 +178,8 @@ int main(int argc, char** argv)
       return out;
     }
   );
-  auto sphereGridFct = Functions::makeAnalyticGridViewFunction(sphere, hostGrid->leafGridView());
+  //auto sphereGridFct = Functions::makeAnalyticGridViewFunction(sphere, hostGrid->leafGridView());
+  auto sphereGridFct = analyticGridFunction<HostGrid>(sphere);
 
   using Grid = CurvedSurfaceGrid<decltype(sphereGridFct)>;
   Grid grid(*hostGrid, sphereGridFct);