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
// ============================================================
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)
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);
// 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.
}
-
+// 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>::
{
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,
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)