From 80c8f38b7527e582c3cbdfadd7fd0da29aeb9e9b Mon Sep 17 00:00:00 2001 From: heltai Date: Wed, 4 Dec 2013 22:27:09 +0000 Subject: [PATCH] Backward compatibility for Cylinder and Ball. Still missing the Half-Thingy family. git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@31888 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/deal.II/grid/tria_boundary_lib.h | 55 ++++++------ deal.II/source/grid/tria_boundary_lib.cc | 90 ++++++++++++++++--- 2 files changed, 105 insertions(+), 40 deletions(-) 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 20e2b5416f..54ad9af979 100644 --- a/deal.II/include/deal.II/grid/tria_boundary_lib.h +++ b/deal.II/include/deal.II/grid/tria_boundary_lib.h @@ -92,6 +92,15 @@ public: // get_new_point (const std::vector > &surrounding_points, // const std::vector &weights) const; + /** + * Since this class is derived from FlatManifold, we need to + * overload only the project_to_manifold function in order for the + * boundary to function properly. + */ + virtual Point + project_to_manifold (const Point candidate) const; + + /** * Refer to the general documentation of * this class and the documentation of the @@ -270,13 +279,12 @@ public: get_new_point_on_line (const typename Triangulation::line_iterator &line) const; /** - * Refer to the general documentation of - * this class and the documentation of the - * base class. + * Since this class is derived from FlatManifold, we need to + * overload only the project_to_manifold function in order for the + * boundary to function properly. */ -// virtual Point -// get_new_point (const std::vector > &surrounding_points, -// const std::vector &weights) const; + virtual Point + project_to_manifold (const Point candidate) const; /** @@ -413,16 +421,17 @@ public: Point get_new_point_on_line (const typename Triangulation::line_iterator &line) const; + /** * Refer to the general documentation of * this class and the documentation of the * base class. */ -// virtual Point -// get_new_point (const std::vector > &surrounding_points, -// const std::vector &weights) const; - + virtual Point + get_new_point(const std::vector > &surrounding_points, + const std::vector &weights) const; + /** * Refer to the general documentation of * this class and the documentation of the @@ -599,15 +608,18 @@ public: virtual Point get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; + /** - * Refer to the general documentation of - * this class and the documentation of the + * Refer to the general + * documentation of this class + * and the documentation of the * base class. */ -// virtual Point -// get_new_point (const std::vector > &surrounding_points, -// const std::vector &weights) const; - + virtual + Point + get_new_point(const std::vector > &surrounding_points, + const std::vector &weights) const; + /** * Refer to the general @@ -719,16 +731,7 @@ public: */ virtual Point get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - /** - * Refer to the general documentation of - * this class and the documentation of the - * base class. - */ -// virtual Point -// get_new_point (const std::vector > &surrounding_points, -// const std::vector &weights) const; - - + /** * Construct a new point on a quad. */ diff --git a/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/source/grid/tria_boundary_lib.cc index dfe6410b08..9bd26fb208 100644 --- a/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/source/grid/tria_boundary_lib.cc @@ -87,6 +87,31 @@ get_new_point_on_line (const typename Triangulation::line_iterator point_on_axis); } +template +Point +CylinderBoundary:: +project_to_manifold (const Point middle) const +{ + // The temptative point is computed from the get_new_point function + // of the FlatManifold class, which in turns calls this class + // internally. We have project this point out to the given radius + // from the axis. to this end, we have to take into account the + // offset point_on_axis and the direction of the axis + const Point vector_from_axis = (middle-point_on_axis) - + ((middle-point_on_axis) * direction) * direction; + // scale it to the desired length + // and put everything back + // together, unless we have a point + // on the axis + if (vector_from_axis.norm() <= 1e-10 * middle.norm()) + return middle; + else + return (vector_from_axis / vector_from_axis.norm() * radius + + ((middle-point_on_axis) * direction) * direction + + point_on_axis); +} + + template<> @@ -343,6 +368,21 @@ get_new_point_on_line (const typename Triangulation::line_iterator &line) c } +template +Point +ConeBoundary:: +project_to_manifold (const Point middle) const +{ + const Point axis = x_1 - x_0; + // To project it on the boundary of the cone we first compute the orthogonal + // projection of the middle point onto the axis of the cone. + const double c = (middle - x_0) * axis / axis.square (); + const Point middle_p = x_0 + c * axis; + // Compute the projection of the middle point on the boundary of the cone. + return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm (); +} + + template <> Point<3> @@ -527,6 +567,28 @@ HyperBallBoundary::get_new_point_on_line (const typename Triangula } +template +Point +HyperBallBoundary::get_new_point(const std::vector > &surrounding_points, + const std::vector &weights) const +{ + Point middle = FlatManifold::get_new_point(surrounding_points, weights); + + 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; +} + + template <> Point<1> @@ -807,28 +869,28 @@ HalfHyperBallBoundary::HalfHyperBallBoundary (const Point center, {} +template +Point +HalfHyperBallBoundary::get_new_point_on_line(const typename Triangulation::line_iterator &line) const +{ + Assert(false, ExcInternalError()); + return Point(); +} + template Point -HalfHyperBallBoundary:: -get_new_point_on_line (const typename Triangulation::line_iterator &line) const +HalfHyperBallBoundary::get_new_point(const std::vector > &surrounding_points, + const std::vector &weights) const { // check whether center of object is at x==x_center, since then it belongs // to the plane part of the boundary. however, this is not the case if it is // at the outer perimeter - const Point line_center = line->center(); - const Point vertices[2] = { line->vertex(0), line->vertex(1) }; - - if ((line_center(0) == this->center(0)) - && - ((std::fabs(vertices[0].distance(this->center)-this->radius) > - 1e-5*this->radius) - || - (std::fabs(vertices[1].distance(this->center)-this->radius) > - 1e-5*this->radius))) - return line_center; + const Point object_center = FlatManifold::get_new_point(surrounding_points, weights); + if (object_center.distance(this->center) < 1e-5) + return object_center; else - return HyperBallBoundary::get_new_point_on_line (line); + return HyperBallBoundary::get_new_point (surrounding_points,weights); } -- 2.39.5