From 4e660d6df8a8a962fa1279a60ebbe9ca0ff91553 Mon Sep 17 00:00:00 2001 From: heltai Date: Mon, 9 Dec 2013 18:49:53 +0000 Subject: [PATCH] Fixed all boundary classes. Removed all references to ancient Boundary class. git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@31959 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-18/step-18.cc | 21 +- deal.II/examples/step-47/step-47.cc | 6 +- deal.II/include/deal.II/grid/manifold.h | 2 +- .../include/deal.II/grid/tria_boundary_lib.h | 607 +---- deal.II/source/grid/manifold.inst.in | 8 + deal.II/source/grid/tria_boundary_lib.cc | 2024 +++++++++-------- 6 files changed, 1128 insertions(+), 1540 deletions(-) diff --git a/deal.II/examples/step-18/step-18.cc b/deal.II/examples/step-18/step-18.cc index 283419e941..f587d740be 100644 --- a/deal.II/examples/step-18/step-18.cc +++ b/deal.II/examples/step-18/step-18.cc @@ -40,7 +40,7 @@ #include #include #include -#include +#include #include #include #include @@ -845,12 +845,12 @@ namespace Step18 // z-axis. Note that the boundary objects need to live as long as the // triangulation does; we can achieve this by making the objects static, // which means that they live as long as the program runs: - static const CylinderManifold inner_cylinder (inner_radius, 2); - static const CylinderManifold outer_cylinder (outer_radius, 2); + static const CylinderBoundary inner_cylinder (inner_radius, 2); + static const CylinderBoundary outer_cylinder (outer_radius, 2); // We then attach these two objects to the triangulation, and make them // correspond to boundary indicators 2 and 3: - triangulation.set_manifold (2, inner_cylinder); - triangulation.set_manifold (3, outer_cylinder); + triangulation.set_boundary (2, inner_cylinder); + triangulation.set_boundary (3, outer_cylinder); // There's one more thing we have to take care of (we should have done so // above already, but for didactic reasons it was more appropriate to @@ -1718,17 +1718,6 @@ namespace Step18 // to the vertex we consider at present, as they may or may not be ordered // in the same order. // - // This inconvenience could be avoided if finite elements have support - // points on vertices (which the one here has; for the concept of support - // points, see @ref GlossSupport "support points"). For such a case, one - // could construct a custom quadrature rule using - // FiniteElement::get_unit_support_points(). The first - // GeometryInfo@::vertices_per_cell*fe.dofs_per_vertex - // quadrature points will then correspond to the vertices of the cell and - // are ordered consistent with cell-@>vertex(i), taking into - // account that support points for vector elements will be duplicated - // fe.dofs_per_vertex times. - // // Another point worth explaining about this short function is the way in // which the triangulation class exports information about its vertices: // through the Triangulation::n_vertices function, it diff --git a/deal.II/examples/step-47/step-47.cc b/deal.II/examples/step-47/step-47.cc index 4dfcf9e10d..87707ba908 100644 --- a/deal.II/examples/step-47/step-47.cc +++ b/deal.II/examples/step-47/step-47.cc @@ -33,7 +33,7 @@ #include #include #include -#include +#include #include #include #include @@ -1065,8 +1065,8 @@ namespace Step47 GridGenerator::hyper_ball (triangulation); //GridGenerator::hyper_cube (triangulation, -1, 1); - static const PolarManifold boundary; - triangulation.set_manifold (0, boundary); + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); triangulation.refine_global (2); } diff --git a/deal.II/include/deal.II/grid/manifold.h b/deal.II/include/deal.II/grid/manifold.h index 36e73a0515..553d9bb0eb 100644 --- a/deal.II/include/deal.II/grid/manifold.h +++ b/deal.II/include/deal.II/grid/manifold.h @@ -280,7 +280,7 @@ public: private: FlatManifold sub_manifold; - std::vector > chart_points; + mutable std::vector > chart_points; }; 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 54ad9af979..a2a6d0ad68 100644 --- a/deal.II/include/deal.II/grid/tria_boundary_lib.h +++ b/deal.II/include/deal.II/grid/tria_boundary_lib.h @@ -48,7 +48,7 @@ DEAL_II_NAMESPACE_OPEN * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007 */ template -class CylinderBoundary : public StraightBoundary +class CylinderBoundary : public FlatManifold { public: /** @@ -60,38 +60,16 @@ public: const unsigned int axis = 0); /** - * Constructor. If constructed - * with this constructor, the - * boundary described is a - * cylinder with an axis that - * points in direction #direction - * and goes through the given - * #point_on_axis. The direction - * may be arbitrarily scaled, and - * the given point may be any - * point on the axis. + * Constructor. If constructed with this constructor, the boundary + * described is a cylinder with an axis that points in direction + * #direction and goes through the given #point_on_axis. The + * direction may be arbitrarily scaled, and the given point may be + * any point on the axis. */ CylinderBoundary (const double radius, const Point &direction, const Point &point_on_axis); - - /** - * Refer to the general documentation of - * this class and the documentation of the - * base class. - */ - 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; - + /** * Since this class is derived from FlatManifold, we need to * overload only the project_to_manifold function in order for the @@ -99,66 +77,21 @@ public: */ virtual Point project_to_manifold (const Point candidate) const; - - - /** - * Refer to the general documentation of - * this class and the documentation of the - * base class. - */ - 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 - * base class. - * - * Calls - * @p get_intermediate_points_between_points. + * Compute the normal to the surface. */ - virtual void - get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, - std::vector > &points) const; - - /** - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - * - * Only implemented for dim=3 - * and for points.size()==1. - */ - virtual void - get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, - std::vector > &points) const; - - /** - * Compute the normals to the - * boundary at the vertices of - * the given face. - * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - */ - virtual void - get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const; - + virtual Point + normal_vector (const Point point) const; + /** * Return the radius of the cylinder. */ double get_radius () const; /** - * Exception. Thrown by the - * @p get_radius if the - * @p compute_radius_automatically, - * see below, flag is set true. + * Exception. Thrown by the @p get_radius if the @p + * compute_radius_automatically, see below, flag is set true. */ DeclException0 (ExcRadiusNotSet); @@ -182,24 +115,8 @@ protected: private: /** - * 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, - std::vector > &points) const; - - /** - * Given a number for the axis, - * return a vector that denotes - * this direction. + * Given a number for the axis, return a vector that denotes this + * direction. */ static Point get_axis_vector (const unsigned int axis); }; @@ -240,21 +157,15 @@ private: * @author Markus Bürg, 2009 */ template -class ConeBoundary : public StraightBoundary +class ConeBoundary : public FlatManifold { public: /** - * Constructor. Here the boundary - * object is constructed. The - * points x_0 and - * x_1 describe the - * starting and ending points of - * the axis of the (truncated) - * cone. radius_0 - * denotes the radius - * corresponding to x_0 - * and radius_1 the one - * corresponding to x_1. + * Constructor. Here the boundary object is constructed. The points + * x_0 and x_1 describe the starting and ending + * points of the axis of the (truncated) cone. radius_0 + * denotes the radius corresponding to x_0 and + * radius_1 the one corresponding to x_1. */ ConeBoundary (const double radius_0, const double radius_1, @@ -262,22 +173,11 @@ public: const Point x_1); /** - * Return the radius of the - * (truncated) cone at given point + * Return the radius of the (truncated) cone at given point * x on the axis. */ double get_radius (const Point x) const; - - /** - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - */ - virtual - Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - + /** * Since this class is derived from FlatManifold, we need to * overload only the project_to_manifold function in order for the @@ -288,58 +188,14 @@ public: /** - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - */ - 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 - * base class. - * - * Calls @p - * get_intermediate_points_between_points. - */ - virtual - void - get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, - std::vector > &points) const; - - /** - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - * - * Only implemented for - * dim=3 and for - * points.size()==1. - */ - virtual - void - get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, - std::vector > &points) const; - - /** - * Compute the normals to the - * boundary at the vertices of - * the given face. + * Compute the normals to the boundary at the given point * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. + * Refer to the general documentation of this class and the + * documentation of the base class. */ virtual - void - get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const; + Point + normal_vector (const Point point) const; protected: /** @@ -363,24 +219,6 @@ protected: * Ending point of the axis. */ const Point x_1; - -private: - /** - * 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, - std::vector > &points) const; }; @@ -394,16 +232,12 @@ private: * The center of the ball and its radius may be given upon construction of * an object of this type. They default to the origin and a radius of 1.0. * - * This class is derived from StraightBoundary rather than from - * Boundary, which would seem natural, since this way we can use the - * StraightBoundary::in_between() function. - * * @ingroup boundary * - * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001 + * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001, Luca Heltai, 2013 */ template -class HyperBallBoundary : public StraightBoundary +class HyperBallBoundary : public FlatManifold { public: /** @@ -411,17 +245,7 @@ public: */ HyperBallBoundary (const Point p = Point(), const double radius = 1.0); - - /** - * Refer to the general documentation of - * this class and the documentation of the - * base class. - */ - 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 @@ -431,72 +255,17 @@ public: 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 - * base class. - */ - 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 - * base class. - * - * Calls - * @p get_intermediate_points_between_points. - */ - virtual - void - get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, - std::vector > &points) const; - - /** - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - * - * Only implemented for dim=3 - * and for points.size()==1. - */ - virtual - void - get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, - std::vector > &points) const; /** - * Implementation of the function - * declared in the base class. + * Implementation of the function declared in the base class. * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. + * Refer to the general documentation of this class and the + * documentation of the base class. */ virtual - Tensor<1,spacedim> - normal_vector (const typename Triangulation::face_iterator &face, - const Point &p) const; + Point + normal_vector (const Point p) const; - /** - * Compute the normals to the - * boundary at the vertices of - * the given face. - * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - */ - virtual - void - get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const; /** * Return the center of the ball. @@ -532,40 +301,15 @@ protected: const double radius; /** - * This flag is @p false for - * this class and for all derived - * classes that set the radius by - * the constructor. For example - * this flag is @p false for the - * HalfHyperBallBoundary - * class but it is @p true for - * the HyperShellBoundary - * class, for example. The - * latter class doesn't get its - * radii by the constructor but - * need to compute the radii - * automatically each time one of - * the virtual functions is + * This flag is @p false for this class and for all derived classes + * that set the radius by the constructor. For example this flag is + * @p false for the HalfHyperBallBoundary class but it is @p true + * for the HyperShellBoundary class, for example. The latter class + * doesn't get its radii by the constructor but need to compute the + * radii automatically each time one of the virtual functions is * called. */ bool compute_radius_automatically; - -private: - - /** - * 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, - std::vector > &points) const; }; @@ -592,74 +336,26 @@ public: HalfHyperBallBoundary (const Point p = Point(), const double radius = 1.0); - /** - * Check if on the line x==0, - * otherwise pass to the base - * class. - */ - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - - /** - * Check if on the line x==0, - * otherwise pass to the base - * class. - */ - 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 - * base class. + * 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; - - - /** - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - * - * Calls - * @p get_intermediate_points_between_points. - */ - virtual void - get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, - std::vector > &points) const; - - /** - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - * - * Only implemented for dim=3 - * and for points.size()==1. - */ - virtual void - get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, - std::vector > &points) const; + /** - * Compute the normals to the - * boundary at the vertices of - * the given face. + * Implementation of the function declared in the base class. * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. + * Refer to the general documentation of this class and the + * documentation of the base class. */ - virtual void - get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const; + virtual + Point + normal_vector (const Point p) const; }; @@ -726,58 +422,13 @@ public: const double inner_radius = -1, const double outer_radius = -1); - /** - * Construct a new point on a line. - */ - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - - /** - * Construct a new point on a quad. - */ - 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 - * base class. - * - * Calls - * @p get_intermediate_points_between_points. - */ - virtual void - get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, - std::vector > &points) const; - - /** - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. - * - * Only implemented for dim=3 - * and for points.size()==1. - */ - virtual void - get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, - std::vector > &points) const; /** - * Compute the normals to the - * boundary at the vertices of - * the given face. - * - * Refer to the general - * documentation of this class - * and the documentation of the - * base class. + * Compute the normal to the surface. */ - virtual void - get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const; - + virtual Point + normal_vector (const Point point) const; + private: /** * Inner and outer radii of the shell. @@ -796,7 +447,7 @@ private: * dim=2,spacedim=3, that is, just the surface. */ template -class TorusBoundary : public Boundary +class TorusBoundary : public ManifoldChart { public: /** @@ -805,93 +456,42 @@ public: */ TorusBoundary (const double R, const double r); -//Boundary Refinenment Functions - /** - * Construct a new point on a line. - */ - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - /** - * Construct a new point on a quad. - */ - 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 - * base class. - */ -// virtual Point -// get_new_point (const std::vector > &surrounding_points, -// const std::vector &weights) const; - - - /** - * Construct a new points on a line. - */ - virtual void get_intermediate_points_on_line ( - const typename Triangulation< dim, spacedim >::line_iterator &line, - std::vector< Point< spacedim > > &points) const ; - - /** - * Construct a new points on a quad. + * Pull back the given point in spacedim to the euclidean dim + * dimensional space. Get the surface coordinates of the Torus, + * i.e., from (x,y,z) to (theta,phi). + * + * Refer to the general documentation of this class for more + * information. */ - virtual void get_intermediate_points_on_quad ( - const typename Triangulation< dim, spacedim >::quad_iterator &quad, - std::vector< Point< spacedim > > &points ) const ; + virtual Point + pull_back(const Point space_point) const; + + /** + * Given a point in the dim dimensianal Euclidean space, this + * method returns a point on the manifold embedded in the spacedim + * Euclidean space. Get the cartesian coordinates of the Torus, + * i.e., from (theta,phi) to (x,y,z). + */ + virtual Point + push_forward(const Point chart_point) const; /** - * Get the normal from cartesian - * coordinates. This normal does not have - * unit length. + * Get the normal from cartesian coordinates. This normal does not + * have unit length. */ - virtual void get_normals_at_vertices ( - const typename Triangulation< dim, spacedim >::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const ; + Point normal_vector(const Point p) const; private: - //Handy functions - /** - * Function that corrects the value and - * sign of angle, that is, given - * angle=tan(abs(y/x)); if - * (y > 0) && (x < 0) then - * correct_angle = Pi - angle, - * etc. - */ - double get_correct_angle(const double angle,const double x,const double y) const; - /** - * Get the cartesian coordinates of the - * Torus, i.e., from (theta,phi) - * to (x,y,z). - */ - Point get_real_coord(const Point &surfP) const ; - - /** - * Get the surface coordinates of the - * Torus, i.e., from (x,y,z) to - * (theta,phi). - */ - Point get_surf_coord(const Point &p) const ; - /** * Get the normal from surface * coordinates. This normal does not have * unit length. */ Point get_surf_norm_from_sp(const Point &surfP) const ; - - /** - * Get the normal from cartesian - * coordinates. This normal does not have - * unit length. - */ - Point get_surf_norm(const Point &p) const ; - + /** * Inner and outer radii of the shell. */ @@ -903,59 +503,6 @@ private: /* -------------- declaration of explicit specializations ------------- */ -#ifndef DOXYGEN - -template <> -Point<1> -HyperBallBoundary<1>:: -get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const; -template <> -void -HyperBallBoundary<1>::get_intermediate_points_on_line ( - const Triangulation<1>::line_iterator &, - std::vector > &) const; -template <> -void -HyperBallBoundary<3>::get_intermediate_points_on_quad ( - const Triangulation<3>::quad_iterator &quad, - std::vector > &points) const; -template <> -void -HyperBallBoundary<1>:: -get_normals_at_vertices (const Triangulation<1,1>::face_iterator &, - Boundary<1,1>::FaceVertexNormals &) const; -template <> -Point<1> -HalfHyperBallBoundary<1>:: -get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const; -template <> -void -HalfHyperBallBoundary<1>:: -get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &, - std::vector > &) const; -template <> -void -HalfHyperBallBoundary<1>:: -get_normals_at_vertices (const Triangulation<1,1>::face_iterator &, - Boundary<1,1>::FaceVertexNormals &) const; -template <> -Point<1> -HalfHyperShellBoundary<1>:: -get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const; -template <> -void -HalfHyperShellBoundary<1>:: -get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &, - std::vector > &) const; -template <> -void -HalfHyperShellBoundary<1>:: -get_normals_at_vertices (const Triangulation<1,1>::face_iterator &, - Boundary<1,1>::FaceVertexNormals &) const; - - -#endif // DOXYGEN - DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/source/grid/manifold.inst.in b/deal.II/source/grid/manifold.inst.in index 2e8b1f8bfd..7f99d819b1 100644 --- a/deal.II/source/grid/manifold.inst.in +++ b/deal.II/source/grid/manifold.inst.in @@ -20,6 +20,14 @@ for (deal_II_space_dimension : SPACE_DIMENSIONS) { template class Manifold; template class FlatManifold; + template class ManifoldChart; +#if deal_II_space_dimension > 1 + template class ManifoldChart<1,deal_II_space_dimension>; +#endif +#if deal_II_space_dimension > 2 + template class ManifoldChart<2,deal_II_space_dimension>; +#endif + } diff --git a/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/source/grid/tria_boundary_lib.cc index 9bd26fb208..02975a3162 100644 --- a/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/source/grid/tria_boundary_lib.cc @@ -59,33 +59,33 @@ CylinderBoundary::get_axis_vector (const unsigned int axis) -template -Point -CylinderBoundary:: -get_new_point_on_line (const typename Triangulation::line_iterator &line) const -{ - // compute a proposed new point - const Point middle = StraightBoundary::get_new_point_on_line (line); - - // we then have to 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 +// Point +// CylinderBoundary:: +// get_new_point_on_line (const typename Triangulation::line_iterator &line) const +// { +// // compute a proposed new point +// const Point middle = StraightBoundary::get_new_point_on_line (line); + +// // we then have to 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 Point @@ -114,173 +114,168 @@ project_to_manifold (const Point middle) const -template<> -Point<3> -CylinderBoundary<3>:: -get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const -{ - const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad); - - // same algorithm as above - const unsigned int spacedim = 3; - - const Point vector_from_axis = (middle-point_on_axis) - - ((middle-point_on_axis) * direction) * direction; - 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<> -Point<3> -CylinderBoundary<2,3>:: -get_new_point_on_quad (const Triangulation<2,3>::quad_iterator &quad) const -{ - const Point<3> middle = StraightBoundary<2,3>::get_new_point_on_quad (quad); - - // same algorithm as above - const unsigned int spacedim = 3; - const Point vector_from_axis = (middle-point_on_axis) - - ((middle-point_on_axis) * direction) * direction; - 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 -Point -CylinderBoundary:: -get_new_point_on_quad (const typename Triangulation::quad_iterator &) const -{ - Assert (false, ExcImpossibleInDim(dim)); - return Point(); -} - - - -template -void -CylinderBoundary::get_intermediate_points_on_line ( - const typename Triangulation::line_iterator &line, - 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<> +// Point<3> +// CylinderBoundary<3>:: +// get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const +// { +// const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad); + +// // same algorithm as above +// const unsigned int spacedim = 3; + +// const Point vector_from_axis = (middle-point_on_axis) - +// ((middle-point_on_axis) * direction) * direction; +// 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 -void -CylinderBoundary::get_intermediate_points_between_points ( - const Point &v0, - const Point &v1, - std::vector > &points) const -{ - const unsigned int n=points.size(); - Assert(n>0, ExcInternalError()); +// template<> +// Point<3> +// CylinderBoundary<2,3>:: +// get_new_point_on_quad (const Triangulation<2,3>::quad_iterator &quad) const +// { +// const Point<3> middle = StraightBoundary<2,3>::get_new_point_on_quad (quad); + +// // same algorithm as above +// const unsigned int spacedim = 3; +// const Point vector_from_axis = (middle-point_on_axis) - +// ((middle-point_on_axis) * direction) * direction; +// 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); +// } - // Do a simple linear interpolation followed by projection, using the same - // algorithm as above - const std::vector > &line_points = this->get_line_support_points(n); - for (unsigned int i=0; i middle = (1-x)*v0 + x*v1; +// template +// Point +// CylinderBoundary:: +// get_new_point_on_quad (const typename Triangulation::quad_iterator &) const +// { +// Assert (false, ExcImpossibleInDim(dim)); +// return Point(); +// } - const Point vector_from_axis = (middle-point_on_axis) - - ((middle-point_on_axis) * direction) * direction; - if (vector_from_axis.norm() <= 1e-10 * middle.norm()) - points[i] = middle; - else - points[i] = (vector_from_axis / vector_from_axis.norm() * radius + - ((middle-point_on_axis) * direction) * direction + - point_on_axis); - } -} +// template +// void +// CylinderBoundary::get_intermediate_points_on_line ( +// const typename Triangulation::line_iterator &line, +// 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 -CylinderBoundary<3>::get_intermediate_points_on_quad ( - const Triangulation<3>::quad_iterator &quad, - std::vector > &points) const -{ - if (points.size()==1) - points[0]=get_new_point_on_quad(quad); - else - { - unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); - Assert(points.size()==m*m, ExcInternalError()); - std::vector > lp0(m); - std::vector > lp1(m); +// template +// void +// CylinderBoundary::get_intermediate_points_between_points ( +// const Point &v0, +// const Point &v1, +// std::vector > &points) const +// { +// const unsigned int n=points.size(); +// Assert(n>0, ExcInternalError()); + +// // Do a simple linear interpolation followed by projection, using the same +// // algorithm as above +// const std::vector > &line_points = this->get_line_support_points(n); + +// for (unsigned int i=0; i middle = (1-x)*v0 + x*v1; + +// const Point vector_from_axis = (middle-point_on_axis) - +// ((middle-point_on_axis) * direction) * direction; +// if (vector_from_axis.norm() <= 1e-10 * middle.norm()) +// points[i] = middle; +// else +// points[i] = (vector_from_axis / vector_from_axis.norm() * radius + +// ((middle-point_on_axis) * direction) * direction + +// point_on_axis); +// } +// } - get_intermediate_points_on_line(quad->line(0), lp0); - get_intermediate_points_on_line(quad->line(1), lp1); - std::vector > lps(m); - for (unsigned int i=0; i +// void +// CylinderBoundary<3>::get_intermediate_points_on_quad ( +// const Triangulation<3>::quad_iterator &quad, +// std::vector > &points) const +// { +// if (points.size()==1) +// points[0]=get_new_point_on_quad(quad); +// else +// { +// unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); +// Assert(points.size()==m*m, ExcInternalError()); + +// std::vector > lp0(m); +// std::vector > lp1(m); + +// get_intermediate_points_on_line(quad->line(0), lp0); +// get_intermediate_points_on_line(quad->line(1), lp1); + +// std::vector > lps(m); +// for (unsigned int i=0; i -void -CylinderBoundary::get_intermediate_points_on_quad ( - const typename Triangulation::quad_iterator &, - std::vector > &) const -{ - Assert (false, ExcImpossibleInDim(dim)); -} +// template +// void +// CylinderBoundary::get_intermediate_points_on_quad ( +// const typename Triangulation::quad_iterator &, +// std::vector > &) const +// { +// Assert (false, ExcImpossibleInDim(dim)); +// } -template <> -void -CylinderBoundary<1>:: -get_normals_at_vertices (const Triangulation<1>::face_iterator &, - Boundary<1,1>::FaceVertexNormals &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} +// template <> +// void +// CylinderBoundary<1>:: +// get_normals_at_vertices (const Triangulation<1>::face_iterator &, +// Boundary<1,1>::FaceVertexNormals &) const +// { +// Assert (false, ExcImpossibleInDim(1)); +// } template -void +Point CylinderBoundary:: -get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const +normal_vector (const Point vertex) const { - for (unsigned int v=0; v::vertices_per_face; ++v) - { - const Point vertex = face->vertex(v); - - const Point vector_from_axis = (vertex-point_on_axis) - - ((vertex-point_on_axis) * direction) * direction; + const Point vector_from_axis = (vertex-point_on_axis) - + ((vertex-point_on_axis) * direction) * direction; - face_vertex_normals[v] = (vector_from_axis / vector_from_axis.norm()); - } + Point normal = (vector_from_axis / vector_from_axis.norm()); + return normal; } @@ -321,51 +316,51 @@ double ConeBoundary::get_radius (Point x) const -template -void -ConeBoundary:: -get_intermediate_points_between_points (const Point &p0, - const Point &p1, - std::vector > &points) const -{ - const unsigned int n = points.size (); - const Point axis = x_1 - x_0; - - Assert (n > 0, ExcInternalError ()); - - const std::vector > &line_points = this->get_line_support_points(n); - - for (unsigned int i=0; i x_i = (1-x)*p0 + x*p1; - // To project this point on the boundary of the cone we first compute - // the orthogonal projection of this point onto the axis of the cone. - const double c = (x_i - x_0) * axis / axis.square (); - const Point x_ip = x_0 + c * axis; - // Compute the projection of the middle point on the boundary of the - // cone. - points[i] = x_ip + get_radius (x_ip) * (x_i - x_ip) / (x_i - x_ip).norm (); - } -} +// template +// void +// ConeBoundary:: +// get_intermediate_points_between_points (const Point &p0, +// const Point &p1, +// std::vector > &points) const +// { +// const unsigned int n = points.size (); +// const Point axis = x_1 - x_0; + +// Assert (n > 0, ExcInternalError ()); + +// const std::vector > &line_points = this->get_line_support_points(n); + +// for (unsigned int i=0; i x_i = (1-x)*p0 + x*p1; +// // To project this point on the boundary of the cone we first compute +// // the orthogonal projection of this point onto the axis of the cone. +// const double c = (x_i - x_0) * axis / axis.square (); +// const Point x_ip = x_0 + c * axis; +// // Compute the projection of the middle point on the boundary of the +// // cone. +// points[i] = x_ip + get_radius (x_ip) * (x_i - x_ip) / (x_i - x_ip).norm (); +// } +// } -template -Point -ConeBoundary:: -get_new_point_on_line (const typename Triangulation::line_iterator &line) const -{ - const Point axis = x_1 - x_0; - // Compute the middle point of the line. - const Point middle = StraightBoundary::get_new_point_on_line (line); - // 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 +// ConeBoundary:: +// get_new_point_on_line (const typename Triangulation::line_iterator &line) const +// { +// const Point axis = x_1 - x_0; +// // Compute the middle point of the line. +// const Point middle = StraightBoundary::get_new_point_on_line (line); +// // 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 @@ -384,133 +379,130 @@ project_to_manifold (const Point middle) const -template <> -Point<3> -ConeBoundary<3>:: -get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const -{ - const int dim = 3; - - const Point axis = x_1 - x_0; - // Compute the middle point of the quad. - const Point middle = StraightBoundary<3>::get_new_point_on_quad (quad); - // Same algorithm as above: 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> +// ConeBoundary<3>:: +// get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const +// { +// const int dim = 3; + +// const Point axis = x_1 - x_0; +// // Compute the middle point of the quad. +// const Point middle = StraightBoundary<3>::get_new_point_on_quad (quad); +// // Same algorithm as above: 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 -ConeBoundary:: -get_new_point_on_quad (const typename Triangulation::quad_iterator &) const -{ - Assert (false, ExcImpossibleInDim (dim)); +// template +// Point +// ConeBoundary:: +// get_new_point_on_quad (const typename Triangulation::quad_iterator &) const +// { +// Assert (false, ExcImpossibleInDim (dim)); - return Point(); -} +// return Point(); +// } -template -void -ConeBoundary:: -get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, - 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 +// ConeBoundary:: +// get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, +// 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 -ConeBoundary<3>:: -get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, - std::vector > &points) const -{ - if (points.size () == 1) - points[0] = get_new_point_on_quad (quad); - else - { - unsigned int n = static_cast (std::sqrt (static_cast (points.size ()))); +// template<> +// void +// ConeBoundary<3>:: +// get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, +// std::vector > &points) const +// { +// if (points.size () == 1) +// points[0] = get_new_point_on_quad (quad); +// else +// { +// unsigned int n = static_cast (std::sqrt (static_cast (points.size ()))); - Assert (points.size () == n * n, ExcInternalError ()); +// Assert (points.size () == n * n, ExcInternalError ()); - std::vector > points_line_0 (n); - std::vector > points_line_1 (n); +// std::vector > points_line_0 (n); +// std::vector > points_line_1 (n); - get_intermediate_points_on_line (quad->line (0), points_line_0); - get_intermediate_points_on_line (quad->line (1), points_line_1); +// get_intermediate_points_on_line (quad->line (0), points_line_0); +// get_intermediate_points_on_line (quad->line (1), points_line_1); - std::vector > points_line_segment (n); +// std::vector > points_line_segment (n); - for (unsigned int i = 0; i < n; ++i) - { - get_intermediate_points_between_points (points_line_0[i], - points_line_1[i], - points_line_segment); +// for (unsigned int i = 0; i < n; ++i) +// { +// get_intermediate_points_between_points (points_line_0[i], +// points_line_1[i], +// points_line_segment); - for (unsigned int j = 0; j < n; ++j) - points[i * n + j] = points_line_segment[j]; - } - } -} +// for (unsigned int j = 0; j < n; ++j) +// points[i * n + j] = points_line_segment[j]; +// } +// } +// } -template -void -ConeBoundary:: -get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &, - std::vector > &) const -{ - Assert (false, ExcImpossibleInDim (dim)); -} +// template +// void +// ConeBoundary:: +// get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &, +// std::vector > &) const +// { +// Assert (false, ExcImpossibleInDim (dim)); +// } -template<> -void -ConeBoundary<1>:: -get_normals_at_vertices (const Triangulation<1>::face_iterator &, - Boundary<1,1>::FaceVertexNormals &) const -{ - Assert (false, ExcImpossibleInDim (1)); -} +// template<> +// void +// ConeBoundary<1>:: +// get_normals_at_vertices (const Triangulation<1>::face_iterator &, +// Boundary<1,1>::FaceVertexNormals &) const +// { +// Assert (false, ExcImpossibleInDim (1)); +// } template -void +Point ConeBoundary:: -get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const +normal_vector (const Point vertex) const { const Point axis = x_1 - x_0; - for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_cell; ++vertex) - { - // Compute the orthogonal projection of the vertex onto the axis of the - // cone. - const double c = (face->vertex (vertex) - x_0) * axis / axis.square (); - const Point vertex_p = x_0 + c * axis; + // Compute the orthogonal projection of the vertex onto the axis of + // the cone. + const double c = (vertex - x_0) * axis / axis.square (); + const Point vertex_p = x_0 + c * axis; // Then compute the vector pointing from the point vertex_p on // the axis to the vertex. - const Point axis_to_vertex = face->vertex (vertex) - vertex_p; + const Point axis_to_vertex = vertex - vertex_p; - face_vertex_normals[vertex] = axis_to_vertex / axis_to_vertex.norm (); - } + Point normal = axis_to_vertex / axis_to_vertex.norm (); + return normal; } @@ -544,27 +536,27 @@ HyperBallBoundary::HyperBallBoundary (const Point p, // } -template -Point -HyperBallBoundary::get_new_point_on_line (const typename Triangulation::line_iterator &line) const -{ - Point middle = StraightBoundary::get_new_point_on_line (line); - - middle -= center; - - double r=0; - if (compute_radius_automatically) - { - const Point vertex_relative = line->vertex(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 +// HyperBallBoundary::get_new_point_on_line (const typename Triangulation::line_iterator &line) const +// { +// Point middle = StraightBoundary::get_new_point_on_line (line); + +// middle -= center; + +// double r=0; +// if (compute_radius_automatically) +// { +// const Point vertex_relative = line->vertex(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 @@ -590,257 +582,253 @@ HyperBallBoundary::get_new_point(const std::vector -template <> -Point<1> -HyperBallBoundary<1,1>:: -get_new_point_on_quad (const Triangulation<1,1>::quad_iterator &) const -{ - Assert (false, ExcInternalError()); - return Point<1>(); -} - - -template <> -Point<2> -HyperBallBoundary<1,2>:: -get_new_point_on_quad (const Triangulation<1,2>::quad_iterator &) const -{ - Assert (false, ExcInternalError()); - return Point<2>(); -} - - - -template -Point -HyperBallBoundary:: -get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const -{ - Point middle = StraightBoundary::get_new_point_on_quad (quad); - - middle -= center; - - double r=0; - if (compute_radius_automatically) - { - const Point vertex_relative = quad->vertex(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 <> -void -HyperBallBoundary<1>::get_intermediate_points_on_line ( - const Triangulation<1>::line_iterator &, - std::vector > &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} - - - -template -void -HyperBallBoundary::get_intermediate_points_on_line ( - const typename Triangulation::line_iterator &line, - 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, - std::vector > &points) const -{ - const unsigned int n=points.size(); - Assert(n>0, ExcInternalError()); - - 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 = p0 - center; - r = std::sqrt(vertex_relative.square()); - } - else - r=radius; - - - const double r2=r*r; - Assert(std::fabs(v0.square()-r2) +// Point<1> +// HyperBallBoundary<1,1>:: +// get_new_point_on_quad (const Triangulation<1,1>::quad_iterator &) const +// { +// Assert (false, ExcInternalError()); +// return Point<1>(); +// } - const double alpha=std::acos((v0*v1)/std::sqrt(v0.square()*v1.square())); - const Point pm=0.5*(v0+v1); - const double h=std::sqrt(pm.square()); +// template <> +// Point<2> +// HyperBallBoundary<1,2>:: +// get_new_point_on_quad (const Triangulation<1,2>::quad_iterator &) const +// { +// Assert (false, ExcInternalError()); +// return Point<2>(); +// } - // n even: m=n/2, - // n odd: m=(n-1)/2 - const std::vector > &line_points = this->get_line_support_points(n); - const unsigned int m=n/2; - for (unsigned int i=0; i +// Point +// HyperBallBoundary:: +// get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +// { +// Point middle = StraightBoundary::get_new_point_on_quad (quad); + +// middle -= center; + +// double r=0; +// if (compute_radius_automatically) +// { +// const Point vertex_relative = quad->vertex(0) - center; +// r = std::sqrt(vertex_relative.square()); +// } +// else +// r=radius; +// // project to boundary +// middle *= r / std::sqrt(middle.square()); + +// middle += center; +// return middle; +// } - // project the points from the straight line to the HyperBallBoundary - for (unsigned int i=0; i +// void +// HyperBallBoundary<1>::get_intermediate_points_on_line ( +// const Triangulation<1>::line_iterator &, +// std::vector > &) const +// { +// Assert (false, ExcImpossibleInDim(1)); +// } -template <> -void -HyperBallBoundary<3>::get_intermediate_points_on_quad ( - const Triangulation<3>::quad_iterator &quad, - std::vector > &points) const -{ - if (points.size()==1) - points[0]=get_new_point_on_quad(quad); - else - { - unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); - Assert(points.size()==m*m, ExcInternalError()); - std::vector > lp0(m); - std::vector > lp1(m); - get_intermediate_points_on_line(quad->line(0), lp0); - get_intermediate_points_on_line(quad->line(1), lp1); +// template +// void +// HyperBallBoundary::get_intermediate_points_on_line ( +// const typename Triangulation::line_iterator &line, +// 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); +// } - std::vector > lps(m); - for (unsigned int i=0; i +// void +// HyperBallBoundary::get_intermediate_points_between_points ( +// const Point &p0, const Point &p1, +// std::vector > &points) const +// { +// const unsigned int n=points.size(); +// Assert(n>0, ExcInternalError()); + +// 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 = p0 - center; +// r = std::sqrt(vertex_relative.square()); +// } +// else +// r=radius; + + +// const double r2=r*r; +// Assert(std::fabs(v0.square()-r2) pm=0.5*(v0+v1); + +// const double h=std::sqrt(pm.square()); + +// // n even: m=n/2, +// // n odd: m=(n-1)/2 +// const std::vector > &line_points = this->get_line_support_points(n); +// const unsigned int m=n/2; +// for (unsigned int i=0; i -void -HyperBallBoundary<2,3>::get_intermediate_points_on_quad ( - const Triangulation<2,3>::quad_iterator &quad, - std::vector > &points) const -{ - if (points.size()==1) - points[0]=get_new_point_on_quad(quad); - else - { - unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); - Assert(points.size()==m*m, ExcInternalError()); - std::vector > lp0(m); - std::vector > lp1(m); +// template <> +// void +// HyperBallBoundary<3>::get_intermediate_points_on_quad ( +// const Triangulation<3>::quad_iterator &quad, +// std::vector > &points) const +// { +// if (points.size()==1) +// points[0]=get_new_point_on_quad(quad); +// else +// { +// unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); +// Assert(points.size()==m*m, ExcInternalError()); + +// std::vector > lp0(m); +// std::vector > lp1(m); + +// get_intermediate_points_on_line(quad->line(0), lp0); +// get_intermediate_points_on_line(quad->line(1), lp1); + +// std::vector > lps(m); +// for (unsigned int i=0; iline(0), lp0); - get_intermediate_points_on_line(quad->line(1), lp1); - std::vector > lps(m); - for (unsigned int i=0; i +// void +// HyperBallBoundary<2,3>::get_intermediate_points_on_quad ( +// const Triangulation<2,3>::quad_iterator &quad, +// std::vector > &points) const +// { +// if (points.size()==1) +// points[0]=get_new_point_on_quad(quad); +// else +// { +// unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); +// Assert(points.size()==m*m, ExcInternalError()); + +// std::vector > lp0(m); +// std::vector > lp1(m); + +// get_intermediate_points_on_line(quad->line(0), lp0); +// get_intermediate_points_on_line(quad->line(1), lp1); + +// std::vector > lps(m); +// for (unsigned int i=0; i -void -HyperBallBoundary::get_intermediate_points_on_quad ( - const typename Triangulation::quad_iterator &, - std::vector > &) const -{ - Assert(false, ExcImpossibleInDim(dim)); -} +// template +// void +// HyperBallBoundary::get_intermediate_points_on_quad ( +// const typename Triangulation::quad_iterator &, +// std::vector > &) const +// { +// Assert(false, ExcImpossibleInDim(dim)); +// } -template -Tensor<1,spacedim> -HyperBallBoundary:: -normal_vector (const typename Triangulation::face_iterator &, - const Point &p) const -{ - const Tensor<1,spacedim> unnormalized_normal = p-center; - return unnormalized_normal/unnormalized_normal.norm(); -} +// template +// Tensor<1,spacedim> +// HyperBallBoundary:: +// normal_vector (const typename Triangulation::face_iterator &, +// const Point &p) const +// { +// const Tensor<1,spacedim> unnormalized_normal = p-center; +// return unnormalized_normal/unnormalized_normal.norm(); +// } -template <> -void -HyperBallBoundary<1>:: -get_normals_at_vertices (const Triangulation<1>::face_iterator &, - Boundary<1,1>::FaceVertexNormals &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} +// template <> +// void +// HyperBallBoundary<1>:: +// get_normals_at_vertices (const Triangulation<1>::face_iterator &, +// Boundary<1,1>::FaceVertexNormals &) const +// { +// Assert (false, ExcImpossibleInDim(1)); +// } -template <> -void -HyperBallBoundary<1,2>:: -get_normals_at_vertices (const Triangulation<1,2>::face_iterator &, - Boundary<1,2>::FaceVertexNormals &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} +// template <> +// void +// HyperBallBoundary<1,2>:: +// get_normals_at_vertices (const Triangulation<1,2>::face_iterator &, +// Boundary<1,2>::FaceVertexNormals &) const +// { +// Assert (false, ExcImpossibleInDim(1)); +// } template -void +Point HyperBallBoundary:: -get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const +normal_vector (const Point vertex) const { - for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) - face_vertex_normals[vertex] = face->vertex(vertex)-center; + return vertex-center; } - - template Point HyperBallBoundary::get_center () const @@ -869,13 +857,13 @@ 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 +// { +// Assert(false, ExcInternalError()); +// return Point(); +// } template @@ -883,117 +871,157 @@ Point 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 + // check whether every point is at distance R + bool all_at_distance_R = true; + for(unsigned int i=0; icenter)-this->radius) >1e-5) + { + all_at_distance_R = false; + break; + } + const Point object_center = FlatManifold::get_new_point(surrounding_points, weights); - if (object_center.distance(this->center) < 1e-5) + if (std::abs(object_center(0)-this->center(0)) < 1e-5 && !all_at_distance_R) return object_center; else - return HyperBallBoundary::get_new_point (surrounding_points,weights); + if(all_at_distance_R) + return HyperBallBoundary::get_new_point (surrounding_points,weights); + else + { + Assert(false, ExcInternalError()); + return Point(); + } } -template <> -Point<1> -HalfHyperBallBoundary<1>:: -get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const -{ - Assert (false, ExcInternalError()); - return Point<1>(); -} +// template <> +// Point<1> +// HalfHyperBallBoundary<1>:: +// get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const +// { +// Assert (false, ExcInternalError()); +// return Point<1>(); +// } -template -Point -HalfHyperBallBoundary:: -get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const -{ - const Point quad_center = quad->center(); - if (quad_center(0) == this->center(0)) - return quad_center; - else - return HyperBallBoundary::get_new_point_on_quad (quad); -} +// template +// Point +// HalfHyperBallBoundary:: +// get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +// { +// const Point quad_center = quad->center(); +// if (quad_center(0) == this->center(0)) +// return quad_center; +// else +// return HyperBallBoundary::get_new_point_on_quad (quad); +// } -template -void -HalfHyperBallBoundary:: -get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, - std::vector > &points) const -{ - // check whether center of object is at x==0, since then it belongs to the - // plane part of the boundary - const Point line_center = line->center(); - if (line_center(0) == this->center(0)) - return StraightBoundary::get_intermediate_points_on_line (line, points); - else - return HyperBallBoundary::get_intermediate_points_on_line (line, points); -} +// template +// void +// HalfHyperBallBoundary:: +// get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, +// std::vector > &points) const +// { +// // check whether center of object is at x==0, since then it belongs to the +// // plane part of the boundary +// const Point line_center = line->center(); +// if (line_center(0) == this->center(0)) +// return StraightBoundary::get_intermediate_points_on_line (line, points); +// else +// return HyperBallBoundary::get_intermediate_points_on_line (line, points); +// } -template -void -HalfHyperBallBoundary:: -get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, - std::vector > &points) const -{ - if (points.size()==1) - points[0]=get_new_point_on_quad(quad); - else - { - // check whether center of object is at x==0, since then it belongs to - // the plane part of the boundary - const Point quad_center = quad->center(); - if (quad_center(0) == this->center(0)) - StraightBoundary::get_intermediate_points_on_quad (quad, points); - else - HyperBallBoundary::get_intermediate_points_on_quad (quad, points); - } -} +// template +// void +// HalfHyperBallBoundary:: +// get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, +// std::vector > &points) const +// { +// if (points.size()==1) +// points[0]=get_new_point_on_quad(quad); +// else +// { +// // check whether center of object is at x==0, since then it belongs to +// // the plane part of the boundary +// const Point quad_center = quad->center(); +// if (quad_center(0) == this->center(0)) +// StraightBoundary::get_intermediate_points_on_quad (quad, points); +// else +// HyperBallBoundary::get_intermediate_points_on_quad (quad, points); +// } +// } -template <> -void -HalfHyperBallBoundary<1>:: -get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &, - std::vector > &) const -{ - Assert (false, ExcInternalError()); -} +// template <> +// void +// HalfHyperBallBoundary<1>:: +// get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &, +// std::vector > &) const +// { +// Assert (false, ExcInternalError()); +// } -template <> -void -HalfHyperBallBoundary<1>:: -get_normals_at_vertices (const Triangulation<1>::face_iterator &, - Boundary<1,1>::FaceVertexNormals &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} +// template <> +// void +// HalfHyperBallBoundary<1>:: +// get_normals_at_vertices (const Triangulation<1>::face_iterator &, +// Boundary<1,1>::FaceVertexNormals &) const +// { +// Assert (false, ExcImpossibleInDim(1)); +// } +// template +// void +// HalfHyperBallBoundary:: +// get_normals_at_vertices (const typename Triangulation::face_iterator &face, +// typename Boundary::FaceVertexNormals &face_vertex_normals) const +// { +// // check whether center of object is at x==0, since then it belongs to the +// // plane part of the boundary +// const Point quad_center = face->center(); +// if (quad_center(0) == this->center(0)) +// StraightBoundary::get_normals_at_vertices (face, face_vertex_normals); +// else +// HyperBallBoundary::get_normals_at_vertices (face, face_vertex_normals); +// } + + template -void +Point HalfHyperBallBoundary:: -get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const +normal_vector (const Point vertex) const { - // check whether center of object is at x==0, since then it belongs to the - // plane part of the boundary - const Point quad_center = face->center(); - if (quad_center(0) == this->center(0)) - StraightBoundary::get_normals_at_vertices (face, face_vertex_normals); + // Check whether the point is inside the flat part, or outside. If + // it is both at distance R from the center and with zero first + // component, throw an exception, since we are on a edge, and normal + // is not well defined. + bool at_distance_R = (std::abs(vertex.distance(this->center)-this->radius)<1e-10); + + if (abs(vertex(0)) < 1e-10 && !at_distance_R) + { + Point normal; + normal[0] = -1; + return normal; + } + else if(at_distance_R) + return vertex-this->center; else - HyperBallBoundary::get_normals_at_vertices (face, face_vertex_normals); + { + Assert(false, + ExcMessage("Normals at edge is not well defined. Use two distinct manifolds!")); + return Point(); + } } @@ -1033,241 +1061,257 @@ HalfHyperShellBoundary::HalfHyperShellBoundary (const Point ¢er, -template -Point -HalfHyperShellBoundary:: -get_new_point_on_line (const typename Triangulation::line_iterator &line) const -{ - switch (dim) - { - // in 2d, first check whether the two end points of the line are on the - // axis of symmetry. if so, then return the mid point - case 2: - { - if ((line->vertex(0)(0) == this->center(0)) - && - (line->vertex(1)(0) == this->center(0))) - return (line->vertex(0) + line->vertex(1))/2; - else - // otherwise we are on the outer or inner part of the shell. proceed - // as in the base class - return HyperShellBoundary::get_new_point_on_line (line); - } - - // in 3d, a line is a straight line if it is on the symmetry plane and if - // not both of its end points are on either the inner or outer sphere - case 3: - { - - if (((line->vertex(0)(0) == this->center(0)) - && - (line->vertex(1)(0) == this->center(0))) - && - !(((std::fabs (line->vertex(0).distance (this->center) - - inner_radius) < 1e-12 * outer_radius) - && - (std::fabs (line->vertex(1).distance (this->center) - - inner_radius) < 1e-12 * outer_radius)) - || - ((std::fabs (line->vertex(0).distance (this->center) - - outer_radius) < 1e-12 * outer_radius) - && - (std::fabs (line->vertex(1).distance (this->center) - - outer_radius) < 1e-12 * outer_radius)))) - return (line->vertex(0) + line->vertex(1))/2; - else - // otherwise we are on the outer or inner part of the shell. proceed - // as in the base class - return HyperShellBoundary::get_new_point_on_line (line); - } - - default: - Assert (false, ExcNotImplemented()); - } - - return Point(); -} - - - -template <> -Point<1> -HalfHyperShellBoundary<1>:: -get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const -{ - Assert (false, ExcInternalError()); - return Point<1>(); -} - +// template +// Point +// HalfHyperShellBoundary:: +// get_new_point_on_line (const typename Triangulation::line_iterator &line) const +// { +// switch (dim) +// { +// // in 2d, first check whether the two end points of the line are on the +// // axis of symmetry. if so, then return the mid point +// case 2: +// { +// if ((line->vertex(0)(0) == this->center(0)) +// && +// (line->vertex(1)(0) == this->center(0))) +// return (line->vertex(0) + line->vertex(1))/2; +// else +// // otherwise we are on the outer or inner part of the shell. proceed +// // as in the base class +// return HyperShellBoundary::get_new_point_on_line (line); +// } + +// // in 3d, a line is a straight line if it is on the symmetry plane and if +// // not both of its end points are on either the inner or outer sphere +// case 3: +// { + +// if (((line->vertex(0)(0) == this->center(0)) +// && +// (line->vertex(1)(0) == this->center(0))) +// && +// !(((std::fabs (line->vertex(0).distance (this->center) +// - inner_radius) < 1e-12 * outer_radius) +// && +// (std::fabs (line->vertex(1).distance (this->center) +// - inner_radius) < 1e-12 * outer_radius)) +// || +// ((std::fabs (line->vertex(0).distance (this->center) +// - outer_radius) < 1e-12 * outer_radius) +// && +// (std::fabs (line->vertex(1).distance (this->center) +// - outer_radius) < 1e-12 * outer_radius)))) +// return (line->vertex(0) + line->vertex(1))/2; +// else +// // otherwise we are on the outer or inner part of the shell. proceed +// // as in the base class +// return HyperShellBoundary::get_new_point_on_line (line); +// } + +// default: +// Assert (false, ExcNotImplemented()); +// } + +// return Point(); +// } -template -Point -HalfHyperShellBoundary:: -get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const -{ - // if this quad is on the symmetry plane, take the center point and project - // it outward to the same radius as the centers of the two radial lines - if ((quad->vertex(0)(0) == this->center(0)) && - (quad->vertex(1)(0) == this->center(0)) && - (quad->vertex(2)(0) == this->center(0)) && - (quad->vertex(3)(0) == this->center(0))) - { - const Point quad_center = (quad->vertex(0) + quad->vertex(1) + - quad->vertex(2) + quad->vertex(3) )/4; - const Point quad_center_offset = quad_center - this->center; - - - if (std::fabs (quad->line(0)->center().distance(this->center) - - quad->line(1)->center().distance(this->center)) - < 1e-12 * outer_radius) - { - // lines 0 and 1 are radial - const double needed_radius - = quad->line(0)->center().distance(this->center); - - return (this->center + - quad_center_offset/quad_center_offset.norm() * needed_radius); - } - else if (std::fabs (quad->line(2)->center().distance(this->center) - - quad->line(3)->center().distance(this->center)) - < 1e-12 * outer_radius) - { - // lines 2 and 3 are radial - const double needed_radius - = quad->line(2)->center().distance(this->center); - - return (this->center + - quad_center_offset/quad_center_offset.norm() * needed_radius); - } - else - Assert (false, ExcInternalError()); - } +// template <> +// Point<1> +// HalfHyperShellBoundary<1>:: +// get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const +// { +// Assert (false, ExcInternalError()); +// return Point<1>(); +// } - // otherwise we are on the outer or inner part of the shell. proceed as in - // the base class - return HyperShellBoundary::get_new_point_on_quad (quad); -} -template -void -HalfHyperShellBoundary:: -get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, - std::vector > &points) const -{ - switch (dim) - { - // in 2d, first check whether the two end points of the line are on the - // axis of symmetry. if so, then return the mid point - case 2: - { - if ((line->vertex(0)(0) == this->center(0)) - && - (line->vertex(1)(0) == this->center(0))) - StraightBoundary::get_intermediate_points_on_line (line, points); - else - // otherwise we are on the outer or inner part of the shell. proceed - // as in the base class - HyperShellBoundary::get_intermediate_points_on_line (line, points); - break; - } +// template +// Point +// HalfHyperShellBoundary:: +// get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +// { +// // if this quad is on the symmetry plane, take the center point and project +// // it outward to the same radius as the centers of the two radial lines +// if ((quad->vertex(0)(0) == this->center(0)) && +// (quad->vertex(1)(0) == this->center(0)) && +// (quad->vertex(2)(0) == this->center(0)) && +// (quad->vertex(3)(0) == this->center(0))) +// { +// const Point quad_center = (quad->vertex(0) + quad->vertex(1) + +// quad->vertex(2) + quad->vertex(3) )/4; +// const Point quad_center_offset = quad_center - this->center; + + +// if (std::fabs (quad->line(0)->center().distance(this->center) - +// quad->line(1)->center().distance(this->center)) +// < 1e-12 * outer_radius) +// { +// // lines 0 and 1 are radial +// const double needed_radius +// = quad->line(0)->center().distance(this->center); + +// return (this->center + +// quad_center_offset/quad_center_offset.norm() * needed_radius); +// } +// else if (std::fabs (quad->line(2)->center().distance(this->center) - +// quad->line(3)->center().distance(this->center)) +// < 1e-12 * outer_radius) +// { +// // lines 2 and 3 are radial +// const double needed_radius +// = quad->line(2)->center().distance(this->center); + +// return (this->center + +// quad_center_offset/quad_center_offset.norm() * needed_radius); +// } +// else +// Assert (false, ExcInternalError()); +// } + +// // otherwise we are on the outer or inner part of the shell. proceed as in +// // the base class +// return HyperShellBoundary::get_new_point_on_quad (quad); +// } - // in 3d, a line is a straight line if it is on the symmetry plane and if - // not both of its end points are on either the inner or outer sphere - case 3: - { - if (((line->vertex(0)(0) == this->center(0)) - && - (line->vertex(1)(0) == this->center(0))) - && - !(((std::fabs (line->vertex(0).distance (this->center) - - inner_radius) < 1e-12 * outer_radius) - && - (std::fabs (line->vertex(1).distance (this->center) - - inner_radius) < 1e-12 * outer_radius)) - || - ((std::fabs (line->vertex(0).distance (this->center) - - outer_radius) < 1e-12 * outer_radius) - && - (std::fabs (line->vertex(1).distance (this->center) - - outer_radius) < 1e-12 * outer_radius)))) - StraightBoundary::get_intermediate_points_on_line (line, points); - else - // otherwise we are on the outer or inner part of the shell. proceed - // as in the base class - HyperShellBoundary::get_intermediate_points_on_line (line, points); - break; - } - default: - Assert (false, ExcNotImplemented()); - } -} +// template +// void +// HalfHyperShellBoundary:: +// get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, +// std::vector > &points) const +// { +// switch (dim) +// { +// // in 2d, first check whether the two end points of the line are on the +// // axis of symmetry. if so, then return the mid point +// case 2: +// { +// if ((line->vertex(0)(0) == this->center(0)) +// && +// (line->vertex(1)(0) == this->center(0))) +// StraightBoundary::get_intermediate_points_on_line (line, points); +// else +// // otherwise we are on the outer or inner part of the shell. proceed +// // as in the base class +// HyperShellBoundary::get_intermediate_points_on_line (line, points); +// break; +// } + +// // in 3d, a line is a straight line if it is on the symmetry plane and if +// // not both of its end points are on either the inner or outer sphere +// case 3: +// { +// if (((line->vertex(0)(0) == this->center(0)) +// && +// (line->vertex(1)(0) == this->center(0))) +// && +// !(((std::fabs (line->vertex(0).distance (this->center) +// - inner_radius) < 1e-12 * outer_radius) +// && +// (std::fabs (line->vertex(1).distance (this->center) +// - inner_radius) < 1e-12 * outer_radius)) +// || +// ((std::fabs (line->vertex(0).distance (this->center) +// - outer_radius) < 1e-12 * outer_radius) +// && +// (std::fabs (line->vertex(1).distance (this->center) +// - outer_radius) < 1e-12 * outer_radius)))) +// StraightBoundary::get_intermediate_points_on_line (line, points); +// else +// // otherwise we are on the outer or inner part of the shell. proceed +// // as in the base class +// HyperShellBoundary::get_intermediate_points_on_line (line, points); + +// break; +// } + +// default: +// Assert (false, ExcNotImplemented()); +// } +// } -template -void -HalfHyperShellBoundary:: -get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, - std::vector > &points) const -{ - Assert (dim < 3, ExcNotImplemented()); +// template +// void +// HalfHyperShellBoundary:: +// get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, +// std::vector > &points) const +// { +// Assert (dim < 3, ExcNotImplemented()); + +// // check whether center of object is at x==0, since then it belongs to the +// // plane part of the boundary +// const Point quad_center = quad->center(); +// if (quad_center(0) == this->center(0)) +// StraightBoundary::get_intermediate_points_on_quad (quad, points); +// else +// HyperShellBoundary::get_intermediate_points_on_quad (quad, points); +// } - // check whether center of object is at x==0, since then it belongs to the - // plane part of the boundary - const Point quad_center = quad->center(); - if (quad_center(0) == this->center(0)) - StraightBoundary::get_intermediate_points_on_quad (quad, points); - else - HyperShellBoundary::get_intermediate_points_on_quad (quad, points); -} +// template <> +// void +// HalfHyperShellBoundary<1>:: +// get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &, +// std::vector > &) const +// { +// Assert (false, ExcInternalError()); +// } -template <> -void -HalfHyperShellBoundary<1>:: -get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &, - std::vector > &) const -{ - Assert (false, ExcInternalError()); -} +// template <> +// void +// HalfHyperShellBoundary<1>:: +// get_normals_at_vertices (const Triangulation<1>::face_iterator &, +// Boundary<1,1>::FaceVertexNormals &) const +// { +// Assert (false, ExcImpossibleInDim(1)); +// } -template <> -void -HalfHyperShellBoundary<1>:: -get_normals_at_vertices (const Triangulation<1>::face_iterator &, - Boundary<1,1>::FaceVertexNormals &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} +// template +// void +// HalfHyperShellBoundary:: +// get_normals_at_vertices (const typename Triangulation::face_iterator &face, +// typename Boundary::FaceVertexNormals &face_vertex_normals) const +// { +// if (face->center()(0) == this->center(0)) +// StraightBoundary::get_normals_at_vertices (face, face_vertex_normals); +// else +// HyperShellBoundary::get_normals_at_vertices (face, face_vertex_normals); +// } template -void +Point HalfHyperShellBoundary:: -get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const +normal_vector (const Point vertex) const { - if (face->center()(0) == this->center(0)) - StraightBoundary::get_normals_at_vertices (face, face_vertex_normals); + if (vertex(0) == this->center(0)) + { + Point n; + n[0] = -1; + return n; + } else - HyperShellBoundary::get_normals_at_vertices (face, face_vertex_normals); + return HyperShellBoundary::normal_vector (vertex); } + template TorusBoundary::TorusBoundary (const double R__, const double r__) @@ -1316,7 +1360,7 @@ TorusBoundary::get_correct_angle(const double angle, template <> Point<3> -TorusBoundary<2,3>::get_real_coord (const Point<2> &surfP) const +TorusBoundary<2,3>::push_forward (const Point<2> surfP) const { const double theta=surfP(0); const double phi=surfP(1); @@ -1330,7 +1374,7 @@ TorusBoundary<2,3>::get_real_coord (const Point<2> &surfP) const template <> Point<2> -TorusBoundary<2,3>::get_surf_coord(const Point<3> &p) const +TorusBoundary<2,3>::pull_back(const Point<3> p) const { const double phi=std::asin(std::abs(p(1))/r); const double Rr_2=p(0)*p(0)+p(2)*p(2); @@ -1356,63 +1400,63 @@ TorusBoundary<2,3>::get_surf_coord(const Point<3> &p) const -template <> -Point<3> -TorusBoundary<2,3>::get_new_point_on_line (const Triangulation<2,3>::line_iterator &line) const -{ - //Just get the average - Point<2> p0=get_surf_coord(line->vertex(0)); - Point<2> p1=get_surf_coord(line->vertex(1)); +// template <> +// Point<3> +// TorusBoundary<2,3>::get_new_point_on_line (const Triangulation<2,3>::line_iterator &line) const +// { +// //Just get the average +// Point<2> p0=get_surf_coord(line->vertex(0)); +// Point<2> p1=get_surf_coord(line->vertex(1)); - Point<2> middle(0,0); +// Point<2> middle(0,0); - //Take care for periodic conditions, For instance phi0= 0, phi1= 3/2*Pi - //middle has to be 7/4*Pi not 3/4*Pi. This also works for -Pi/2 + Pi, middle - //is 5/4*Pi - for (unsigned int i=0; i<2; i++) - if (std::abs(p0(i)-p1(i))> numbers::PI) - middle(i)=2*numbers::PI; +// //Take care for periodic conditions, For instance phi0= 0, phi1= 3/2*Pi +// //middle has to be 7/4*Pi not 3/4*Pi. This also works for -Pi/2 + Pi, middle +// //is 5/4*Pi +// for (unsigned int i=0; i<2; i++) +// if (std::abs(p0(i)-p1(i))> numbers::PI) +// middle(i)=2*numbers::PI; - middle+= p0 + p1; - middle*=0.5; +// middle+= p0 + p1; +// middle*=0.5; - Point<3> midReal=get_real_coord(middle); - return midReal; -} +// Point<3> midReal=get_real_coord(middle); +// return midReal; +// } -template <> -Point<3> -TorusBoundary<2,3>::get_new_point_on_quad (const Triangulation<2,3>::quad_iterator &quad) const -{ - //Just get the average - Point<2> p[4]; +// template <> +// Point<3> +// TorusBoundary<2,3>::get_new_point_on_quad (const Triangulation<2,3>::quad_iterator &quad) const +// { +// //Just get the average +// Point<2> p[4]; - for (unsigned int i=0; i<4; i++) - p[i]=get_surf_coord(quad->vertex(i)); +// for (unsigned int i=0; i<4; i++) +// p[i]=get_surf_coord(quad->vertex(i)); - Point<2> middle(0,0); +// Point<2> middle(0,0); - //Take care for periodic conditions, see get_new_point_on_line() above - //For instance phi0= 0, phi1= 3/2*Pi middle has to be 7/4*Pi not 3/4*Pi - //This also works for -Pi/2 + Pi + Pi- Pi/2, middle is 5/4*Pi - for (unsigned int i=0; i<2; i++) - for (unsigned int j=1; j<4; j++) - { - if (std::abs(p[0](i)-p[j](i))> numbers::PI) - { - middle(i)+=2*numbers::PI; - } - } +// //Take care for periodic conditions, see get_new_point_on_line() above +// //For instance phi0= 0, phi1= 3/2*Pi middle has to be 7/4*Pi not 3/4*Pi +// //This also works for -Pi/2 + Pi + Pi- Pi/2, middle is 5/4*Pi +// for (unsigned int i=0; i<2; i++) +// for (unsigned int j=1; j<4; j++) +// { +// if (std::abs(p[0](i)-p[j](i))> numbers::PI) +// { +// middle(i)+=2*numbers::PI; +// } +// } - for (unsigned int i=0; i<4; i++) - middle+=p[i]; +// for (unsigned int i=0; i<4; i++) +// middle+=p[i]; - middle*= 0.25; +// middle*= 0.25; - return get_real_coord(middle); -} +// return get_real_coord(middle); +// } @@ -1440,136 +1484,136 @@ TorusBoundary<2,3>:: get_surf_norm_from_sp(const Point<2> &surfP) const //Normal field without unit length template <> Point<3> -TorusBoundary<2,3>::get_surf_norm(const Point<3> &p) const +TorusBoundary<2,3>::normal_vector(const Point<3> p) const { - Point<2> surfP=get_surf_coord(p); + Point<2> surfP=pull_back(p); return get_surf_norm_from_sp(surfP); } -template<> -void -TorusBoundary<2,3>:: -get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator &line, - std::vector< Point< 3 > > &points) const -{ - //Almost the same implementation as StraightBoundary<2,3> - unsigned int npoints=points.size(); - if (npoints==0) return; - - Point<2> p[2]; - - for (unsigned int i=0; i<2; i++) - p[i]=get_surf_coord(line->vertex(i)); - - unsigned int offset[2]; - offset[0]=0; - offset[1]=0; - - //Take care for periodic conditions & negative angles, see - //get_new_point_on_line() above. Because we dont have a symmetric - //interpolation (just the middle) we need to add 2*Pi to each almost zero - //and negative angles. - for (unsigned int i=0; i<2; i++) - for (unsigned int j=1; j<2; j++) - { - if (std::abs(p[0](i)-p[j](i))> numbers::PI) - { - offset[i]++; - break; - } - } - - for (unsigned int i=0; i<2; i++) - for (unsigned int j=0; j<2; j++) - if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles - p[j](i)+=2*numbers::PI*offset[i]; - - - Point<2> target; - const std::vector > &line_points = this->get_line_support_points(npoints); - for (unsigned int i=0; i +// void +// TorusBoundary<2,3>:: +// get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator &line, +// std::vector< Point< 3 > > &points) const +// { +// //Almost the same implementation as StraightBoundary<2,3> +// unsigned int npoints=points.size(); +// if (npoints==0) return; + +// Point<2> p[2]; + +// for (unsigned int i=0; i<2; i++) +// p[i]=get_surf_coord(line->vertex(i)); + +// unsigned int offset[2]; +// offset[0]=0; +// offset[1]=0; + +// //Take care for periodic conditions & negative angles, see +// //get_new_point_on_line() above. Because we dont have a symmetric +// //interpolation (just the middle) we need to add 2*Pi to each almost zero +// //and negative angles. +// for (unsigned int i=0; i<2; i++) +// for (unsigned int j=1; j<2; j++) +// { +// if (std::abs(p[0](i)-p[j](i))> numbers::PI) +// { +// offset[i]++; +// break; +// } +// } + +// for (unsigned int i=0; i<2; i++) +// for (unsigned int j=0; j<2; j++) +// if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles +// p[j](i)+=2*numbers::PI*offset[i]; + + +// Point<2> target; +// const std::vector > &line_points = this->get_line_support_points(npoints); +// for (unsigned int i=0; i -void -TorusBoundary<2,3>:: -get_intermediate_points_on_quad (const Triangulation< 2, 3 >::quad_iterator &quad, - std::vector< Point< 3 > > &points )const -{ - //Almost the same implementation as StraightBoundary<2,3> - const unsigned int n=points.size(), - m=static_cast(std::sqrt(static_cast(n))); - // is n a square number - Assert(m*m==n, ExcInternalError()); - - Point<2> p[4]; - - for (unsigned int i=0; i<4; i++) - p[i]=get_surf_coord(quad->vertex(i)); - - Point<2> target; - unsigned int offset[2]; - offset[0]=0; - offset[1]=0; - - //Take care for periodic conditions & negative angles, see - //get_new_point_on_line() above. Because we dont have a symmetric - //interpolation (just the middle) we need to add 2*Pi to each almost zero - //and negative angles. - for (unsigned int i=0; i<2; i++) - for (unsigned int j=1; j<4; j++) - { - if (std::abs(p[0](i)-p[j](i))> numbers::PI) - { - offset[i]++; - break; - } - } - for (unsigned int i=0; i<2; i++) - for (unsigned int j=0; j<4; j++) - if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles - p[j](i)+=2*numbers::PI*offset[i]; - const std::vector > &line_points = this->get_line_support_points(m); - for (unsigned int i=0; i +// void +// TorusBoundary<2,3>:: +// get_intermediate_points_on_quad (const Triangulation< 2, 3 >::quad_iterator &quad, +// std::vector< Point< 3 > > &points )const +// { +// //Almost the same implementation as StraightBoundary<2,3> +// const unsigned int n=points.size(), +// m=static_cast(std::sqrt(static_cast(n))); +// // is n a square number +// Assert(m*m==n, ExcInternalError()); + +// Point<2> p[4]; + +// for (unsigned int i=0; i<4; i++) +// p[i]=get_surf_coord(quad->vertex(i)); + +// Point<2> target; +// unsigned int offset[2]; +// offset[0]=0; +// offset[1]=0; + +// //Take care for periodic conditions & negative angles, see +// //get_new_point_on_line() above. Because we dont have a symmetric +// //interpolation (just the middle) we need to add 2*Pi to each almost zero +// //and negative angles. +// for (unsigned int i=0; i<2; i++) +// for (unsigned int j=1; j<4; j++) +// { +// if (std::abs(p[0](i)-p[j](i))> numbers::PI) +// { +// offset[i]++; +// break; +// } +// } + +// for (unsigned int i=0; i<2; i++) +// for (unsigned int j=0; j<4; j++) +// if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles +// p[j](i)+=2*numbers::PI*offset[i]; + +// const std::vector > &line_points = this->get_line_support_points(m); +// for (unsigned int i=0; i -void -TorusBoundary<2,3>:: -get_normals_at_vertices (const Triangulation<2,3 >::face_iterator &face, - Boundary<2,3>::FaceVertexNormals &face_vertex_normals) const -{ - for (unsigned int i=0; i::vertices_per_face; i++) - face_vertex_normals[i]=get_surf_norm(face->vertex(i)); -} +// template<> +// void +// TorusBoundary<2,3>:: +// get_normals_at_vertices (const Triangulation<2,3 >::face_iterator &face, +// Boundary<2,3>::FaceVertexNormals &face_vertex_normals) const +// { +// for (unsigned int i=0; i::vertices_per_face; i++) +// face_vertex_normals[i]=get_surf_norm(face->vertex(i)); +// } -- 2.39.5