]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove all remaining usages of the Boundary class.
authorDavid Wells <wellsd2@rpi.edu>
Sat, 12 May 2018 14:44:10 +0000 (10:44 -0400)
committerDavid Wells <wellsd2@rpi.edu>
Sun, 13 May 2018 18:21:59 +0000 (14:21 -0400)
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 <int dim>
class Boundary {
  public:
     /**
      *  Typedef an array of the needed number
      *  of old points.
      */
    typedef const Point<dim>* PointArray[1<<(dim-1)];

     /**
      *  This function calculates the position
      *  of the new vertex.
      */
    virtual Point<dim> in_between (const PointArray &neighbors) const = 0;
};

/**
    Specialisation of \Ref{Boundary}<dim>, 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 <int dim>
class StraightBoundary : public Boundary<dim> {
  public:
     /**
      *  This function calculates the position
      *  of the new vertex.
      */
    virtual Point<dim> in_between (const PointArray &neighbors) const {
      Point<dim> 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.

36 files changed:
contrib/python-bindings/include/triangulation_wrapper.h
contrib/python-bindings/source/export_triangulation.cc
doc/doxygen/headers/boundary.h [deleted file]
examples/step-32/doc/intro.dox
examples/step-59/step-59.cc
include/deal.II/grid/tria.h
include/deal.II/grid/tria_accessor.h
include/deal.II/grid/tria_accessor.templates.h
include/deal.II/grid/tria_boundary.h [deleted file]
include/deal.II/grid/tria_boundary_lib.h [deleted file]
source/dofs/dof_accessor.cc
source/dofs/dof_accessor_get.cc
source/dofs/dof_accessor_set.cc
source/fe/fe.cc
source/fe/fe_values.cc
source/fe/mapping_c1.cc
source/fe/mapping_fe_field.cc
source/fe/mapping_manifold.cc
source/fe/mapping_q_generic.cc
source/grid/CMakeLists.txt
source/grid/tria.cc
source/grid/tria_accessor.cc
source/grid/tria_boundary.cc [deleted file]
source/grid/tria_boundary.inst.in [deleted file]
source/grid/tria_boundary_lib.cc [deleted file]
source/grid/tria_boundary_lib.inst.in [deleted file]
tests/fe/fe_nedelec_sz_non_rect_face.cc
tests/grid/normal_vector_01.cc
tests/matrix_free/faces_value_optimization.cc
tests/matrix_free/interpolate_functions_common.h
tests/matrix_free/matrix_vector_faces_04.cc
tests/matrix_free/matrix_vector_faces_common.h
tests/matrix_free/multigrid_dg_periodic.cc
tests/matrix_free/multigrid_dg_sip_01.cc
tests/matrix_free/multigrid_dg_sip_02.cc
tests/mpi/step-37.cc

index a36f8c62ccb90ea64828dca07eb436561b278974..440ad6d799af36e936154eab9d04e27aa10639c2 100644 (file)
@@ -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 &center,
                                      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 &center,
                                   const double  radius = 1.);
index e4d9402df3e951109deda27e55ce95d30ae1076f..83c94df619dffdd203186c9ead53bd51075b8115 100644 (file)
@@ -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 (file)
index 3945fd1..0000000
+++ /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:
- * <ul>
- *
- *   <li> 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.
- *
- *   <li> 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).
- *
- *   <li> 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.
- *
- * </ul>
- *
- * In the context of triangulations, each object stores a number
- * called <tt>manifold_id</tt>, and each face of a cell that is
- * located at the boundary of the domain stores a number called
- * <tt>boundary_id</tt> 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
- */
index 7b125c0b1cde2c412046a558533ac7c9014a7bee..55e60a1c2e6240eee39b1bd432d64b9657ecb23b 100644 (file)
@@ -885,8 +885,7 @@ the following quantities:
   <a target="_top" href="http://en.wikipedia.org/wiki/Outer_core">outer earth
   core</a> 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.
 
   <li>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
index 248cc96ec9aa6912a767d7bc8981947c2c5434ab..98f6170ca6f8d18e1ba61d5aefa46f5008c355d1 100644 (file)
@@ -39,7 +39,6 @@
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_tools.h>
 
index 7ea247c60562301ec730a783a1a9c4c8a37019c3..c4e1e60769f883116d7ee6591d17a85f2b57ab40 100644 (file)
@@ -42,8 +42,6 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-template <int dim, int spacedim> class Boundary;
-template <int dim, int spacedim> class StraightBoundary;
 template <int dim, int spacedim> class Manifold;
 
 namespace GridTools
@@ -1277,16 +1275,6 @@ private:
   typedef dealii::internal::TriangulationImplementation::Iterators<dim, spacedim> 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<dim,spacedim> 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 <tt>set_boundary(number)</tt>, 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<dim,spacedim> &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<dim,spacedim> &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
index 89e5199380d9aae5d5ae02c5513f87cbdddc9998..112d8ca654cfb084b068a20f141a42ec770ef8e0 100644 (file)
@@ -36,7 +36,6 @@ template <typename Accessor> class TriaRawIterator;
 template <typename Accessor> class TriaIterator;
 template <typename Accessor> class TriaActiveIterator;
 
-template <int dim, int spacedim> class Boundary;
 template <int dim, int spacedim> 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<dim,spacedim> &get_boundary () const;
