From c6d7084a72a8546b25665d8783b8c3dc9cfbd50f Mon Sep 17 00:00:00 2001 From: heltai Date: Mon, 13 Jan 2014 18:11:15 +0000 Subject: [PATCH] Removed all deprecated methods from Boundary and StraightBoundary. git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32202 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/grid/tria_boundary.h | 489 +------- deal.II/source/fe/mapping_c1.cc | 10 +- deal.II/source/grid/tria_boundary.cc | 971 --------------- deal.II/source/grid/tria_boundary_lib.cc | 1151 ------------------ 4 files changed, 54 insertions(+), 2567 deletions(-) diff --git a/deal.II/include/deal.II/grid/tria_boundary.h b/deal.II/include/deal.II/grid/tria_boundary.h index 85abb09c76..b80a2bcd2f 100644 --- a/deal.II/include/deal.II/grid/tria_boundary.h +++ b/deal.II/include/deal.II/grid/tria_boundary.h @@ -34,290 +34,93 @@ template class Triangulation; /** - * This class is used to represent a boundary to a triangulation. + * This class is used to represent a boundary to a triangulation, + * and as such, it is a specialization of the Manifold + * class, where the topological dimension of the Manifold is fixed + * to be (dim-1), i.e., the boundary of a triangulation with + * topological dimension dim. + * * When a triangulation creates a new vertex on the boundary of the * domain, it determines the new vertex' coordinates through the * following code (here in two dimensions): * @code * ... - * Point<2> new_vertex = boundary.get_new_point_on_line (line); + * Point<2> new_vertex = boundary.get_new_point (points, weights); * ... * @endcode - * @p line denotes the line at the boundary that shall be refined - * and for which we seek the common point of the two child lines. + * @p points is a vector containing the points which define the object + * that shall be refined, while @p weights are used to decide at what + * location the child point should be created with respect to its + * parents. * * In 3D, a new vertex may be placed on the middle of a line or on * the middle of a side. Respectively, the library calls * @code * ... * Point<3> new_line_vertices[4] - * = { boundary.get_new_point_on_line (face->line(0)), - * boundary.get_new_point_on_line (face->line(1)), - * boundary.get_new_point_on_line (face->line(2)), - * boundary.get_new_point_on_line (face->line(3)) }; + * = { boundary.get_new_point (vertices_of_line_0, w_line), + * boundary.get_new_point (vertices_of_line_1, w_line), + * boundary.get_new_point (vertices_of_line_2, w_line), + * boundary.get_new_point (vertices_of_line_3, w_line) }; * ... * @endcode - * to get the four midpoints of the lines bounding the quad at the + * where @p w_line is a vector containing the values (1/2, 1/2). + * This return the four midpoints of the lines bounding the quad at the * boundary, and after that * @code * ... - * Point<3> new_quad_vertex = boundary.get_new_point_on_quad (face); + * Point<3> new_quad_vertex = boundary.get_new_point (vertices_of_face, w_quad); * ... * @endcode * to get the midpoint of the face. It is guaranteed that this order * (first lines, then faces) holds, so you can use information from - * the children of the four lines of a face, since these already exist - * at the time the midpoint of the face is to be computed. + * the children of the four lines of a face, since these already + * exist at the time the midpoint of the face is to be + * computed. This in fact is exploited in the library by passing a + * list of 8 points (the actual vertices, and the four points + * computed above) together with 8 weights which minimize cell + * distortion. * - * Since iterators are passed to the functions, you may use information - * about boundary indicators and the like, as well as all other information - * provided by these objects. + * This class is derived from FlatManifold. A minimal + * implementation for derived classes is given by overloading the + * method Manifold::project_to_manifold, which is called + * internally by the FlatManifold class (from which + * Boundary is derived) with a guess computed by + * considering the Manifold "flat", i.e., by putting the new point + * in the weighted average of the surrounding points. * * There are specialisations, StraightBoundary, which places * the new point right into the middle of the given points, and * HyperBallBoundary creating a hyperball with given radius * around a given center point. * + * @deprecated A new Manifold class was introduced which + * generalises the functionality of this + * class. Boundary is no longer necessary, and should + * be replaced by any equivalent Manifold description. + * * @ingroup boundary - * @author Wolfgang Bangerth, 1999, 2001, 2009, Ralf Hartmann, 2001, 2008 + * @author Wolfgang Bangerth, 1999, 2001, 2009, Ralf Hartmann, 2001, 2008, + * Luca Heltai 2013, 2014 */ template class Boundary : public FlatManifold { public: - /** - * Type keeping information about the normals at the vertices of a face of a - * cell. Thus, there are GeometryInfo::vertices_per_face - * normal vectors, that define the tangent spaces of the boundary at the - * vertices. Note that the vectors stored in this object are not required to - * be normalized, nor to actually point outward, as one often will only want - * to check for orthogonality to define the tangent plane; if a function - * requires the normals to be normalized, then it must do so itself. - * - * For obvious reasons, this type is not useful in 1d. - */ - typedef Tensor<1,spacedim> FaceVertexNormals[GeometryInfo::vertices_per_face]; - /** * Destructor. Does nothing here, but needs to be declared to make it * virtual. */ virtual ~Boundary (); - - /** - * Return the point which shall become the new middle vertex of the two - * children of a regular line. In 2D, this line is a line at the boundary, - * while in 3d, it is bounding a face at the boundary (the lines therefore - * is also on the boundary). - */ - virtual - Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const = 0; - - /** - * Return the point which shall become the common point of the four children - * of a quad at the boundary in three or more spatial dimensions. This - * function therefore is only useful in at least three dimensions and should - * not be called for lower dimensions. - * - * This function is called after the four lines bounding the given @p quad - * are refined, so you may want to use the information provided by - * quad->line(i)->child(j), i=0...3, j=0,1. - * - * Because in 2D, this function is not needed, it is not made pure virtual, - * to avoid the need to overload it. The default implementation throws an - * error in any case, however. - */ - virtual - Point - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; - - /** - * Depending on dim=2 or dim=3 this function calls the - * get_new_point_on_line or the get_new_point_on_quad function. It throws an - * exception for dim=1. This wrapper allows dimension independent - * programming. - */ - Point - get_new_point_on_face (const typename Triangulation::face_iterator &face) const; - - /** - * Return intermediate points on a line spaced according to the interior - * support points of the 1D Gauss-Lobatto quadrature formula. - * - * The number of points requested is given by the size of the vector @p - * points. It is the task of the derived classes to arrange the points in - * approximately equal distances. - * - * This function is called by the @p MappingQ class. This happens on each - * face line of a cells that has got at least one boundary line. - * - * As this function is not needed for @p MappingQ1, it is not made pure - * virtual, to avoid the need to overload it. The default implementation - * throws an error in any case, however. - */ - virtual - void - get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, - std::vector > &points) const; - - /** - * Return intermediate points on a line spaced according to the tensor - * product of the interior support points of the 1D Gauss-Lobatto quadrature - * formula. - * - * The number of points requested is given by the size of the vector @p - * points. It is required that this number is a square of another integer, - * i.e. n=points.size()=m*m. It is the task of the derived classes - * to arrange the points such they split the quad into (m+1)(m+1) - * approximately equal-sized subquads. - * - * This function is called by the MappingQ<3> class. This happens - * each face quad of cells in 3d that has got at least one boundary face - * quad. - * - * As this function is not needed for @p MappingQ1, it is not made pure - * virtual, to avoid the need to overload it. The default implementation - * throws an error in any case, however. - */ - virtual - void - get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, - std::vector > &points) const; - - /** - * Depending on dim=2 or dim=3 this function calls the - * get_intermediate_points_on_line or the get_intermediate_points_on_quad - * function. It throws an exception for dim=1. This wrapper allows - * dimension independent programming. - */ - void - get_intermediate_points_on_face (const typename Triangulation::face_iterator &face, - std::vector > &points) const; - - /** - * Return the normal vector to the surface at the point p. If p is not in - * fact on the surface, but only close-by, try to return something - * reasonable, for example the normal vector at the surface point closest to - * p. (The point p will in fact not normally lie on the actual surface, but - * rather be a quadrature point mapped by some polynomial mapping; the - * mapped surface, however, will not usually coincide with the actual - * surface.) - * - * The face iterator gives an indication which face this function is - * supposed to compute the normal vector for. This is useful if the - * boundary of the domain is composed of different nondifferential pieces - * (for example when using the StraightBoundary class to approximate a - * geometry that is completely described by the coarse mesh, with piecewise - * (bi-)linear components between the vertices, but where the boundary may - * have a kink at the vertices itself). - * - * @note Implementations of this function should be able to assume that the - * point p lies within or close to the face described by the first - * argument. In turn, callers of this function should ensure that this is in - * fact the case. - */ - virtual - Tensor<1,spacedim> - normal_vector (const typename Triangulation::face_iterator &face, - const Point &p) const; - - /** - * Compute the normal vectors to the boundary at each vertex of the given - * face. It is not required that the normal vectors be normed - * somehow. Neither is it required that the normals actually point outward. - * - * This function is needed to compute data for C1 mappings. The default - * implementation is to throw an error, so you need not overload this - * function in case you do not intend to use C1 mappings. - * - * Note that when computing normal vectors at a vertex where the boundary is - * not differentiable, you have to make sure that you compute the one-sided - * limits, i.e. limit with respect to points inside the given face. - */ - virtual - void - get_normals_at_vertices (const typename Triangulation::face_iterator &face, - FaceVertexNormals &face_vertex_normals) const; - - /** - * Given a candidate point and a line segment characterized by the iterator, - * return a point that lies on the surface described by this object. This - * function is used in some mesh smoothing algorithms that try to move - * around points in order to improve the mesh quality but need to ensure - * that points that were on the boundary remain on the boundary. - * - * If spacedim==1, then the line represented by the line iterator is the - * entire space (i.e. it is a cell, not a part of the boundary), and the - * returned point equals the given input point. - * - * Derived classes do not need to implement this function unless mesh - * smoothing algorithms are used with a particular boundary object. The - * default implementation of this function throws an exception of type - * ExcPureFunctionCalled. - */ - virtual - Point - project_to_surface (const typename Triangulation::line_iterator &line, - const Point &candidate) const; - - /** - * Same function as above but for a point that is to be projected onto the - * area characterized by the given quad. - * - * If spacedim<=2, then the surface represented by the quad iterator is the - * entire space (i.e. it is a cell, not a part of the boundary), and the - * returned point equals the given input point. - */ - virtual - Point - project_to_surface (const typename Triangulation::quad_iterator &quad, - const Point &candidate) const; - - /** - * Same function as above but for a point that is to be projected onto the - * area characterized by the given quad. - * - * If spacedim<=3, then the manifold represented by the hex iterator is the - * entire space (i.e. it is a cell, not a part of the boundary), and the - * returned point equals the given input point. - */ - virtual - Point - project_to_surface (const typename Triangulation::hex_iterator &hex, - const Point &candidate) const; - -protected: - /** - * Returns the support points of the Gauss-Lobatto quadrature formula used - * for intermediate points. - * - * @note Since the boundary description is closely tied to the unit cell - * support points of MappingQ, new boundary descriptions need to explicitly - * use these Gauss-Lobatto points and not equidistant points. - */ - const std::vector > & - get_line_support_points (const unsigned int n_intermediate_points) const; - -private: - /** - * Point generator for the intermediate points on a boundary. - */ - mutable std::vector > > points; - - /** - * Mutex for protecting the points array. - */ - mutable Threads::Mutex mutex; + } DEAL_II_DEPRECATED; - /** * Specialization of Boundary, which places the new point * right into the middle of the given points. The middle is defined - * as the arithmetic mean of the points. + * as the weighted average of the surrounding points. * * This class does not really describe a boundary in the usual * sense. By placing new points in the middle of old ones, it rather @@ -325,9 +128,16 @@ private: * polygon/polyhedron defined by the boundary of the initial coarse * triangulation. * + * @deprecated The functionality of this class is equivalent, but + * less general, to that of FlatManifold. Please use + * FlatManifold instead of this + * class. StraightBoundary will be removed in future + * releases. + * * @ingroup boundary * - * @author Wolfgang Bangerth, 1998, 2001, Ralf Hartmann, 2001 + * @author Wolfgang Bangerth, 1998, 2001, Ralf Hartmann, 2001, Luca + * Heltai 2013, 2014 */ template class StraightBoundary : public Boundary @@ -337,208 +147,7 @@ public: * Default constructor. Some compilers require this for some reasons. */ StraightBoundary (); - - /** - * Let the new point be the arithmetic mean of the two vertices of the line. - * - * Refer to the general documentation of this class and the documentation of - * the base class for more information. - */ - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - - /** - * Let the new point be the arithmetic mean of the four vertices of this - * quad and the four midpoints of the lines, which are already created at - * the time of calling this function. - * - * Refer to the general documentation of this class and the documentation of - * the base class for more information. - */ - virtual - Point - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; - - /** - * Gives n=points.size() points that splits the StraightBoundary - * line into $n+1$ partitions of equal lengths. - * - * Refer to the general documentation of this class and the documentation of - * the base class. - */ - virtual - void - get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, - std::vector > &points) const; - - /** - * Gives n=points.size()=m*m points that splits the - * StraightBoundary quad into $(m+1)(m+1)$ subquads of equal size. - * - * Refer to the general documentation of this class and the documentation of - * the base class. - */ - 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; - - /** - * Given a candidate point and a line segment characterized by the iterator, - * return a point that lies on the surface described by this object. This - * function is used in some mesh smoothing algorithms that try to move - * around points in order to improve the mesh quality but need to ensure - * that points that were on the boundary remain on the boundary. - * - * The point returned is the projection of the candidate point onto the line - * through the two vertices of the given line iterator. - * - * If spacedim==1, then the line represented by the line iterator is the - * entire space (i.e. it is a cell, not a part of the boundary), and the - * returned point equals the given input point. - */ - virtual - Point - project_to_surface (const typename Triangulation::line_iterator &line, - const Point &candidate) const; - - /** - * Same function as above but for a point that is to be projected onto the - * area characterized by the given quad. - * - * The point returned is the projection of the candidate point onto the - * bilinear surface spanned by the four vertices of the given quad iterator. - * - * If spacedim<=2, then the surface represented by the quad iterator is the - * entire space (i.e. it is a cell, not a part of the boundary), and the - * returned point equals the given input point. - */ - virtual - Point - project_to_surface (const typename Triangulation::quad_iterator &quad, - const Point &candidate) const; - - /** - * Same function as above but for a point that is to be projected onto the - * area characterized by the given quad. - * - * The point returned is the projection of the candidate point onto the - * trilinear manifold spanned by the eight vertices of the given hex - * iterator. - * - * If spacedim<=3, then the manifold represented by the hex iterator is the - * entire space (i.e. it is a cell, not a part of the boundary), and the - * returned point equals the given input point. - */ - virtual - Point - project_to_surface (const typename Triangulation::hex_iterator &hex, - const Point &candidate) const; -}; - - - -/* -------------- declaration of explicit specializations ------------- */ - -#ifndef DOXYGEN - -template <> -Point<1> -Boundary<1,1>:: -get_new_point_on_face (const Triangulation<1,1>::face_iterator &) const; - -template <> -void -Boundary<1,1>:: -get_intermediate_points_on_face (const Triangulation<1,1>::face_iterator &, - std::vector > &) const; - -template <> -Point<2> -Boundary<1,2>:: -get_new_point_on_face (const Triangulation<1,2>::face_iterator &) const; - -template <> -void -Boundary<1,2>:: -get_intermediate_points_on_face (const Triangulation<1,2>::face_iterator &, - std::vector > &) const; - - - -template <> -Point<3> -Boundary<1,3>:: -get_new_point_on_face (const Triangulation<1,3>::face_iterator &) const; - -template <> -void -Boundary<1,3>:: -get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &, - std::vector > &) const; -template <> -void -StraightBoundary<1,1>:: -get_normals_at_vertices (const Triangulation<1,1>::face_iterator &, - Boundary<1,1>::FaceVertexNormals &) const; -template <> -void -StraightBoundary<2,2>:: -get_normals_at_vertices (const Triangulation<2,2>::face_iterator &face, - Boundary<2,2>::FaceVertexNormals &face_vertex_normals) const; -template <> -void -StraightBoundary<3,3>:: -get_normals_at_vertices (const Triangulation<3,3>::face_iterator &face, - Boundary<3,3>::FaceVertexNormals &face_vertex_normals) const; - -template <> -Point<3> -StraightBoundary<3,3>:: -get_new_point_on_quad (const Triangulation<3,3>::quad_iterator &quad) const; - -template <> -void -StraightBoundary<1,1>:: -get_intermediate_points_on_line (const Triangulation<1,1>::line_iterator &, - std::vector > &) const; - -template <> -void -StraightBoundary<3,3>:: -get_intermediate_points_on_quad (const Triangulation<3,3>::quad_iterator &quad, - std::vector > &points) const; - -template <> -Point<3> -StraightBoundary<1,3>:: -project_to_surface (const Triangulation<1, 3>::quad_iterator &quad, - const Point<3> &y) const; - - -#endif // DOXYGEN +} DEAL_II_DEPRECATED; DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/mapping_c1.cc b/deal.II/source/fe/mapping_c1.cc index caa97ba586..3196ececa5 100644 --- a/deal.II/source/fe/mapping_c1.cc +++ b/deal.II/source/fe/mapping_c1.cc @@ -59,6 +59,9 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) { const Triangulation::line_iterator line = cell->line(line_no); + std::vector > points(2); + points[0] = line->vertex(0); + points[1] = line->vertex(1); if (line->at_boundary()) { @@ -66,11 +69,8 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce // from the boundary description const Manifold &boundary = line->get_boundary(); - - std::vector > points(2); + std::vector > face_vertex_normals(2); - points[0] = line->vertex(0); - points[1] = line->vertex(1); boundary.get_normals_at_points (face_vertex_normals, points); // then transform them into interpolation points for a cubic @@ -140,7 +140,7 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce // not at boundary { static const StraightBoundary straight_boundary; - straight_boundary.get_intermediate_points_on_line (line, line_points); + straight_boundary.get_intermediate_points (line_points, points); a.insert (a.end(), line_points.begin(), line_points.end()); }; }; diff --git a/deal.II/source/grid/tria_boundary.cc b/deal.II/source/grid/tria_boundary.cc index d28ab418b9..3ed27b97b2 100644 --- a/deal.II/source/grid/tria_boundary.cc +++ b/deal.II/source/grid/tria_boundary.cc @@ -14,18 +14,9 @@ // // --------------------------------------------------------------------- -#include #include -#include -#include -#include -#include -#include - DEAL_II_NAMESPACE_OPEN - - /* -------------------------- Boundary --------------------- */ @@ -33,245 +24,6 @@ template Boundary::~Boundary () {} - - -template -Point -Boundary:: -get_new_point_on_quad (const typename Triangulation::quad_iterator &) const -{ - Assert (false, ExcPureFunctionCalled()); - return Point(); -} - -template -void -Boundary:: -get_intermediate_points_on_line (const typename Triangulation::line_iterator &, - std::vector > &) const -{ - Assert (false, ExcPureFunctionCalled()); -} - - - -template -void -Boundary:: -get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &, - std::vector > &) const -{ - Assert (false, ExcPureFunctionCalled()); -} - - - -template -Point -Boundary:: -get_new_point_on_face (const typename Triangulation::face_iterator &face) const -{ - Assert (dim>1, ExcImpossibleInDim(dim)); - - switch (dim) - { - case 2: - return get_new_point_on_line (face); - case 3: - return get_new_point_on_quad (face); - } - - return Point(); -} - - -template -void -Boundary:: -get_intermediate_points_on_face (const typename Triangulation::face_iterator &face, - std::vector > &points) const -{ - Assert (dim>1, ExcImpossibleInDim(dim)); - - switch (dim) - { - case 2: - get_intermediate_points_on_line (face, points); - break; - case 3: - get_intermediate_points_on_quad (face, points); - break; - default: - Assert (false, ExcNotImplemented()); - } -} - - - -template <> -Point<1> -Boundary<1,1>:: -get_new_point_on_face (const Triangulation<1,1>::face_iterator &) const -{ - Assert (false, ExcImpossibleInDim(1)); - return Point<1>(); -} - - -template <> -void -Boundary<1,1>:: -get_intermediate_points_on_face (const Triangulation<1,1>::face_iterator &, - std::vector > &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} - - - -template <> -Point<2> -Boundary<1,2>:: -get_new_point_on_face (const Triangulation<1,2>::face_iterator &) const -{ - Assert (false, ExcImpossibleInDim(1)); - return Point<2>(); -} - - -template <> -void -Boundary<1,2>:: -get_intermediate_points_on_face (const Triangulation<1,2>::face_iterator &, - std::vector > &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} - - - -template <> -Point<3> -Boundary<1,3>:: -get_new_point_on_face (const Triangulation<1,3>::face_iterator &) const -{ - Assert (false, ExcImpossibleInDim(1)); - return Point<3>(); -} - - -template <> -void -Boundary<1,3>:: -get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &, - std::vector > &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} - - - - -template -Tensor<1,spacedim> -Boundary:: -normal_vector (const typename Triangulation::face_iterator &, - const Point &) const -{ - Assert (false, ExcPureFunctionCalled()); - return Tensor<1,spacedim>(); -} - - - -template -void -Boundary:: -get_normals_at_vertices (const typename Triangulation::face_iterator &, - FaceVertexNormals &) const -{ - Assert (false, ExcPureFunctionCalled()); -} - - - -template -Point -Boundary:: -project_to_surface (const typename Triangulation::line_iterator &, - const Point &trial_point) const -{ - if (spacedim <= 1) - return trial_point; - else - { - Assert (false, ExcPureFunctionCalled()); - return Point(); - } -} - - - -template -Point -Boundary:: -project_to_surface (const typename Triangulation::quad_iterator &, - const Point &trial_point) const -{ - if (spacedim <= 2) - return trial_point; - else - { - Assert (false, ExcPureFunctionCalled()); - return Point(); - } -} - - - -template -Point -Boundary:: -project_to_surface (const typename Triangulation::hex_iterator &, - const Point &trial_point) const -{ - if (spacedim <= 3) - return trial_point; - else - { - Assert (false, ExcPureFunctionCalled()); - return Point(); - } -} - - - -template -const std::vector > & -Boundary:: -get_line_support_points (const unsigned int n_intermediate_points) const -{ - if (points.size() <= n_intermediate_points || - points[n_intermediate_points].get() == 0) - { - Threads::Mutex::ScopedLock lock(mutex); - if (points.size() <= n_intermediate_points) - points.resize(n_intermediate_points+1); - - // another thread might have created points in the meantime - if (points[n_intermediate_points].get() == 0) - { - std_cxx1x::shared_ptr > - quadrature (new QGaussLobatto<1>(n_intermediate_points+2)); - points[n_intermediate_points] = quadrature; - } - } - return points[n_intermediate_points]->get_points(); -} - - - - /* -------------------------- StraightBoundary --------------------- */ @@ -279,729 +31,6 @@ template StraightBoundary::StraightBoundary () {} - -template -Point -StraightBoundary:: -get_new_point_on_line (const typename Triangulation::line_iterator &line) const -{ - return (line->vertex(0) + line->vertex(1)) / 2; -} - - -namespace -{ - // compute the new midpoint of a quad -- - // either of a 2d cell on a manifold in 3d - // or of a face of a 3d triangulation in 3d - template - Point<3> - compute_new_point_on_quad (const typename Triangulation::quad_iterator &quad) - { - // generate a new point in the middle of - // the face based on the points on the - // edges and the vertices. - // - // there is a pathological situation when - // this face is on a straight boundary, but - // one of its edges and the face behind it - // are not; if that face is refined first, - // the new point in the middle of that edge - // may not be at the same position as - // quad->line(.)->center() would have been, - // but would have been moved to the - // non-straight boundary. We cater to that - // situation by using existing edge - // midpoints if available, or center() if - // not - // - // note that this situation can not happen - // during mesh refinement, as there the - // edges are refined first and only then - // the face. thus, the check whether a line - // has children does not lead to the - // situation where the new face midpoints - // have different positions depending on - // which of the two cells is refined first. - // - // the situation where the edges aren't - // refined happens when a higher order - // MappingQ requests the midpoint of a - // face, though, and it is for these cases - // that we need to have the check available - // - // note that the factor of 1/8 for each - // of the 8 surrounding points isn't - // chosen arbitrarily. rather, we may ask - // where the harmonic map would place the - // point (0,0) if we map the square - // [-1,1]^2 onto the domain that is - // described using the 4 vertices and 4 - // edge point points of this quad. we can - // then discretize the harmonic map using - // four cells and Q1 elements on each of - // the quadrants of the square [-1,1]^2 - // and see where the midpoint would land - // (this is the procedure we choose, for - // example, in - // GridGenerator::laplace_solve) and it - // turns out that it will land at the - // mean of the 8 surrounding - // points. whether a discretization of - // the harmonic map with only 4 cells is - // adequate is a different question - // altogether, of course. - return (quad->vertex(0) + quad->vertex(1) + - quad->vertex(2) + quad->vertex(3) + - (quad->line(0)->has_children() ? - quad->line(0)->child(0)->vertex(1) : - quad->line(0)->center()) + - (quad->line(1)->has_children() ? - quad->line(1)->child(0)->vertex(1) : - quad->line(1)->center()) + - (quad->line(2)->has_children() ? - quad->line(2)->child(0)->vertex(1) : - quad->line(2)->center()) + - (quad->line(3)->has_children() ? - quad->line(3)->child(0)->vertex(1) : - quad->line(3)->center()) ) / 8; - } -} - - - -template -Point -StraightBoundary:: -get_new_point_on_quad (const typename Triangulation::quad_iterator &) const -{ - Assert (false, ExcImpossibleInDim(dim)); - return Point(); -} - - -template <> -Point<3> -StraightBoundary<2,3>:: -get_new_point_on_quad (const Triangulation<2,3>::quad_iterator &quad) const -{ - return compute_new_point_on_quad<2> (quad); -} - - - -template <> -Point<3> -StraightBoundary<3>:: -get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const -{ - return compute_new_point_on_quad<3> (quad); -} - - - -template <> -void -StraightBoundary<1>:: -get_intermediate_points_on_line (const Triangulation<1>::line_iterator &, - std::vector > &) const -{ - Assert(false, ExcImpossibleInDim(1)); -} - -template <> -void -StraightBoundary<1, 2>:: -get_intermediate_points_on_line (const Triangulation<1, 2>::line_iterator &line, - std::vector > &points) const -{ - const unsigned int spacedim = 2; - const unsigned int n=points.size(); - Assert(n>0, ExcInternalError()); - - // Use interior points of QGaussLobatto quadrature formula support points - // for consistency with MappingQ - const std::vector > &line_points = this->get_line_support_points(n); - const Point vertices[2] = { line->vertex(0), - line->vertex(1) - }; - - for (unsigned int i=0; i -void -StraightBoundary:: -get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, - std::vector > &points) const -{ - const unsigned int n=points.size(); - Assert(n>0, ExcInternalError()); - - // Use interior points of QGaussLobatto quadrature formula support points - // for consistency with MappingQ - const std::vector > &line_points = this->get_line_support_points(n); - - const Point vertices[2] = { line->vertex(0), - line->vertex(1) - }; - - for (unsigned int i=0; i -void -StraightBoundary:: -get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &, - std::vector > &) const -{ - Assert(false, ExcImpossibleInDim(dim)); -} - - - -template <> -void -StraightBoundary<3>:: -get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, - std::vector > &points) const -{ - const unsigned int spacedim = 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()); - - const std::vector > &line_points = this->get_line_support_points(m); - - const Point vertices[4] = { quad->vertex(0), - quad->vertex(1), - quad->vertex(2), - quad->vertex(3) - }; - - for (unsigned int i=0; i -void -StraightBoundary<2,3>:: -get_intermediate_points_on_quad (const Triangulation<2,3>::quad_iterator &quad, - std::vector > &points) const -{ - const unsigned int spacedim = 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()); - - const std::vector > &line_points = this->get_line_support_points(m); - - const Point vertices[4] = { quad->vertex(0), - quad->vertex(1), - quad->vertex(2), - quad->vertex(3) - }; - - for (unsigned int i=0; i -Tensor<1,1> -StraightBoundary<1,1>:: -normal_vector (const Triangulation<1,1>::face_iterator &, - const Point<1> &) const -{ - Assert (false, ExcNotImplemented()); - return Tensor<1,1>(); -} - - -template <> -Tensor<1,2> -StraightBoundary<1,2>:: -normal_vector (const Triangulation<1,2>::face_iterator &, - const Point<2> &) const -{ - Assert (false, ExcNotImplemented()); - return Tensor<1,2>(); -} - - -template <> -Tensor<1,3> -StraightBoundary<1,3>:: -normal_vector (const Triangulation<1,3>::face_iterator &, - const Point<3> &) const -{ - Assert (false, ExcNotImplemented()); - return Tensor<1,3>(); -} - - -namespace internal -{ - namespace - { - /** - * Compute the normalized cross product of a set of dim-1 basis - * vectors. - */ - Tensor<1,2> - normalized_alternating_product (const Tensor<1,2> (&basis_vectors)[1]) - { - Tensor<1,2> tmp; - cross_product (tmp, basis_vectors[0]); - return tmp/tmp.norm(); - } - - - - Tensor<1,3> - normalized_alternating_product (const Tensor<1,3> ( &)[1]) - { - // we get here from StraightBoundary<2,3>::normal_vector, but - // the implementation below is bogus for this case anyway - // (see the assert at the beginning of that function). - Assert (false, ExcNotImplemented()); - return Tensor<1,3>(); - } - - - - Tensor<1,3> - normalized_alternating_product (const Tensor<1,3> (&basis_vectors)[2]) - { - Tensor<1,3> tmp; - cross_product (tmp, basis_vectors[0], basis_vectors[1]); - return tmp/tmp.norm(); - } - - } -} - - -template -Tensor<1,spacedim> -StraightBoundary:: -normal_vector (const typename Triangulation::face_iterator &face, - const Point &p) const -{ - // I don't think the implementation below will work when dim!=spacedim; - // in fact, I believe that we don't even have enough information here, - // because we would need to know not only about the tangent vectors - // of the face, but also of the cell, to compute the normal vector. - // Someone will have to think about this some more. - Assert (dim == spacedim, ExcNotImplemented()); - - // in order to find out what the normal vector is, we first need to - // find the reference coordinates of the point p on the given face, - // or at least the reference coordinates of the closest point on the - // face - // - // in other words, we need to find a point xi so that f(xi)=||F(xi)-p||^2->min - // where F(xi) is the mapping. this algorithm is implemented in - // MappingQ1::transform_real_to_unit_cell but only for cells, - // while we need it for faces here. it's also implemented in somewhat - // more generality there using the machinery of the MappingQ1 class - // while we really only need it for a specific case here - // - // in any case, the iteration we use here is a Gauss-Newton's iteration with - // xi^{n+1} = xi^n - H(xi^n)^{-1} J(xi^n) - // where - // J(xi) = (grad F(xi))^T (F(xi)-p) - // and - // H(xi) = [grad F(xi)]^T [grad F(xi)] - // In all this, - // F(xi) = sum_v vertex[v] phi_v(xi) - // We get the shape functions phi_v from an object of type FE_Q(1) - - // we start with the point xi=1/2, xi=(1/2,1/2), ... - const unsigned int facedim = dim-1; - - Point xi; - for (unsigned int i=0; i linear_fe(1); - - const double eps = 1e-12; - Tensor<1,spacedim> grad_F[facedim]; - while (true) - { - Point F; - for (unsigned int v=0; v::vertices_per_cell; ++v) - F += face->vertex(v) * linear_fe.shape_value(v, xi); - - for (unsigned int i=0; i::vertices_per_cell; ++v) - grad_F[i] += face->vertex(v) * linear_fe.shape_grad(v, xi)[i]; - } - - Tensor<1,facedim> J; - for (unsigned int i=0; i H; - for (unsigned int i=0; i delta_xi = -invert(H) * J; - xi += delta_xi; - - if (delta_xi.norm() < eps) - break; - } - - // so now we have the reference coordinates xi of the point p. - // we then have to compute the normal vector, which we can do - // by taking the (normalize) alternating product of all the tangent - // vectors given by grad_F - return internal::normalized_alternating_product(grad_F); -} - - - -template <> -void -StraightBoundary<1>:: -get_normals_at_vertices (const Triangulation<1>::face_iterator &, - Boundary<1,1>::FaceVertexNormals &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} - -template <> -void -StraightBoundary<1,2>:: -get_normals_at_vertices (const Triangulation<1,2>::face_iterator &, - Boundary<1,2>::FaceVertexNormals &) const -{ - Assert (false, ExcNotImplemented()); -} - - -template <> -void -StraightBoundary<1,3>:: -get_normals_at_vertices (const Triangulation<1,3>::face_iterator &, - Boundary<1,3>::FaceVertexNormals &) const -{ - Assert (false, ExcNotImplemented()); -} - - - -template <> -void -StraightBoundary<2>:: -get_normals_at_vertices (const Triangulation<2>::face_iterator &face, - Boundary<2,2>::FaceVertexNormals &face_vertex_normals) const -{ - const Tensor<1,2> tangent = face->vertex(1) - face->vertex(0); - for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) - // compute normals from tangent - face_vertex_normals[vertex] = Point<2>(tangent[1], - -tangent[0]); -} - -template <> -void -StraightBoundary<2,3>:: -get_normals_at_vertices (const Triangulation<2,3>::face_iterator &face, - Boundary<2,3>::FaceVertexNormals &face_vertex_normals) const -{ - const Tensor<1,3> tangent = face->vertex(1) - face->vertex(0); - for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) - // compute normals from tangent - face_vertex_normals[vertex] = Point<3>(tangent[1], - -tangent[0],0); - Assert(false, ExcNotImplemented()); -} - - - - -template <> -void -StraightBoundary<3>:: -get_normals_at_vertices (const Triangulation<3>::face_iterator &face, - Boundary<3,3>::FaceVertexNormals &face_vertex_normals) const -{ - const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face; - - static const unsigned int neighboring_vertices[4][2]= - { {1,2},{3,0},{0,3},{2,1}}; - for (unsigned int vertex=0; vertex tangents[2] - = { face->vertex(neighboring_vertices[vertex][0]) - - face->vertex(vertex), - face->vertex(neighboring_vertices[vertex][1]) - - face->vertex(vertex) - }; - - // then compute the normal by - // taking the cross - // product. since the normal is - // not required to be - // normalized, no problem here - cross_product (face_vertex_normals[vertex], - tangents[0], tangents[1]); - }; -} - - - -template -Point -StraightBoundary:: -project_to_surface (const typename Triangulation::line_iterator &line, - const Point &trial_point) const -{ - if (spacedim <= 1) - return trial_point; - else - { - // find the point that lies on - // the line p1--p2. the - // formulas pan out to - // something rather simple - // because the mapping to the - // line is linear - const Point p1 = line->vertex(0), - p2 = line->vertex(1); - const double s = (trial_point-p1)*(p2-p1) / ((p2-p1)*(p2-p1)); - return p1 + s*(p2-p1); - } -} - - - -namespace internal -{ - template - Point - compute_projection (const Iterator &object, - const Point &y, - internal::int2type) - { - // let's look at this for - // simplicity for a quad (dim==2) - // in a space with spacedim>2: - - // all points on the surface are given by - // x(\xi) = sum_i v_i phi_x(\xi) - // where v_i are the vertices of the quad, - // and \xi=(\xi_1,\xi_2) are the reference - // coordinates of the quad. so what we are - // trying to do is find a point x on - // the surface that is closest to the point - // y. there are different ways - // to solve this problem, but in the end - // it's a nonlinear problem and we have to - // find reference coordinates \xi so that - // J(\xi) = 1/2 || x(\xi)-y ||^2 - // is minimal. x(\xi) is a function that - // is dim-linear in \xi, so J(\xi) is - // a polynomial of degree 2*dim that - // we'd like to minimize. unless dim==1, - // we'll have to use a Newton - // method to find the - // answer. This leads to the - // following formulation of - // Newton steps: - // - // Given \xi_k, find \delta\xi_k so that - // H_k \delta\xi_k = - F_k - // where H_k is an approximation to the - // second derivatives of J at \xi_k, and - // F_k is the first derivative of J. - // We'll iterate this a number of times - // until the right hand side is small - // enough. As a stopping criterion, we - // terminate if ||\delta\xi|| xi; - for (unsigned int d=0; d x_k; - for (unsigned int i=0; i::vertices_per_cell; ++i) - x_k += object->vertex(i) * - GeometryInfo::d_linear_shape_function (xi, i); - - do - { - Tensor<1,dim> F_k; - for (unsigned int i=0; i::vertices_per_cell; ++i) - F_k += (x_k-y)*object->vertex(i) * - GeometryInfo::d_linear_shape_function_gradient (xi, i); - - Tensor<2,dim> H_k; - for (unsigned int i=0; i::vertices_per_cell; ++i) - for (unsigned int j=0; j::vertices_per_cell; ++j) - { - Tensor<2,dim> tmp; - outer_product (tmp, - GeometryInfo::d_linear_shape_function_gradient (xi, i), - GeometryInfo::d_linear_shape_function_gradient (xi, j)); - H_k += (object->vertex(i) * object->vertex(j)) * tmp; - } - - const Point delta_xi = - invert(H_k) * F_k; - xi += delta_xi; - - x_k = Point(); - for (unsigned int i=0; i::vertices_per_cell; ++i) - x_k += object->vertex(i) * - GeometryInfo::d_linear_shape_function (xi, i); - - if (delta_xi.norm() < 1e-5) - break; - } - while (true); - - return x_k; - } - - - // specialization for a quad in 1d - template - Point<1> - compute_projection (const Iterator &, - const Point<1> &y, - /* it's a quad: */internal::int2type<2>) - { - return y; - } - - // specialization for a quad in 2d - template - Point<2> - compute_projection (const Iterator &, - const Point<2> &y, - /* it's a quad: */internal::int2type<2>) - { - return y; - } -} - - - - - -template <> -Point<3> -StraightBoundary<1,3>:: -project_to_surface (const Triangulation<1, 3>::quad_iterator &, - const Point<3> &y) const -{ - return y; -} - -//TODO[SP]: This is just a horrible way out to make it compile in codim 2. -template -Point -StraightBoundary:: -project_to_surface (const typename Triangulation::quad_iterator &quad, - const Point &y) const -{ - if (spacedim <= 2) - return y; - else - return internal::compute_projection (quad, y, - /* it's a quad */internal::int2type<2>()); -} - - - -template -Point -StraightBoundary:: -project_to_surface (const typename Triangulation::hex_iterator &, - const Point &trial_point) const -{ - if (spacedim <= 3) - return trial_point; - else - { - // we can presumably call the - // same function as above (it's - // written in a generic way) - // but someone needs to check - // whether that actually yields - // the correct result - Assert (false, ExcNotImplemented()); - return Point(); - } -} - - - // explicit instantiations #include "tria_boundary.inst" diff --git a/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/source/grid/tria_boundary_lib.cc index 61bd284c53..025a58b879 100644 --- a/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/source/grid/tria_boundary_lib.cc @@ -58,35 +58,6 @@ 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:: @@ -111,160 +82,6 @@ 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 -// 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); -// } -// } - - - -// 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); - -// 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<1>:: -// get_normals_at_vertices (const Triangulation<1>::face_iterator &, -// Boundary<1,1>::FaceVertexNormals &) const -// { -// Assert (false, ExcImpossibleInDim(1)); -// } - - - - template Point CylinderBoundary:: @@ -314,54 +131,6 @@ 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 -// 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:: @@ -377,114 +146,6 @@ 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 -// ConeBoundary:: -// get_new_point_on_quad (const typename Triangulation::quad_iterator &) const -// { -// Assert (false, ExcImpossibleInDim (dim)); - -// 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<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 ()); - -// 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); - -// 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 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<1>:: -// get_normals_at_vertices (const Triangulation<1>::face_iterator &, -// Boundary<1,1>::FaceVertexNormals &) const -// { -// Assert (false, ExcImpossibleInDim (1)); -// } - - - template Point ConeBoundary:: @@ -516,47 +177,6 @@ HyperBallBoundary::HyperBallBoundary (const Point p, compute_radius_automatically(false) {} -// template -// Point -// HyperBallBoundary::get_new_point (const std::vector > &surrounding_points, -// const std::vector &weights) const -// { -// Assert(surrounding_points.size() == weights.size(), -// ExcDimensionMismatch(surrounding_points.size(), weights.size())); - -// double radius = 0; -// for(unsigned int i=0; i p = FlatManifold::get_new_point(surrounding_points, weights); -// p = p-center; -// p = p/p.norm()*radius+center; -// return 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 @@ -581,246 +201,6 @@ 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) 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<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 -// 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 -// 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,2>:: -// get_normals_at_vertices (const Triangulation<1,2>::face_iterator &, -// Boundary<1,2>::FaceVertexNormals &) const -// { -// Assert (false, ExcImpossibleInDim(1)); -// } - - - template Point HyperBallBoundary:: @@ -857,15 +237,6 @@ 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(const std::vector > &surrounding_points, @@ -887,109 +258,6 @@ HalfHyperBallBoundary::get_new_point(const std::vector > &surrou } - -// 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 -// 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<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:: -// 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 Point HalfHyperBallBoundary:: @@ -1051,239 +319,6 @@ 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_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); -// } - - - -// 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()); - -// // 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_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 Point HalfHyperShellBoundary:: @@ -1390,67 +425,6 @@ TorusBoundary<2,3>::pull_back(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)); - -// 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; - -// middle+= p0 + p1; -// middle*=0.5; - -// 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]; - -// for (unsigned int i=0; i<4; i++) -// p[i]=get_surf_coord(quad->vertex(i)); - -// 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; -// } -// } - -// for (unsigned int i=0; i<4; i++) -// middle+=p[i]; - -// middle*= 0.25; - -// return get_real_coord(middle); -// } - - - //Normal field without unit length template <> Point<3> @@ -1483,131 +457,6 @@ TorusBoundary<2,3>::normal_vector(const Point<3> p) const } - - -// 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_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)); -// } - - - // explicit instantiations #include "tria_boundary_lib.inst" -- 2.39.5