]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Improve SphericalManifold::get_new_point
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Thu, 2 Nov 2017 14:54:12 +0000 (15:54 +0100)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Thu, 2 Nov 2017 14:54:12 +0000 (15:54 +0100)
source/grid/manifold_lib.cc

index 1fccfef375ebaa5fb516f0a052bd2075029fadde..fddb143ab00eece3fa37a5850dfb6a34d9c0fc23 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
+
+namespace internal
+{
+  // Rotate a given unit vector u around the axis dir
+  // where the angle is given by the length of dir.
+  // This is the exponential map for a sphere.
+  Tensor<1,3>
+  apply_exponential_map (const Tensor<1,3> &u,
+                         const Tensor<1,3> &dir)
+  {
+    const double theta = dir.norm();
+    if (theta<1.e-10)
+      {
+        return u;
+      }
+    else
+      {
+        const Tensor<1,3> dirUnit = dir/theta;
+        const Tensor<1,3> tmp = cos(theta)*u + sin(theta)*dirUnit;
+        return tmp/tmp.norm();
+      }
+  }
+
+  // Returns the direction to go from v to u
+  // projected to the plane perpendicular to the unit vector v.
+  // This one is more stable when u and v are nearly equal.
+  Tensor<1,3>
+  projected_direction (const Tensor<1,3> &u,
+                       const Tensor<1,3> &v)
+  {
+    Tensor<1,3> ans = u-v;
+    ans -= (ans*v)*v;
+    return ans;             // ans = (u-v) - ((u-v)*v)*v
+  }
+
+  // helper function to compute a vector orthogonal to a given one.
+  Point<3>
+  compute_normal(const Tensor<1,3> &vector)
+  {
+    Assert(vector.norm_square() != 0.,
+           ExcMessage("The direction parameter must not be zero!"));
+    Point<3> normal;
+    if (std::abs(vector[0]) >= std::abs(vector[1])
+        && std::abs(vector[0]) >= std::abs(vector[2]))
+      {
+        normal[1]=-1.;
+        normal[2]=-1.;
+        normal[0]=(vector[1]+vector[2])/vector[0];
+      }
+    else if (std::abs(vector[1]) >= std::abs(vector[0])
+             && std::abs(vector[1]) >= std::abs(vector[2]))
+      {
+        normal[0]=-1.;
+        normal[2]=-1.;
+        normal[1]=(vector[0]+vector[2])/vector[1];
+      }
+    else
+      {
+        normal[0]=-1.;
+        normal[1]=-1.;
+        normal[2]=(vector[0]+vector[1])/vector[2];
+      }
+    return normal;
+  }
+}
+
+
 // ============================================================
 // PolarManifold
 // ============================================================
