* Generate a hyper-ball intersected with the positive orthant relate to @p
* center, which contains three elements in 2d and four in 3d. The boundary
* indicators for the final triangulations are 0 for the curved boundary
- * and 1 for the cut plane. The appropriate boundary class is
- * HyperBallBoundary.
+ * and 1 for the cut plane. The appropriate manifold class is
+ * SphericalManifold.
*/
void generate_quarter_hyper_ball(PointWrapper ¢er,
const double radius = 1.);
* Generate a half hyper-ball around @p center, which contains four elements
* in 2d and 6 in 3d. The cut plane is perpendicular to the x-axis. The
* boundary indicators for the final triangulation are 0 for the curved
- * boundary and 1 for the cut plane. The appropriate boundary class is
- * HalfHyperBallBoundary, or HyperBallBoundary.
+ * boundary and 1 for the cut plane. The appropriate manifold class is
+ * SphericalManifold.
*/
void generate_half_hyper_ball(PointWrapper ¢er,
const double radius = 1.);
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"
;
"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"
;
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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
- */
<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
#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>
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
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.
*/
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
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
template <typename Accessor> class TriaIterator;
template <typename Accessor> class TriaActiveIterator;
-template <int dim, int spacedim> class Boundary;
template <int dim, int spacedim> class Manifold;
*/
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.
*
-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
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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ü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> ¢er = 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> ¢er = 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
#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>
#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>
#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>
#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>
#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>
#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>
#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>
#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>
#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>
}
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
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);
}
}
}
{
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;
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);
}
}
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);
}
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
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
)
#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>
}
-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
-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,
}
-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
}
-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> &
#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>
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);
make_array_view(points_and_weights.second.begin(),
points_and_weights.second.end()));
}
+ else
+ {
+ return get_new_point_on_object(obj);
+ }
}
}
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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
-}
-
-
-
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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> ¢er)
- :
- 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> ¢er,
- 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
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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
-}
-
-
#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>
-// 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
#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>
{
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;
#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>
{
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();
#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>
#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>
#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>
#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>
#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>
projects / dealii.git / commitdiff