-
   /**
    * Return a constant reference to the manifold object used for this object.
    *
index eba6866591eb3afd6d67b459f1a079862714ffd0..6abfff53120c5e1d3edbb54437737371fb126f0e 100644 (file)
@@ -1998,29 +1998,6 @@ TriaAccessor<structdim, dim, spacedim>::at_boundary () const
 
 
 
-template <int structdim, int dim, int spacedim>
-const Boundary<dim,spacedim> &
-TriaAccessor<structdim, dim, spacedim>::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 <int structdim, int dim, int spacedim>
 const Manifold<dim,spacedim> &
 TriaAccessor<structdim, dim, spacedim>::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 (file)
index f3c3d3e..0000000
+++ /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 <deal.II/base/config.h>
-#include <deal.II/base/subscriptor.h>
-#include <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/thread_management.h>
-#include <deal.II/base/point.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/manifold.h>
-
-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 <int dim, int spacedim=dim>
-class DEAL_II_DEPRECATED Boundary : public FlatManifold<dim, spacedim>
-{
-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<Manifold<dim,spacedim> > 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<dim,spacedim>::line_iterator &line,
-                                   std::vector<Point<spacedim> > &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. <tt>n=points.size()=m*m</tt>. It is the task of the derived classes
-   * to arrange the points such they split the quad into <tt>(m+1)(m+1)</tt>
-   * 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<dim,spacedim>::quad_iterator &quad,
-                                   std::vector<Point<spacedim> > &points) const;
-
-  /**
-   * Depending on <tt>dim=2</tt> or <tt>dim=3</tt> this function calls the
-   * get_intermediate_points_on_line or the get_intermediate_points_on_quad
-   * function. It throws an exception for <tt>dim=1</tt>. This wrapper allows
-   * dimension independent programming.
-   */
-  void
-  get_intermediate_points_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face,
-                                   std::vector<Point<spacedim> > &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<spacedim>
-  project_to_surface (const typename Triangulation<dim,spacedim>::line_iterator &line,
-                      const Point<spacedim> &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<spacedim>
-  project_to_surface (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
-                      const Point<spacedim> &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<spacedim>
-  project_to_surface (const typename Triangulation<dim,spacedim>::hex_iterator &hex,
-                      const Point<spacedim> &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<Point<1> > &
-  get_line_support_points (const unsigned int n_intermediate_points) const;
-
-private:
-  /**
-   * Point generator for the intermediate points on a boundary.
-   */
-  mutable std::vector<std::shared_ptr<QGaussLobatto<1> > > points;
-
-  /**
-   * Mutex for protecting the points array.
-   */
-  mutable Threads::Mutex mutex;
-};
-
-
-
-/**
- * Specialization of Boundary<dim,spacedim>, 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 <int dim, int spacedim=dim>
-class DEAL_II_DEPRECATED StraightBoundary : public Boundary<dim,spacedim>
-{
-public:
-  /**
-   * Default constructor. Some compilers require this for some reasons.
-   */
-  StraightBoundary ();
-
-  /**
-   * Clone this Boundary.
-   */
-  virtual std::unique_ptr<Manifold<dim,spacedim> > 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<spacedim>
-  get_new_point_on_line (const typename Triangulation<dim,spacedim>::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<spacedim>
-  get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const override;
-
-  /**
-   * Return <tt>n=points.size()</tt> 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<dim,spacedim>::line_iterator &line,
-                                   std::vector<Point<spacedim> > &points) const override;
-
-  /**
-   * Return <tt>n=points.size()=m*m</tt> 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<dim,spacedim>::quad_iterator &quad,
-                                   std::vector<Point<spacedim> > &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<dim,spacedim>::face_iterator &face,
-                 const Point<spacedim> &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<dim,spacedim>::face_iterator &face,
-                           typename Boundary<dim,spacedim>::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<spacedim>
-  project_to_surface (const typename Triangulation<dim,spacedim>::line_iterator &line,
-                      const Point<spacedim> &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<spacedim>
-  project_to_surface (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
-                      const Point<spacedim> &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<spacedim>
-  project_to_surface (const typename Triangulation<dim,spacedim>::hex_iterator &hex,
-                      const Point<spacedim> &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<Point<1> > &) const;
-
-template <>
-void
-Boundary<1,2>::
-get_intermediate_points_on_face (const Triangulation<1,2>::face_iterator &,
-                                 std::vector<Point<2> > &) const;
-
-template <>
-void
-Boundary<1,3>::
-get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &,
-                                 std::vector<Point<3> > &) 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<Point<3> > &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 (file)
index ff25666..0000000
+++ /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 <deal.II/base/config.h>
-#include <deal.II/grid/tria_boundary.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-/**
- * Boundary object for the hull of a cylinder.  In three dimensions, points
- * are projected on a circular tube along the <tt>x-</tt>, <tt>y-</tt> or
- * <tt>z</tt>-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 <int dim, int spacedim = dim>
-class DEAL_II_DEPRECATED CylinderBoundary : public StraightBoundary<dim,spacedim>
-{
-public:
-  /**
-   * Constructor. Using default values for the constructor arguments yields a
-   * circular tube along the x-axis (<tt>axis=0</tt>). Choose <tt>axis=1</tt>
-   * or <tt>axis=2</tt> 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<spacedim> &direction,
-                    const Point<spacedim> &point_on_axis);
-
-  /**
-   * Clone this Boundary object.
-   */
-  virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
-
-  /**
-   * Refer to the general documentation of this class and the documentation of
-   * the base class.
-   */
-  virtual Point<spacedim>
-  get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const override;
-
-  /**
-   * Refer to the general documentation of this class and the documentation of
-   * the base class.
-   */
-  virtual Point<spacedim>
-  get_new_point_on_quad (const typename Triangulation<dim,spacedim>::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<dim,spacedim>::line_iterator &line,
-                                   std::vector<Point<spacedim> > &points) const override;
-
-  /**
-   * Refer to the general documentation of this class and the documentation of
-   * the base class.
-   *
-   * Only implemented for <tt>dim=3</tt> and for <tt>points.size()==1</tt>.
-   */
-  virtual void
-  get_intermediate_points_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
-                                   std::vector<Point<spacedim> > &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<dim,spacedim>::face_iterator &face,
-                           typename Boundary<dim,spacedim>::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<spacedim> direction;
-
-  /**
-   * An arbitrary point on the axis.
-   */
-  const Point<spacedim> 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<spacedim> &p0, const Point<spacedim> &p1,
-                                               std::vector<Point<spacedim> > &points) const;
-
-  /**
-   * Given a number for the axis, return a vector that denotes this direction.
-   */
-  static Point<spacedim> 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<dim> triangulation;
- *  GridGenerator::truncated_cone (triangulation);
- *  Point<dim> p1, p2;
- *  p1[0] = -1;
- *  p2[0] = 1;
- *  const ConeBoundary<dim> 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&uuml;rg, 2009
- */
-template <int dim>
-class DEAL_II_DEPRECATED ConeBoundary : public StraightBoundary<dim>
-{
-public:
-  /**
-   * Constructor. Here the boundary object is constructed. The points
-   * <tt>x_0</tt> and <tt>x_1</tt> describe the starting and ending points of
-   * the axis of the (truncated) cone. <tt>radius_0</tt> denotes the radius
-   * corresponding to <tt>x_0</tt> and <tt>radius_1</tt> the one corresponding
-   * to <tt>x_1</tt>.
-   */
-  ConeBoundary (const double radius_0,
-                const double radius_1,
-                const Point<dim> x_0,
-                const Point<dim> x_1);
-
-  /**
-   * Clone this Boundary object.
-   */
-  virtual std::unique_ptr<Manifold<dim,dim> > clone() const override;
-
-  /**
-   * Return the radius of the (truncated) cone at given point <tt>x</tt> on
-   * the axis.
-   */
-  double get_radius (const Point<dim> x) const;
-
-  /**
-   * Refer to the general documentation of this class and the documentation of
-   * the base class.
-   */
-  virtual
-  Point<dim>
-  get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const override;
-
-  /**
-   * Refer to the general documentation of this class and the documentation of
-   * the base class.
-   */
-  virtual
-  Point<dim>
-  get_new_point_on_quad (const typename Triangulation<dim>::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<dim>::line_iterator &line,
-                                   std::vector<Point<dim> > &points) const override;
-
-  /**
-   * Refer to the general documentation of this class and the documentation of
-   * the base class.
-   *
-   * Only implemented for <tt>dim=3</tt> and for <tt>points.size()==1</tt>.
-   */
-  virtual
-  void
-  get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
-                                   std::vector<Point<dim> > &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<dim>::face_iterator &face,
-                           typename Boundary<dim>::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<dim>::face_iterator &face,
-                 const Point<dim> &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<dim> x_0;
-
-  /**
-   * Ending point of the axis.
-   */
-  const Point<dim> 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<dim> &p0,
-                                          const Point<dim> &p1,
-                                          std::vector<Point<dim> > &points) const;
-};
-
-
-
-/**
- * Specialization of Boundary<dim>, 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 <int dim, int spacedim=dim>
-class DEAL_II_DEPRECATED HyperBallBoundary : public StraightBoundary<dim,spacedim>
-{
-public:
-  /**
-   * Constructor
-   */
-  HyperBallBoundary (const Point<spacedim> p      = Point<spacedim>(),
-                     const double     radius = 1.0);
-
-  /**
-   * Clone this Boundary object.
-   */
-  virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
-
-  /**
-   * Refer to the general documentation of this class and the documentation of
-   * the base class.
-   */
-  virtual
-  Point<spacedim>
-  get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const override;
-
-  /**
-   * Refer to the general documentation of this class and the documentation of
-   * the base class.
-   */
-  virtual
-  Point<spacedim>
-  get_new_point_on_quad (const typename Triangulation<dim,spacedim>::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<dim,spacedim>::line_iterator &line,
-                                   std::vector<Point<spacedim> > &points) const override;
-
-  /**
-   * Refer to the general documentation of this class and the documentation of
-   * the base class.
-   *
-   * Only implemented for <tt>dim=3</tt> and for <tt>points.size()==1</tt>.
-   */
-  virtual
-  void
-  get_intermediate_points_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
-                                   std::vector<Point<spacedim> > &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<dim,spacedim>::face_iterator &face,
-                 const Point<spacedim> &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<dim,spacedim>::face_iterator &face,
-                           typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const override;
-
-  /**
-   * Return the center of the ball.
-   */
-  Point<spacedim>
-  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<spacedim> 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<spacedim> &p0, const Point<spacedim> &p1,
-                                               std::vector<Point<spacedim> > &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 <int dim>
-class DEAL_II_DEPRECATED HalfHyperBallBoundary : public HyperBallBoundary<dim>
-{
-public:
-  /**
-   * Constructor
-   */
-  HalfHyperBallBoundary (const Point<dim> p      = Point<dim>(),
-                         const double     radius = 1.0);
-
-  /**
-   * Clone this Boundary object.
-   */
-  virtual std::unique_ptr<Manifold<dim,dim> > clone() const override;
-
-  /**
-   * Check if on the line <tt>x==0</tt>, otherwise pass to the base class.
-   */
-  virtual Point<dim>
-  get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const override;
-
-  /**
-   * Check if on the line <tt>x==0</tt>, otherwise pass to the base class.
-   */
-  virtual Point<dim>
-  get_new_point_on_quad (const typename Triangulation<dim>::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<dim>::line_iterator &line,
-                                   std::vector<Point<dim> > &points) const override;
-
-  /**
-   * Refer to the general documentation of this class and the documentation of
-   * the base class.
-   *
-   * Only implemented for <tt>dim=3</tt> and for <tt>points.size()==1</tt>.
-   */
-  virtual void
-  get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
-                                   std::vector<Point<dim> > &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<dim>::face_iterator &face,
-                           typename Boundary<dim>::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 <int dim>
-class DEAL_II_DEPRECATED HyperShellBoundary : public HyperBallBoundary<dim>
-{
-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<dim> &center = Point<dim>());
-
-  /**
-   * Clone this Boundary object.
-   */
-  virtual std::unique_ptr<Manifold<dim,dim> > 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 <int dim>
-class DEAL_II_DEPRECATED HalfHyperShellBoundary : public HyperShellBoundary<dim>
-{
-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<dim> &center = Point<dim>(),
-                          const double inner_radius = -1,
-                          const double outer_radius = -1);
-
-  /**
-   * Clone this Boundary object.
-   */
-  virtual std::unique_ptr<Manifold<dim,dim> > clone() const override;
-
-  /**
-   * Construct a new point on a line.
-   */
-  virtual Point<dim>
-  get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const override;
-
-  /**
-   * Construct a new point on a quad.
-   */
-  virtual Point<dim>
-  get_new_point_on_quad (const typename Triangulation<dim>::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<dim>::line_iterator &line,
-                                   std::vector<Point<dim> > &points) const override;
-
-  /**
-   * Refer to the general documentation of this class and the documentation of
-   * the base class.
-   *
-   * Only implemented for <tt>dim=3</tt> and for <tt>points.size()==1</tt>.
-   */
-  virtual void
-  get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
-                                   std::vector<Point<dim> > &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<dim>::face_iterator &face,
-                           typename Boundary<dim>::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 <tt>dim=2</tt>,<tt>spacedim=3</tt>, that
- * is, just the surface.
- *
- * @deprecated This class has been deprecated in favor of the more general
- * TorusManifold class.
- */
-template <int dim, int spacedim>
-class DEAL_II_DEPRECATED TorusBoundary : public Boundary<dim,spacedim>
-{
-public:
-  /**
-   * Constructor.<tt>R</tt> has to be greater than <tt>r</tt>.
-   */
-  TorusBoundary (const double R, const double r);
-
-  /**
-   * Clone this Boundary object.
-   */
-  virtual std::unique_ptr<Manifold<dim,spacedim> > clone() const override;
-
-//Boundary Refinement Functions
-  /**
-   * Construct a new point on a line.
-   */
-  virtual Point<spacedim>
-  get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const override;
-
-  /**
-   * Construct a new point on a quad.
-   */
-  virtual Point<spacedim>
-  get_new_point_on_quad (const typename Triangulation<dim,spacedim>::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<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const override;
-
-private:
-  //Handy functions
-  /**
-   * Function that corrects the value and sign of angle, that is, given
-   * <tt>angle=tan(abs(y/x))</tt>; if <tt> (y > 0) && (x < 0) </tt> then
-   * <tt>correct_angle = Pi - angle</tt>, 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
-   * <tt>(theta,phi)</tt> to <tt>(x,y,z)</tt>.
-   */
-  Point<spacedim>  get_real_coord(const Point<dim> &surfP) const;
-
-  /**
-   * Get the surface coordinates of the Torus, i.e., from <tt>(x,y,z)</tt> to
-   * <tt>(theta,phi)</tt>.
-   */
-  Point<dim>       get_surf_coord(const Point<spacedim> &p) const;
-
-  /**
-   * Get the normal from surface coordinates. This normal does not have unit
-   * length.
-   */
-  Point<spacedim>  get_surf_norm_from_sp(const Point<dim> &surfP)      const;
-
-  /**
-   * Get the normal from cartesian coordinates. This normal does not have unit
-   * length.
-   */
-  Point<spacedim>  get_surf_norm(const Point<spacedim> &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<Point<1> > &) const;
-template <>
-void
-HyperBallBoundary<3>::get_intermediate_points_on_quad (
-  const Triangulation<3>::quad_iterator &quad,
-  std::vector<Point<3> > &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<Point<1> > &) 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<Point<1> > &) 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
index 42e8c7e85e0dd42f9a5875b50ac2522235c261ea..ec61cb7ac92171251dfe82e146bf40043bc7abf4 100644 (file)
@@ -17,7 +17,6 @@
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_levels.h>
 #include <deal.II/hp/dof_handler.h>
-#include <deal.II/grid/tria_boundary.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/grid/tria_iterator.templates.h>
 #include <deal.II/fe/fe.h>
index d823552b1947ee001c258297c6c6c789f73547f4..83087b67eb9363a6b5229321ef3d7c87f6fd4860 100644 (file)
@@ -28,7 +28,6 @@
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_levels.h>
 #include <deal.II/hp/dof_handler.h>
-#include <deal.II/grid/tria_boundary.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/grid/tria_iterator.templates.h>
 #include <deal.II/fe/fe.h>
index e4eb96f89643e98be4eab7750b1d46987cb785ca..3ab878e32c23645c42791398d22d69c11b8a3925 100644 (file)
@@ -28,7 +28,6 @@
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_levels.h>
 #include <deal.II/hp/dof_handler.h>
-#include <deal.II/grid/tria_boundary.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/grid/tria_iterator.templates.h>
 #include <deal.II/fe/fe.h>
index 6110359bc786d15c51696935795cb212331fab09..138e40671c6c9aac297db1e3a02ac564cf7fb8dd 100644 (file)
@@ -22,7 +22,6 @@
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/grid/tria_boundary.h>
 
 #include <algorithm>
 #include <functional>
index 2b2d0bfcf317506a8bf36e18c3f8f86d7d5d0ee4..1218751bdf6b331a0a2ecb72eca7082e691e4e70 100644 (file)
@@ -33,7 +33,6 @@
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/grid/tria_accessor.h>
-#include <deal.II/grid/tria_boundary.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/fe/mapping_q1.h>
 #include <deal.II/fe/fe_values.h>
index 73d089152f1be04cd7bf427e0d70981d0fce029a..4ae4bba504fd8d3765492dfd07759ce3541ebeff 100644 (file)
@@ -15,9 +15,9 @@
 
 #include <deal.II/base/std_cxx14/memory.h>
 
+#include <deal.II/grid/manifold.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary.h>
 
 #include <deal.II/fe/mapping_c1.h>
 
index 3009f8cc8f26790c230ccf5697ffcd8afb444cc0..3f5a57172def2c9ba053fb59e19645f8b20c22e9 100644 (file)
@@ -34,7 +34,6 @@
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/fe/fe_tools.h>
 #include <deal.II/fe/fe_values.h>
index 01fdcffeae1f437e7a00c918ae7a983639fb2552..1c6edd9c20abdf02115203f9984c24cd44d554b6 100644 (file)
@@ -25,7 +25,6 @@
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/grid/manifold.h>
 #include <deal.II/fe/fe_tools.h>
index bff949b15b5f6d38701b083a4f779552d0d827fb..c74c284523b32e8999adfcd98e72025a9cc71a10 100644 (file)
@@ -27,7 +27,6 @@
 #include <deal.II/lac/tensor_product_matrix.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary.h>
 #include <deal.II/grid/manifold_lib.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/fe/fe_tools.h>
@@ -3612,11 +3611,9 @@ add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterat
     }
   else
     // otherwise call the more complicated functions and ask for inner points
-    // from the boundary description
+    // from the manifold description
     {
       std::vector<Point<spacedim> > 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<GeometryInfo<dim>::lines_per_cell; ++line_no)
         {
           const typename Triangulation<dim,spacedim>::line_iterator
@@ -3633,34 +3630,20 @@ add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterat
               cell->get_manifold() :
               line->get_manifold() );
 
-          if (const Boundary<dim,spacedim> *boundary
-              = dynamic_cast<const Boundary<dim,spacedim> *>(&manifold))
+          const std::array<Point<spacedim>, 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<dim>::line_to_cell_vertices(line_no, 0)),
+              cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 1))
             }
