From: heltai Date: Mon, 20 Jan 2014 10:05:24 +0000 (+0000) Subject: Reverted back the way HyperBallBoundary was creating intermediate points. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=aa6eea5e68bd7372422175251dff85741f0b9312;p=dealii-svn.git Reverted back the way HyperBallBoundary was creating intermediate points. git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32261 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/grid/tria_boundary_lib.h b/deal.II/include/deal.II/grid/tria_boundary_lib.h index 971558bcef..ee1a6319b4 100644 --- a/deal.II/include/deal.II/grid/tria_boundary_lib.h +++ b/deal.II/include/deal.II/grid/tria_boundary_lib.h @@ -312,6 +312,19 @@ protected: * called. */ bool compute_radius_automatically; + + +private: + + /** + * Called by @p get_new_point for special cases where surrounding + * points are two. Returns the point which is at xi between + * p0 and p1, where xi is between 0 and 1. + */ + Point get_intermediate_point (const Point &p0, + const Point &p1, + const double xi) const; + }; diff --git a/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/source/grid/tria_boundary_lib.cc index 42c35d8b0c..e0f2017682 100644 --- a/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/source/grid/tria_boundary_lib.cc @@ -183,21 +183,87 @@ Point HyperBallBoundary::get_new_point(const std::vector > &surrounding_points, const std::vector &weights) const { - Point middle = FlatManifold::get_new_point(surrounding_points, weights); - middle -= center; + // Special cases where the points are vertices in one or two dimensions + switch(surrounding_points.size()) { + case 2: + return get_intermediate_point(surrounding_points[0], surrounding_points[1], weights[1]); + break; + case 4: + { + Point v[2]; + // Coordinates of the point in the reference element + double xi1 = weights[1]+weights[3]; // = x + double xi2 = weights[2]+weights[3]; // = y + + v[0] = get_intermediate_point(surrounding_points[0], + surrounding_points[1], + xi1); + v[1] = get_intermediate_point(surrounding_points[2], + surrounding_points[3], + xi1); + return get_intermediate_point(v[0], v[1], xi2); + } + break; + default: + { + Point middle = FlatManifold::get_new_point(surrounding_points, weights); + middle -= center; + + double r=0; + if (compute_radius_automatically) + { + const Point vertex_relative = surrounding_points[0] - center; + r = std::sqrt(vertex_relative.square()); + } + else + r=radius; + // project to boundary + middle *= r / std::sqrt(middle.square()); + middle += center; + return middle; + } + break; + } +} + + +template +Point +HyperBallBoundary::get_intermediate_point(const Point &p0, + const Point &p1, + const double xi) const { + const Point v0=p0-center, + v1=p1-center; + const double length=std::sqrt((v1-v0).square()); + + double eps=1e-12; double r=0; if (compute_radius_automatically) { - const Point vertex_relative = surrounding_points[0] - center; + const Point vertex_relative = p0 - center; r = std::sqrt(vertex_relative.square()); } else r=radius; - // project to boundary - middle *= r / std::sqrt(middle.square()); - middle += center; - return middle; + + const double r2=r*r; + Assert(std::fabs(v0.square()-r2) pm=0.5*(v0+v1); + + const double h=std::sqrt(pm.square()); + + const double beta = alpha * (xi-0.5); + const double d=h*std::tan(beta); + Point p=pm+d/length*(v1-v0); + + p *= r / std::sqrt(p.square()); + p += center; + + return p; }