Update the ellipsoid and torus parametrization

doc/geometries/ellipsoid.py

@@ -66,6 +66,16 @@ for i in range(3):
 print("};")
 print("")
 
+
+phi = x**2/a**2 + y**2/b**2 + z**2/c**2 - 1
+Jphi = [diff(phi,X[i]) for i in range(3)]
+print("Jphi = ")
+print("return {")
+for i in range(3):
+  print("    ",ccode(simplify(Jphi[i])),",")
+print("};")
+print("")
+
 # P(F)_j
 def project1(F):
   return simplify_all([

doc/geometries/torus.py

+from sympy import *
+import numpy as np
+import collections
+
+x, y, z = symbols("x y z", real=True)
+R, r = symbols("R r", positive=True)
+
+X = [x,y,z]
+
+phi0 = sqrt(x**2 + y**2) - R
+phi = phi0**2 + z**2 - r**2
+
+Jphi = [diff(phi,X[i]) for i in range(3)]
+print("Jphi = ")
+print("return {")
+for i in range(3):
+  print("    ",ccode(simplify(Jphi[i])),",")
+print("};")
+print("")
+
+# simplify entries of an array
+def simplify_all(A):
+  if isinstance(A, collections.Iterable):
+    return list(map(lambda a: simplify_all(a), A))
+  else:
+    return simplify(A) #.subs(phi0, r**2-z**2)
+
+def add(A,B):
+  if isinstance(A, collections.Iterable):
+    return list(map(lambda a,b: add(a,b), A,B))
+  else:
+    return A + B
+
+def negate(A):
+  if isinstance(A, collections.Iterable):
+    return list(map(lambda a: negate(a), A))
+  else:
+    return -A
+
+# Euclidean gradient
+def Grad0(f):
+  return simplify_all([
+    diff(f, X[i])
+  for i in range(3)])
+
+
+
+# grad_phi = Grad0(phi)
+# div = sqrt(grad_phi[0]**2 + grad_phi[1]**2 + grad_phi[2]**2)
+# N = simplify_all([ grad_phi[i]/div for i in range(3) ])
+
+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("")
+
+#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("")
+
+# P(F)_j
+def project1(F):
+  return [
+    np.sum([
+      P[i][j] * F[i]
+    for i in range(3)])
+    for j in range(3)]
+
+# P(F)_kl
+def project2(F):
+  return [[
+    np.sum([np.sum([
+      P[i][k] * P[j][l] * F[i][j]
+    for j in range(3)]) for i in range(3)])
+    for l in range(3)]  for k in range(3)]
+
+# P(F)_lmn
+def project3(F):
+  return [[[
+    np.sum([np.sum([np.sum([
+      P[i][l] * P[j][m] * P[k][n] * F[i][j][k]
+    for k in range(3)]) for j in range(3)]) for i in range(3)])
+    for n in range(3)]  for m in range(3)]  for l in range(3)]
+
+# P(F)_mnop
+def project4(F):
+  return [[[[
+    np.sum([np.sum([np.sum([np.sum([
+      P[i][m] * P[j][n] * P[k][o] * P[l][p] * F[i][j][k][l]
+    for l in range(3)]) for k in range(3)]) for j in range(3)]) for i in range(3)])
+    for p in range(3)]  for o in range(3)]  for n in range(3)]  for m in range(3)]
+
+
+
+
+# surface gradient (covariant derivative of scalars)
+def grad0(f):
+  return project1(Grad0(f))
+
+def Grad1(T):
+  return [[
+    diff(T[i], X[j])
+  for j in range(3)] for i in range(3)]
+
+# 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("")
+
+# covariant derivative of vector field
+def grad1(T):
+  return add(project2(Grad1(T)), [[
+    np.sum([
+      B[i][j] * T[k] * N[k]
+    for k in range(3)])
+    for j in range(3)] for i in range(3)])
+
+def Grad2(T):
+  return [[[
+    diff(T[i][j], X[k])
+  for k in range(3)] for j in range(3)] for i in range(3)]
+
+def grad2(T):
+  return [[[
+    np.sum([np.sum([np.sum([
+      P[L1][I1] * P[L2][I2] * P[L3][K] * diff(T[L1][L2], X[L3])
+    for L3 in range(3)]) for L2 in range(3)]) for L1 in range(3)]) +
+    np.sum([np.sum([
+      B[K][I1] * P[J2][I2] * T[L][J2] * N[L]
+    for J2 in range(3)]) for L in range(3)]) +
+    np.sum([np.sum([
+      B[K][I2] * P[J1][I1] * T[J1][L] * N[L]
+    for J1 in range(3)]) for L in range(3)])
+    for K in range(3)]  for I2 in range(3)] for I1 in range(3)]
+
+def Grad3(T):
+  return [[[[
+    diff(T[i][j][k], X[l])
+  for l in range(3)] for k in range(3)] for j in range(3)] for i in range(3)]
+
+def grad3(T):
+  return add(project4(Grad3(T)), [[[[
+    np.sum([np.sum([np.sum([
+      B[K][I1] * P[J2][I2] * P[J3][I3] * T[L][J2][J3] * N[L]
+    for J2 in range(3)]) for J3 in range(3)]) for L in range(3)]) +
+    np.sum([np.sum([np.sum([
+      B[K][I2] * P[J1][I1] * P[J3][I3] * T[J1][L][J3] * N[L]
+    for J1 in range(3)]) for J3 in range(3)]) for L in range(3)]) +
+    np.sum([np.sum([np.sum([
+      B[K][I3] * P[J1][I1] * P[J2][I2] * T[J1][J2][L] * N[L]
+    for J1 in range(3)]) for J2 in range(3)]) for L in range(3)])
+  for K in range(3)] for I3 in range(3)] for I2 in range(3)] for I1 in range(3)])
+
+
+
+# normal-rotation of scalar field
+def rot0(f):
+  return [diff(f*N[2],y) - diff(f*N[1],z),
+          diff(f*N[0],z) - diff(f*N[2],x),
+          diff(f*N[1],x) - diff(f*N[0],y)]
+
+def Div1(F):
+  return diff(F[0],x) + diff(F[1],y) + diff(F[2],z)
+
+def div1(F):
+  return Div1(project1(F))
+
+# def div2(t):
+#   F = Matrix([div1(t.row(0).T), div1(t.row(1).T), div1(t.row(2).T)])
+#   return P*F
+
+# div(T)_I1,I2
+def div3(T):
+  return [[
+    np.sum([np.sum([np.sum([np.sum([np.sum([
+      P[L1][I1] * P[L2][I2] * P[L3][K] * P[L4][K] * diff(T[L1][L2][L3],X[L4])
+    for K in range(3)]) for L4 in range(3)]) for L3 in range(3)]) for L2 in range(3)]) for L1 in range(3)]) +
+    np.sum([np.sum([
+      B[I3][I1] * T[L][I2][I3] * N[L] +
+      B[I3][I2] * T[I1][L][I3] * N[L] +
+      B[I3][I3] * T[I1][I2][L] * N[L]
+    for I3 in range(3)]) for L in range(3)])
+    for I2 in range(3)]  for I1 in range(3)]
+
+
+
+#p0 = simplify_all( rot0(z) ) # => vector   # rot0(x*y*z)
+#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("};")
+
+
+#F = simplify( -div2(grad1(p0)) + p0 )
+#print("F = ", F)
+#print("F*N = ", simplify(F.dot(N)))
\ No newline at end of file

dune/curvedsurfacegrid/gridfunctions/ellipsoidgridfunction.hh

@@ -9,19 +9,37 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/math.hh>
 #include "analyticgridfunction.hh"
+#include "implicitsurfaceprojection.hh"
 
 namespace Dune
 {
   // Ellipsoid functor
-  template< class T >
+  template< class T = double >
   class EllipsoidProjection
   {
     using Domain = FieldVector<T,3>;
     using Jacobian = FieldMatrix<T,3,3>;
 
-    T a_;
-    T b_;
-    T c_;
+    // Implicit function representation of the ellipsoid surface
+    struct Phi
+    {
+      T a_, b_, c_;
+
+      // phi(x,y,z) = (x/a)^2 + (y/b)^2 + (z/c)^2 = 1
+      T operator() (const FieldVector<T,3>& x) const
+      {
+        return x[0]*x[0]/(a_*a_) + x[1]*x[1]/(b_*b_) + x[2]*x[2]/(c_*c_) - 1;
+      }
+
+      // grad(phi)
+      friend auto derivative (Phi phi)
+      {
+        return [a=phi.a_,b=phi.b_,c=phi.c_](const Domain& x) -> Domain
+        {
+          return { T(2*x[0]/(a*a)), T(2*x[1]/(b*b)), T(2*x[2]/(c*c)) };
+        };
+      }
+    };
 
   public:
     //! Constructor of ellipsoid by major axes
@@ -29,35 +47,13 @@ namespace Dune
       : a_(a)
       , b_(b)
       , c_(c)
+      , implicit_(Phi{a,b,c},100)
     {}
 
     //! project the coordinate to the ellipsoid
-    // NOTE: This is not a closest-point projection, but a spherical-coordinate projection
-    Domain operator() (const Domain& X) const
-    {
-      using std::sqrt;
-      T x = X[0], y = X[1], z = X[2];
-      T nrm = sqrt(x*x + y*y + z*z);
-      return { a_*x/nrm , b_*y/nrm , c_*z/nrm };
-    }
-
-    //! derivative of the projection
-    friend auto derivative (const EllipsoidProjection& ellipsoid)
+    Domain operator() (const Domain& x) const
     {
-      return [a=ellipsoid.a_,b=ellipsoid.b_,c=ellipsoid.c_](Domain const& X)
-      {
-        using std::sqrt;
-        T x = X[0], y = X[1], z = X[2];
-        T x2 = x*x, y2 = y*y, z2 = z*z;
-
-        T nrm2 = x2 + y2 + z2;
-        T nrm3 = sqrt(nrm2)*nrm2;
-        return Jacobian{
-          { a*(y2 + z2)/nrm3 ,  -b*y*x/nrm3 ,       -c*z*x/nrm3 },
-          {-a*x*y/nrm3 ,         b*(x2 + z2)/nrm3 , -c*z*y/nrm3 },
-          {-a*x*z/nrm3 ,        -b*y*z/nrm3 ,        c*(x2 + y2)/nrm3 }
-        };
-      };
+      return implicit_(x);
     }
 
     //! Normal vector
@@ -95,13 +91,8 @@ namespace Dune
     }
 
   private:
-    std::array<T,2> angles (Domain x) const
-    {
-      using std::acos; using std::atan2;
-      x /= x.two_norm();
-
-      return { atan2(x[1], x[0]), acos(x[2]) };
-    }
+    T a_, b_, c_;
+    ImplicitSurfaceProjection<Phi> implicit_;
   };
 
   //! construct a grid function representing a sphere parametrization

dune/curvedsurfacegrid/gridfunctions/implicitsurfaceprojection.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_IMPLICIT_SURFACE_PROJECTION_HH
+#define DUNE_CURVED_SURFACE_GRID_IMPLICIT_SURFACE_PROJECTION_HH
+
+#include <cmath>
+#include <type_traits>
+
+#include <dune/common/parametertree.hh>
+
+namespace Dune {
+
+//! Closest-point projection to surface given by implicit function
+/**
+ * Surface S is given by zero-levelset of function F.
+ * We assume that F is differentiable in order to evaluate normals and to
+ * do an iterative projection.
+ **/
+template< class F >
+class ImplicitSurfaceProjection
+{
+  using Functor = F;
+  using DFunctor = std::decay_t<decltype(derivative(std::declval<F>()))>;
+
+public:
+  //! Constructor for a given differentiable function with surface = { x : phi(x) == 0 }
+  ImplicitSurfaceProjection (const F& phi, int maxIter = 10)
+    : maxIter_(maxIter)
+    , phi_(phi)
+    , dphi_(derivative(phi_))
+  {}
+
+  //! Evaluation of the closest-point projection
+  /**
+   * Simple iterative algorithm proposed by Demlow and Dziuk in
+   * AN ADAPTIVE FINITE ELEMENT METHOD FOR THE LAPLACE-BELTRAMI OPERATOR ON IMPLICITLY DEFINED SURFACES
+   **/
+  template< class Domain >
+  Domain operator() (const Domain& x0) const
+  {
+    using std::sqrt;
+    using T = typename Domain::value_type;
+
+    Domain x = x0;
+    T tol = sqrt(std::numeric_limits<T>::epsilon());
+
+    auto phi = phi_(x);
+    auto sign = phi > 0 ? 1 : phi < 0 ? -1 : 0;
+    auto grad_phi = dphi_(x);
+    auto grad_phi_norm2 = grad_phi.two_norm2();
+
+    T err = 1;
+    for (int i = 0; i < maxIter_ && err > tol; ++i) {
+      auto y = x - grad_phi * (phi/grad_phi_norm2);
+
+      auto dist = sign * (y - x0).two_norm();
+      auto grad_phi_y = dphi_(y);
+
+      auto& dx = grad_phi_y;
+      dx *= -dist/grad_phi_y.two_norm();
+      x = x0 + dx;
+
+      phi = phi_(x);
+      grad_phi = dphi_(x);
+      grad_phi_norm2 = grad_phi.two_norm2();
+      err = sqrt(phi*phi / grad_phi_norm2 + (grad_phi / grad_phi_norm2 - dx / dx.two_norm()).two_norm2());
+    }
+
+    return x;
+  }
+
+  //! Normal vector given by grad(F)/|grad(F)|
+  template< class Domain >
+  Domain normal (const Domain& x0) const
+  {
+    auto grad_phi = dphi_(x0);
+    return grad_phi / grad_phi.two_norm();
+  }
+
+private:
+  int maxIter_;
+  Functor phi_;
+  DFunctor dphi_;
+};
+
+//! Construct a grid function representing a projection to an implicit surface
+template< class Grid, class F >
+auto implicitSurfaceGridFunction (const F& phi, int maxIter = 10)
+{
+  return analyticGridFunction<Grid>(ImplicitSurfaceProjection<F>{phi, maxIter});
+}
+
+} // end namespace Dune
+
+#endif // DUNE_CURVED_SURFACE_GRID_IMPLICIT_SURFACE_PROJECTION_HH

dune/curvedsurfacegrid/gridfunctions/torusgridfunction.hh

@@ -9,6 +9,7 @@
 #include <type_traits>
 
 #include <dune/common/math.hh>
+
 #include "analyticgridfunction.hh"
 
 namespace Dune
@@ -20,93 +21,74 @@ namespace Dune
     using Domain = FieldVector<T,3>;
     using Jacobian = FieldMatrix<T,3,3>;
 
-    T R_;
-    T r_;
+    // Implicit function representation of the torus surface
+    struct Phi
+    {
+      T R_, r_;
+
+      // phi(x,y,z) = (sqrt(x^2 + y^2) - R)+2 + z^2 = r^2
+      T operator() (const Domain& x) const
+      {
+        T phi0 = sqrt(x[0]*x[0] + x[1]*x[1]);
+        return (phi0 - R_)*(phi0 - R_) + x[2]*x[2] - r_*r_;
+      }
+
+      // grad(phi)
+      friend auto derivative (Phi phi)
+      {
+        return [R=phi.R_,r=phi.r_](const Domain& x) -> Domain
+        {
+          T phi0 = sqrt(x[0]*x[0] + x[1]*x[1]);
+          return {
+            -2*R*x[0]/phi0 + 2*x[0] ,
+            -2*R*x[1]/phi0 + 2*x[1] ,
+             2*x[2]
+          };
+        };
+      }
+    };
 
   public:
+    //! Construction of the torus function by major and minor radius
     TorusProjection (T R, T r)
       : R_(R)
       , r_(r)
+      , implicit_(Phi{R,r},100)
     {}
 
-    // closest point projection
-    Domain operator() (Domain x) const
-    {
-      using std::abs;
-
-      // on center of inner ring with z=0
-      T factor = R_ / (x[0]*x[0] + x[1]*x[1]);
-      Domain c{x[0] * factor, x[1] * factor, T(0)};
-      x -= c;
-
-      factor = r_ / x.two_norm();
-      x *= factor;
-      x += c;
-
-      // x is in the zero-levelset of the implicit function phi
-      assert(abs(phi(x)) < 10*std::numeric_limits<T>::epsilon());
-      return x;
-    }
-
-    //! Implicit representaion of the torus as phi(x) = 0
-    T phi (const Domain& x) const
+    //! Closest point projection
+    Domain operator() (const Domain& x) const
     {
-      using std::sqrt;
-      T phi0 = sqrt(x[0]*x[0] + x[1]*x[1]) - R_;
-      return phi0*phi0 + x[2]*x[2] - r_*r_;
+      return implicit_(x);
     }
 
-    Domain normal (Domain X) const
+    //! Normal vector
+    Domain normal (const Domain& 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);
-      // assert(x2y2_1_2 > std::numeric_limits<T>::epsilon());
-
-      // return { x - 2*x/x2y2_1_2, y - 2*y/x2y2_1_2, z };
-
-      T x2 = x*x, y2 = y*y, z2 = z*z;
-      auto factor1 = x2 + y2 + z2 - 5;
-      auto factor2 = x2 + y2 + z2 + 3;
-      auto factor1_2 = factor1*factor1;
-      auto factor2_2 = factor2*factor2;
-      auto div = sqrt(x2*factor1_2 + y2*factor1_2 + z2*factor2_2);
-      return {
-        x*factor1/div,
-        y*factor1/div,
-        z*factor2/div
-      };
-
+      return implicit_.normal(x);
     }
 
+    //! Mean curvature
     T mean_curvature (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;
-
-      T x2 = x*x, y2 = y*y, z2 = z*z;
-      T x3 = x*x2,y3 = y*y2,z3 = z*z2;
-      auto factor1 = x2 + y2 + z2 - 5;
-      auto factor2 = x2 + y2 + z2 + 3;
-      auto factor1_2 = factor1*factor1;
-      auto factor2_2 = factor2*factor2;
-      auto div = sqrt(x2*factor1_2 + y2*factor1_2 + z2*factor2_2);
-      auto div2= power(div, 3);
-      return -(2*x2/div + x*factor1*(-2*x3*factor1 - 2*x*y2*factor1 - 2*x*z2*factor2 - x*factor1_2)/div2 + 2*y2/div + y*factor1*(-2*x2*y*factor1 - 2*y3*factor1 - 2*y*z2*factor2 - y*factor1_2)/div2 + 2*z2/div + z*factor2*(-2*x2*z*factor1 - 2*y2*z*factor1 - 2*z3*factor2 - z*factor2_2)/div2 + 2*factor1/div + factor2/div) ;
+      T x2y2_1_2 = sqrt(x*x + y*y);
+
+      return -(3 - 2/x2y2_1_2)/2;
     }
 
-    //! surface area of the torus = 4*pi^2*r1*r2
+    //! surface area of the torus = 4*pi^2*R*r