From 7ea03960a4ee8a3150760fe406c7efefc20ac4ed Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 5 Nov 2017 22:57:17 +0100 Subject: [PATCH] Merge points in spherical manifold --- source/grid/manifold_lib.cc | 56 +++++++++++++++++++++++++++---------- 1 file changed, 41 insertions(+), 15 deletions(-) diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index d6af4b50c5..7c7325e231 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -421,20 +421,46 @@ get_new_point (const ArrayView> &vertices, // If the candidate happens to coincide with a normalized // direction, we return it. Otherwise, the Hessian would be singular. - boost::container::small_vector, 100> directions(n_points); + boost::container::small_vector, 100> directions; + boost::container::small_vector merged_weights; for (unsigned int i=0; i direction(vertices[i]-center); + const double norm = direction.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) + direction /= norm; + if ((candidate - direction).norm_square() < tolerance*tolerance) return center + rho * candidate; + + // append direction. check if the normalized candidate direction is + // the same as a previous direction (to a tighter tolerance (1e-14)^2 + // than the outer ones to really not miss anything) -> in that case we + // can simply add the weights. Since the trigonometric functions used + // below are quite expensive, it makes sense to merge the points here, + // even if this search loop is of quadratic complexity loop (but we + // rarely have more than 9 points) + Tensor<1,3> direction_3d; + for (unsigned int c=0; c vPerp; Tensor<2,2> Hessian; Tensor<1,2> gradient; @@ -456,15 +482,15 @@ get_new_point (const ArrayView> &vertices, // Then compute its contribution to the Hessian. gradient = 0.; Hessian = 0.; - for (unsigned int i=0; i1.e-15) + for (unsigned int i=0; i1.e-15) { vPerp = internal::projected_direction(directions[i], xVec); const double sintheta = vPerp.norm(); if (sintheta> &vertices, gradlocal[0] = cosphi; gradlocal[1] = sinphi; - gradient += (weights[i]*theta)*gradlocal; + gradient += (merged_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; + const double tt = (theta*sinthetaInv)*costheta; + const double offdiag = cosphi*sinphi*merged_weights[i]*(1.0-tt); + Hessian[0][0] += merged_weights[i]*(cosphiSq+tt*sinphiSq); Hessian[0][1] += offdiag; Hessian[1][0] += offdiag; - Hessian[1][1] += weights[i]*sinphiSq+tt*cosphiSq; + Hessian[1][1] += merged_weights[i]*(sinphiSq+tt*cosphiSq); } } -- 2.39.5