From 96d8e7a4bd961c6e78dfab10edfd94812a376689 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 2 Nov 2017 15:54:12 +0100 Subject: [PATCH] Improve SphericalManifold::get_new_point --- source/grid/manifold_lib.cc | 280 ++++++++++++++++++++++++++++-------- 1 file changed, 217 insertions(+), 63 deletions(-) diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 1fccfef375..fddb143ab0 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -25,6 +25,73 @@ 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 &p1, const Point &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 &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(center + w*v2 + (1-w)*v1); @@ -215,11 +283,12 @@ get_intermediate_point (const Point &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 &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 Point SphericalManifold:: @@ -281,71 +357,149 @@ get_new_point (const ArrayView> &vertices, { const unsigned int n_points = vertices.size(); - double rho = 0.0; - Tensor<1,spacedim> candidate; - for (unsigned int i = 0; i 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]) 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 - Point - compute_normal(const Tensor<1,spacedim> &vector) + Tensor<1, 3> xVec; + for (unsigned int c=0; c, 100> directions(n_points); + for (unsigned int i=0; i 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 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; i1.e-15) + { + vPerp = internal::projected_direction(directions[i], xVec); + const double sintheta = vPerp.norm(); + if (sinthetatolerance, 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 CylindricalManifold::CylindricalManifold(const unsigned int axis, @@ -362,7 +516,7 @@ CylindricalManifold::CylindricalManifold(const Point &d const Point &point_on_axis_, const double tolerance) : ChartManifold(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) -- 2.39.5