-          else
-            {
-              const std::array<Point<spacedim>, 2> vertices
-              {
-                {
-                  cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0)),
-                  cell->vertex(GeometryInfo<dim>::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<const Boundary<3,3> *>(&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; i<tmp_points.size(); ++i)
-            a.push_back(tmp_points[fe_q->adjust_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<Point<3>, 200>
-          tmp_points(GeometryInfo<2>::vertices_per_cell
-                     + GeometryInfo<2>::lines_per_cell*(polynomial_degree-1));
-          for (unsigned int v=0; v<GeometryInfo<2>::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<GeometryInfo<2>::lines_per_cell; ++line)
-              for (unsigned int i=0; i<polynomial_degree-1; ++i)
-                tmp_points[4+line*(polynomial_degree-1)+i] =
-                  a[GeometryInfo<3>::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<Point<3>, 200>
+      tmp_points(GeometryInfo<2>::vertices_per_cell
+                 + GeometryInfo<2>::lines_per_cell*(polynomial_degree-1));
+      for (unsigned int v=0; v<GeometryInfo<2>::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<GeometryInfo<2>::lines_per_cell; ++line)
+          for (unsigned int i=0; i<polynomial_degree-1; ++i)
+            tmp_points[4+line*(polynomial_degree-1)+i] =
+              a[GeometryInfo<3>::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<Point<3> >                  &a) const
 {
-  if (const Boundary<2,3> *boundary =
-        dynamic_cast<const Boundary<2,3> *>(&cell->get_manifold()))
-    {
-      std::vector<Point<3> > 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<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
-      for (unsigned int i=0; i<GeometryInfo<2>::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<polynomial_degree-1; ++q2)
-        for (unsigned int q1=0; q1<polynomial_degree-1; ++q1, ++q)
-          {
-            Point<2> point(line_support_points.point(q1+1)[0],
-                           line_support_points.point(q2+1)[0]);
-            for (unsigned int i=0; i<GeometryInfo<2>::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<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
+  for (unsigned int i=0; i<GeometryInfo<2>::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<polynomial_degree-1; ++q2)
+    for (unsigned int q1=0; q1<polynomial_degree-1; ++q1, ++q)
+      {
+        Point<2> point(line_support_points.point(q1+1)[0],
+                       line_support_points.point(q2+1)[0]);
+        for (unsigned int i=0; i<GeometryInfo<2>::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);
 }
 
 
index c73798600ee4d20ab505e1e6404d25779d7d300a..6de7024fd7b3c60f7b3540e37c7b07bcf60aafdb 100644 (file)
@@ -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
   )
index 693d00d854ecaa061089123e5724b4235b050840..44029d51a6551ce19d109ead3dc910064ae4900f 100644 (file)
@@ -23,7 +23,6 @@
 #include <deal.II/grid/tria_levels.h>
 #include <deal.II/grid/tria_faces.h>
 #include <deal.II/grid/manifold.h>
-#include <deal.II/grid/tria_boundary.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/grid/grid_tools.h>
@@ -8791,11 +8790,6 @@ namespace internal
 }
 
 
-template <int dim, int spacedim>
-const StraightBoundary<dim,spacedim>
-Triangulation<dim, spacedim>::straight_boundary = StraightBoundary<dim,spacedim>();
-
-
 
 template <int dim, int spacedim>
 const unsigned int
@@ -8942,14 +8936,6 @@ Triangulation<dim,spacedim>::get_mesh_smoothing() const
 
 
 
-template <int dim, int spacedim>
-void
-Triangulation<dim, spacedim>::set_boundary (const types::manifold_id m_number,
-                                            const Boundary<dim, spacedim> &boundary_object)
-{
-  set_manifold(m_number, boundary_object);
-}
-
 template <int dim, int spacedim>
 void
 Triangulation<dim, spacedim>::set_manifold (const types::manifold_id m_number,
@@ -8962,13 +8948,6 @@ Triangulation<dim, spacedim>::set_manifold (const types::manifold_id m_number,
 }
 
 
-template <int dim, int spacedim>
-void
-Triangulation<dim, spacedim>::set_boundary (const types::manifold_id m_number)
-{
-  set_manifold(m_number);
-}
-
 
 template <int dim, int spacedim>
 void
@@ -9065,18 +9044,6 @@ Triangulation<dim, spacedim>::set_all_manifold_ids_on_boundary (const types::bou
 }
 
 
-template <int dim, int spacedim>
-const Boundary<dim,spacedim> &
-Triangulation<dim, spacedim>::get_boundary (const types::manifold_id m_number) const
-{
-  const Boundary<dim, spacedim> *man =
-    dynamic_cast<const Boundary<dim, spacedim> *>(&get_manifold(m_number));
-  Assert(man != nullptr,
-         ExcMessage("You tried to get a Boundary, but I only have a Manifold."));
-
-  return *man;
-}
-
 
 template <int dim, int spacedim>
 const Manifold<dim,spacedim> &
index 7fc45b7ffb242c233b638f61c9a59cd9191d74f3..f6339a562660d43d63a900a81f456c56d64382c7 100644 (file)
@@ -19,7 +19,6 @@
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_accessor.templates.h>
 #include <deal.II/grid/tria_iterator.templates.h>
-#include <deal.II/grid/tria_boundary.h>
 #include <deal.II/base/geometry_info.h>
 #include <deal.II/base/quadrature.h>
 #include <deal.II/grid/grid_tools.h>
@@ -1047,13 +1046,7 @@ namespace
   Point<spacedim> get_new_point_on_object(const TriaAccessor<structdim, dim, spacedim> &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<const Boundary<dim,spacedim> *>(&obj.get_manifold()) != nullptr)
-      return get_new_point_on_object(obj);
-    else
+    if (use_interpolation)
       {
         TriaRawIterator<TriaAccessor<structdim, dim, spacedim> > 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 (file)
index 6cbab1f..0000000
+++ /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 <deal.II/base/tensor.h>
-#include <deal.II/base/std_cxx14/memory.h>
-#include <deal.II/grid/tria_boundary.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/fe/fe_q.h>
-#include <cmath>
-
-DEAL_II_NAMESPACE_OPEN
-
-
-
-/* -------------------------- Boundary --------------------- */
-
-
-template <int dim, int spacedim>
-std::unique_ptr<Manifold<dim, spacedim> >
-Boundary<dim,spacedim>::clone() const
-{
-  AssertThrow(false, ExcPureFunctionCalled());
-  return std::unique_ptr<Manifold<dim,spacedim> >();
-}
-
-
-
-template <int dim, int spacedim>
-void
-Boundary<dim, spacedim>::
-get_intermediate_points_on_line (const typename Triangulation<dim, spacedim>::line_iterator &,
-                                 std::vector<Point<spacedim> > &) const
-{
-  Assert (false, ExcPureFunctionCalled());
-}
-
-
-
-template <int dim, int spacedim>
-void
-Boundary<dim, spacedim>::
-get_intermediate_points_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &,
-                                 std::vector<Point<spacedim> > &) const
-{
-  Assert (false, ExcPureFunctionCalled());
-}
-
-
-template <int dim, int spacedim>
-void
-Boundary<dim,spacedim>::
-get_intermediate_points_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face,
-                                 std::vector<Point<spacedim> > &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<Point<1> > &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
-
-
-template <>
-void
-Boundary<1,2>::
-get_intermediate_points_on_face (const Triangulation<1,2>::face_iterator &,
-                                 std::vector<Point<2> > &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
-
-
-template <>
-void
-Boundary<1,3>::
-get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &,
-                                 std::vector<Point<3> > &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-Boundary<dim, spacedim>::
-project_to_surface (const typename Triangulation<dim, spacedim>::line_iterator &,
-                    const Point<spacedim>                                &trial_point) const
-{
-  if (spacedim <= 1)
-    return trial_point;
-  else
-    {
-      Assert (false, ExcPureFunctionCalled());
-      return Point<spacedim>();
-    }
-}
-
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-Boundary<dim, spacedim>::
-project_to_surface (const typename Triangulation<dim, spacedim>::quad_iterator &,
-                    const Point<spacedim>                                &trial_point) const
-{
-  if (spacedim <= 2)
-    return trial_point;
-  else
-    {
-      Assert (false, ExcPureFunctionCalled());
-      return Point<spacedim>();
-    }
-}
-
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-Boundary<dim, spacedim>::
-project_to_surface (const typename Triangulation<dim, spacedim>::hex_iterator &,
-                    const Point<spacedim>                                &trial_point) const
-{
-  if (spacedim <= 3)
-    return trial_point;
-  else
-    {
-      Assert (false, ExcPureFunctionCalled());
-      return Point<spacedim>();
-    }
-}
-
-
-
-template <int dim, int spacedim>
-const std::vector<Point<1> > &
-Boundary<dim,spacedim>::
-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<QGaussLobatto<1> >
-                                        (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 <int dim, int spacedim>
-StraightBoundary<dim, spacedim>::StraightBoundary () = default;
-
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-StraightBoundary<dim, spacedim>::
-get_new_point_on_line (const typename Triangulation<dim, spacedim>::line_iterator &line) const
-{
-  return (line->vertex(0) + line->vertex(1)) / 2;
-}
-
-
-template <int dim, int spacedim>
-std::unique_ptr<Manifold<dim, spacedim> >
-StraightBoundary<dim,spacedim>::clone() const
-{
-  return std_cxx14::make_unique<StraightBoundary<dim,spacedim> >();
-}
-
-
-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 <int dim>
-  Point<3>
-  compute_new_point_on_quad (const typename Triangulation<dim, 3>::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 <int dim, int spacedim>
-Point<spacedim>
-StraightBoundary<dim, spacedim>::
-get_new_point_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
-{
-  return FlatManifold<dim,spacedim>::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 <int dim, int spacedim>
-void
-StraightBoundary<dim, spacedim>::
-get_intermediate_points_on_line (const typename Triangulation<dim, spacedim>::line_iterator &line,
-                                 std::vector<Point<spacedim> > &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<Point<1> > &line_points = this->get_line_support_points(n);
-
-  const Point<spacedim> vertices[2] = { line->vertex(0),
-                                        line->vertex(1)
-                                      };
-
-  for (unsigned int i=0; i<n; ++i)
-    {
-      const double x = line_points[1+i][0];
-      points[i] = (1-x)*vertices[0] + x*vertices[1];
-    }
-}
-
-
-
-
-template <int dim, int spacedim>
-void
-StraightBoundary<dim, spacedim>::
-get_intermediate_points_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &,
-                                 std::vector<Point<spacedim> > &) const
-{
-  Assert(false, ExcImpossibleInDim(dim));
-}
-
-
-
-template <>
-void
-StraightBoundary<3>::
-get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad,
-                                 std::vector<Point<3> > &points) const
-{
-  const unsigned int spacedim = 3;
-
-  const unsigned int n=points.size(),
-                     m=static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
-  // is n a square number
-  Assert(m*m==n, ExcInternalError());
-
-  const std::vector<Point<1> > &line_points = this->get_line_support_points(m);
-
-  const Point<spacedim> vertices[4] = { quad->vertex(0),
-                                        quad->vertex(1),
-                                        quad->vertex(2),
-                                        quad->vertex(3)
-                                      };
-
-  for (unsigned int i=0; i<m; ++i)
-    {
-      const double y=line_points[1+i][0];
-      for (unsigned int j=0; j<m; ++j)
-        {
-          const double x=line_points[1+j][0];
-          points[i*m+j]=((1-x) * vertices[0] +
-                         x     * vertices[1]) * (1-y) +
-                        ((1-x) * vertices[2] +
-                         x     * vertices[3]) * y;
-        }
-    }
-}
-
-
-
-template <>
-void
-StraightBoundary<2,3>::
-get_intermediate_points_on_quad (const Triangulation<2,3>::quad_iterator &quad,
-                                 std::vector<Point<3> > &points) const
-{
-  const unsigned int spacedim = 3;
-
-  const unsigned int n=points.size(),
-                     m=static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
-  // is n a square number
-  Assert(m*m==n, ExcInternalError());
-
-  const std::vector<Point<1> > &line_points = this->get_line_support_points(m);
-
-  const Point<spacedim> vertices[4] = { quad->vertex(0),
-                                        quad->vertex(1),
-                                        quad->vertex(2),
-                                        quad->vertex(3)
-                                      };
-
-  for (unsigned int i=0; i<m; ++i)
-    {
-      const double y=line_points[1+i][0];
-      for (unsigned int j=0; j<m; ++j)
-        {
-          const double x=line_points[1+j][0];
-          points[i*m+j]=((1-x) * vertices[0] +
-                         x     * vertices[1]) * (1-y) +
-                        ((1-x) * vertices[2] +
-                         x     * vertices[3]) * y;
-        }
-    }
-}
-
-
-
-template <int dim, int spacedim>
-Tensor<1,spacedim>
-StraightBoundary<dim,spacedim>::
-normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &face,
-               const Point<spacedim>                                     &p) const
-{
-  return FlatManifold<dim, spacedim>::normal_vector(face, p);
-}
-
-
-
-template <int dim, int spacedim>
-void
-StraightBoundary<dim, spacedim>::
-get_normals_at_vertices (const typename Triangulation<dim, spacedim>::face_iterator &face,
-                         typename Boundary<dim,spacedim>::FaceVertexNormals         &face_vertex_normals) const
-{
-  FlatManifold<dim, spacedim>::get_normals_at_vertices(face, face_vertex_normals);
-}
-
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-StraightBoundary<dim, spacedim>::
-project_to_surface (const typename Triangulation<dim, spacedim>::line_iterator &line,
-                    const Point<spacedim>                                &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<spacedim> 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 <typename Iterator, int spacedim, int dim>
-  Point<spacedim>
-  compute_projection (const Iterator        &object,
-                      const Point<spacedim> &y,
-                      std::integral_constant<int, dim>)
-  {
-    // 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||<eps.
-    //
-    // As for the Hessian, the best choice
-    // would be
-    //   H_k = J''(\xi_k)
-    // but we'll opt for the simpler
-    // Gauss-Newton form
-    //   H_k = A^T A
-    // i.e.
-    //   (H_k)_{nm} = \sum_{i,j} v_i*v_j *
-    //                   \partial_n phi_i *
-    //                   \partial_m phi_j
-    // we start at xi=(0.5,0.5).
-    Point<dim> xi;
-    for (unsigned int d=0; d<dim; ++d)
-      xi[d] = 0.5;
-
-    Point<spacedim> x_k;
-    for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-      x_k += object->vertex(i) *
-             GeometryInfo<dim>::d_linear_shape_function (xi, i);
-
-    do
-      {
-        Tensor<1,dim> F_k;
-        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-          F_k += (x_k-y)*object->vertex(i) *
-                 GeometryInfo<dim>::d_linear_shape_function_gradient (xi, i);
-
-        Tensor<2,dim> H_k;
-        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-          for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
-            {
-              Tensor<2, dim> tmp = outer_product(
-                                     GeometryInfo<dim>::d_linear_shape_function_gradient(xi, i),
-                                     GeometryInfo<dim>::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<spacedim>();
-        for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-          x_k += object->vertex(i) *
-                 GeometryInfo<dim>::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 <typename Iterator>
-  Point<1>
-  compute_projection (const Iterator &,
-                      const Point<1> &y,
-                      /* it's a quad: */std::integral_constant<int, 2>)
-  {
-    return y;
-  }
-
-  // specialization for a quad in 2d
-  template <typename Iterator>
-  Point<2>
-  compute_projection (const Iterator &,
-                      const Point<2> &y,
-                      /* it's a quad: */std::integral_constant<int, 2>)
-  {
-    return y;
-  }
-}
-
-
-
-
-
-template <>
-Point<3>
-StraightBoundary<1,3>::
-project_to_surface (const Triangulation<1, 3>::quad_iterator &,
-                    const Point<3>  &y) const
-{
-  return y;
-}
-
-//TODO[SP]: This is just a horrible way out to make it compile in codim 2.
-template <int dim, int spacedim>
-Point<spacedim>
-StraightBoundary<dim, spacedim>::
-project_to_surface (const typename Triangulation<dim, spacedim>::quad_iterator &quad,
-                    const Point<spacedim>  &y) const
-{
-  if (spacedim <= 2)
-    return y;
-  else
-    return internal::compute_projection (quad, y,
-                                         /* it's a quad */std::integral_constant<int, 2>());
-}
-
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-StraightBoundary<dim, spacedim>::
-project_to_surface (const typename Triangulation<dim, spacedim>::hex_iterator &,
-                    const Point<spacedim>                                &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<spacedim>();
-    }
-}
-
-
-
-// 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 (file)
index a896f63..0000000
+++ /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<deal_II_dimension, deal_II_space_dimension>;
-    template class StraightBoundary<deal_II_dimension, deal_II_space_dimension>;
-#endif
-}
-
-
-
diff --git a/source/grid/tria_boundary_lib.cc b/source/grid/tria_boundary_lib.cc
deleted file mode 100644 (file)
index dc87dcd..0000000
+++ /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 <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_accessor.h>
-#include <deal.II/base/tensor.h>
-#include <deal.II/base/std_cxx14/memory.h>
-#include <cmath>
-
-
-
-DEAL_II_NAMESPACE_OPEN
-
-
-
-template <int dim, int spacedim>
-CylinderBoundary<dim,spacedim>::CylinderBoundary (const double radius,
-                                                  const unsigned int axis)
-  :
-  radius(radius),
-  direction (get_axis_vector (axis)),
-  point_on_axis (Point<spacedim>())
-{}
-
-
-template <int dim, int spacedim>
-CylinderBoundary<dim,spacedim>::CylinderBoundary (const double           radius,
-                                                  const Point<spacedim> &direction,
-                                                  const Point<spacedim> &point_on_axis)
-  :
-  radius(radius),
-  direction (direction / direction.norm()),
-  point_on_axis (point_on_axis)
-{}
-
-
-
-template<int dim, int spacedim>
-std::unique_ptr<Manifold<dim, spacedim> >
-CylinderBoundary<dim,spacedim>::clone() const
-{
-  return std_cxx14::make_unique<CylinderBoundary<dim, spacedim> >(radius, direction, point_on_axis);
-}
-
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-CylinderBoundary<dim,spacedim>::get_axis_vector (const unsigned int axis)
-{
-  Assert (axis < spacedim, ExcIndexRange (axis, 0, spacedim));
-
-  Point<spacedim> axis_vector;
-  axis_vector[axis] = 1;
-  return axis_vector;
-}
-
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-CylinderBoundary<dim,spacedim>::
-get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const
-{
-  // compute a proposed new point
-  const Point<spacedim> middle = StraightBoundary<dim,spacedim>::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<spacedim>(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 <int dim, int spacedim>
-Point<spacedim>
-CylinderBoundary<dim,spacedim>::
-get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &) const
-{
-  Assert (false, ExcImpossibleInDim(dim));
-  return Point<spacedim>();
-}
-
-
-
-template <int dim, int spacedim>
-void
-CylinderBoundary<dim,spacedim>::get_intermediate_points_on_line (
-  const typename Triangulation<dim,spacedim>::line_iterator &line,
-  std::vector<Point<spacedim> > &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 <int dim, int spacedim>
-void
-CylinderBoundary<dim,spacedim>::get_intermediate_points_between_points (
-  const Point<spacedim> &v0,
-  const Point<spacedim> &v1,
-  std::vector<Point<spacedim> > &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<Point<1> > &line_points = this->get_line_support_points(n);
-
-  for (unsigned int i=0; i<n; ++i)
-    {
-      const double x = line_points[i+1][0];
-      const Point<spacedim> 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<spacedim>(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<Point<3> > &points) const
-{
-  if (points.size()==1)
-    points[0]=get_new_point_on_quad(quad);
-  else
-    {
-      unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
-      Assert(points.size()==m*m, ExcInternalError());
-
-      std::vector<Point<3> > lp0(m);
-      std::vector<Point<3> > lp1(m);
-
-      get_intermediate_points_on_line(quad->line(0), lp0);
-      get_intermediate_points_on_line(quad->line(1), lp1);
-
-      std::vector<Point<3> > lps(m);
-      for (unsigned int i=0; i<m; ++i)
-        {
-          get_intermediate_points_between_points(lp0[i], lp1[i], lps);
-
-          for (unsigned int j=0; j<m; ++j)
-            points[i*m+j]=lps[j];
-        }
-    }
-}
-
-
-
-template <int dim, int spacedim>
-void
-CylinderBoundary<dim,spacedim>::get_intermediate_points_on_quad (
-  const typename Triangulation<dim,spacedim>::quad_iterator &,
-  std::vector<Point<spacedim> > &) 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 <int dim, int spacedim>
-void
-CylinderBoundary<dim,spacedim>::
-get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
-                         typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const
-{
-  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
-    {
-      const Point<spacedim> 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 <int dim, int spacedim>
-double
-CylinderBoundary<dim,spacedim>::get_radius () const
-{
-  return radius;
-}
-
-
-//======================================================================//
-
-template <int dim>
-ConeBoundary<dim>::ConeBoundary (const double radius_0,
-                                 const double radius_1,
-                                 const Point<dim> x_0,
-                                 const Point<dim> x_1)
-  :
-  radius_0 (radius_0),
-  radius_1 (radius_1),
-  x_0 (x_0),
-  x_1 (x_1)
-{}
-
-
-
-template<int dim>
-std::unique_ptr<Manifold<dim,dim> > ConeBoundary<dim>::clone() const
-{
-  return std_cxx14::make_unique<ConeBoundary<dim> >(radius_0,radius_1,x_0,x_1);
-}
-
-
-
-template <int dim>
-double ConeBoundary<dim>::get_radius (Point<dim> 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 <int dim>
-void
-ConeBoundary<dim>::
-get_intermediate_points_between_points (const Point<dim> &p0,
-                                        const Point<dim> &p1,
-                                        std::vector<Point<dim> > &points) const
-{
-  const unsigned int n = points.size ();
-  const Tensor<1,dim> axis = x_1 - x_0;
-
-  Assert (n > 0, ExcInternalError ());
-
-  const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
-
-  for (unsigned int i=0; i<n; ++i)
-    {
-      const double x = line_points[i+1][0];
-
-      // Compute the current point.
-      const Point<dim> 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<dim> 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 <int dim>
-Point<dim>
-ConeBoundary<dim>::
-get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
-{
-  const Tensor<1,dim> axis = x_1 - x_0;
-  // Compute the middle point of the line.
-  const Point<dim> middle = StraightBoundary<dim>::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<dim> 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<dim> 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<dim> 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 <int dim>
-Point<dim>
-ConeBoundary<dim>::
-get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &) const
-{
-  Assert (false, ExcImpossibleInDim (dim));
-
-  return Point<dim>();
-}
-
-
-
-template <int dim>
-void
-ConeBoundary<dim>::
-get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
-                                 std::vector<Point<dim> > &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<Point<3> > &points) const
-{
-  if (points.size () == 1)
-    points[0] = get_new_point_on_quad (quad);
-  else
-    {
-      unsigned int n = static_cast<unsigned int> (std::sqrt (static_cast<double> (points.size ())));
-
-      Assert (points.size () == n * n, ExcInternalError ());
-
-      std::vector<Point<3> > points_line_0 (n);
-      std::vector<Point<3> > 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<Point<3> > 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 <int dim>
-void
-ConeBoundary<dim>::
-get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &,
-                                 std::vector<Point<dim> > &) 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 <int dim>
-void
-ConeBoundary<dim>::
-get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
-                         typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
-{
-  const Tensor<1,dim> axis = x_1 - x_0;
-
-  for (unsigned int vertex = 0; vertex < GeometryInfo<dim>::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<dim> vertex_p = x_0 + c * axis;
-      // Then compute the vector pointing from the point <tt>vertex_p</tt> 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 <int dim>
-inline
-Tensor<1,dim>
-ConeBoundary<dim>::
-normal_vector (const typename Triangulation<dim>::face_iterator &,
-               const Point<dim> &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<dim>(), ExcNotImplemented());
-  for (unsigned int d=0; d<dim; ++d)
-    if (d != dim-1)  // don't test the last component of the vector
-      Assert (this->x_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 <int dim, int spacedim>
-HyperBallBoundary<dim,spacedim>::HyperBallBoundary (const Point<spacedim> p,
-                                                    const double     radius)
-  :
-  center(p),
-  radius(radius),
-  compute_radius_automatically(false)
-{}
-
-
-
-template<int dim, int spacedim>
-std::unique_ptr<Manifold<dim, spacedim> >
-HyperBallBoundary<dim,spacedim>::clone() const
-{
-  return std_cxx14::make_unique<HyperBallBoundary<dim,spacedim> >(center, radius);
-}
-
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-HyperBallBoundary<dim,spacedim>::get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const
-{
-  Point<spacedim> middle = StraightBoundary<dim,spacedim>::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 <int dim, int spacedim>
-Point<spacedim>
-HyperBallBoundary<dim,spacedim>::
-get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const
-{
-  Point<spacedim> middle = StraightBoundary<dim,spacedim>::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<Point<1> > &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
-
-
-
-template <int dim, int spacedim>
-void
-HyperBallBoundary<dim,spacedim>::get_intermediate_points_on_line (
-  const typename Triangulation<dim,spacedim>::line_iterator &line,
-  std::vector<Point<spacedim> > &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 <int dim, int spacedim>
-void
-HyperBallBoundary<dim,spacedim>::get_intermediate_points_between_points (
-  const Point<spacedim> &p0, const Point<spacedim> &p1,
-  std::vector<Point<spacedim> > &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)<eps*r*r)
-         &&
-         (std::fabs(v1*v1-r*r)<eps*r*r),
-         ExcMessage("In computing the location of midpoints of an edge of a "
-                    "HyperBallBoundary, one of the vertices of the edge "
-                    "does not have the expected distance from the center "
-                    "of the sphere. This may happen if the geometry has "
-                    "been deformed, or if the HyperBallBoundary object has "
-                    "been assigned to a part of the boundary that is not "
-                    "in fact a sphere (e.g., the sides of a quarter shell "
-                    "or one octant of a ball)."));
-
-  const double alpha=std::acos((v0*v1)/std::sqrt((v0*v0)*(v1*v1)));
-  const Tensor<1,spacedim> pm=0.5*(v0+v1);
-
-  const double h=pm.norm();
-
-  // n even:  m=n/2,
-  // n odd:   m=(n-1)/2
-  const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
-  const unsigned int m=n/2;
-  for (unsigned int i=0; i<m ; ++i)
-    {
-      const double beta = alpha * (line_points[i+1][0]-0.5);
-      const double d = h*std::tan(beta);
-      points[i]      = Point<spacedim>(pm+d/length*(v1-v0));
-      points[n-1-i]  = Point<spacedim>(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<spacedim>(pm);
-
-
-  // project the points from the straight line to the HyperBallBoundary
-  for (unsigned int i=0; i<n; ++i)
-    {
-      points[i] *= r / std::sqrt(points[i].square());
-      points[i] += center;
-    }
-}
-
-
-
-template <>
-void
-HyperBallBoundary<3>::get_intermediate_points_on_quad (
-  const Triangulation<3>::quad_iterator &quad,
-  std::vector<Point<3> > &points) const
-{
-  if (points.size()==1)
-    points[0]=get_new_point_on_quad(quad);
-  else
-    {
-      unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
-      Assert(points.size()==m*m, ExcInternalError());
-
-      std::vector<Point<3> > lp0(m);
-      std::vector<Point<3> > lp1(m);
-
-      get_intermediate_points_on_line(quad->line(0), lp0);
-      get_intermediate_points_on_line(quad->line(1), lp1);
-
-      std::vector<Point<3> > lps(m);
-      for (unsigned int i=0; i<m; ++i)
-        {
-          get_intermediate_points_between_points(lp0[i], lp1[i], lps);
-
-          for (unsigned int j=0; j<m; ++j)
-            points[i*m+j]=lps[j];
-        }
-    }
-}
-
-
-
-template <>
-void
-HyperBallBoundary<2,3>::get_intermediate_points_on_quad (
-  const Triangulation<2,3>::quad_iterator &quad,
-  std::vector<Point<3> > &points) const
-{
-  if (points.size()==1)
-    points[0]=get_new_point_on_quad(quad);
-  else
-    {
-      unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
-      Assert(points.size()==m*m, ExcInternalError());
-
-      std::vector<Point<3> > lp0(m);
-      std::vector<Point<3> > lp1(m);
-
-      get_intermediate_points_on_line(quad->line(0), lp0);
-      get_intermediate_points_on_line(quad->line(1), lp1);
-
-      std::vector<Point<3> > lps(m);
-      for (unsigned int i=0; i<m; ++i)
-        {
-          get_intermediate_points_between_points(lp0[i], lp1[i], lps);
-
-          for (unsigned int j=0; j<m; ++j)
-            points[i*m+j]=lps[j];
-        }
-    }
-}
-
-
-
-template <int dim, int spacedim>
-void
-HyperBallBoundary<dim,spacedim>::get_intermediate_points_on_quad (
-  const typename Triangulation<dim,spacedim>::quad_iterator &,
-  std::vector<Point<spacedim> > &) const
-{
-  Assert(false, ExcImpossibleInDim(dim));
-}
-
-
-
-template <int dim, int spacedim>
-Tensor<1,spacedim>
-HyperBallBoundary<dim,spacedim>::
-normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &,
-               const Point<spacedim> &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 <int dim, int spacedim>
-void
-HyperBallBoundary<dim,spacedim>::
-get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
-                         typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const
-{
-  for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
-    face_vertex_normals[vertex] = face->vertex(vertex)-center;
-}
-
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-HyperBallBoundary<dim,spacedim>::get_center () const
-{
-  return center;
-}
-
-
-
-template <int dim, int spacedim>
-double
-HyperBallBoundary<dim,spacedim>::get_radius () const
-{
-  Assert(!compute_radius_automatically, ExcRadiusNotSet());
-  return radius;
-}
-
-
-/* ---------------------------------------------------------------------- */
-
-
-template <int dim>
-HalfHyperBallBoundary<dim>::HalfHyperBallBoundary (const Point<dim> center,
-                                                   const double     radius) :
-  HyperBallBoundary<dim> (center, radius)
-{}
-
-
-
-template<int dim>
-std::unique_ptr<Manifold<dim,dim> > HalfHyperBallBoundary<dim>::clone() const
-{
-  return std_cxx14::make_unique<HalfHyperBallBoundary<dim> >(this->get_center(),
-                                                             this->get_radius());
-}
-
-
-
-template <int dim>
-Point<dim>
-HalfHyperBallBoundary<dim>::
-get_new_point_on_line (const typename Triangulation<dim>::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<dim> line_center = line->center();
-  const Point<dim> 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<dim>::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 <int dim>
-Point<dim>
-HalfHyperBallBoundary<dim>::
-get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
-{
-  const Point<dim> quad_center = quad->center();
-  if (quad_center(0) == this->center(0))
-    return quad_center;
-  else
-    return HyperBallBoundary<dim>::get_new_point_on_quad (quad);
-}
-
-
-
-template <int dim>
-void
-HalfHyperBallBoundary<dim>::
-get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
-                                 std::vector<Point<dim> > &points) const
-{
-  // check whether center of object is at x==0, since then it belongs to the
-  // plane part of the boundary
-  const Point<dim> line_center = line->center();
-  if (line_center(0) == this->center(0))
-    return StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
-  else
-    return HyperBallBoundary<dim>::get_intermediate_points_on_line (line, points);
-}
-
-
-
-template <int dim>
-void
-HalfHyperBallBoundary<dim>::
-get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
-                                 std::vector<Point<dim> > &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<dim> quad_center = quad->center();
-      if (quad_center(0) == this->center(0))
-        StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
-      else
-        HyperBallBoundary<dim>::get_intermediate_points_on_quad (quad, points);
-    }
-}
-
-
-
-template <>
-void
-HalfHyperBallBoundary<1>::
-get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &,
-                                 std::vector<Point<1> > &) 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 <int dim>
-void
-HalfHyperBallBoundary<dim>::
-get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
-                         typename Boundary<dim>::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<dim> quad_center = face->center();
-  if (quad_center(0) == this->center(0))
-    StraightBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
-  else
-    HyperBallBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
-}
-
-
-/* ---------------------------------------------------------------------- */
-
-
-
-template <int dim>
-HyperShellBoundary<dim>::HyperShellBoundary (const Point<dim> &center)
-  :
-  HyperBallBoundary<dim>(center, 0.)
-{
-  this->compute_radius_automatically=true;
-}
-
-
-
-template<int dim>
-std::unique_ptr<Manifold<dim,dim> >
-HyperShellBoundary<dim>::clone() const
-{
-  return std_cxx14::make_unique<HyperShellBoundary<dim> >(this->get_center());
-}
-
-
-
-/* ---------------------------------------------------------------------- */
-
-
-
-
-template <int dim>
-HalfHyperShellBoundary<dim>::HalfHyperShellBoundary (const Point<dim> &center,
-                                                     const double inner_radius,
-                                                     const double outer_radius)
-  :
-  HyperShellBoundary<dim> (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<int dim>
-std::unique_ptr<Manifold<dim,dim> >
-HalfHyperShellBoundary<dim>::clone() const
-{
-  return std_cxx14::make_unique<HalfHyperShellBoundary<dim> >(this->get_center(),
-                                                              inner_radius,
-                                                              outer_radius);
-}
-
-
-
-template <int dim>
-Point<dim>
-HalfHyperShellBoundary<dim>::
-get_new_point_on_line (const typename Triangulation<dim>::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<dim>::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<dim>::get_new_point_on_line (line);
-    }
-
-    default:
-      Assert (false, ExcNotImplemented());
-    }
-
-  return Point<dim>();
-}
-
-
-
-template <>
-Point<1>
-HalfHyperShellBoundary<1>::
-get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
-{
-  Assert (false, ExcInternalError());
-  return Point<1>();
-}
-
-
-
-
-template <int dim>
-Point<dim>
-HalfHyperShellBoundary<dim>::
-get_new_point_on_quad (const typename Triangulation<dim>::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<dim> 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<dim>::get_new_point_on_quad (quad);
-}
-
-
-
-template <int dim>
-void
-HalfHyperShellBoundary<dim>::
-get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
-                                 std::vector<Point<dim> > &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<dim>::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<dim>::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<dim>::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<dim>::get_intermediate_points_on_line (line, points);
-
-      break;
-    }
-
-    default:
-      Assert (false, ExcNotImplemented());
-    }
-}
-
-
-
-template <int dim>
-void
-HalfHyperShellBoundary<dim>::
-get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
-                                 std::vector<Point<dim> > &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<dim> quad_center = quad->center();
-  if (quad_center(0) == this->center(0))
-    StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
-  else
-    HyperShellBoundary<dim>::get_intermediate_points_on_quad (quad, points);
-}
-
-
-
-template <>
-void
-HalfHyperShellBoundary<1>::
-get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &,
-                                 std::vector<Point<1> > &) 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 <int dim>
-void
-HalfHyperShellBoundary<dim>::
-get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
-                         typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
-{
-  if (face->center()(0) == this->center(0))
-    StraightBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
-  else
-    HyperShellBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
-}
-
-
-
-
-template <int dim, int spacedim>
-TorusBoundary<dim,spacedim>::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<int dim, int spacedim>
-std::unique_ptr<Manifold<dim, spacedim> >
-TorusBoundary<dim,spacedim>::clone() const
-{
-  return std_cxx14::make_unique<TorusBoundary<dim,spacedim> >(R,r);
-}
-
-
-
-template <int dim, int spacedim>
-double
-TorusBoundary<dim,spacedim>::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<Point<1> > &line_points = this->get_line_support_points(npoints);
-  for (unsigned int i=0; i<npoints; i++)
-    {
-      const double x = line_points[i+1][0];
-      target=  (1-x)*p[0] + x*p[1];
-      points[i]=get_real_coord(target);
-    }
-}
-
-
-
-template <>
-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<unsigned int>(std::sqrt(static_cast<double>(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<Point<1> > &line_points = this->get_line_support_points(m);
-  for (unsigned int i=0; i<m; ++i)
-    {
-      const double y=line_points[i+1][0];
-      for (unsigned int j=0; j<m; ++j)
-        {
-          const double x=line_points[j+1][0];
-          target=((1-x) * p[0] +
-                  x     * p[1]) * (1-y) +
-                 ((1-x) * p[2] +
-                  x     * p[3]) * y;
-
-          points[i*m+j]=get_real_coord(target);
-        }
-    }
-}
-
-
-
-template <>
-void
-TorusBoundary<2,3>::
-get_normals_at_vertices (const Triangulation<2,3 >::face_iterator &face,
-                         Boundary<2,3>::FaceVertexNormals &face_vertex_normals) const
-{
-  for (unsigned int i=0; i<GeometryInfo<2>::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 (file)
index 6933a50..0000000
+++ /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<deal_II_dimension>;
-    template class ConeBoundary<deal_II_dimension>;
-    template class HyperBallBoundary<deal_II_dimension>;
-    template class HalfHyperBallBoundary<deal_II_dimension>;
-    template class HyperShellBoundary<deal_II_dimension>;
-    template class HalfHyperShellBoundary<deal_II_dimension>;
-
-#if deal_II_dimension != 3
-    template class HyperBallBoundary<deal_II_dimension,deal_II_dimension+1>;
-#endif
-#if deal_II_dimension == 3
-    template class CylinderBoundary<deal_II_dimension-1,deal_II_dimension>;
-#endif
-}
-
-
index d7317b04e5e10a52f65711703f13d49f98b4fc60..25fdc8afb66a304b71dcc3b2839522e15c31da2e 100644 (file)
@@ -57,7 +57,6 @@
 #include <deal.II/grid/grid_refinement.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
 
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/dofs/dof_tools.h>
index 8d4a31994d5eb84c6d4dc78aacfd2bd740e7db06..dca99113200b8c90ba44daab452859295b09de07 100644 (file)
@@ -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
 
 
 
index 125e375919668ce9445c2c782a7bbbb0acef3d96..a159142555de42bba2bde126b7d32980552e8891 100644 (file)
@@ -20,7 +20,6 @@
 #include "../tests.h"
 #include <deal.II/base/function.h>
 #include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/matrix_free/matrix_free.h>
 #include <deal.II/matrix_free/fe_evaluation.h>
@@ -161,8 +160,6 @@ void test ()
 {
   Triangulation<dim> tria;
   GridGenerator::hyper_ball (tria);
-  static const HyperBallBoundary<dim> boundary;
-  tria.set_boundary (0, boundary);
   tria.begin_active ()->set_refine_flag();
   tria.execute_coarsening_and_refinement();
   typename Triangulation<dim>::active_cell_iterator cell, endc;
index 13dff441e6b1e92f014dc185789754644d8ad955..b9a90ab0e7308ab73348ab824378dcadcf14e688 100644 (file)
@@ -28,7 +28,6 @@
 #include <deal.II/lac/vector.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/dofs/dof_tools.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/lac/constraint_matrix.h>
index e5515b4dc737f8adcb2ed42806221a10e38b5d27..02baa02c7ca286ae8aacfa0a1c64ad7bd73ce039 100644 (file)
@@ -32,8 +32,6 @@ void test ()
 {
   Triangulation<dim> tria;
   GridGenerator::hyper_ball (tria);
-  static const HyperBallBoundary<dim> boundary;
-  tria.set_boundary (0, boundary);
   tria.refine_global(4-dim);
   tria.begin_active()->set_refine_flag();
   tria.last()->set_refine_flag();
index b5ce5d95b30542364a43179b4196ba31f8acb1c5..068dbc2cd36f35301fb51d662029b7db90082499 100644 (file)
@@ -26,7 +26,6 @@
 #include <deal.II/lac/vector.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/dofs/dof_tools.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/lac/constraint_matrix.h>
index 92adece7d9281c32a5bcdf509d61b8b44e31f22b..c87826aa8767f40d4cb5530de9362edcfe150bfc 100644 (file)
@@ -27,7 +27,6 @@
 #include <deal.II/lac/precondition.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/dofs/dof_tools.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/lac/constraint_matrix.h>
index 5fe468aec506db6b45eda6d30129133c73008d2b..0adc32168fde9dc49f129c5bc18637bb046d8978 100644 (file)
@@ -28,7 +28,6 @@
 #include <deal.II/lac/precondition.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/dofs/dof_tools.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/lac/constraint_matrix.h>
index d4cd524236c42a4db34a9155ab5062ba871141fb..227e2815aadadff6f2dd459f5ddb862f0971ee68 100644 (file)
@@ -29,7 +29,6 @@
 #include <deal.II/lac/precondition.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/dofs/dof_tools.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/lac/constraint_matrix.h>
index 1d5308071631a25e4b1b4d35e0b31e199240664f..f16aec4c77dabde0788799459b99316fa59f8152 100644 (file)
@@ -35,7 +35,6 @@
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/grid/grid_generator.h>
 
 #include <deal.II/multigrid/multigrid.h>

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.