From 105e1da69e02bd72d1826bc518aba682b6d8da51 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Mon, 10 Feb 2014 15:28:44 +0000 Subject: [PATCH] Compiling intermediate branch manifold id git-svn-id: https://svn.dealii.org/branches/branch_manifold_id_intermediate@32447 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/grid/tria_accessor.templates.h | 4 +- .../include/deal.II/grid/tria_boundary_lib.h | 614 ++++++++++++++---- .../deal.II/numerics/vector_tools.templates.h | 6 +- deal.II/source/grid/grid_tools.cc | 30 +- deal.II/source/grid/tria.cc | 14 +- deal.II/source/grid/tria_accessor.cc | 1 - 6 files changed, 511 insertions(+), 158 deletions(-) diff --git a/deal.II/include/deal.II/grid/tria_accessor.templates.h b/deal.II/include/deal.II/grid/tria_accessor.templates.h index 7e9cffa7e3..6dcb43c052 100644 --- a/deal.II/include/deal.II/grid/tria_accessor.templates.h +++ b/deal.II/include/deal.II/grid/tria_accessor.templates.h @@ -1936,7 +1936,7 @@ TriaAccessor::at_boundary () const template -const Manifold & +const Boundary & TriaAccessor::get_boundary () const { return get_manifold(); @@ -1944,7 +1944,7 @@ TriaAccessor::get_boundary () const template -const Manifold & +const Boundary & TriaAccessor::get_manifold () const { Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed()); 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 8daf006889..2f683b964e 100644 --- a/deal.II/include/deal.II/grid/tria_boundary_lib.h +++ b/deal.II/include/deal.II/grid/tria_boundary_lib.h @@ -20,8 +20,6 @@ #include #include -#include -#include DEAL_II_NAMESPACE_OPEN @@ -50,7 +48,7 @@ DEAL_II_NAMESPACE_OPEN * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007 */ template -class CylinderBoundary : public FlatManifold +class CylinderBoundary : public StraightBoundary { public: /** @@ -62,41 +60,87 @@ 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); /** - * 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. + * Refer to the general documentation of + * this class and the documentation of the + * base class. */ - virtual Point - project_to_manifold (const std::vector > & vertices, - const Point &candidate) const; - + virtual Point + get_new_point_on_line (const typename Triangulation::line_iterator &line) const; /** - * Compute the normal to the surface. + * Refer to the general documentation of + * this class and the documentation of the + * base class. */ virtual Point - normal_vector (const std::vector > & vertices, - const Point &point) const; - + 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. + */ + virtual void + get_normals_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) 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); @@ -120,8 +164,24 @@ protected: private: /** - * Given a number for the axis, return a vector that denotes this - * direction. + * 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. */ static Point get_axis_vector (const unsigned int axis); }; @@ -162,15 +222,21 @@ private: * @author Markus Bürg, 2009 */ template -class ConeBoundary : public FlatManifold +class ConeBoundary : public StraightBoundary { 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, @@ -178,31 +244,75 @@ 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; + + /** + * 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; + /** - * 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. + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + * + * Calls @p + * get_intermediate_points_between_points. */ - virtual Point - project_to_manifold (const std::vector > & vertices, - const Point &candidate) const; + 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 given point + * 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. + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. */ virtual - Point - normal_vector (const std::vector > & vertices, - const Point &point) const; + void + get_normals_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) const; protected: /** @@ -226,6 +336,24 @@ 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; }; @@ -239,12 +367,16 @@ protected: * 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, Luca Heltai, 2013 + * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001 */ template -class HyperBallBoundary : public FlatManifold +class HyperBallBoundary : public StraightBoundary { public: /** @@ -252,28 +384,81 @@ 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(const std::vector > &surrounding_points, - const std::vector &weights) const; + 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_on_quad (const typename Triangulation::quad_iterator &quad) const; /** - * 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. + * Calls + * @p get_intermediate_points_between_points. */ virtual - Point - normal_vector (const std::vector > & vertices, - const Point &p) const; + 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. + * + * 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; + + /** + * 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. @@ -309,28 +494,40 @@ 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_new_point for special cases where surrounding - * points are two. Returns the point which is at xi between - * p0 and p1, where xi is between 0 and 1. + /** + * 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. */ - Point get_intermediate_point (const Point &p0, - const Point &p1, - const double xi) const; - + void get_intermediate_points_between_points (const Point &p0, const Point &p1, + std::vector > &points) const; }; @@ -357,27 +554,61 @@ public: HalfHyperBallBoundary (const Point p = Point(), const double radius = 1.0); - /** - * Refer to the general documentation of this class and the - * documentation of the base class. + * 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(const std::vector > &surrounding_points, - const std::vector &weights) const; + virtual Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; - /** - * 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. + * Calls + * @p get_intermediate_points_between_points. */ - virtual - Point - normal_vector (const std::vector > & vertices, - const Point &p) const; + 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; }; @@ -393,21 +624,21 @@ public: * @author Wolfgang Bangerth, 1999 */ template -class HyperShellBoundary : public PolarManifold +class HyperShellBoundary : public HyperBallBoundary { public: /** - * Constructor. The center of the spheres defaults to the origin. + * Constructor. The center of the + * spheres defaults to the + * origin. + * + * Calls the constructor of its + * base @p HyperBallBoundary + * class with a dummy radius as + * argument. This radius will be + * ignored */ HyperShellBoundary (const Point ¢er = Point()); - - /** - * Project to manifold. The vertices are used to compute - * automatically the radius. - */ - virtual Point - project_to_manifold(const std::vector > &vertices, - const Point &p) const; }; @@ -444,14 +675,58 @@ 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; /** - * Compute the normal to the surface. + * Construct a new point on a quad. */ - virtual Point - normal_vector (const std::vector > & vertices, - const Point &point) const; - + 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. + */ + virtual void + get_normals_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) const; + private: /** * Inner and outer radii of the shell. @@ -470,7 +745,7 @@ private: * dim=2,spacedim=3, that is, just the surface. */ template -class TorusBoundary : public ManifoldChart +class TorusBoundary : public Boundary { public: /** @@ -479,43 +754,83 @@ public: */ TorusBoundary (const double R, const double r); +//Boundary Refinenment Functions /** - * 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. + * 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; + + /** + * Construct a new points on a line. */ - 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; + virtual void get_intermediate_points_on_line ( + const typename Triangulation< dim, spacedim >::line_iterator &line, + std::vector< Point< spacedim > > &points) const ; /** - * Get the normal from cartesian coordinates. This normal does not - * have unit length. + * Construct a new points on a quad. */ - Point normal_vector(const std::vector > & vertices, - const Point &p) const; + virtual void get_intermediate_points_on_quad ( + const typename Triangulation< dim, spacedim >::quad_iterator &quad, + std::vector< Point< spacedim > > &points ) const ; + + /** + * 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 ; 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. */ @@ -527,6 +842,59 @@ 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/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index abb0ee58f6..6b57aea315 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -4424,12 +4424,10 @@ namespace VectorTools // sign of the normal vector provided by the boundary // if they should point in different directions. this is the // case in tests/deal.II/no_flux_11. - std::vector > vertices(GeometryInfo::vertices_per_face); - for(unsigned int v=0; v::vertices_per_face; ++v) - vertices[v] = cell->face(face_no)->vertex(v); Point normal_vector = (cell->face(face_no)->get_boundary() - .normal_vector (vertices, fe_values.quadrature_point(i))); + .normal_vector (cell->face(face_no), + fe_values.quadrature_point(i))); if (normal_vector * fe_values.normal_vector(i) < 0) normal_vector *= -1; Assert (std::fabs(normal_vector.norm() - 1) < 1e-14, diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index c03cf1acc7..44cdca00d7 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -1744,16 +1744,14 @@ next_cell: fix_up_object (const Iterator &object, const bool respect_manifold) { - const Manifold + const Boundary *manifold = (respect_manifold ? &object->get_boundary() : 0); const unsigned int structdim = Iterator::AccessorType::structure_dimension; const unsigned int spacedim = Iterator::AccessorType::space_dimension; - std::vector > vertices(GeometryInfo::vertices_per_cell); - for(unsigned int i=0; ivertex(i); // right now we can only deal // with cells that have been @@ -1817,16 +1815,16 @@ next_cell: / eps); else - gradient[d] - = ((objective_function (object, - manifold->project_to_manifold(vertices, - object_mid_point + h)) - - - objective_function (object, - manifold->project_to_manifold(vertices, - object_mid_point - h))) - / - eps); + gradient[d] + = ((objective_function (object, + manifold->project_to_surface(object, + object_mid_point + h)) + - + objective_function (object, + manifold->project_to_surface(object, + object_mid_point - h))) + / + eps); } // sometimes, the @@ -1858,8 +1856,8 @@ next_cell: gradient; if (respect_manifold == true) - object_mid_point = manifold->project_to_manifold(vertices, - object_mid_point); + object_mid_point = manifold->project_to_surface(object, + object_mid_point); // compute current value of the // objective function diff --git a/deal.II/source/grid/tria.cc b/deal.II/source/grid/tria.cc index be24961e2f..a3701f8c80 100644 --- a/deal.II/source/grid/tria.cc +++ b/deal.II/source/grid/tria.cc @@ -9585,23 +9585,13 @@ Triangulation::set_mesh_smoothing(const MeshSmoothing mesh_smooth template void -Triangulation::set_boundary (const types::boundary_id number, +Triangulation::set_boundary (const types::boundary_id m_number, const Boundary &boundary_object) -{ - set_manifold(m_number, boundary_object); -} - - - -template -void -Triangulation::set_manifold (const types::manifold_id m_number, - const Boundary &manifold_object) { Assert(m_number < numbers::invalid_manifold_id, ExcIndexRange(m_number,0,numbers::invalid_manifold_id)); - manifold[m_number] = &manifold_object; + manifold[m_number] = &boundary_object; } diff --git a/deal.II/source/grid/tria_accessor.cc b/deal.II/source/grid/tria_accessor.cc index d092696095..c0df81098c 100644 --- a/deal.II/source/grid/tria_accessor.cc +++ b/deal.II/source/grid/tria_accessor.cc @@ -22,7 +22,6 @@ #include #include #include -#include #include #include -- 2.39.5