From 2aecaec577ca1541f124462ff56be3403d31271e Mon Sep 17 00:00:00 2001 From: Ralf Hartmann Date: Tue, 30 Jan 2001 16:21:18 +0000 Subject: [PATCH] New get_intermediate_points_between_points function. git-svn-id: https://svn.dealii.org/trunk@3842 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/grid/tria_boundary_lib.h | 20 ++- .../deal.II/source/grid/tria_boundary_lib.cc | 147 ++++++++++-------- 2 files changed, 102 insertions(+), 65 deletions(-) diff --git a/deal.II/deal.II/include/grid/tria_boundary_lib.h b/deal.II/deal.II/include/grid/tria_boundary_lib.h index a0e52fda18..859187a253 100644 --- a/deal.II/deal.II/include/grid/tria_boundary_lib.h +++ b/deal.II/deal.II/include/grid/tria_boundary_lib.h @@ -63,7 +63,8 @@ class HyperBallBoundary : public StraightBoundary * this class and the documentation of the * base class. * - * Only implemented for @p{dim=2}. + * Calls + * @p{get_intermediate_points_between_points}. */ virtual void get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, @@ -92,6 +93,23 @@ class HyperBallBoundary : public StraightBoundary double get_radius () const; protected: + + /** + * Called by + * @p{get_intermediate_points_on_line} + * and by + * @p{get_intermediate_points_on_quad}. + * + * Refer to the general + * documentation of + * @p{get_intermediate_points_on_line} + * in the documentation of the + * base class. + */ + void get_intermediate_points_between_points (const Point &p0, const Point &p1, + typename std::vector > &points) const; + + /** * Center point of the hyperball. */ diff --git a/deal.II/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/deal.II/source/grid/tria_boundary_lib.cc index dc11bc36b3..2049f93fff 100644 --- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc @@ -80,86 +80,97 @@ get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) c template <> void HyperBallBoundary<1>::get_intermediate_points_on_line ( - const Triangulation<1>::line_iterator &line, - std::vector > &points) const + const Triangulation<1>::line_iterator &, + std::vector > &) const { Assert(false, ExcInternalError()); } #else + template void HyperBallBoundary::get_intermediate_points_on_line ( const typename Triangulation::line_iterator &line, typename std::vector > &points) const +{ + if (points.size()==1) + points[0]=get_new_point_on_line(line); + else + get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points); +} + + +template +void +HyperBallBoundary::get_intermediate_points_between_points ( + const Point &p0, const Point &p1, + typename std::vector > &points) const { const unsigned int n=points.size(); Assert(n>0, ExcInternalError()); - if (n==1) - points[0]=get_new_point_on_line(line); + + const Point v0=p0-center, + v1=p1-center; + const double length=sqrt((v1-v0).square()); + + double eps=1e-14; + Assert(fabs(v0.square()-radius*radius) pm=0.5*(v0+v1); + + const double h=sqrt(pm.square()); + + double beta=0; + unsigned int left_index=0, right_index=0; + if ((n+1)%2==0) + { + // if the number of + // parts is even insert + // the midpoint + left_index=(n-1)/2; + right_index=left_index; + points[left_index]=pm; + --left_index; + ++right_index; + beta=d_alpha; + } else { - Point v0=line->vertex(0)-center, - v1=line->vertex(1)-center; - const double length=sqrt((v1-v0).square()); - - double eps=1e-14; - Assert(fabs(v0.square()-radius*radius) pm=StraightBoundary::get_new_point_on_line (line)-center; - - const double h=sqrt(pm.square()); - - double beta=0; - unsigned int left_index=0, right_index=0; - if ((n+1)%2==0) - { - // if the number of - // parts is even insert - // the midpoint - left_index=(n-1)/2; - right_index=left_index; - points[left_index]=pm; - --left_index; - ++right_index; - beta=d_alpha; - } - else - { - beta=(n/2+1)*d_alpha-alpha/2.; - right_index=n/2; - left_index=n/2-1; - } + beta=(n/2+1)*d_alpha-alpha/2.; + right_index=n/2; + left_index=n/2-1; + } - // n even: m=n/2, - // n odd: m=(n-1)/2 - const unsigned int m=n/2; - for (unsigned int i=0; i=1, ExcInternalError()); + // n even: m=n/2, + // n odd: m=(n-1)/2 + const unsigned int m=n/2; + for (unsigned int i=0; i=1, ExcInternalError()); - double d=h*tan(beta); - points[right_index]=pm+d/length*(v1-v0); - points[left_index]=pm-d/length*(v1-v0); - } + double d=h*tan(beta); + points[right_index]=pm+d/length*(v1-v0); + points[left_index]=pm-d/length*(v1-v0); + } - // project the points from the - // straight line to the - // HyperBallBoundary - for (unsigned int i=0; i::get_intermediate_points_on_quad ( const Triangulation<3>::quad_iterator &quad, std::vector > &points) const { - Assert(points.size()==1, ExcNotImplemented()); - StraightBoundary<3>::get_intermediate_points_on_quad (quad, points); - Point<3> &middle=points[0]; + unsigned int m=static_cast (sqrt(points.size())); + Assert(points.size()==m*m, ExcInternalError()); + + vector > lp3(m); + vector > lp1(m); - middle -= center; - // project to boundary - middle *= radius / sqrt(middle.square()); + get_intermediate_points_on_line(quad->line(3), lp3); + get_intermediate_points_on_line(quad->line(1), lp1); - middle += center; + vector > lps(m); + for (unsigned int i=0; i