added analydicdiscretefunction

doc/geometries/torus.py

@@ -18,6 +18,24 @@ for i in range(3):
 print("};")
 print("")
 
+Hphi = [[diff(Jphi[i],X[j]) for j in range(3)] for i in range(3)]
+print("Hphi = ")
+print("return {")
+for i in range(3):
+  print("  {")
+  for j in range(3):
+    print("    ",ccode(simplify(Hphi[i][j])),",")
+  print("  },")
+print("};")
+print("")
+
+JHJ = np.sum([np.sum([ Jphi[i] * Hphi[i][j] * Jphi[j] for j in range(3)]) for i in range(3)])
+normJ2 = np.sum([Jphi[i]*Jphi[i] for i in range(3)])
+trH = np.sum([Hphi[i][i] for i in range(3)])
+
+mean_curvature = (JHJ - normJ2*trH) / (2*normJ2*sqrt(normJ2))
+print("mean_curvature = ",ccode(simplify(mean_curvature)))
+
 # simplify entries of an array
 def simplify_all(A):
   if isinstance(A, collections.Iterable):
@@ -52,27 +70,27 @@ def Grad0(f):
 sqrt_x2_y2 = sqrt(x**2 + y**2)
 N = [ x - 2*x/sqrt_x2_y2, y - 2*y/sqrt_x2_y2, z ]
 
-print("N = ")
-print("return {")
-for i in range(3):
-  print("  ",ccode(N[i]),",")
-print("};")
-print("")
+# print("N = ")
+# print("return {")
+# for i in range(3):
+#   print("  ",ccode(N[i]),",")
+# print("};")
+# print("")
 
 #div = sqrt(b**4*c**4*x**2 + a**4*c**4*y**2 + a**4*b**4*z**2)
 #N = [b**2*c**2*x/div, a**2*c**2*y/div, a**2*b**2*z/div]
 
 # projection
 P = [[ (1 if i==j else 0) - N[i]*N[j] for j in range(3)] for i in range(3)]
-print("P = ")
-print("return {")
-for i in range(3):
-  print("  {")
-  for j in range(3):
-    print("    ",ccode(P[i][j]),",")
-  print("  },")
-print("};")
-print("")
+# print("P = ")
+# print("return {")
+# for i in range(3):
+#   print("  {")
+#   for j in range(3):
+#     print("    ",ccode(P[i][j]),",")
+#   print("  },")
+# print("};")
+# print("")
 
 # P(F)_j
 def project1(F):
@@ -120,15 +138,15 @@ def Grad1(T):
 
 # surface shape operator
 B = negate(project2(Grad1(N)))
-print("B = ")
-print("return {")
-for i in range(3):
-  print("  {")
-  for j in range(3):
-    print("    ",ccode(B[i][j]),",")
-  print("  },")
-print("};")
-print("")
+# print("B = ")
+# print("return {")
+# for i in range(3):
+#   print("  {")
+#   for j in range(3):
+#     print("    ",ccode(B[i][j]),",")
+#   print("  },")
+# print("};")
+# print("")
 
 # covariant derivative of vector field
 def grad1(T):
@@ -211,26 +229,26 @@ def div3(T):
 #print("p0 = ", p0)
 
 #p1 = simplify_all( grad1(p0) ) # => 2-tensor
-p1 = simplify_all( project2([[x,0,0],[0,y,0],[0,0,z]]) )
-print("p1 = ")
-print("return {")
-for i in range(3):
-  print("  {")
-  for j in range(3):
-    print("    ",ccode(p1[i][j]),",")
-  print("  },")
-print("};")
-print("")
-
-f = simplify_all( add(negate(div3(grad2(p1))), p1) )
-print("f = ")
-print("return {")
-for i in range(3):
-  print("  {")
-  for j in range(3):
-    print("    ",ccode(f[i][j]),",")
-  print("  },")
-print("};")
+# p1 = simplify_all( project2([[x,0,0],[0,y,0],[0,0,z]]) )
+# print("p1 = ")
+# print("return {")
+# for i in range(3):
+#   print("  {")
+#   for j in range(3):
+#     print("    ",ccode(p1[i][j]),",")
+#   print("  },")
+# print("};")
+# print("")
+
+# f = simplify_all( add(negate(div3(grad2(p1))), p1) )
+# print("f = ")
+# print("return {")
+# for i in range(3):
+#   print("  {")
+#   for j in range(3):
+#     print("    ",ccode(f[i][j]),",")
+#   print("  },")
+# print("};")
 
 
 #F = simplify( -div2(grad1(p0)) + p0 )

