/**
- * This class is used to represent a boundary to a triangulation.
+ * This class is used to represent a boundary to a triangulation,
+ * and as such, it is a specialization of the Manifold<spacedim>
+ * class, where the topological dimension of the Manifold is fixed
+ * to be (dim-1), i.e., the boundary of a triangulation with
+ * topological dimension dim.
+ *
* When a triangulation creates a new vertex on the boundary of the
* domain, it determines the new vertex' coordinates through the
* following code (here in two dimensions):
* @code
* ...
- * Point<2> new_vertex = boundary.get_new_point_on_line (line);
+ * Point<2> new_vertex = boundary.get_new_point (points, weights);
* ...
* @endcode
- * @p line denotes the line at the boundary that shall be refined
- * and for which we seek the common point of the two child lines.
+ * @p points is a vector containing the points which define the object
+ * that shall be refined, while @p weights are used to decide at what
+ * location the child point should be created with respect to its
+ * parents.
*
* In 3D, a new vertex may be placed on the middle of a line or on
* the middle of a side. Respectively, the library calls
* @code
* ...
* Point<3> new_line_vertices[4]
- * = { boundary.get_new_point_on_line (face->line(0)),
- * boundary.get_new_point_on_line (face->line(1)),
- * boundary.get_new_point_on_line (face->line(2)),
- * boundary.get_new_point_on_line (face->line(3)) };
+ * = { boundary.get_new_point (vertices_of_line_0, w_line),
+ * boundary.get_new_point (vertices_of_line_1, w_line),
+ * boundary.get_new_point (vertices_of_line_2, w_line),
+ * boundary.get_new_point (vertices_of_line_3, w_line) };
* ...
* @endcode
- * to get the four midpoints of the lines bounding the quad at the
+ * where @p w_line is a vector containing the values (1/2, 1/2).
+ * This return the four midpoints of the lines bounding the quad at the
* boundary, and after that
* @code
* ...
- * Point<3> new_quad_vertex = boundary.get_new_point_on_quad (face);
+ * Point<3> new_quad_vertex = boundary.get_new_point (vertices_of_face, w_quad);
* ...
* @endcode
* to get the midpoint of the face. It is guaranteed that this order
* (first lines, then faces) holds, so you can use information from
- * the children of the four lines of a face, since these already exist
- * at the time the midpoint of the face is to be computed.
+ * the children of the four lines of a face, since these already
+ * exist at the time the midpoint of the face is to be
+ * computed. This in fact is exploited in the library by passing a
+ * list of 8 points (the actual vertices, and the four points
+ * computed above) together with 8 weights which minimize cell
+ * distortion.
*
- * Since iterators are passed to the functions, you may use information
- * about boundary indicators and the like, as well as all other information
- * provided by these objects.
+ * This class is derived from FlatManifold<spacedim>. A minimal
+ * implementation for derived classes is given by overloading the
+ * method Manifold<spacedim>::project_to_manifold, which is called
+ * internally by the FlatManifold<spacedim> class (from which
+ * Boundary<dim,spacedim> is derived) with a guess computed by
+ * considering the Manifold "flat", i.e., by putting the new point
+ * in the weighted average of the surrounding points.
*
* There are specialisations, StraightBoundary, which places
* the new point right into the middle of the given points, and
* HyperBallBoundary creating a hyperball with given radius
* around a given center point.
*
+ * @deprecated A new Manifold<spacedim> class was introduced which
+ * generalises the functionality of this
+ * class. Boundary<dim,spacedim> is no longer necessary, and should
+ * be replaced by any equivalent Manifold<spacedim> description.
+ *
* @ingroup boundary
- * @author Wolfgang Bangerth, 1999, 2001, 2009, Ralf Hartmann, 2001, 2008
+ * @author Wolfgang Bangerth, 1999, 2001, 2009, Ralf Hartmann, 2001, 2008,
+ * Luca Heltai 2013, 2014
*/
template <int dim, int spacedim=dim>
class Boundary : public FlatManifold<spacedim>
{
public:
- /**
- * Type keeping information about the normals at the vertices of a face of a
- * cell. Thus, there are <tt>GeometryInfo<dim>::vertices_per_face</tt>
- * normal vectors, that define the tangent spaces of the boundary at the
- * vertices. Note that the vectors stored in this object are not required to
- * be normalized, nor to actually point outward, as one often will only want
- * to check for orthogonality to define the tangent plane; if a function
- * requires the normals to be normalized, then it must do so itself.
- *
- * For obvious reasons, this type is not useful in 1d.
- */
- typedef Tensor<1,spacedim> FaceVertexNormals[GeometryInfo<dim>::vertices_per_face];
-
/**
* Destructor. Does nothing here, but needs to be declared to make it
* virtual.
*/
virtual ~Boundary ();
-
- /**
- * Return the point which shall become the new middle vertex of the two
- * children of a regular line. In 2D, this line is a line at the boundary,
- * while in 3d, it is bounding a face at the boundary (the lines therefore
- * is also on the boundary).
- */
- virtual
- Point<spacedim>
- get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const = 0;
-
- /**
- * Return the point which shall become the common point of the four children
- * of a quad at the boundary in three or more spatial dimensions. This
- * function therefore is only useful in at least three dimensions and should
- * not be called for lower dimensions.
- *
- * This function is called after the four lines bounding the given @p quad
- * are refined, so you may want to use the information provided by
- * <tt>quad->line(i)->child(j)</tt>, <tt>i=0...3</tt>, <tt>j=0,1</tt>.
- *
- * Because in 2D, this function is not needed, it is not made pure virtual,
- * to avoid the need to overload it. The default implementation throws an
- * error in any case, however.
- */
- virtual
- Point<spacedim>
- get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
-
- /**
- * Depending on <tt>dim=2</tt> or <tt>dim=3</tt> this function calls the
- * get_new_point_on_line or the get_new_point_on_quad function. It throws an
- * exception for <tt>dim=1</tt>. This wrapper allows dimension independent
- * programming.
- */
- Point<spacedim>
- get_new_point_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face) const;
-
- /**
- * Return intermediate points on a line spaced according to the interior
- * support points of the 1D Gauss-Lobatto quadrature formula.
- *
- * The number of points requested is given by the size of the vector @p
- * points. It is the task of the derived classes to arrange the points in
- * approximately equal distances.
- *
- * This function is called by the @p MappingQ class. This happens on each
- * face line of a cells that has got at least one boundary line.
- *
- * As this function is not needed for @p MappingQ1, it is not made pure
- * virtual, to avoid the need to overload it. The default implementation
- * throws an error in any case, however.
- */
- virtual
- void
- get_intermediate_points_on_line (const typename Triangulation<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.
- *
- * This function is called by the <tt>MappingQ<3></tt> class. This happens
- * each face quad of cells in 3d that has got at least one boundary face
- * quad.
- *
- * As this function is not needed for @p MappingQ1, it is not made pure
- * virtual, to avoid the need to overload it. The default implementation
- * throws an error in any case, however.
- */
- virtual
- void
- get_intermediate_points_on_quad (const typename Triangulation<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;
-
- /**
- * Return the normal vector to the surface at the point p. If p is not in
- * fact on the surface, but only close-by, try to return something
- * reasonable, for example the normal vector at the surface point closest to
- * p. (The point p will in fact not normally lie on the actual surface, but
- * rather be a quadrature point mapped by some polynomial mapping; the
- * mapped surface, however, will not usually coincide with the actual
- * surface.)
- *
- * The face iterator gives an indication which face this function is
- * supposed to compute the normal vector for. This is useful if the
- * boundary of the domain is composed of different nondifferential pieces
- * (for example when using the StraightBoundary class to approximate a
- * geometry that is completely described by the coarse mesh, with piecewise
- * (bi-)linear components between the vertices, but where the boundary may
- * have a kink at the vertices itself).
- *
- * @note Implementations of this function should be able to assume that the
- * point p lies within or close to the face described by the first
- * argument. In turn, callers of this function should ensure that this is in
- * fact the case.
- */
- virtual
- Tensor<1,spacedim>
- normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &face,
- const Point<spacedim> &p) const;
-
- /**
- * Compute the normal vectors to the boundary at each vertex of the given
- * face. It is not required that the normal vectors be normed
- * somehow. Neither is it required that the normals actually point outward.
- *
- * This function is needed to compute data for C1 mappings. The default
- * implementation is to throw an error, so you need not overload this
- * function in case you do not intend to use C1 mappings.
- *
- * Note that when computing normal vectors at a vertex where the boundary is
- * not differentiable, you have to make sure that you compute the one-sided
- * limits, i.e. limit with respect to points inside the given face.
- */
- virtual
- void
- get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
- FaceVertexNormals &face_vertex_normals) const;
-
- /**
- * Given a candidate point and a line segment characterized by the iterator,
- * return a point that lies on the surface described by this object. This
- * function is used in some mesh smoothing algorithms that try to move
- * around points in order to improve the mesh quality but need to ensure
- * that points that were on the boundary remain on the boundary.
- *
- * If spacedim==1, then the line represented by the line iterator is the
- * entire space (i.e. it is a cell, not a part of the boundary), and the
- * returned point equals the given input point.
- *
- * Derived classes do not need to implement this function unless mesh
- * smoothing algorithms are used with a particular boundary object. The
- * default implementation of this function throws an exception of type
- * ExcPureFunctionCalled.
- */
- virtual
- Point<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:
- /**
- * Returns the support points of the Gauss-Lobatto quadrature formula used
- * for intermediate points.
- *
- * @note Since the boundary description is closely tied to the unit cell
- * support points of MappingQ, new boundary descriptions need to explicitly
- * use these Gauss-Lobatto points and not equidistant points.
- */
- const std::vector<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_cxx1x::shared_ptr<QGaussLobatto<1> > > points;
-
- /**
- * Mutex for protecting the points array.
- */
- mutable Threads::Mutex mutex;
+
} DEAL_II_DEPRECATED;
-
/**
* 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.
+ * as the weighted average of the surrounding points.
*
* This class does not really describe a boundary in the usual
* sense. By placing new points in the middle of old ones, it rather
* polygon/polyhedron defined by the boundary of the initial coarse
* triangulation.
*
+ * @deprecated The functionality of this class is equivalent, but
+ * less general, to that of FlatManifold<spacedim>. Please use
+ * FlatManifold<spacedim> instead of this
+ * class. StraightBoundary<dim,spacedim> will be removed in future
+ * releases.
+ *
* @ingroup boundary
*
- * @author Wolfgang Bangerth, 1998, 2001, Ralf Hartmann, 2001
+ * @author Wolfgang Bangerth, 1998, 2001, Ralf Hartmann, 2001, Luca
+ * Heltai 2013, 2014
*/
template <int dim, int spacedim=dim>
class StraightBoundary : public Boundary<dim,spacedim>
* Default constructor. Some compilers require this for some reasons.
*/
StraightBoundary ();
-
- /**
- * Let the new point be the arithmetic mean of the two vertices of the line.
- *
- * Refer to the general documentation of this class and the documentation of
- * the base class for more information.
- */
- virtual Point<spacedim>
- get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
-
- /**
- * Let the new point be the arithmetic mean of the four vertices of this
- * quad and the four midpoints of the lines, which are already created at
- * the time of calling this function.
- *
- * Refer to the general documentation of this class and the documentation of
- * the base class for more information.
- */
- virtual
- Point<spacedim>
- get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
-
- /**
- * Gives <tt>n=points.size()</tt> points that splits the StraightBoundary
- * line into $n+1$ partitions of equal lengths.
- *
- * Refer to the general documentation of this class and the documentation of
- * the base class.
- */
- virtual
- void
- get_intermediate_points_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line,
- std::vector<Point<spacedim> > &points) const;
-
- /**
- * Gives <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;
-
- /**
- * 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;
-
- /**
- * 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;
-
- /**
- * 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;
-
- /**
- * 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;
-
- /**
- * 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;
-};
-
-
-
-/* -------------- declaration of explicit specializations ------------- */
-
-#ifndef DOXYGEN
-
-template <>
-Point<1>
-Boundary<1,1>::
-get_new_point_on_face (const Triangulation<1,1>::face_iterator &) const;
-
-template <>
-void
-Boundary<1,1>::
-get_intermediate_points_on_face (const Triangulation<1,1>::face_iterator &,
- std::vector<Point<1> > &) const;
-
-template <>
-Point<2>
-Boundary<1,2>::
-get_new_point_on_face (const Triangulation<1,2>::face_iterator &) const;
-
-template <>
-void
-Boundary<1,2>::
-get_intermediate_points_on_face (const Triangulation<1,2>::face_iterator &,
- std::vector<Point<2> > &) const;
-
-
-
-template <>
-Point<3>
-Boundary<1,3>::
-get_new_point_on_face (const Triangulation<1,3>::face_iterator &) const;
-
-template <>
-void
-Boundary<1,3>::
-get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &,
- std::vector<Point<3> > &) const;
-template <>
-void
-StraightBoundary<1,1>::
-get_normals_at_vertices (const Triangulation<1,1>::face_iterator &,
- Boundary<1,1>::FaceVertexNormals &) const;
-template <>
-void
-StraightBoundary<2,2>::
-get_normals_at_vertices (const Triangulation<2,2>::face_iterator &face,
- Boundary<2,2>::FaceVertexNormals &face_vertex_normals) const;
-template <>
-void
-StraightBoundary<3,3>::
-get_normals_at_vertices (const Triangulation<3,3>::face_iterator &face,
- Boundary<3,3>::FaceVertexNormals &face_vertex_normals) const;
-
-template <>
-Point<3>
-StraightBoundary<3,3>::
-get_new_point_on_quad (const Triangulation<3,3>::quad_iterator &quad) const;
-
-template <>
-void
-StraightBoundary<1,1>::
-get_intermediate_points_on_line (const Triangulation<1,1>::line_iterator &,
- std::vector<Point<1> > &) 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_DEPRECATED;
DEAL_II_NAMESPACE_CLOSE
for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
{
const Triangulation<dim>::line_iterator line = cell->line(line_no);
+ std::vector<Point<dim> > points(2);
+ points[0] = line->vertex(0);
+ points[1] = line->vertex(1);
if (line->at_boundary())
{
// from the boundary description
const Manifold<dim> &boundary
= line->get_boundary();
-
- std::vector<Point<dim> > points(2);
+
std::vector<Point<dim> > face_vertex_normals(2);
- points[0] = line->vertex(0);
- points[1] = line->vertex(1);
boundary.get_normals_at_points (face_vertex_normals, points);
// then transform them into interpolation points for a cubic
// not at boundary
{
static const StraightBoundary<dim> straight_boundary;
- straight_boundary.get_intermediate_points_on_line (line, line_points);
+ straight_boundary.get_intermediate_points (line_points, points);
a.insert (a.end(), line_points.begin(), line_points.end());
};
};
//
// ---------------------------------------------------------------------
-#include <deal.II/base/tensor.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 --------------------- */
Boundary<dim, spacedim>::~Boundary ()
{}
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-Boundary<dim, spacedim>::
-get_new_point_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &) const
-{
- Assert (false, ExcPureFunctionCalled());
- return Point<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>
-Point<spacedim>
-Boundary<dim,spacedim>::
-get_new_point_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face) const
-{
- Assert (dim>1, ExcImpossibleInDim(dim));
-
- switch (dim)
- {
- case 2:
- return get_new_point_on_line (face);
- case 3:
- return get_new_point_on_quad (face);
- }
-
- return Point<spacedim>();
-}
-
-
-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 <>
-Point<1>
-Boundary<1,1>::
-get_new_point_on_face (const Triangulation<1,1>::face_iterator &) const
-{
- Assert (false, ExcImpossibleInDim(1));
- return Point<1>();
-}
-
-
-template <>
-void
-Boundary<1,1>::
-get_intermediate_points_on_face (const Triangulation<1,1>::face_iterator &,
- std::vector<Point<1> > &) const
-{
- Assert (false, ExcImpossibleInDim(1));
-}
-
-
-
-template <>
-Point<2>
-Boundary<1,2>::
-get_new_point_on_face (const Triangulation<1,2>::face_iterator &) const
-{
- Assert (false, ExcImpossibleInDim(1));
- return Point<2>();
-}
-
-
-template <>
-void
-Boundary<1,2>::
-get_intermediate_points_on_face (const Triangulation<1,2>::face_iterator &,
- std::vector<Point<2> > &) const
-{
- Assert (false, ExcImpossibleInDim(1));
-}
-
-
-
-template <>
-Point<3>
-Boundary<1,3>::
-get_new_point_on_face (const Triangulation<1,3>::face_iterator &) const
-{
- Assert (false, ExcImpossibleInDim(1));
- return Point<3>();
-}
-
-
-template <>
-void
-Boundary<1,3>::
-get_intermediate_points_on_face (const Triangulation<1,3>::face_iterator &,
- std::vector<Point<3> > &) const
-{
- Assert (false, ExcImpossibleInDim(1));
-}
-
-
-
-
-template <int dim, int spacedim>
-Tensor<1,spacedim>
-Boundary<dim, spacedim>::
-normal_vector (const typename Triangulation<dim, spacedim>::face_iterator &,
- const Point<spacedim> &) const
-{
- Assert (false, ExcPureFunctionCalled());
- return Tensor<1,spacedim>();
-}
-
-
-
-template <int dim, int spacedim>
-void
-Boundary<dim, spacedim>::
-get_normals_at_vertices (const typename Triangulation<dim, spacedim>::face_iterator &,
- FaceVertexNormals &) const
-{
- Assert (false, ExcPureFunctionCalled());
-}
-
-
-
-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() == 0)
- {
- Threads::Mutex::ScopedLock lock(mutex);
- if (points.size() <= n_intermediate_points)
- points.resize(n_intermediate_points+1);
-
- // another thread might have created points in the meantime
- if (points[n_intermediate_points].get() == 0)
- {
- std_cxx1x::shared_ptr<QGaussLobatto<1> >
- quadrature (new QGaussLobatto<1>(n_intermediate_points+2));
- points[n_intermediate_points] = quadrature;
- }
- }
- return points[n_intermediate_points]->get_points();
-}
-
-
-
-
/* -------------------------- StraightBoundary --------------------- */
StraightBoundary<dim, spacedim>::StraightBoundary ()
{}
-
-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;
-}
-
-
-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 &) const
-{
- Assert (false, ExcImpossibleInDim(dim));
- return Point<spacedim>();
-}
-
-
-template <>
-Point<3>
-StraightBoundary<2,3>::
-get_new_point_on_quad (const Triangulation<2,3>::quad_iterator &quad) const
-{
- return compute_new_point_on_quad<2> (quad);
-}
-
-
-
-template <>
-Point<3>
-StraightBoundary<3>::
-get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
-{
- return compute_new_point_on_quad<3> (quad);
-}
-
-
-
-template <>
-void
-StraightBoundary<1>::
-get_intermediate_points_on_line (const Triangulation<1>::line_iterator &,
- std::vector<Point<1> > &) const
-{
- Assert(false, ExcImpossibleInDim(1));
-}
-
-template <>
-void
-StraightBoundary<1, 2>::
-get_intermediate_points_on_line (const Triangulation<1, 2>::line_iterator &line,
- std::vector<Point<2> > &points) const
-{
- const unsigned int spacedim = 2;
- const unsigned int n=points.size();
- Assert(n>0, ExcInternalError());
-
- // Use interior points of QGaussLobatto quadrature formula support points
- // for consistency with MappingQ
- const std::vector<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[i+1][0];
- points[i] = (1-x)*vertices[0] + x*vertices[1];
- }
-}
-
-
-
-
-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 <>
-Tensor<1,1>
-StraightBoundary<1,1>::
-normal_vector (const Triangulation<1,1>::face_iterator &,
- const Point<1> &) const
-{
- Assert (false, ExcNotImplemented());
- return Tensor<1,1>();
-}
-
-
-template <>
-Tensor<1,2>
-StraightBoundary<1,2>::
-normal_vector (const Triangulation<1,2>::face_iterator &,
- const Point<2> &) const
-{
- Assert (false, ExcNotImplemented());
- return Tensor<1,2>();
-}
-
-
-template <>
-Tensor<1,3>
-StraightBoundary<1,3>::
-normal_vector (const Triangulation<1,3>::face_iterator &,
- const Point<3> &) const
-{
- Assert (false, ExcNotImplemented());
- return Tensor<1,3>();
-}
-
-
-namespace internal
-{
- namespace
- {
- /**
- * Compute the normalized cross product of a set of dim-1 basis
- * vectors.
- */
- Tensor<1,2>
- normalized_alternating_product (const Tensor<1,2> (&basis_vectors)[1])
- {
- Tensor<1,2> tmp;
- cross_product (tmp, basis_vectors[0]);
- return tmp/tmp.norm();
- }
-
-
-
- Tensor<1,3>
- normalized_alternating_product (const Tensor<1,3> ( &)[1])
- {
- // we get here from StraightBoundary<2,3>::normal_vector, but
- // the implementation below is bogus for this case anyway
- // (see the assert at the beginning of that function).
- Assert (false, ExcNotImplemented());
- return Tensor<1,3>();
- }
-
-
-
- Tensor<1,3>
- normalized_alternating_product (const Tensor<1,3> (&basis_vectors)[2])
- {
- Tensor<1,3> tmp;
- cross_product (tmp, basis_vectors[0], basis_vectors[1]);
- return tmp/tmp.norm();
- }
-
- }
-}
-
-
-template <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
-{
- // I don't think the implementation below will work when dim!=spacedim;
- // in fact, I believe that we don't even have enough information here,
- // because we would need to know not only about the tangent vectors
- // of the face, but also of the cell, to compute the normal vector.
- // Someone will have to think about this some more.
- Assert (dim == spacedim, ExcNotImplemented());
-
- // in order to find out what the normal vector is, we first need to
- // find the reference coordinates of the point p on the given face,
- // or at least the reference coordinates of the closest point on the
- // face
- //
- // in other words, we need to find a point xi so that f(xi)=||F(xi)-p||^2->min
- // where F(xi) is the mapping. this algorithm is implemented in
- // MappingQ1<dim,spacedim>::transform_real_to_unit_cell but only for cells,
- // while we need it for faces here. it's also implemented in somewhat
- // more generality there using the machinery of the MappingQ1 class
- // while we really only need it for a specific case here
- //
- // in any case, the iteration we use here is a Gauss-Newton's iteration with
- // xi^{n+1} = xi^n - H(xi^n)^{-1} J(xi^n)
- // where
- // J(xi) = (grad F(xi))^T (F(xi)-p)
- // and
- // H(xi) = [grad F(xi)]^T [grad F(xi)]
- // In all this,
- // F(xi) = sum_v vertex[v] phi_v(xi)
- // We get the shape functions phi_v from an object of type FE_Q<dim-1>(1)
-
- // we start with the point xi=1/2, xi=(1/2,1/2), ...
- const unsigned int facedim = dim-1;
-
- Point<facedim> xi;
- for (unsigned int i=0; i<facedim; ++i)
- xi[i] = 1./2;
-
- FE_Q<facedim> linear_fe(1);
-
- const double eps = 1e-12;
- Tensor<1,spacedim> grad_F[facedim];
- while (true)
- {
- Point<spacedim> F;
- for (unsigned int v=0; v<GeometryInfo<facedim>::vertices_per_cell; ++v)
- F += face->vertex(v) * linear_fe.shape_value(v, xi);
-
- for (unsigned int i=0; i<facedim; ++i)
- {
- grad_F[i] = 0;
- for (unsigned int v=0; v<GeometryInfo<facedim>::vertices_per_cell; ++v)
- grad_F[i] += face->vertex(v) * linear_fe.shape_grad(v, xi)[i];
- }
-
- Tensor<1,facedim> J;
- for (unsigned int i=0; i<facedim; ++i)
- for (unsigned int j=0; j<spacedim; ++j)
- J[i] += grad_F[i][j] * (F-p)[j];
-
- Tensor<2,facedim> H;
- for (unsigned int i=0; i<facedim; ++i)
- for (unsigned int j=0; j<facedim; ++j)
- for (unsigned int k=0; k<spacedim; ++k)
- H[i][j] += grad_F[i][k] * grad_F[j][k];
-
- const Point<facedim> delta_xi = -invert(H) * J;
- xi += delta_xi;
-
- if (delta_xi.norm() < eps)
- break;
- }
-
- // so now we have the reference coordinates xi of the point p.
- // we then have to compute the normal vector, which we can do
- // by taking the (normalize) alternating product of all the tangent
- // vectors given by grad_F
- return internal::normalized_alternating_product(grad_F);
-}
-
-
-
-template <>
-void
-StraightBoundary<1>::
-get_normals_at_vertices (const Triangulation<1>::face_iterator &,
- Boundary<1,1>::FaceVertexNormals &) const
-{
- Assert (false, ExcImpossibleInDim(1));
-}
-
-template <>
-void
-StraightBoundary<1,2>::
-get_normals_at_vertices (const Triangulation<1,2>::face_iterator &,
- Boundary<1,2>::FaceVertexNormals &) const
-{
- Assert (false, ExcNotImplemented());
-}
-
-
-template <>
-void
-StraightBoundary<1,3>::
-get_normals_at_vertices (const Triangulation<1,3>::face_iterator &,
- Boundary<1,3>::FaceVertexNormals &) const
-{
- Assert (false, ExcNotImplemented());
-}
-
-
-
-template <>
-void
-StraightBoundary<2>::
-get_normals_at_vertices (const Triangulation<2>::face_iterator &face,
- Boundary<2,2>::FaceVertexNormals &face_vertex_normals) const
-{
- const Tensor<1,2> tangent = face->vertex(1) - face->vertex(0);
- for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_face; ++vertex)
- // compute normals from tangent
- face_vertex_normals[vertex] = Point<2>(tangent[1],
- -tangent[0]);
-}
-
-template <>
-void
-StraightBoundary<2,3>::
-get_normals_at_vertices (const Triangulation<2,3>::face_iterator &face,
- Boundary<2,3>::FaceVertexNormals &face_vertex_normals) const
-{
- const Tensor<1,3> tangent = face->vertex(1) - face->vertex(0);
- for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_face; ++vertex)
- // compute normals from tangent
- face_vertex_normals[vertex] = Point<3>(tangent[1],
- -tangent[0],0);
- Assert(false, ExcNotImplemented());
-}
-
-
-
-
-template <>
-void
-StraightBoundary<3>::
-get_normals_at_vertices (const Triangulation<3>::face_iterator &face,
- Boundary<3,3>::FaceVertexNormals &face_vertex_normals) const
-{
- const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face;
-
- static const unsigned int neighboring_vertices[4][2]=
- { {1,2},{3,0},{0,3},{2,1}};
- for (unsigned int vertex=0; vertex<vertices_per_face; ++vertex)
- {
- // first define the two tangent
- // vectors at the vertex by
- // using the two lines
- // radiating away from this
- // vertex
- const Tensor<1,3> tangents[2]
- = { face->vertex(neighboring_vertices[vertex][0])
- - face->vertex(vertex),
- face->vertex(neighboring_vertices[vertex][1])
- - face->vertex(vertex)
- };
-
- // then compute the normal by
- // taking the cross
- // product. since the normal is
- // not required to be
- // normalized, no problem here
- cross_product (face_vertex_normals[vertex],
- tangents[0], tangents[1]);
- };
-}
-
-
-
-template <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,
- internal::int2type<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 (tmp,
- 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 Point<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: */internal::int2type<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: */internal::int2type<2>)
- {
- return y;
- }
-}
-
-
-
-
-
-template <>
-Point<3>
-StraightBoundary<1,3>::
-project_to_surface (const Triangulation<1, 3>::quad_iterator &,
- const Point<3> &y) const
-{
- return y;
-}
-
-//TODO[SP]: This is just a horrible way out to make it compile in codim 2.
-template <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 */internal::int2type<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"
}
-
-// 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 Point<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 (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>::
}
-
-
-// template<>
-// Point<3>
-// CylinderBoundary<3>::
-// get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
-// {
-// const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
-
-// // same algorithm as above
-// const unsigned int spacedim = 3;
-
-// const Point<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 (vector_from_axis / vector_from_axis.norm() * radius +
-// ((middle-point_on_axis) * direction) * direction +
-// point_on_axis);
-// }
-
-// template<>
-// Point<3>
-// CylinderBoundary<2,3>::
-// get_new_point_on_quad (const Triangulation<2,3>::quad_iterator &quad) const
-// {
-// const Point<3> middle = StraightBoundary<2,3>::get_new_point_on_quad (quad);
-
-// // same algorithm as above
-// const unsigned int spacedim = 3;
-// const Point<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 (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 Point<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] = (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>
Point<spacedim>
CylinderBoundary<dim,spacedim>::
}
-
-// 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 Point<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.square ();
-// 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 Point<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.square ();
-// 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>::
}
-
-// template <>
-// Point<3>
-// ConeBoundary<3>::
-// get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
-// {
-// const int dim = 3;
-
-// const Point<dim> axis = x_1 - x_0;
-// // Compute the middle point of the quad.
-// const Point<dim> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
-// // Same algorithm as above: To project it on the boundary of the cone we
-// // first compute the orthogonal projection of the middle point onto the axis
-// // of the cone.
-// const double c = (middle - x_0) * axis / axis.square ();
-// const Point<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>
Point<dim>
ConeBoundary<dim>::
compute_radius_automatically(false)
{}
-// template <int dim, int spacedim>
-// Point<spacedim>
-// HyperBallBoundary<dim,spacedim>::get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
-// const std::vector<double> &weights) const
-// {
-// Assert(surrounding_points.size() == weights.size(),
-// ExcDimensionMismatch(surrounding_points.size(), weights.size()));
-
-// double radius = 0;
-// for(unsigned int i=0; i<surrounding_points.size(); ++i)
-// radius += weights[i]*(surrounding_points[i]-center).norm();
-
-// Point<spacedim> p = FlatManifold<spacedim>::get_new_point(surrounding_points, weights);
-// p = p-center;
-// p = p/p.norm()*radius+center;
-// return p;
-// }
-
-
-// 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)
-// {
-// const Point<spacedim> vertex_relative = line->vertex(0) - center;
-// r = std::sqrt(vertex_relative.square());
-// }
-// else
-// r=radius;
-// // project to boundary
-// middle *= r / std::sqrt(middle.square());
-// middle += center;
-// return middle;
-// }
-
template <int dim, int spacedim>
Point<spacedim>
}
-
-// 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)
-// {
-// const Point<spacedim> vertex_relative = quad->vertex(0) - center;
-// r = std::sqrt(vertex_relative.square());
-// }
-// else
-// r=radius;
-// // project to boundary
-// middle *= r / std::sqrt(middle.square());
-
-// middle += center;
-// return middle;
-// }
-
-
-
-// template <>
-// void
-// HyperBallBoundary<1>::get_intermediate_points_on_line (
-// const Triangulation<1>::line_iterator &,
-// std::vector<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 Point<spacedim> v0=p0-center,
-// v1=p1-center;
-// const double length=std::sqrt((v1-v0).square());
-
-// double eps=1e-12;
-// double r=0;
-// if (compute_radius_automatically)
-// {
-// const Point<spacedim> vertex_relative = p0 - center;
-// r = std::sqrt(vertex_relative.square());
-// }
-// else
-// r=radius;
-
-
-// const double r2=r*r;
-// Assert(std::fabs(v0.square()-r2)<eps*r2, ExcInternalError());
-// Assert(std::fabs(v1.square()-r2)<eps*r2, ExcInternalError());
-
-// const double alpha=std::acos((v0*v1)/std::sqrt(v0.square()*v1.square()));
-// const Point<spacedim> pm=0.5*(v0+v1);
-
-// const double h=std::sqrt(pm.square());
-
-// // n even: m=n/2,
-// // n odd: m=(n-1)/2
-// const std::vector<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]=pm+d/length*(v1-v0);
-// points[n-1-i]=pm-d/length*(v1-v0);
-// }
-
-// if ((n+1)%2==0)
-// // if the number of parts is even insert the midpoint
-// points[(n-1)/2]=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>
Point<spacedim>
HyperBallBoundary<dim,spacedim>::
{}
-// template <int dim>
-// Point<dim>
-// HalfHyperBallBoundary<dim>::get_new_point_on_line(const typename Triangulation<dim>::line_iterator &line) const
-// {
-// Assert(false, ExcInternalError());
-// return Point<dim>();
-// }
-
-
template <int dim>
Point<dim>
HalfHyperBallBoundary<dim>::get_new_point(const std::vector<Point<dim> > &surrounding_points,
}
-
-// 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>
Point<dim>
HalfHyperBallBoundary<dim>::
}
-
-// 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 Point<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>
Point<dim>
HalfHyperShellBoundary<dim>::
}
-
-// 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>
}
-
-
-// template<>
-// void
-// TorusBoundary<2,3>::
-// get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator &line,
-// std::vector< Point< 3 > > &points) const
-// {
-// //Almost the same implementation as StraightBoundary<2,3>
-// unsigned int npoints=points.size();
-// if (npoints==0) return;
-
-// Point<2> p[2];
-
-// for (unsigned int i=0; i<2; i++)
-// p[i]=get_surf_coord(line->vertex(i));
-
-// unsigned int offset[2];
-// offset[0]=0;
-// offset[1]=0;
-
-// //Take care for periodic conditions & negative angles, see
-// //get_new_point_on_line() above. Because we dont have a symmetric
-// //interpolation (just the middle) we need to add 2*Pi to each almost zero
-// //and negative angles.
-// for (unsigned int i=0; i<2; i++)
-// for (unsigned int j=1; j<2; j++)
-// {
-// if (std::abs(p[0](i)-p[j](i))> numbers::PI)
-// {
-// offset[i]++;
-// break;
-// }
-// }
-
-// for (unsigned int i=0; i<2; i++)
-// for (unsigned int j=0; j<2; j++)
-// if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
-// p[j](i)+=2*numbers::PI*offset[i];
-
-
-// Point<2> target;
-// const std::vector<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 dont have a symmetric
-// //interpolation (just the middle) we need to add 2*Pi to each almost zero
-// //and negative angles.
-// for (unsigned int i=0; i<2; i++)
-// for (unsigned int j=1; j<4; j++)
-// {
-// if (std::abs(p[0](i)-p[j](i))> numbers::PI)
-// {
-// offset[i]++;
-// break;
-// }
-// }
-
-// for (unsigned int i=0; i<2; i++)
-// for (unsigned int j=0; j<4; j++)
-// if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
-// p[j](i)+=2*numbers::PI*offset[i];
-
-// const std::vector<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"
projects / dealii-svn.git / commitdiff