@@ -177,9 +244,6 @@ get_intermediate_point (const Point<spacedim> &p1,
                         const Point<spacedim> &p2,
                         const double w) const
 {
-  Assert(w >=0.0 && w <= 1.0,
-         ExcMessage("w should be in the range [0.0,1.0]."));
-
   const double tol = 1e-10;
 
   if ( p1.distance(p2) < tol || w < tol)
@@ -203,6 +267,10 @@ get_intermediate_point (const Point<spacedim> &p1,
   const Tensor<1,spacedim> e1 = v1/r1;
   const Tensor<1,spacedim> e2 = v2/r2;
 
+  // Treat points that are collinear with the center special.
+  if ((e1 + e2).norm_square() == 0.)
+    return center;
+
   if ((e1 - e2).norm_square() < tol*tol)
     return Point<spacedim>(center + w*v2 + (1-w)*v1);
 
@@ -215,11 +283,12 @@ get_intermediate_point (const Point<spacedim> &p1,
   // Normal to v1 in the plane described by v1,v2,and the origin.
   // Since p1 and p2 do not coincide n is not zero and well defined.
   Tensor<1,spacedim> n = v2 - (v2*e1)*e1;
-  Assert( n.norm() > 0,
-          ExcInternalError("n should be different from the null vector."
-                           "Probably this means v1==v2 or v2==0."));
+  const double n_norm = n.norm();
+  Assert( n_norm > 0,
+          ExcInternalError("n should be different from the null vector. "
+                           "Probably, this means v1==v2 or v2==0."));
 
-  n /= n.norm();
+  n /= n_norm;
 
   // Find the point Q along O,v1 such that
   // P1,V,P2 has measure sigma.
@@ -272,7 +341,14 @@ get_tangent_vector (const Point<spacedim> &p1,
 }
 
 
-
+// The main part of the implementation uses the ideas in the publication
+//
+// Buss, Samuel R., and Jay P. Fillmore.
+// "Spherical averages and applications to spherical splines and interpolation."
+// ACM Transactions on Graphics (TOG) 20.2 (2001): 95-126.
+//
+// and in particular the implementation provided at
+// http://math.ucsd.edu/~sbuss/ResearchWeb/spheremean/
 template <int dim, int spacedim>
 Point<spacedim>
 SphericalManifold<dim,spacedim>::
@@ -281,71 +357,149 @@ get_new_point (const ArrayView<const Point<spacedim>> &vertices,
 {
   const unsigned int n_points = vertices.size();
 
-  double rho = 0.0;
-  Tensor<1,spacedim> candidate;
-  for (unsigned int i = 0; i<n_points; i++)
-    {
-      const Tensor<1,spacedim> direction (vertices[i]-center);
-      rho += direction.norm()*weights[i];
-      candidate += direction*weights[i];
-    }
+  if (n_points == 2)
+    return get_intermediate_point(vertices[0], vertices[1], weights[1]);
 
-  // Unit norm direction.
-  const double norm = candidate.norm();
-
-  if (norm == 0)
-    return center;
+  const double tolerance = 1e-10;
+  const int max_iterations = 10;
 
-  candidate /= norm;
+  double rho = 0.;
+  Tensor<1,spacedim> candidate;
 
-  return center+rho*candidate;
-}
+  // Step 1:
+  // Perform a simple average ...
+  {
+    double total_weights = 0.;
+    for (unsigned int i = 0; i < n_points; i++)
+      {
+        if (std::abs(1-weights[i])<tolerance)
+          return vertices[i];
+        const Tensor<1, spacedim> direction(vertices[i] - center);
+        rho += direction.norm() * weights[i];
+        candidate += direction * weights[i];
+        total_weights += weights[i];
+      }
+    // ... and normalize if the candidate is different from the origin.
+    const double norm = candidate.norm();
+    if (norm == 0.)
+      return center;
+    candidate /= norm;
+    rho /= total_weights;
+  }
 
-// ============================================================
-// CylindricalManifold
-// ============================================================
+  if (spacedim<2)
+    return center + rho*candidate;
 
-namespace internal
-{
-  namespace CylindricalManifold
+  // Step 2:
+  // Do Newton-style iterations to improve the estimate.
+  //
+  // In this step, we consider all points and directions to be embedded
+  // in a three-dimensional space.
   {
-    namespace
-    {
-      // helper function to compute a vector orthogonal to a given one.
-      template <int spacedim>
-      Point<spacedim>
-      compute_normal(const Tensor<1,spacedim> &vector)
+    Tensor<1, 3> xVec;
+    for (unsigned int c=0; c<spacedim; ++c)
+      xVec[c]=candidate[c];
+
+    // If the candidate happens to coincide with a normalized
+    // direction, we return it. Otherwise, the Hessian would be singular.
+    boost::container::small_vector<Tensor<1, 3>, 100> directions(n_points);
+    for (unsigned int i=0; i<n_points; ++i)
       {
-        AssertThrow(vector.norm() != 0.,
-                    ExcMessage("The direction parameter must not be zero!"));
-        Point<3> normal;
-        if (std::abs(vector[0]) >= std::abs(vector[1])
-            && std::abs(vector[0]) >= std::abs(vector[2]))
-          {
-            normal[1]=-1.;
-            normal[2]=-1.;
-            normal[0]=(vector[1]+vector[2])/vector[0];
-          }
-        else if (std::abs(vector[1]) >= std::abs(vector[0])
-                 && std::abs(vector[1]) >= std::abs(vector[2]))
-          {
-            normal[0]=-1.;
-            normal[2]=-1.;
-            normal[1]=(vector[0]+vector[2])/vector[1];
-          }
-        else
-          {
-            normal[0]=-1.;
-            normal[1]=-1.;
-            normal[2]=(vector[0]+vector[1])/vector[2];
-          }
-        return normal;
+        for (unsigned int c = 0; c < spacedim; ++c)
+          directions[i][c] = vertices[i][c] - center[c];
+        const double norm = directions[i].norm();
+        Assert(norm != 0.,
+               ExcMessage("One of the vertices coincides with the center. "
+                          "This is not allowed!"));
+        directions[i] /= norm;
+        if ((xVec - directions[i]).norm_square() < tolerance*tolerance)
+          return center + rho * candidate;
       }
-    }
+
+    Tensor<1,3> vPerp;
+    Tensor<2,2> Hessian;
+    Tensor<1,2> gradient;
+    Tensor<1,2> gradlocal;
+
+    // On success we exit the loop early.
+    // Otherwise, we just take the result after max_iterations steps.
+    for (unsigned int i=0; i<max_iterations; ++i)
+      {
+        // Step 2a: Find new descent direction
+
+        // Get local basis for the estimate xVec
+        const Tensor<1,3> Clocalx = internal::compute_normal(xVec);
+        const Tensor<1,3> Clocaly = cross_product_3d(xVec, Clocalx);
+
+        // For each vertices vector, compute the tangent vector from xVec
+        // towards the vertices vector -- its length is the spherical length
+        // from xVec to the vertices vector.
+        // Then compute its contribution to the Hessian.
+        gradient = 0.;
+        Hessian = 0.;
+        for (unsigned int i=0; i<n_points; ++i)
+          if (std::abs(weights[i])>1.e-15)
+            {
+              vPerp = internal::projected_direction(directions[i], xVec);
+              const double sintheta = vPerp.norm();
+              if (sintheta<tolerance)
+                {
+                  Hessian[0][0]+=weights[i];
+                  Hessian[1][1]+=weights[i];
+                }
+              else
+                {
+                  const double costheta = (directions[i])*xVec;
+                  const double theta = atan2(sintheta, costheta);
+                  const double sinthetaInv = 1.0/sintheta;
+
+                  vPerp *= sinthetaInv;
+                  const double cosphi = vPerp*Clocalx;
+                  const double sinphi = vPerp*Clocaly;
+
+                  gradlocal[0] = cosphi;
+                  gradlocal[1] = sinphi;
+                  gradient += (weights[i]*theta)*gradlocal;
+
+                  const double sinphiSq = sinphi*sinphi;
+                  const double cosphiSq = cosphi*cosphi;
+                  const double tt = weights[i]*(theta*sinthetaInv)*costheta;
+                  const double offdiag = cosphi*sinphi*(weights[i]-tt);
+                  Hessian[0][0] += weights[i]*cosphiSq+tt*sinphiSq;
+                  Hessian[0][1] += offdiag;
+                  Hessian[1][0] += offdiag;
+                  Hessian[1][1] += weights[i]*sinphiSq+tt*cosphiSq;
+                }
+            }
+
+        Assert(determinant(Hessian)>tolerance, ExcInternalError());
+
+        const Tensor<2,2> inverse_Hessian = invert(Hessian);
+
+        const Tensor<1,2> xDisplocal = inverse_Hessian*gradient;
+        const Tensor<1,3> xDisp = xDisplocal[0]*Clocalx + xDisplocal[1]*Clocaly;
+
+        // Step 2b: rotate xVec in direction xDisp for a new candidate.
+        const Tensor<1,3> xVecOld = xVec;
+        xVec = internal::apply_exponential_map(xVec, xDisp);
+
+        // Step 3c: return the new candidate if we didn't move
+        if ((xVec-xVecOld).norm_square() < tolerance*tolerance)
+          break;
+      }
+
+    for (unsigned int c=0; c<spacedim; ++c)
+      candidate[c] = xVec[c];
+
+    Assert (spacedim == 3 || std::abs(xVec[2]) < tolerance,
+            ExcInternalError());
   }
+  return center + rho*candidate;
 }
 
-
+// ============================================================
+// CylindricalManifold
+// ============================================================
 
 template <int dim, int spacedim>
 CylindricalManifold<dim, spacedim>::CylindricalManifold(const unsigned int axis,
@@ -362,7 +516,7 @@ CylindricalManifold<dim, spacedim>::CylindricalManifold(const Point<spacedim> &d
                                                         const Point<spacedim> &point_on_axis_,
                                                         const double tolerance) :
   ChartManifold<dim,3,3>(Tensor<1,3>({0,2.*numbers::PI,0})),
-              normal_direction(internal::CylindricalManifold::compute_normal(direction_)),
+              normal_direction(internal::compute_normal(direction_)),
               direction (direction_/direction_.norm()),
               point_on_axis (point_on_axis_),
               tolerance(tolerance)

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.