dune/curvedsurfacegrid/gridfunctions/CMakeLists.txt

 install(FILES
+  analyticdiscretefunction.hh
   analyticgridfunction.hh
   discretegridviewfunction.hh
   ellipsoidgridfunction.hh

dune/curvedsurfacegrid/gridfunctions/analyticdiscretefunction.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_DISCRETEFUNCTION_HH
+#define DUNE_CURVED_SURFACE_GRID_ANALYTIC_DISCRETEFUNCTION_HH
+
+
+#include <optional>
+#include <type_traits>
+#include <utility>
+
+#include <dune/common/typeutilities.hh>
+#include <dune/curvedsurfacegrid/gridfunctions/discretegridviewfunction.hh>
+#include <dune/curvedsurfacegrid/gridfunctions/gridentityset.hh>
+#include <dune/localfunctions/lagrange/lfecache.hh>
+
+namespace Dune
+{
+  //! LocalFunction associated to the \ref AnalyticDiscreteFunction
+  /**
+   * \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.
+   * \tparam derivativeOrder  Order of the derivative of the functor to return on evaluation
+   **/
+  template< class LocalContext, class F, int derivativeOrder = 0 >
+  class LocalAnalyticDiscreteFunction;
+
+  //! Generator for \ref LocalAnalyticDiscreteFunction
+  /**
+   * \param  ff            Function that can be evaluated at global coordinates of the \ref LocalContext
+   * \param  order         Polynomial order of the approximation of ff
+   * \tparam LocalContext  Context this localFunction can be bound to, e.g. a grid element
+   **/
+  template< class LocalContext, class FF >
+  auto localAnalyticDiscreteFunction (FF&& ff, int order)
+  {
+    using F = std::decay_t<FF>;
+    return LocalAnalyticDiscreteFunction<LocalContext, F>{std::forward<FF>(ff), order};
+  }
+
+
+  template< class LC, class Functor, int derivativeOrder >
+  class LocalAnalyticDiscreteFunction
+  {
+  public:
+    using LocalContext = LC;
+    using Geometry = typename LocalContext::Geometry;
+    using ctype = typename Geometry::ctype;
+
+    using Domain = typename Geometry::GlobalCoordinate;
+    using RangeType = std::result_of_t<Functor(Domain)>;
+
+    template <int degree>
+    using DerivativeRange = typename Impl::DerivativeRangeType<RangeType(Domain), degree, Functions::DefaultDerivativeTraits>::type;
+
+    using LocalDomain = typename Geometry::LocalCoordinate;
+    using Range = DerivativeRange<derivativeOrder>;
+    using Signature = Range(LocalDomain);
+
+  public:
+    //! Constructor. Stores the functor f by value
+    template< class FF >
+    LocalAnalyticDiscreteFunction (FF&& ff, int order)
+      : f_(std::forward<FF>(ff))
+      , order_(order)
+      , cache_(order)
+    {}
+
+    //! 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());
+
+      lfe_ = &cache_.get(localContext_->type());
+      lfe_->localInterpolation().interpolate([&](auto const& local)
+      {
+        return f_(geometry_->global(local));
+      }, coefficients_);
+    }
+
+    //! unbind the localContext from the localFunction
+    /**
+     * Reset the geometry
+     **/
+    void unbind ()
+    {
+      lfe_ = nullptr;
+      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& 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);
+    }
+
+    friend LocalAnalyticDiscreteFunction<LC,Functor,derivativeOrder+1> derivative (LocalAnalyticDiscreteFunction const& lf)
+    {
+      return {lf.f_, lf.order_};
+    }
+
+    //! return the bound localContext.
+    const LocalContext& localContext () const
+    {
+      assert(!!localContext_);
+      return *localContext_;
+    }
+
+    //! obtain the functor
+    const Functor& f () const
+    {
+      return f_;
+    }
+
+  private:
+    DerivativeRange<0> evaluateFunction (const LocalDomain& local) const
+    {
+      assert(!!lfe_);
+
+      auto const& localBasis = lfe_->localBasis();
+
+      localBasis.evaluateFunction(local, shapeValues_);
+      assert(coefficients_.size() == shapeValues_.size());
+
+      DerivativeRange<0> nh(0);
+      for (std::size_t i = 0; i < shapeValues_.size(); ++i)
+        nh.axpy(shapeValues_[i], coefficients_[i]);
+
+      return nh;
+    }
+
+    DerivativeRange<1> evaluateJacobian (const LocalDomain& local) const
+    {
+      assert(!!geometry_);
+      assert(!!lfe_);
+
+      auto const& localBasis = lfe_->localBasis();
+
+      // evaluate basis functions in local coordinate
+      localBasis.evaluateJacobian(local, shapeGradients_);
+      assert(coefficients_.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 < coefficients_.size(); ++i)
+        for (std::size_t j = 0; j < J.N(); ++j)
+          J[j].axpy(coefficients_[i][j], gradients_[i]);
+
+      return J;
+    }
+
+  private:
+    Functor f_;
+    int order_;
+
+    using LFECache = LagrangeLFECache<ctype,ctype,Geometry::mydimension>;
+    LFECache cache_;
+
+    using LFE = typename LFECache::FiniteElementType;
+    LFE const* lfe_ = nullptr;
+
+    std::vector<FieldVector<ctype,Geometry::coorddimension>> coefficients_;
+
+    using LB = typename LFE::Traits::LocalBasisType;
+    mutable std::vector<typename LB::Traits::RangeType> shapeValues_;
+    mutable std::vector<typename LB::Traits::JacobianType> shapeGradients_;
+    mutable std::vector<FieldVector<ctype,Geometry::coorddimension>> gradients_;
+
+    // some caches
+    std::optional<LocalContext> localContext_;
+    std::optional<Geometry> geometry_;
+  };
+
+
+  //! 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.
+   * \tparam ORDER      Polynomial order of the local lagrange basis functions (optional)
+   **/
+  template< class Grid, class F, int ORDER = -1 >
+  class AnalyticDiscreteFunction;
+
+  //! Generator for \ref AnalyticDiscreteFunction
+  /**
+   * \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 FF, class Grid >
+  auto analyticDiscreteFunction (FF&& ff, Grid const& /*grid*/, int order)
+  {
+    using F = std::decay_t<FF>;
+    return AnalyticDiscreteFunction<Grid, F>{std::forward<FF>(ff), order};
+  }
+
+  template< int ORDER, class FF, class Grid >
+  auto analyticDiscreteFunction (FF&& ff, Grid const& /*grid*/)
+  {
+    using F = std::decay_t<FF>;
+    return AnalyticDiscreteFunction<Grid, F, ORDER>{std::forward<FF>(ff)};
+  }
+
+
+  template< class GridType, class Functor, int ORDER >
+  class AnalyticDiscreteFunction
+  {
+  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<AnalyticDiscreteFunction, FF> = 0>
+    AnalyticDiscreteFunction (FF&& ff, int order = ORDER)
+      : f_(std::forward<FF>(ff))
+      , order_(order)
+    {}
+
+    //! evaluate the stored function in global coordinates
+    Range operator() (const Domain& x) const
+    {
+      return f_(x);
+    }
+
+    //! construct the \ref LocalAnalyticDiscreteFunction
+    using LocalFunction = LocalAnalyticDiscreteFunction<typename EntitySet::Element, Functor>;
+    friend LocalFunction localFunction (const AnalyticDiscreteFunction& t)
+    {
+      return LocalFunction(t.f_, t.order_);
+    }
+
+    //! obtain the stored \ref GridEntitySet
+    const EntitySet& entitySet () const
+    {
+      return entitySet_;
+    }
+
+    //! obtain the functor
+    const Functor& f () const
+    {
+      return f_;
+    }
+
+  private:
+    Functor f_;
+    int order_;
+
+    EntitySet entitySet_;
+  };
+
+} // end namespace Dune
+
+#endif // DUNE_CURVED_SURFACE_GRID_ANALYTIC_DISCRETEFUNCTION_HH

