add ellipoid gridfunction and python generator

doc/geometries/ellipsoid.py

+from sympy import *
+import numpy as np
+import collections
+
+x, y, z = symbols("x y z")
+a, b, c = symbols("a b c")
+u, v = symbols("u v")
+
+# 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)
+
+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
+
+
+X = [x,y,z]
+
+# normal vector
+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]
+
+print("N = ")
+print("return {")
+for i in range(3):
+  print("  ",ccode(N[i]),",")
+print("};")
+print("")
+
+# projection
+P = [[ (1 if i==j else 0) - N[i]*N[j] for j in range(3)] for i in range(3)]
+
+# parametrization
+nrm_X = sqrt(x**2 + y**2 + z**2)
+X0 = [x/nrm_X, y/nrm_X, z/nrm_X]
+u = atan(X0[1]/X0[0])
+v = acos(X0[2])
+X1 = [a*cos(u)*sin(v), b*sin(u)*sin(v), c*cos(v)]
+
+# jacobian of parametrization
+J = simplify_all([[diff(X1[i],X[j]) for i in range(3)] for j in range(3)])
+print("J = ")
+print("return {")
+for i in range(3):
+  print("  {")
+  for j in range(3):
+    print("    ",ccode(J[i][j]),",")
+  print("  },")
+print("};")
+print("")
+
+print("equal = ",simplify((x**2 + y**2 + z**2)**2 - (x**4 + 2*x**2*y**2 + 2*x**2*z**2 + y**4 + 2*y**2*z**2 + z**4)))
+
+# P(F)_j
+def project1(F):
+  return simplify_all([
+    np.sum([
+      P[i][j] * F[i]
+    for i in range(3)])
+    for j in range(3)])
+
+# P(F)_kl
+def project2(F):
+  return simplify_all([[
+    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 simplify_all([[[
+    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 simplify_all([[[[
+    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)])
+
+
+# Euclidean gradient
+def Grad0(f):
+  return simplify_all([
+    diff(f, X[i])
+  for i in range(3)])
+
+# surface gradient (covariant derivative of scalars)
+def grad0(f):
+  return project1(Grad0(f))
+
+def Grad1(T):
+  return simplify_all([[
+    diff(T[i], X[j])
+  for j in range(3)] for i in range(3)])
+
+# surface shape operator
+#B = negate(project2(Grad1(N)))
+
+# covariant derivative of vector field
+def grad1(T):
+  return simplify_all(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 simplify_all([[[
+    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 simplify_all([[[
+    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 simplify_all([[[[
+    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 simplify_all(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 simplify_all([[
+    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(x*y*z) ) # => vector
+#print("p0 = ", p0)
+
+# p1 = simplify_all( grad1(p0) ) # => 2-tensor
+
+# print("p1 = ", p1)
+
+# f = simplify_all( add(negate(div3(grad2(p1))), p1) )
+# print("f = ", f)
+
+#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

+// -*- 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/gridfunctions/analyticgridfunction.hh>
+#include <dune/functions/common/defaultderivativetraits.hh>
+
+namespace Dune
+{
+  // Ellipsoid functor
+  template< class T >
+  class EllipsoidProjection
+  {
+    T a_;
+    T b_;
+    T c_;
+
+  public:
+    //! Constructor of ellipsoid by major axes
+    EllipsoidProjection (T a, T b, T c)
+      : a_(a)
+      , b_(b)
+      , c_(c)
+    {}
+
+    //! project the coordinate to the ellipsoid
+    template< class Domain >
+    Domain operator() (const Domain& X) const
+    {
+      using std::sin; using std::cos;
+      auto [phi,theta] = angles(X);
+      return {a_*cos(phi)*sin(theta), b_*sin(phi)*sin(theta), c_*cos(theta)};
+    }
+
+    //! derivative of the projection
+#if0
+    friend auto derivative (const EllipsoidProjection& 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;
+
+        T x = X[0], y = X[1], z = X[2];
+        T x2 = x*x, y2 = y*y, z2 = z*z;
+        T x5 = x2*x2*x;
+
+        T nrm0 = x2 + y2;
+        T nrm1 = x2 + y2 + z2;
+        T nrm2 = sqrt(nrm0/nrm1);
+        T nrm3 = sqrt(nrm0/x2);
+        T nrm4 = power(sqrt(nrm1), 3);
+        T nrm5 = sqrt(nrm0)*nrm1*nrm1/sqrt(nrm1);
+
+        return {
+          {
+             a*x*nrm3*(y2 + z2)/nrm5 ,
+            -b*y*nrm0/(nrm3*nrm5) ,
+            -c*z*x/nrm4
+          },
+          {
+            -a*x2*y*nrm3/nrm5 ,
+             b*nrm0*nrm0*nrm0*(x2 + z2)/(x5*power(nrm3, 5)*nrm5) ,
+            -c*y*z/nrm4
+          },
+          {
+            -a*z*nrm0/(nrm3*nrm5) ,
+            -b*y*z*nrm0/(x*nrm3*nrm5) ,
+             c*(x2 + y2)/nrm4
+          }
+        };
+
+
+      };
+    }
+#endif
+
+    //! Normal vector
+    template< class Domain >
+    Domain normal (const Domain& X) const
+    {
+      using std::sqrt;
+      T x = X[0], y = X[1], z = X[2];
+      T a2 = a_*a_, b2 = b_*b_, c2 = c_*c_;
+
+      auto div = sqrt(b2*b2*c2*c2*x*x + a2*a2*c2*c2*y*y + a2*a2*b2*b2*z*z);
+      return {b2*c2*x/div, a2*c2*y/div, a2*b2*z/div};
+    }
+    
+    //! Mean curvature
+    template< class Domain >
+    T mean_curvature (const Domain& X) const
+    {
+      using std::sqrt; using std::abs;
+      T x = X[0], y = X[1], z = X[2];
+      T a2 = a_*a_, b2 = b_*b_, c2 = c_*c_;
+
+      auto div = 2*a2*b2*c2*power(sqrt(x*x/(a2*a2) + y*y/(b2*b2) + z*z/(c2*c2)), 3);
+      return abs(x*x + y*y + z*z - a2 - b2 - c2)/div;
+    }
+
+    //! Gaussian curvature
+    template< class Domain >
+    T gauss_curvature (const Domain& X) const 
+    {
+      T x = X[0], y = X[1], z = X[2];
+      T a2 = a_*a_, b2 = b_*b_, c2 = c_*c_;
+
+      auto div = a2*b2*c2*power(x*x/(a2*a2) + y*y/(b2*b2) + z*z/(c2*c2), 2);
+      return T(1)/div;
+    }
+
+  private:
+    FieldVector<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])};
+    }
+  };
+
+  //! construct a grid function representing a sphere parametrization
+  template< class Grid, class T >
+  auto ellipsoidGridFunction (T a, T b, T c)
+  {
+    return analyticGridFunction<Grid>(EllipsoidProjection<T>{a,b,c});
+  }
+
+} // end namespace Dune
+
+#endif // DUNE_CURVED_SURFACE_GRID_SPHERE_GRIDFUNCTION_HH