From: David Wells Date: Sat, 12 May 2018 14:44:10 +0000 (-0400) Subject: Remove all remaining usages of the Boundary class. X-Git-Tag: v9.1.0-rc1~1139^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1621ff316ead9849f950e534c40c03ddf4d91066;p=dealii.git Remove all remaining usages of the Boundary class. Requiescat in pace March 4, 1998 (commit 8: 184639502b4) - May 12, 2018 (41,058 commits and counting). Manifold has not gone too far from that original design: in the beginning we had the in_between function: template class Boundary { public: /** * Typedef an array of the needed number * of old points. */ typedef const Point* PointArray[1<<(dim-1)]; /** * This function calculates the position * of the new vertex. */ virtual Point in_between (const PointArray &neighbors) const = 0; }; /** Specialisation of \Ref{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. */ template class StraightBoundary : public Boundary { public: /** * This function calculates the position * of the new vertex. */ virtual Point in_between (const PointArray &neighbors) const { Point p; for (int i=0; i<(1<<(dim-1)); ++i) p += *neighbors[i]; p /= (1<<(dim-1))*1.0; return p; }; }; which is essentially Manifold::get_intermediate_point without the weight. --- diff --git a/contrib/python-bindings/include/triangulation_wrapper.h b/contrib/python-bindings/include/triangulation_wrapper.h index a36f8c62cc..440ad6d799 100644 --- a/contrib/python-bindings/include/triangulation_wrapper.h +++ b/contrib/python-bindings/include/triangulation_wrapper.h @@ -231,8 +231,8 @@ namespace python * Generate a hyper-ball intersected with the positive orthant relate to @p * center, which contains three elements in 2d and four in 3d. The boundary * indicators for the final triangulations are 0 for the curved boundary - * and 1 for the cut plane. The appropriate boundary class is - * HyperBallBoundary. + * and 1 for the cut plane. The appropriate manifold class is + * SphericalManifold. */ void generate_quarter_hyper_ball(PointWrapper ¢er, const double radius = 1.); @@ -241,8 +241,8 @@ namespace python * Generate a half hyper-ball around @p center, which contains four elements * in 2d and 6 in 3d. The cut plane is perpendicular to the x-axis. The * boundary indicators for the final triangulation are 0 for the curved - * boundary and 1 for the cut plane. The appropriate boundary class is - * HalfHyperBallBoundary, or HyperBallBoundary. + * boundary and 1 for the cut plane. The appropriate manifold class is + * SphericalManifold. */ void generate_half_hyper_ball(PointWrapper ¢er, const double radius = 1.); diff --git a/contrib/python-bindings/source/export_triangulation.cc b/contrib/python-bindings/source/export_triangulation.cc index e4d9402df3..83c94df619 100644 --- a/contrib/python-bindings/source/export_triangulation.cc +++ b/contrib/python-bindings/source/export_triangulation.cc @@ -225,9 +225,9 @@ namespace python const char generate_quarter_hyper_ball_docstring [] = "Generate a hyper-ball intersected with the positive orthant relate to \n" "center, which contains three elements in 2d and four in 3d. The \n" - "boundary indicators for the final triangulations are 0 for the curved \n" - "boundary and 1 for the cut plane. The appropriate boundary class is \n" - "HyperBallBoundary. \n" + "boundary indicators for the final triangulations are 0 for the curved \n" + "boundary and 1 for the cut plane. The appropriate Manifold class is \n" + "SphericalManifold. \n" ; @@ -236,8 +236,8 @@ namespace python "Generate a half hyper-ball around center, which contains four \n" "elements in 2d and 6 in 3d. The cut plane is perpendicular to the \n" "x-axis. The boundary indicators for the final triangulation are 0 for \n" - "the curved boundary and 1 for the cut plane. The appropriate boundary \n" - "class is HalfHyperBallBoundary, or HyperBallBoundary. \n" + "the curved boundary and 1 for the cut plane. The appropriate manifold \n" + "class is SphericalManifold. \n" ; diff --git a/doc/doxygen/headers/boundary.h b/doc/doxygen/headers/boundary.h deleted file mode 100644 index 3945fd131b..0000000000 --- a/doc/doxygen/headers/boundary.h +++ /dev/null @@ -1,131 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2003 - 2018 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - -/** - * @defgroup boundary Boundary and manifold description for triangulations - * - * @warning This module describes usage of the Boundary classes, which are - * deprecated: see the - * @ref manifold - * module for detailed information on how to describe curved boundaries with - * the Manifold classes, which are the replacements for the Boundary classes. - * - * The classes in this module are concerned with the description of the - * geometry of a domain in which a Triangulation lives. This geometry - * description is necessary in three contexts: - *
    - * - *
  • Mesh refinement: Whenever a cell is refined, it is necessary - * to introduce at least one new vertex. In the simplest case, one - * assumes that the cells and their faces consists of straight line - * segments, bilinear surface or trilinear volumes between the - * vertices of the original, coarsest mesh, and the next vertex is - * simply put into the middle of the old ones. This is the default - * behavior of the Triangulation class, and is described by the - * StraightBoundary and FlatManifold classes. - * - * On the other hand, if one deals with curved geometries and - * boundaries, this is not the appropriate thing to do. The classes - * derived from the Manifold and Boundary base classes describe the - * geometry of a domain. One can then attach an object of a class - * derived from this base classes to the Triangulation object using - * the Triangulation::set_boundary() or - * Triangulation::set_manifold() functions, and the Triangulation - * will ask the manifold object where a new vertex should be located - * upon mesh refinement. Several classes already exist to support - * the boundary of the most common geometries, e.g., - * CylinderBoundary, HyperBallBoundary, or HyperShellBoundary. - * - *
  • Integration: When using higher order finite element methods, - * it is often necessary to compute cell terms (like cell - * contributions to the matrix and right hand side of the linear - * system) using curved approximations of the geometry, rather than - * the straight line approximation. The actual implementation of - * such curved elements happens in the Mapping class (see the @ref - * mapping module), which however obtains its information about the - * manifold description from the classes described here. The same - * is, of course, true when integrating boundary terms (e.g., - * inhomogenous Neumann boundary conditions). - * - *
  • In cases where a Triangulation is embedded into a higher - * dimensional space, i.e., whenever the second template argument of - * the Triangulation class is explicitly specified and larger than - * the first (for an example, see step-34), the manifold description - * objects serve as a tool to describe the geometry not only of the - * boundary of the domain but of the domain itself, in case the - * domain is a manifold that is in fact curved. In these cases, one - * can use the Triangulation::set_manifold() function to indicate - * what manifold description to use when refining the curve, or when - * computing integrals using high order mappings. - * - *
- * - * In the context of triangulations, each object stores a number - * called manifold_id, and each face of a cell that is - * located at the boundary of the domain stores a number called - * boundary_id that uniquely identifies which part of the - * boundary this face is on. If nothing is specified at creation time, - * each boundary face has a zero boundary id and each triangulation - * object has a flat_manifold_id. On the other hand, the boundary - * id of faces and the manifold_id of objects can be set either at - * creation time or later by looping over all cells and querying their - * faces. - * - * It is then possible to associate objects describing the geometry to certain - * manifold_id values. - * - * Before version 8.2, the library allowed only boundary faces to follow a - * curved geometric description. Since version 8.2 this has been introduced also - * for interior faces and cells, and the boundary_id has been separated from the - * manifold_id. The former is used to identify the type of boundary conditions - * to apply, while the latter is used to identify the geometry and describe how - * new vertices should be created upon refinement of the mesh, or where high - * order Mapping objects should place their support points on the exact - * geometry. - * - * Since version 9.0 of the library, the boundary_id associated to the boundary - * faces is ignored by Manifold objects, and Manifold descriptors can only be - * attached to manifold ids. - * - * The behavior of the Triangulation class with respect to geometry descriptions - * is the following: Triangulation::set_manifold() attaches a manifold - * descriptor to the specified manifold_id. The function expects a Manifold - * descriptor, and you could describe both the interior and the boundary of the - * domain using the same object. - * - * Whenever a new vertex is needed in an object, the Triangulation queries the - * manifold_id of the object which needs refinement. If the query resulted in a - * number different from numbers::flat_manifold_id, then the Triangulation looks - * whether a previous call to Triangulation::set_manifold() was performed with - * the given id. If it was, then the triangulation uses the stored object to - * obtain new vertices, otherwise it uses a FlatManifold object. - * - * @note This behavior is **not** backward compatible to that of deal.II - * versions prior to 9.0. If one ignores the manifold_id of an object (i.e., if - * it has never been set), by default it is and remains set to - * numbers::flat_manifold_id. In previous versions of the library, the - * boundary_id of the boundary faces would be queried and used in place of the - * manifold_id. If you have old programs that only set boundary ids, you should - * modify them to use manifold ids instead (or you could use - * GridTools::copy_boundary_to_manifold_ids or - * GridTools::map_boundary_to_manifold_ids) - * - * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" - * @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators" - * - * @ingroup grid - * @author Wolfgang Bangerth, Luca Heltai, 1998-2018 - */ diff --git a/examples/step-32/doc/intro.dox b/examples/step-32/doc/intro.dox index 7b125c0b1c..55e60a1c2e 100644 --- a/examples/step-32/doc/intro.dox +++ b/examples/step-32/doc/intro.dox @@ -885,8 +885,7 @@ the following quantities: outer earth core starts). The radii are therefore $R_0=(6371-2890)\text{km}, R_1=(6371-35)\text{km}$. This domain is conveniently generated using the - GridGenerator::hyper_shell() function, and we use a HyperShellBoundary - objects for the inner and outer boundary. + GridGenerator::hyper_shell() function.
  • At the interface between crust and mantle, the temperature is between 500 and 900 degrees Celsius, whereas at its bottom it is around 4000 degrees diff --git a/examples/step-59/step-59.cc b/examples/step-59/step-59.cc index 248cc96ec9..98f6170ca6 100644 --- a/examples/step-59/step-59.cc +++ b/examples/step-59/step-59.cc @@ -39,7 +39,6 @@ #include #include #include -#include #include #include diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h index 7ea247c605..c4e1e60769 100644 --- a/include/deal.II/grid/tria.h +++ b/include/deal.II/grid/tria.h @@ -42,8 +42,6 @@ DEAL_II_NAMESPACE_OPEN -template class Boundary; -template class StraightBoundary; template class Manifold; namespace GridTools @@ -1277,16 +1275,6 @@ private: typedef dealii::internal::TriangulationImplementation::Iterators IteratorSelector; public: - /** - * Default manifold object. This is used for those objects for which no - * boundary description has been explicitly set using set_manifold(). - * - * @deprecated This member variable has been deprecated in favor of creating - * an independent FlatManifold. - */ - DEAL_II_DEPRECATED - static const StraightBoundary straight_boundary; - /** * Declare some symbolic names for mesh smoothing algorithms. The meaning of * these flags is documented in the Triangulation class. @@ -1705,68 +1693,6 @@ public: */ virtual const MeshSmoothing &get_mesh_smoothing() const; - /** - * If @p dim==spacedim, assign a boundary object to a certain part of the - * boundary of a the triangulation. If a face with boundary number @p number - * is refined, this object is used to find the location of new vertices on - * the boundary (see the results section of step-49 for a more in-depth - * discussion of this, with examples). It is also used for non-linear - * (i.e.: non-Q1) transformations of cells to the unit cell in shape - * function calculations. - * - * If @p dim!=spacedim the boundary object is in fact the exact manifold - * that the triangulation is approximating (for example a circle - * approximated by a polygon triangulation). As above, the refinement is - * made in such a way that the new points are located on the exact manifold. - * - * Numbers of boundary objects correspond to material numbers of faces at - * the boundary, for instance the material id in a UCD input file. They are - * not necessarily consecutive but must be in the range - * 0-(types::boundary_id-1). Material IDs on boundaries are also called - * boundary indicators and are accessed with accessor functions of that - * name. - * - * The @p boundary_object is not copied and MUST persist until the - * triangulation is destroyed. This is also true for triangulations - * generated from this one by @p copy_triangulation. - * - * It is possible to remove or replace the boundary object during the - * lifetime of a non-empty triangulation. Usually, this is done before the - * first refinement and is dangerous afterwards. Removal of a boundary - * object is done by set_boundary(number), i.e. the function of - * same name but only one argument. This operation then replaces the - * boundary object given before by a straight boundary approximation. - * - * @ingroup boundary - * - * @deprecated This method has been deprecated. Use - * Triangulation::set_manifold() instead. - * - * @see - * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" - */ - DEAL_II_DEPRECATED - void set_boundary (const types::manifold_id number, - const Boundary &boundary_object); - - - /** - * Reset those parts of the boundary with the given number to use a straight - * boundary approximation. This is the default state of a triangulation, and - * undoes assignment of a different boundary object by the function of same - * name and two arguments. - * - * @ingroup boundary - * - * @deprecated This method has been deprecated. Use - * Triangulation::reset_manifold() instead. - * - * @see - * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" - */ - DEAL_II_DEPRECATED - void set_boundary (const types::manifold_id number); - /** * Assign a manifold object to a certain part of the triangulation. If * an object with manifold number @p number is refined, this object is used @@ -1870,22 +1796,6 @@ public: void set_all_manifold_ids_on_boundary (const types::boundary_id b_id, const types::manifold_id number); - - /** - * Return a constant reference to a boundary object used for this - * triangulation. Number is the same as in @p set_boundary - * - * @ingroup boundary - * - * @deprecated This method has been deprecated. Use - * Triangulation::get_manifold() instead. - * - * @see - * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" - */ - DEAL_II_DEPRECATED - const Boundary &get_boundary (const types::manifold_id number) const; - /** * Return a constant reference to a Manifold object used for this * triangulation. Number is the same as in @p set_manifold diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index 89e5199380..112d8ca654 100644 --- a/include/deal.II/grid/tria_accessor.h +++ b/include/deal.II/grid/tria_accessor.h @@ -36,7 +36,6 @@ template class TriaRawIterator; template class TriaIterator; template class TriaActiveIterator; -template class Boundary; template class Manifold; @@ -1018,16 +1017,6 @@ public: */ bool at_boundary () const; - /** - * Return a constant reference to the manifold object used for this object. - * This function exists for backward compatibility and calls get_manifold() - * internally. - * - * @deprecated The classes derived from Boundary have been deprecated in - * favor of the Manifold classes. - */ - DEAL_II_DEPRECATED const Boundary &get_boundary () const; - /** * Return a constant reference to the manifold object used for this object. * diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h index eba6866591..6abfff5312 100644 --- a/include/deal.II/grid/tria_accessor.templates.h +++ b/include/deal.II/grid/tria_accessor.templates.h @@ -1998,29 +1998,6 @@ TriaAccessor::at_boundary () const -template -const Boundary & -TriaAccessor::get_boundary () const -{ - Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed()); - - // Get the default (manifold_id) - const types::manifold_id mi = this->objects().manifold_id[this->present_index]; - - // In case this is not valid, check - // the boundary id, after having - // casted it to a manifold id - if (mi == numbers::invalid_manifold_id) - return this->tria->get_boundary(structdim < dim ? - this->objects().boundary_or_material_id[this->present_index].boundary_id: - dim < spacedim ? - this->objects().boundary_or_material_id[this->present_index].material_id: - numbers::invalid_manifold_id); - else - return this->tria->get_boundary(mi); -} - - template const Manifold & TriaAccessor::get_manifold () const diff --git a/include/deal.II/grid/tria_boundary.h b/include/deal.II/grid/tria_boundary.h deleted file mode 100644 index f3c3d3e4b5..0000000000 --- a/include/deal.II/grid/tria_boundary.h +++ /dev/null @@ -1,427 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 1998 - 2018 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - -#ifndef dealii_tria_boundary_h -#define dealii_tria_boundary_h - - -/*---------------------------- boundary-function.h ---------------------------*/ - -#include -#include -#include -#include -#include -#include -#include - -DEAL_II_NAMESPACE_OPEN - -/** - * This class is used to represent a boundary to a triangulation. 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); - * ... - * @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. - * - * 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)) }; - * ... - * @endcode - * to get 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); - * ... - * @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. - * - * 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. - * - * There are specializations, 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 This class has been deprecated in favor of the more general - * Manifold class. - * - * @ingroup boundary - * @author Wolfgang Bangerth, 1999, 2001, 2009, Ralf Hartmann, 2001, 2008, - * Luca Heltai, 2014 - */ -template -class DEAL_II_DEPRECATED Boundary : public FlatManifold -{ -public: - - /** - * Destructor. Does nothing here, but needs to be declared to make it - * virtual. - */ - virtual ~Boundary () override = default; - - /** - * Clone this Boundary. - */ - virtual std::unique_ptr > clone() const override; - - /** - * 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 derived classes to arrange the points in - * approximately equal distances along the length of the line segment on the - * boundary bounded by the vertices of the first argument. - * - * Among other places in the library, this function is called by the Mapping - * classes, for example the @p MappingQGeneric class. On the other hand, not - * all mapping classes actually require intermediate points on lines (for - * example, $Q_1$ mappings do not). Consequently this function is not made - * pure virtual, to allow users to define their own boundary classes without - * having to overload this function. However, the default implementation - * throws an error in any case and can, consequently, not be used if you use - * a mapping that does need the information provided by this function. - */ - 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. - * - * Among other places in the library, this function is called by the Mapping - * classes, for example the @p MappingQGeneric class. On the other hand, not - * all mapping classes actually require intermediate points on quads (for - * example, $Q_1$ mappings do not). Consequently this function is not made - * pure virtual, to allow users to define their own boundary classes without - * having to overload this function. However, the default implementation - * throws an error in any case and can, consequently, not be used if you use - * a mapping that does need the information provided by this function. - */ - 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; - - /** - * 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: - /** - * Return 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; -}; - - - -/** - * 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. - * - * This class does not really describe a boundary in the usual sense. By - * placing new points in the middle of old ones, it rather assumes that the - * boundary of the domain is given by the polygon/polyhedron defined by the - * boundary of the initial coarse triangulation. - * - * @deprecated This class has been deprecated in favor of the more general - * FlatManifold class. - * - * @ingroup boundary - * - * @author Wolfgang Bangerth, 1998, 2001, Ralf Hartmann, 2001 - */ -template -class DEAL_II_DEPRECATED StraightBoundary : public Boundary -{ -public: - /** - * Default constructor. Some compilers require this for some reasons. - */ - StraightBoundary (); - - /** - * Clone this Boundary. - */ - virtual std::unique_ptr > clone() const override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * Return n=points.size() points that split 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 override; - - /** - * Return 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; -}; - - - -/* -------------- declaration of explicit specializations ------------- */ - -#ifndef DOXYGEN - -template <> -void -Boundary<1,1>:: -get_intermediate_points_on_face (const Triangulation<1,1>::face_iterator &, - std::vector > &) const; - -template <> -void -Boundary<1,2>:: -get_intermediate_points_on_face (const Triangulation<1,2>::face_iterator &, - std::vector > &) const; - -template <> -void -Boundary<1,3>:: -get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &, - std::vector > &) const; - -template <> -Point<3> -StraightBoundary<3,3>:: -get_new_point_on_quad (const Triangulation<3,3>::quad_iterator &quad) 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_NAMESPACE_CLOSE - -#endif diff --git a/include/deal.II/grid/tria_boundary_lib.h b/include/deal.II/grid/tria_boundary_lib.h deleted file mode 100644 index ff256667a2..0000000000 --- a/include/deal.II/grid/tria_boundary_lib.h +++ /dev/null @@ -1,841 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 1999 - 2018 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - -#ifndef dealii_tria_boundary_lib_h -#define dealii_tria_boundary_lib_h - - -#include -#include - -DEAL_II_NAMESPACE_OPEN - -/** - * Boundary object for the hull of a cylinder. In three dimensions, points - * are projected on a circular tube along the x-, y- or - * z-axis (when using the first constructor of this class), or an - * arbitrarily oriented cylinder described by the direction of its axis and a - * point located on the axis. The radius of the tube can be given - * independently. Similar to HyperBallBoundary, new points are projected by - * dividing the straight line between the old two points and adjusting the - * radius from the axis. - * - * This class was developed to be used in conjunction with the @p cylinder - * function of GridGenerator. It should be used for the hull of the cylinder - * only (boundary indicator 0). Its use is discussed in detail in the results - * section of step-49. - * - * 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. - * - * @deprecated The boundary classes have been deprecated in favor of the more - * general Manifold classes. For this class CylindricalManifold is an appropriate - * replacement which is expected to produce better results under mesh refinement. - * For a two-dimensional cylinder use FlatManifold instead. - * - * @ingroup boundary - * - * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007 - */ -template -class DEAL_II_DEPRECATED CylinderBoundary : public StraightBoundary -{ -public: - /** - * Constructor. Using default values for the constructor arguments yields a - * circular tube along the x-axis (axis=0). Choose axis=1 - * or axis=2 for a tube along the y- or z-axis, respectively. - */ - CylinderBoundary (const double radius = 1.0, - 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. - */ - CylinderBoundary (const double radius, - const Point &direction, - const Point &point_on_axis); - - /** - * Clone this Boundary object. - */ - virtual std::unique_ptr > clone() const override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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. - */ - DeclException0 (ExcRadiusNotSet); - - -protected: - /** - * Radius of the cylinder. - */ - const double radius; - - /** - * The direction vector of the axis. - */ - const Point direction; - - /** - * An arbitrary point on the axis. - */ - const Point point_on_axis; - -private: - - /** - * Called by @p get_intermediate_points_on_line and by @p - * get_intermediate_points_on_quad. - * - * Refer to the general documentation of @p get_intermediate_points_on_line - * in the documentation of the base class. - */ - void get_intermediate_points_between_points (const Point &p0, const Point &p1, - std::vector > &points) const; - - /** - * Given a number for the axis, return a vector that denotes this direction. - */ - static Point get_axis_vector (const unsigned int axis); -}; - - - -/** - * Boundary object for the hull of a (truncated) cone with two different radii - * at the two ends. If one radius is chosen to be 0 the object describes the - * boundary of a cone. In three dimensions, points are projected on an - * arbitrarily oriented (truncated) cone described by the two endpoints and - * the corresponding radii. Similar to HyperBallBoundary, new points are - * projected by dividing the straight line between the old two points and - * adjusting the radius from the axis. - * - * 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. - * - * As an example of use, consider the following code snippet: - * @code - * Triangulation triangulation; - * GridGenerator::truncated_cone (triangulation); - * Point p1, p2; - * p1[0] = -1; - * p2[0] = 1; - * const ConeBoundary boundary (1, 0.5, p1, p2); - * triangulation.set_boundary (0, boundary); - * triangulation.refine_global (2); - * @endcode - * This will produce the following meshes after the two refinements we - * perform, in 2d and 3d, respectively: - * - * @image html cone_2d.png - * @image html cone_3d.png - * - * @deprecated The boundary classes have been deprecated in favor of the more - * general Manifold classes. For this class CylindricalManifold is an appropriate - * replacement for spacedim==3, else StraightManifold can be used. - * - * @author Markus Bürg, 2009 - */ -template -class DEAL_II_DEPRECATED 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. - */ - ConeBoundary (const double radius_0, - const double radius_1, - const Point x_0, - const Point x_1); - - /** - * Clone this Boundary object. - */ - virtual std::unique_ptr > clone() const override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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,dim> - normal_vector (const typename Triangulation::face_iterator &face, - const Point &p) const override; - - -protected: - /** - * First radius of the (truncated) cone. - */ - const double radius_0; - - /** - * Second radius of the (truncated) cone. - */ - const double radius_1; - - /** - * Starting point of the axis. - */ - const Point x_0; - - /** - * 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; -}; - - - -/** - * Specialization of Boundary, which places the new point on the boundary - * of a ball in arbitrary dimension. It works by projecting the point in the - * middle of the old points onto the ball. The middle is defined as the - * arithmetic mean of the points. - * - * 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 - * - * @deprecated The boundary classes have been deprecated in favor of the more - * general Manifold classes. For this class SphericalManifold is an appropriate - * replacement which is expected to produce better results under mesh refinement. - * - * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001 - */ -template -class DEAL_II_DEPRECATED HyperBallBoundary : public StraightBoundary -{ -public: - /** - * Constructor - */ - HyperBallBoundary (const Point p = Point(), - const double radius = 1.0); - - /** - * Clone this Boundary object. - */ - virtual std::unique_ptr > clone() const override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * Return the center of the ball. - */ - Point - get_center () const; - - /** - * Return the radius of the ball. - */ - double - get_radius () const; - - /** - * Exception. Thrown by the @p get_radius if the @p - * compute_radius_automatically, see below, flag is set true. - */ - DeclException0 (ExcRadiusNotSet); - - -protected: - - /** - * Center point of the hyperball. - */ - const Point center; - - /** - * Radius of the hyperball. - */ - 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 called. - */ - bool compute_radius_automatically; - -private: - - /** - * Called by @p get_intermediate_points_on_line and by @p - * get_intermediate_points_on_quad. - * - * Refer to the general documentation of @p get_intermediate_points_on_line - * in the documentation of the base class. - */ - void get_intermediate_points_between_points (const Point &p0, const Point &p1, - std::vector > &points) const; -}; - - - -/** - * Variant of HyperBallBoundary which denotes a half hyper ball where the - * first coordinate is restricted to the range $x>=0$ (or $x>=center(0)$). In - * two dimensions, this equals the right half circle, in three space - * dimensions it is a half ball. This class might be useful for computations - * with rotational symmetry, where one dimension is the radius from the axis - * of rotation. - * - * @ingroup boundary - * - * @deprecated The boundary classes have been deprecated in favor of the more - * general Manifold classes. For this class SphericalManifold is an appropriate - * replacement which is expected to produce better results under mesh refinement. - * - * @author Wolfgang Bangerth, 1999, 2001 - */ -template -class DEAL_II_DEPRECATED HalfHyperBallBoundary : public HyperBallBoundary -{ -public: - /** - * Constructor - */ - HalfHyperBallBoundary (const Point p = Point(), - const double radius = 1.0); - - /** - * Clone this Boundary object. - */ - virtual std::unique_ptr > clone() const override; - - /** - * 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 override; - - /** - * Check if on the line x==0, otherwise pass to the base class. - */ - virtual Point - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; -}; - - - -/** - * Class describing the boundaries of a hyper shell. Only the center of the - * two spheres needs to be given, the radii of inner and outer sphere are - * computed automatically upon calling one of the virtual functions. - * - * @ingroup boundary - * - * @deprecated The boundary classes have been deprecated in favor of the more - * general Manifold classes. For this class SphericalManifold is an appropriate - * replacement which is expected to produce better results under mesh refinement. - * - * @author Wolfgang Bangerth, 1999 - */ -template -class DEAL_II_DEPRECATED HyperShellBoundary : public HyperBallBoundary -{ -public: - /** - * Constructor. The center of the spheres defaults to the origin. - * - * Call the constructor of its base @p HyperBallBoundary class with a dummy - * radius as argument. This radius will be ignored - */ - HyperShellBoundary (const Point ¢er = Point()); - - /** - * Clone this Boundary object. - */ - virtual std::unique_ptr > clone() const override; -}; - - - -/** - * Variant of HyperShellBoundary which denotes a half hyper shell where the - * first coordinate is restricted to the range $x>=0$ (or $x>=center(0)$). In - * two dimensions, this equals the right half arc, in three space dimensions - * it is a half shell. This class might be useful for computations with - * rotational symmetry, where one dimension is the radius from the axis of - * rotation. - * - * @ingroup boundary - * - * @deprecated The boundary classes have been deprecated in favor of the more - * general Manifold classes. For this class SphericalManifold is an appropriate - * replacement which is expected to produce better results under mesh refinement. - * - * @author Wolfgang Bangerth, 2000, 2009 - */ -template -class DEAL_II_DEPRECATED HalfHyperShellBoundary : public HyperShellBoundary -{ -public: - /** - * Constructor. The center of the spheres defaults to the origin. - * - * If the radii are not specified, the class tries to infer them from the - * location of points on the boundary. This works in 2d, but not in 3d. As a - * consequence, in 3d these radii must be given. - */ - HalfHyperShellBoundary (const Point ¢er = Point(), - const double inner_radius = -1, - const double outer_radius = -1); - - /** - * Clone this Boundary object. - */ - virtual std::unique_ptr > clone() const override; - - /** - * Construct a new point on a line. - */ - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const override; - - /** - * Construct a new point on a quad. - */ - virtual Point - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const override; - - /** - * 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 override; - - /** - * 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 override; - - /** - * 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 override; - -private: - /** - * Inner and outer radii of the shell. - */ - const double inner_radius; - const double outer_radius; -}; - - -/** - * Class describing the boundary of the torus. The axis of the torus is the - * $y$-axis while the plane of the torus is the $x$-$z$ plane. A torus of this - * kind can be generated by GridGenerator::torus. - * - * This class is only implemented for dim=2,spacedim=3, that - * is, just the surface. - * - * @deprecated This class has been deprecated in favor of the more general - * TorusManifold class. - */ -template -class DEAL_II_DEPRECATED TorusBoundary : public Boundary -{ -public: - /** - * Constructor.R has to be greater than r. - */ - TorusBoundary (const double R, const double r); - - /** - * Clone this Boundary object. - */ - virtual std::unique_ptr > clone() const override; - -//Boundary Refinement Functions - /** - * Construct a new point on a line. - */ - virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const override; - - /** - * Construct a new point on a quad. - */ - virtual Point - get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const override; - - /** - * Construct a new points on a line. - */ - virtual void get_intermediate_points_on_line ( - const typename Triangulation< dim, spacedim >::line_iterator &line, - std::vector< Point< spacedim > > &points) const override; - - /** - * Construct a new points on a quad. - */ - virtual void get_intermediate_points_on_quad ( - const typename Triangulation< dim, spacedim >::quad_iterator &quad, - std::vector< Point< spacedim > > &points ) const override; - - /** - * 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 override; - -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. - */ - const double R; - const double r; -}; - - - -/* -------------- 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/source/dofs/dof_accessor.cc b/source/dofs/dof_accessor.cc index 42e8c7e85e..ec61cb7ac9 100644 --- a/source/dofs/dof_accessor.cc +++ b/source/dofs/dof_accessor.cc @@ -17,7 +17,6 @@ #include #include #include -#include #include #include #include diff --git a/source/dofs/dof_accessor_get.cc b/source/dofs/dof_accessor_get.cc index d823552b19..83087b67eb 100644 --- a/source/dofs/dof_accessor_get.cc +++ b/source/dofs/dof_accessor_get.cc @@ -28,7 +28,6 @@ #include #include #include -#include #include #include #include diff --git a/source/dofs/dof_accessor_set.cc b/source/dofs/dof_accessor_set.cc index e4eb96f896..3ab878e32c 100644 --- a/source/dofs/dof_accessor_set.cc +++ b/source/dofs/dof_accessor_set.cc @@ -28,7 +28,6 @@ #include #include #include -#include #include #include #include diff --git a/source/fe/fe.cc b/source/fe/fe.cc index 6110359bc7..138e40671c 100644 --- a/source/fe/fe.cc +++ b/source/fe/fe.cc @@ -22,7 +22,6 @@ #include #include #include -#include #include #include diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index 2b2d0bfcf3..1218751bdf 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -33,7 +33,6 @@ #include #include #include -#include #include #include #include diff --git a/source/fe/mapping_c1.cc b/source/fe/mapping_c1.cc index 73d089152f..4ae4bba504 100644 --- a/source/fe/mapping_c1.cc +++ b/source/fe/mapping_c1.cc @@ -15,9 +15,9 @@ #include +#include #include #include -#include #include diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index 3009f8cc8f..3f5a57172d 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -34,7 +34,6 @@ #include #include #include -#include #include #include #include diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 01fdcffeae..1c6edd9c20 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -25,7 +25,6 @@ #include #include #include -#include #include #include #include diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index bff949b15b..c74c284523 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -27,7 +27,6 @@ #include #include #include -#include #include #include #include @@ -3612,11 +3611,9 @@ add_line_support_points (const typename Triangulation::cell_iterat } else // otherwise call the more complicated functions and ask for inner points - // from the boundary description + // from the manifold description { std::vector > tmp_points; - // loop over each of the lines, and if it is at the boundary, then first - // get the boundary description and second compute the points on it for (unsigned int line_no=0; line_no::lines_per_cell; ++line_no) { const typename Triangulation::line_iterator @@ -3633,34 +3630,20 @@ add_line_support_points (const typename Triangulation::cell_iterat cell->get_manifold() : line->get_manifold() ); - if (const Boundary *boundary - = dynamic_cast *>(&manifold)) + const std::array, 2> vertices + { { - tmp_points.resize(this->polynomial_degree-1); - boundary->get_intermediate_points_on_line(line, tmp_points); - if (dim != 3 || cell->line_orientation(line_no)) - a.insert (a.end(), tmp_points.begin(), tmp_points.end()); - else - a.insert (a.end(), tmp_points.rbegin(), tmp_points.rend()); + cell->vertex(GeometryInfo::line_to_cell_vertices(line_no, 0)), + cell->vertex(GeometryInfo::line_to_cell_vertices(line_no, 1)) } - else - { - const std::array, 2> vertices - { - { - cell->vertex(GeometryInfo::line_to_cell_vertices(line_no, 0)), - cell->vertex(GeometryInfo::line_to_cell_vertices(line_no, 1)) - } - }; + }; - const std::size_t n_rows = support_point_weights_perimeter_to_interior[0].size(0); - a.resize(a.size() + n_rows); - auto a_view = make_array_view(a.end() - n_rows, a.end()); - manifold.get_new_points(make_array_view(vertices.begin(), - vertices.end()), - support_point_weights_perimeter_to_interior[0], - a_view); - } + const std::size_t n_rows = support_point_weights_perimeter_to_interior[0].size(0); + a.resize(a.size() + n_rows); + auto a_view = make_array_view(a.end() - n_rows, a.end()); + manifold.get_new_points(make_array_view(vertices.begin(), vertices.end()), + support_point_weights_perimeter_to_interior[0], + a_view); } } } @@ -3683,12 +3666,10 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell, { const Triangulation<3>::face_iterator face = cell->face(face_no); - // select the correct mappings for the present face +#ifdef DEBUG const bool face_orientation = cell->face_orientation(face_no), face_flip = cell->face_flip (face_no), face_rotation = cell->face_rotation (face_no); - -#ifdef DEBUG const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face, lines_per_face = GeometryInfo<3>::lines_per_face; @@ -3708,54 +3689,29 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell, face_no, i, face_orientation, face_flip, face_rotation)), ExcInternalError()); #endif - - // On a quad, we have to check whether the manifold should determine the - // point distribution from all surrounding points (new manifold code) or - // the old-style Boundary code should simply return the intermediate - // points. The second check is to find out whether the Boundary object - // is actually a StraightBoundary (the default flat manifold assigned to - // the triangulation if no manifold is assigned). - const Boundary<3,3> *boundary = - dynamic_cast *>(&face->get_manifold()); - if (boundary != nullptr && - std::string(typeid(*boundary).name()).find("StraightBoundary") == - std::string::npos) - { - // ask the boundary/manifold object to return intermediate points on it - tmp_points.resize((polynomial_degree-1)*(polynomial_degree-1)); - boundary->get_intermediate_points_on_quad(face, tmp_points); - for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation(i, - face_orientation, - face_flip, - face_rotation)]); - } - else - { - // extract the points surrounding a quad from the points - // already computed. First get the 4 vertices and then the points on - // the four lines - boost::container::small_vector, 200> - tmp_points(GeometryInfo<2>::vertices_per_cell - + GeometryInfo<2>::lines_per_cell*(polynomial_degree-1)); - for (unsigned int v=0; v::vertices_per_cell; ++v) - tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no,v)]; - if (polynomial_degree > 1) - for (unsigned int line=0; line::lines_per_cell; ++line) - for (unsigned int i=0; i::vertices_per_cell + - (polynomial_degree-1)* - GeometryInfo<3>::face_to_cell_lines(face_no,line) + i]; - - const std::size_t n_rows = support_point_weights_perimeter_to_interior[1].size(0); - a.resize(a.size() + n_rows); - auto a_view = make_array_view(a.end() - n_rows, a.end()); - face->get_manifold().get_new_points (make_array_view(tmp_points.begin(), - tmp_points.end()), - support_point_weights_perimeter_to_interior[1], - a_view); - } + // extract the points surrounding a quad from the points + // already computed. First get the 4 vertices and then the points on + // the four lines + boost::container::small_vector, 200> + tmp_points(GeometryInfo<2>::vertices_per_cell + + GeometryInfo<2>::lines_per_cell*(polynomial_degree-1)); + for (unsigned int v=0; v::vertices_per_cell; ++v) + tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no,v)]; + if (polynomial_degree > 1) + for (unsigned int line=0; line::lines_per_cell; ++line) + for (unsigned int i=0; i::vertices_per_cell + + (polynomial_degree-1)* + GeometryInfo<3>::face_to_cell_lines(face_no,line) + i]; + + const std::size_t n_rows = support_point_weights_perimeter_to_interior[1].size(0); + a.resize(a.size() + n_rows); + auto a_view = make_array_view(a.end() - n_rows, a.end()); + face->get_manifold().get_new_points (make_array_view(tmp_points.begin(), + tmp_points.end()), + support_point_weights_perimeter_to_interior[1], + a_view); } } @@ -3767,38 +3723,28 @@ MappingQGeneric<2,3>:: add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, std::vector > &a) const { - if (const Boundary<2,3> *boundary = - dynamic_cast *>(&cell->get_manifold())) - { - std::vector > points((polynomial_degree-1)*(polynomial_degree-1)); - boundary->get_intermediate_points_on_quad(cell, points); - a.insert(a.end(), points.begin(), points.end()); - } - else - { - std::array, GeometryInfo<2>::vertices_per_cell> vertices; - for (unsigned int i=0; i::vertices_per_cell; ++i) - vertices[i] = cell->vertex(i); - - Table<2,double> weights(Utilities::fixed_power<2>(polynomial_degree-1), - GeometryInfo<2>::vertices_per_cell); - for (unsigned int q=0, q2=0; q2 point(line_support_points.point(q1+1)[0], - line_support_points.point(q2+1)[0]); - for (unsigned int i=0; i::vertices_per_cell; ++i) - weights(q,i) = GeometryInfo<2>::d_linear_shape_function(point, i); - } - // TODO: use all surrounding points once Boundary path is removed - const std::size_t n_rows = weights.size(0); - a.resize(a.size() + n_rows); - auto a_view = make_array_view(a.end() - n_rows, a.end()); - cell->get_manifold().get_new_points(make_array_view(vertices.begin(), - vertices.end()), - weights, - a_view); - } + std::array, GeometryInfo<2>::vertices_per_cell> vertices; + for (unsigned int i=0; i::vertices_per_cell; ++i) + vertices[i] = cell->vertex(i); + + Table<2,double> weights(Utilities::fixed_power<2>(polynomial_degree-1), + GeometryInfo<2>::vertices_per_cell); + for (unsigned int q=0, q2=0; q2 point(line_support_points.point(q1+1)[0], + line_support_points.point(q2+1)[0]); + for (unsigned int i=0; i::vertices_per_cell; ++i) + weights(q,i) = GeometryInfo<2>::d_linear_shape_function(point, i); + } + + const std::size_t n_rows = weights.size(0); + a.resize(a.size() + n_rows); + auto a_view = make_array_view(a.end() - n_rows, a.end()); + cell->get_manifold().get_new_points(make_array_view(vertices.begin(), + vertices.end()), + weights, + a_view); } diff --git a/source/grid/CMakeLists.txt b/source/grid/CMakeLists.txt index c73798600e..6de7024fd7 100644 --- a/source/grid/CMakeLists.txt +++ b/source/grid/CMakeLists.txt @@ -27,8 +27,6 @@ SET(_unity_include_src manifold_lib.cc persistent_tria.cc tria_accessor.cc - tria_boundary.cc - tria_boundary_lib.cc tria_faces.cc tria_levels.cc tria_objects.cc @@ -63,8 +61,6 @@ SET(_inst manifold.inst.in manifold_lib.inst.in tria_accessor.inst.in - tria_boundary.inst.in - tria_boundary_lib.inst.in tria.inst.in tria_objects.inst.in ) diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 693d00d854..44029d51a6 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -23,7 +23,6 @@ #include #include #include -#include #include #include #include @@ -8791,11 +8790,6 @@ namespace internal } -template -const StraightBoundary -Triangulation::straight_boundary = StraightBoundary(); - - template const unsigned int @@ -8942,14 +8936,6 @@ Triangulation::get_mesh_smoothing() const -template -void -Triangulation::set_boundary (const types::manifold_id m_number, - const Boundary &boundary_object) -{ - set_manifold(m_number, boundary_object); -} - template void Triangulation::set_manifold (const types::manifold_id m_number, @@ -8962,13 +8948,6 @@ Triangulation::set_manifold (const types::manifold_id m_number, } -template -void -Triangulation::set_boundary (const types::manifold_id m_number) -{ - set_manifold(m_number); -} - template void @@ -9065,18 +9044,6 @@ Triangulation::set_all_manifold_ids_on_boundary (const types::bou } -template -const Boundary & -Triangulation::get_boundary (const types::manifold_id m_number) const -{ - const Boundary *man = - dynamic_cast *>(&get_manifold(m_number)); - Assert(man != nullptr, - ExcMessage("You tried to get a Boundary, but I only have a Manifold.")); - - return *man; -} - template const Manifold & diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 7fc45b7ffb..f6339a5626 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -19,7 +19,6 @@ #include #include #include -#include #include #include #include @@ -1047,13 +1046,7 @@ namespace Point get_new_point_on_object(const TriaAccessor &obj, const bool use_interpolation) { - // if we should not interpolate from the surroundings of the current - // object or if the object is of the old Boundary type, only use the - // vertices points - if (use_interpolation == false || - dynamic_cast *>(&obj.get_manifold()) != nullptr) - return get_new_point_on_object(obj); - else + if (use_interpolation) { TriaRawIterator > it(obj); const auto points_and_weights = Manifolds::get_default_points_and_weights(it, use_interpolation); @@ -1062,6 +1055,10 @@ namespace make_array_view(points_and_weights.second.begin(), points_and_weights.second.end())); } + else + { + return get_new_point_on_object(obj); + } } } diff --git a/source/grid/tria_boundary.cc b/source/grid/tria_boundary.cc deleted file mode 100644 index 6cbab1f677..0000000000 --- a/source/grid/tria_boundary.cc +++ /dev/null @@ -1,658 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 1998 - 2018 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - -#include -#include -#include -#include -#include -#include -#include -#include - -DEAL_II_NAMESPACE_OPEN - - - -/* -------------------------- Boundary --------------------- */ - - -template -std::unique_ptr > -Boundary::clone() const -{ - AssertThrow(false, ExcPureFunctionCalled()); - return std::unique_ptr >(); -} - - - -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 -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 <> -void -Boundary<1,1>:: -get_intermediate_points_on_face (const Triangulation<1,1>::face_iterator &, - std::vector > &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} - - -template <> -void -Boundary<1,2>:: -get_intermediate_points_on_face (const Triangulation<1,2>::face_iterator &, - std::vector > &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} - - -template <> -void -Boundary<1,3>:: -get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &, - std::vector > &) const -{ - Assert (false, ExcImpossibleInDim(1)); -} - - -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() == nullptr) - { - 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() == nullptr) - points[n_intermediate_points] = std_cxx14::make_unique > - (n_intermediate_points+2); - } - return points[n_intermediate_points]->get_points(); -} - - - - -/* -------------------------- StraightBoundary --------------------- */ - -// At least clang < 3.9.0 complains if we move this definition to its -// declaration when a 'const StraightBoundary' object is built. -template -StraightBoundary::StraightBoundary () = default; - - - -template -Point -StraightBoundary:: -get_new_point_on_line (const typename Triangulation::line_iterator &line) const -{ - return (line->vertex(0) + line->vertex(1)) / 2; -} - - -template -std::unique_ptr > -StraightBoundary::clone() const -{ - return std_cxx14::make_unique >(); -} - - -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 &quad) const -{ - return FlatManifold::get_new_point_on_quad(quad); -} - - -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:: -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,spacedim> -StraightBoundary:: -normal_vector (const typename Triangulation::face_iterator &face, - const Point &p) const -{ - return FlatManifold::normal_vector(face, p); -} - - - -template -void -StraightBoundary:: -get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const -{ - FlatManifold::get_normals_at_vertices(face, face_vertex_normals); -} - - - -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, - std::integral_constant) - { - // 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( - 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 Tensor<1,dim> 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: */std::integral_constant) - { - return y; - } - - // specialization for a quad in 2d - template - Point<2> - compute_projection (const Iterator &, - const Point<2> &y, - /* it's a quad: */std::integral_constant) - { - 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 */std::integral_constant()); -} - - - -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" - -DEAL_II_NAMESPACE_CLOSE diff --git a/source/grid/tria_boundary.inst.in b/source/grid/tria_boundary.inst.in deleted file mode 100644 index a896f6366c..0000000000 --- a/source/grid/tria_boundary.inst.in +++ /dev/null @@ -1,27 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 1998 - 2016 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - - -for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) -{ -#if deal_II_dimension <= deal_II_space_dimension - template class Boundary; - template class StraightBoundary; -#endif -} - - - diff --git a/source/grid/tria_boundary_lib.cc b/source/grid/tria_boundary_lib.cc deleted file mode 100644 index dc87dcd737..0000000000 --- a/source/grid/tria_boundary_lib.cc +++ /dev/null @@ -1,1594 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 1999 - 2018 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - -#include -#include -#include -#include -#include -#include -#include - - - -DEAL_II_NAMESPACE_OPEN - - - -template -CylinderBoundary::CylinderBoundary (const double radius, - const unsigned int axis) - : - radius(radius), - direction (get_axis_vector (axis)), - point_on_axis (Point()) -{} - - -template -CylinderBoundary::CylinderBoundary (const double radius, - const Point &direction, - const Point &point_on_axis) - : - radius(radius), - direction (direction / direction.norm()), - point_on_axis (point_on_axis) -{} - - - -template -std::unique_ptr > -CylinderBoundary::clone() const -{ - return std_cxx14::make_unique >(radius, direction, point_on_axis); -} - - - -template -Point -CylinderBoundary::get_axis_vector (const unsigned int axis) -{ - Assert (axis < spacedim, ExcIndexRange (axis, 0, spacedim)); - - Point axis_vector; - axis_vector[axis] = 1; - return axis_vector; -} - - - -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 Tensor<1,spacedim> 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 Point(vector_from_axis / vector_from_axis.norm() * radius + - ((middle-point_on_axis) * direction) * direction + - point_on_axis); -} - - - -template <> -Point<3> -CylinderBoundary<3>:: -get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const -{ - const Point<3> middle = StraightBoundary<3,3>::get_new_point_on_quad (quad); - - // same algorithm as above - const unsigned int spacedim = 3; - - const Tensor<1,spacedim> 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 Point<3>(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 Tensor<1,spacedim> 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 Point<3>(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 Tensor<1,spacedim> 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] = Point(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 -void -CylinderBoundary:: -get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const -{ - for (unsigned int v=0; v::vertices_per_face; ++v) - { - const Point vertex = face->vertex(v); - - const Tensor<1,spacedim> vector_from_axis = (vertex-point_on_axis) - - ((vertex-point_on_axis) * direction) * direction; - - face_vertex_normals[v] = (vector_from_axis / vector_from_axis.norm()); - } -} - - - -template -double -CylinderBoundary::get_radius () const -{ - return radius; -} - - -//======================================================================// - -template -ConeBoundary::ConeBoundary (const double radius_0, - const double radius_1, - const Point x_0, - const Point x_1) - : - radius_0 (radius_0), - radius_1 (radius_1), - x_0 (x_0), - x_1 (x_1) -{} - - - -template -std::unique_ptr > ConeBoundary::clone() const -{ - return std_cxx14::make_unique >(radius_0,radius_1,x_0,x_1); -} - - - -template -double ConeBoundary::get_radius (Point x) const -{ - for (unsigned int i = 0; i < dim; ++i) - if ((x_1 (i) - x_0 (i)) != 0) - return (radius_1 - radius_0) * (x (i) - x_0 (i)) / (x_1 (i) - x_0 (i)) + radius_0; - - return 0; -} - - - -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 Tensor<1,dim> 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*axis); - 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 Tensor<1,dim> 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*axis); - const Point middle_p = x_0 + c * axis; - // Compute the projection of the middle point on the boundary of the cone. - return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm (); -} - - - -template <> -Point<3> -ConeBoundary<3>:: -get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const -{ - const int dim = 3; - - const Tensor<1,dim> axis = x_1 - x_0; - // Compute the middle point of the quad. - const Point middle = StraightBoundary<3,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*axis); - 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 -void -ConeBoundary:: -get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const -{ - const Tensor<1,dim> axis = x_1 - x_0; - - for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_face; ++vertex) - { - // Compute the orthogonal projection of the vertex onto the axis of the - // cone. - const double c = (face->vertex (vertex) - x_0) * axis / (axis*axis); - const Point vertex_p = x_0 + c * axis; - // Then compute the vector pointing from the point vertex_p on - // the axis to the vertex. - const Tensor<1,dim> axis_to_vertex = face->vertex (vertex) - vertex_p; - - face_vertex_normals[vertex] = axis_to_vertex / axis_to_vertex.norm (); - } -} - - -template -inline -Tensor<1,dim> -ConeBoundary:: -normal_vector (const typename Triangulation::face_iterator &, - const Point &p) const -{ -// TODO only for cone opening along z-axis - Assert (dim == 3, ExcNotImplemented()); - Assert (this->radius_0 == 0., ExcNotImplemented()); - Assert (this->x_0 == Point(), ExcNotImplemented()); - for (unsigned int d=0; dx_1[d] == 0., ExcNotImplemented()); - Assert (this->x_1[dim-1] > 0, ExcNotImplemented()); - - const double c_squared = (this->radius_1 / this->x_1[dim-1])*(this->radius_1 / this->x_1[dim-1]); - Tensor<1,dim> normal = p; - normal[0] *= -2.0/c_squared; - normal[1] *= -2.0/c_squared; - normal[dim-1] *= 2.0; - - return normal/normal.norm(); -} - - - -//======================================================================// - -template -HyperBallBoundary::HyperBallBoundary (const Point p, - const double radius) - : - center(p), - radius(radius), - compute_radius_automatically(false) -{} - - - -template -std::unique_ptr > -HyperBallBoundary::clone() const -{ - return std_cxx14::make_unique >(center, radius); -} - - - -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) - r = (line->vertex(0) - center).norm(); - else - r = radius; - - // project to boundary - middle *= r / std::sqrt(middle.square()); - middle += center; - return middle; -} - - - -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) - r = (quad->vertex(0) - center).norm(); - 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 Tensor<1,spacedim> v0=p0-center, - v1=p1-center; - const double length=(v1-v0).norm(); - - double eps=1e-12; - (void)eps; - double r=0; - if (compute_radius_automatically) - r = (p0 - center).norm(); - else - r = radius; - - Assert((std::fabs(v0*v0-r*r) pm=0.5*(v0+v1); - - const double h=pm.norm(); - - // 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(pm+d/length*(v1-v0)); - points[n-1-i] = Point(pm-d/length*(v1-v0)); - } - - if ((n+1)%2==0) - // if the number of parts is even insert the midpoint - points[(n-1)/2] = Point(pm); - - - // project the points from the straight line to the HyperBallBoundary - 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 -void -HyperBallBoundary:: -get_normals_at_vertices (const typename Triangulation::face_iterator &face, - typename Boundary::FaceVertexNormals &face_vertex_normals) const -{ - for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) - face_vertex_normals[vertex] = face->vertex(vertex)-center; -} - - - -template -Point -HyperBallBoundary::get_center () const -{ - return center; -} - - - -template -double -HyperBallBoundary::get_radius () const -{ - Assert(!compute_radius_automatically, ExcRadiusNotSet()); - return radius; -} - - -/* ---------------------------------------------------------------------- */ - - -template -HalfHyperBallBoundary::HalfHyperBallBoundary (const Point center, - const double radius) : - HyperBallBoundary (center, radius) -{} - - - -template -std::unique_ptr > HalfHyperBallBoundary::clone() const -{ - return std_cxx14::make_unique >(this->get_center(), - this->get_radius()); -} - - - -template -Point -HalfHyperBallBoundary:: -get_new_point_on_line (const typename Triangulation::line_iterator &line) const -{ - // check whether center of object is at x==x_center, since then it belongs - // to the plane part of the boundary. however, this is not the case if it is - // at the outer perimeter - const Point line_center = line->center(); - const Point vertices[2] = { line->vertex(0), line->vertex(1) }; - - if ((line_center(0) == this->center(0)) - && - ((std::fabs(vertices[0].distance(this->center)-this->radius) > - 1e-5*this->radius) - || - (std::fabs(vertices[1].distance(this->center)-this->radius) > - 1e-5*this->radius))) - return line_center; - else - return HyperBallBoundary::get_new_point_on_line (line); -} - - - -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 -HyperShellBoundary::HyperShellBoundary (const Point ¢er) - : - HyperBallBoundary(center, 0.) -{ - this->compute_radius_automatically=true; -} - - - -template -std::unique_ptr > -HyperShellBoundary::clone() const -{ - return std_cxx14::make_unique >(this->get_center()); -} - - - -/* ---------------------------------------------------------------------- */ - - - - -template -HalfHyperShellBoundary::HalfHyperShellBoundary (const Point ¢er, - const double inner_radius, - const double outer_radius) - : - HyperShellBoundary (center), - inner_radius (inner_radius), - outer_radius (outer_radius) -{ - if (dim > 2) - Assert ((inner_radius >= 0) && - (outer_radius > 0) && - (outer_radius > inner_radius), - ExcMessage ("Inner and outer radii must be specified explicitly in 3d.")); -} - - - -template -std::unique_ptr > -HalfHyperShellBoundary::clone() const -{ - return std_cxx14::make_unique >(this->get_center(), - inner_radius, - outer_radius); -} - - - -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 Tensor<1,dim> 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 -TorusBoundary::TorusBoundary (const double R_, - const double r_) - : - R(R_), - r(r_) -{ - Assert (false, ExcNotImplemented()); -} - - - -template <> -TorusBoundary<2,3>::TorusBoundary (const double R_, - const double r_) - : - R(R_), - r(r_) -{ - Assert (R>r, ExcMessage("Outer radius must be greater than inner radius.")); -} - -template -std::unique_ptr > -TorusBoundary::clone() const -{ - return std_cxx14::make_unique >(R,r); -} - - - -template -double -TorusBoundary::get_correct_angle(const double angle, - const double x, - const double y) const -{ - if (y>=0) - { - if (x >=0) - return angle; - - return numbers::PI-angle; - } - - if (x <=0) - return numbers::PI+angle; - - return 2.0*numbers::PI-angle; -} - - - -template <> -Point<3> -TorusBoundary<2,3>::get_real_coord (const Point<2> &surfP) const -{ - const double theta=surfP(0); - const double phi=surfP(1); - - return Point<3> ((R+r*std::cos(phi))*std::cos(theta), - r*std::sin(phi), - (R+r*std::cos(phi))*std::sin(theta)); -} - - - -template <> -Point<2> -TorusBoundary<2,3>::get_surf_coord(const Point<3> &p) const -{ - const double phi=std::asin(std::abs(p(1))/r); - const double Rr_2=p(0)*p(0)+p(2)*p(2); - - Point<2> surfP; - surfP(1)=get_correct_angle(phi,Rr_2-R*R,p(1));//phi - - if (std::abs(p(0))<1.E-5) - { - if (p(2)>=0) - surfP(0) = numbers::PI*0.5; - else - surfP(0) = -numbers::PI*0.5; - } - else - { - const double theta = std::atan(std::abs(p(2)/p(0))); - surfP(0)= get_correct_angle(theta,p(0),p(2)); - } - - return surfP; -} - - - -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> -TorusBoundary<2,3>::get_surf_norm_from_sp(const Point<2> &surfP) const -{ - - Point<3> n; - double theta=surfP[0]; - double phi=surfP[1]; - - double f=R+r*std::cos(phi); - - n[0]=r*std::cos(phi)*std::cos(theta)*f; - n[1]=r*std::sin(phi)*f; - n[2]=r*std::sin(theta)*std::cos(phi)*f; - - return n; -} - - - -//Normal field without unit length -template <> -Point<3> -TorusBoundary<2,3>::get_surf_norm(const Point<3> &p) const -{ - - Point<2> surfP=get_surf_coord(p); - return get_surf_norm_from_sp(surfP); - -} - - - -template <> -void -TorusBoundary<2,3>:: -get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator &line, - std::vector< Point< 3 > > &points) const -{ - //Almost the same implementation as StraightBoundary<2,3> - unsigned int npoints=points.size(); - if (npoints==0) return; - - Point<2> p[2]; - - for (unsigned int i=0; i<2; i++) - p[i]=get_surf_coord(line->vertex(i)); - - unsigned int offset[2]; - offset[0]=0; - offset[1]=0; - - //Take care for periodic conditions & negative angles, see - //get_new_point_on_line() above. Because we don't 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 don't 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" - -DEAL_II_NAMESPACE_CLOSE diff --git a/source/grid/tria_boundary_lib.inst.in b/source/grid/tria_boundary_lib.inst.in deleted file mode 100644 index 6933a50412..0000000000 --- a/source/grid/tria_boundary_lib.inst.in +++ /dev/null @@ -1,34 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 1999 - 2016 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - -for (deal_II_dimension : DIMENSIONS) -{ - template class CylinderBoundary; - template class ConeBoundary; - template class HyperBallBoundary; - template class HalfHyperBallBoundary; - template class HyperShellBoundary; - template class HalfHyperShellBoundary; - -#if deal_II_dimension != 3 - template class HyperBallBoundary; -#endif -#if deal_II_dimension == 3 - template class CylinderBoundary; -#endif -} - - diff --git a/tests/fe/fe_nedelec_sz_non_rect_face.cc b/tests/fe/fe_nedelec_sz_non_rect_face.cc index d7317b04e5..25fdc8afb6 100644 --- a/tests/fe/fe_nedelec_sz_non_rect_face.cc +++ b/tests/fe/fe_nedelec_sz_non_rect_face.cc @@ -57,7 +57,6 @@ #include #include #include -#include #include #include diff --git a/tests/grid/normal_vector_01.cc b/tests/grid/normal_vector_01.cc index 8d4a31994d..dca9911320 100644 --- a/tests/grid/normal_vector_01.cc +++ b/tests/grid/normal_vector_01.cc @@ -15,8 +15,8 @@ -// test that at the vertices, Boundary::normal_vector returns the same as -// Boundary::get_normals_at_vertices once the latter vectors are normalized +// test that at the vertices, Manifold::normal_vector returns the same as +// Manifold::get_normals_at_vertices once the latter vectors are normalized diff --git a/tests/matrix_free/faces_value_optimization.cc b/tests/matrix_free/faces_value_optimization.cc index 125e375919..a159142555 100644 --- a/tests/matrix_free/faces_value_optimization.cc +++ b/tests/matrix_free/faces_value_optimization.cc @@ -20,7 +20,6 @@ #include "../tests.h" #include #include -#include #include #include #include @@ -161,8 +160,6 @@ void test () { Triangulation tria; GridGenerator::hyper_ball (tria); - static const HyperBallBoundary boundary; - tria.set_boundary (0, boundary); tria.begin_active ()->set_refine_flag(); tria.execute_coarsening_and_refinement(); typename Triangulation::active_cell_iterator cell, endc; diff --git a/tests/matrix_free/interpolate_functions_common.h b/tests/matrix_free/interpolate_functions_common.h index 13dff441e6..b9a90ab0e7 100644 --- a/tests/matrix_free/interpolate_functions_common.h +++ b/tests/matrix_free/interpolate_functions_common.h @@ -28,7 +28,6 @@ #include #include #include -#include #include #include #include diff --git a/tests/matrix_free/matrix_vector_faces_04.cc b/tests/matrix_free/matrix_vector_faces_04.cc index e5515b4dc7..02baa02c7c 100644 --- a/tests/matrix_free/matrix_vector_faces_04.cc +++ b/tests/matrix_free/matrix_vector_faces_04.cc @@ -32,8 +32,6 @@ void test () { Triangulation tria; GridGenerator::hyper_ball (tria); - static const HyperBallBoundary boundary; - tria.set_boundary (0, boundary); tria.refine_global(4-dim); tria.begin_active()->set_refine_flag(); tria.last()->set_refine_flag(); diff --git a/tests/matrix_free/matrix_vector_faces_common.h b/tests/matrix_free/matrix_vector_faces_common.h index b5ce5d95b3..068dbc2cd3 100644 --- a/tests/matrix_free/matrix_vector_faces_common.h +++ b/tests/matrix_free/matrix_vector_faces_common.h @@ -26,7 +26,6 @@ #include #include #include -#include #include #include #include diff --git a/tests/matrix_free/multigrid_dg_periodic.cc b/tests/matrix_free/multigrid_dg_periodic.cc index 92adece7d9..c87826aa87 100644 --- a/tests/matrix_free/multigrid_dg_periodic.cc +++ b/tests/matrix_free/multigrid_dg_periodic.cc @@ -27,7 +27,6 @@ #include #include #include -#include #include #include #include diff --git a/tests/matrix_free/multigrid_dg_sip_01.cc b/tests/matrix_free/multigrid_dg_sip_01.cc index 5fe468aec5..0adc32168f 100644 --- a/tests/matrix_free/multigrid_dg_sip_01.cc +++ b/tests/matrix_free/multigrid_dg_sip_01.cc @@ -28,7 +28,6 @@ #include #include #include -#include #include #include #include diff --git a/tests/matrix_free/multigrid_dg_sip_02.cc b/tests/matrix_free/multigrid_dg_sip_02.cc index d4cd524236..227e2815aa 100644 --- a/tests/matrix_free/multigrid_dg_sip_02.cc +++ b/tests/matrix_free/multigrid_dg_sip_02.cc @@ -29,7 +29,6 @@ #include #include #include -#include #include #include #include diff --git a/tests/mpi/step-37.cc b/tests/mpi/step-37.cc index 1d53080716..f16aec4c77 100644 --- a/tests/mpi/step-37.cc +++ b/tests/mpi/step-37.cc @@ -35,7 +35,6 @@ #include #include #include -#include #include #include