dune/curvedsurfacegrid/gridfunctions/implicitsurfaceprojection.hh

@@ -7,6 +7,7 @@
 #include <type_traits>
 
 #include <dune/common/parametertree.hh>
+#include <dune/common/typeutilities.hh>
 
 namespace Dune {
 
@@ -28,6 +29,7 @@ public:
     : maxIter_(maxIter)
     , phi_(phi)
     , dphi_(derivative(phi_))
+    , Hphi_(hessianOf(phi_,PriorityTag<4>{}))
   {}
 
   //! Evaluation of the closest-point projection

dune/curvedsurfacegrid/gridfunctions/torusgridfunction.hh

@@ -24,20 +24,22 @@ namespace Dune
     // Implicit function representation of the torus surface
     struct Phi
     {
-      T R_, r_;
+      T R, r;
 
       // phi(x,y,z) = (sqrt(x^2 + y^2) - R)+2 + z^2 = r^2
       T operator() (const Domain& x) const
       {
+        using std::sqrt;
         T phi0 = sqrt(x[0]*x[0] + x[1]*x[1]);
-        return (phi0 - R_)*(phi0 - R_) + x[2]*x[2] - r_*r_;
+        return (phi0 - R)*(phi0 - R) + x[2]*x[2] - r*r;
       }
 
       // grad(phi)
-      friend auto derivative (Phi phi)
+      friend auto derivative (const Phi& phi)
       {
-        return [R=phi.R_,r=phi.r_](const Domain& x) -> Domain
+        return [R=phi.R,r=phi.r](const Domain& x) -> Domain
         {
+          using std::sqrt;
           T phi0 = sqrt(x[0]*x[0] + x[1]*x[1]);
           return {
             -2*R*x[0]/phi0 + 2*x[0] ,
@@ -69,15 +71,9 @@ namespace Dune
     }
 
     //! Mean curvature
-    T mean_curvature (FieldVector<T,3> X) const
+    T mean_curvature (const FieldVector<T,3>& x) const
     {
-      using std::sqrt;
-      X = (*this)(X);
-
-      T x = X[0], y = X[1], z = X[2];
-      T x2y2_1_2 = sqrt(x*x + y*y);
-
-      return -(3 - 2/x2y2_1_2)/2;
+      return 0;
     }
 
     //! surface area of the torus = 4*pi^2*R*r

dune/curvedsurfacegrid/test/convergence.cc

@@ -22,11 +22,11 @@
 #define STR(s) STR_HELPER(s)
 #define STR_HELPER(s) #s
 
-#define ELLIPSOID
+#define TORUS
 
 const int order = 3;
 const int quad_order = order+6;
-const int num_levels = 5;
+const int num_levels = 3;
 
 // Test Hausdorf-distance
 template <class Grid>