--- /dev/null
+// ---------------------------------------------------------------------
+// $Id: tria_boundary.h 30130 2013-07-23 13:01:18Z heltai $
+//
+// Copyright (C) 2013 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 __deal2__manifold_h
+#define __deal2__manifold_h
+
+
+/*---------------------------- boundary-function.h ---------------------------*/
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/subscriptor.h>
+#include <deal.II/base/point.h>
+#include <deal.II/grid/tria.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+template <int dim, int space_dim> class Triangulation;
+
+
+
+/**
+ * This class is used to represent a geometry to a triangulation.
+ * When a triangulation creates a new vertex in the
+ * domain, it determines the new vertex' coordinates through the
+ * following code (here in two dimensions):
+ * @code
+ * ...
+ * Point<2> new_vertex = manifold.get_new_point(surrounding_points, weights);
+ * ...
+ * @endcode
+ * @p surrounding_points is a vector containing all points which surround
+ * the vertex to be created, while weights contains a list of the weights to use
+ * when computing the new location.
+ *
+ * For example, for a flat manifold, the new_vertex is computed as:
+ * @code
+ * ...
+ * Point<dim> new_vertex;
+ *
+ * for(unsigned int i=0; i<surrounding_points.size(); ++i)
+ * new_vertex += weights[i]*surrounding_points[i];
+ * ...
+ * @endcode
+ *
+ * Each different manifold description will override the function above,
+ * by assuming that the surrounding points are on the manifold, and computing
+ * whatever is meaninful out of the weighted average of the surrounding_points.
+ *
+ * @ingroup manifold
+ * @author Luca Heltai, 2013
+ */
+template <int spacedim>
+class Manifold : public Subscriptor
+{
+public:
+
+ /**
+ * Destructor. Does nothing here, but
+ * needs to be declared to make it
+ * virtual.
+ */
+ virtual ~Manifold ();
+
+ /**
+ * Return the point which shall
+ * become the new vertex surrounded
+ * by the given points.
+ */
+ virtual
+ Point<spacedim>
+ get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const = 0;
+
+
+ /**
+ * 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_manifold (const Point<spacedim> &candidate) const;
+};
+
+
+
+/**
+ * Specialization of Manifold<spacedim>, which places the new point
+ * right into the middle of the given points. The middle is defined
+ * as the weighted mean of the points.
+ *
+ * @ingroup manifold
+ *
+ * @author Luca Heltai, 2013
+ */
+template <int spacedim>
+class FlatManifold: public Manifold<spacedim>
+{
+public:
+ /**
+ * Default constructor. Some
+ * compilers require this for
+ * some reasons.
+ */
+ FlatManifold ();
+
+ /**
+ * 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(const std::vector<Point<spacedim> > &surrounding_points,
+ const std::vector<double> &weights) const;
+};
+
+
+
+/* -------------- declaration of explicit specializations ------------- */
+
+#ifndef DOXYGEN
+
+#endif // DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id: tria_boundary_lib.h 30130 2013-07-23 13:01:18Z heltai $
+//
+// Copyright (C) 1999 - 2013 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 __deal2__manifold_lib_h
+#define __deal2__manifold_lib_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/grid/manifold.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.
+// *
+// * @ingroup boundary
+// *
+// * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007
+// */
+// template <int dim, int spacedim = dim>
+// class 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);
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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
+// *
+// * @author Markus Bürg, 2009
+// */
+// template <int dim>
+// class 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);
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// 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
+// *
+// * @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001
+// */
+// template <int dim, int spacedim=dim>
+// class HyperBallBoundary : public StraightBoundary<dim,spacedim>
+// {
+// public:
+// /**
+// * Constructor
+// */
+// HyperBallBoundary (const Point<spacedim> p = Point<spacedim>(),
+// const double radius = 1.0);
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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
+// *
+// * @author Wolfgang Bangerth, 1999, 2001
+// */
+// template <int dim>
+// class HalfHyperBallBoundary : public HyperBallBoundary<dim>
+// {
+// public:
+// /**
+// * Constructor
+// */
+// HalfHyperBallBoundary (const Point<dim> p = Point<dim>(),
+// const double radius = 1.0);
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+// };
+//
+//
+//
+// /**
+// * 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
+// *
+// * @author Wolfgang Bangerth, 1999
+// */
+// template <int dim>
+// class HyperShellBoundary : public HyperBallBoundary<dim>
+// {
+// public:
+// /**
+// * Constructor. The center of the
+// * spheres defaults to the
+// * origin.
+// *
+// * Calls 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>());
+// };
+//
+//
+//
+// /**
+// * 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
+// *
+// * @author Wolfgang Bangerth, 2000, 2009
+// */
+// template <int dim>
+// class 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);
+//
+// /**
+// * Construct a new point on a line.
+// */
+// virtual Point<dim>
+// get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
+//
+// /**
+// * Construct a new point on a quad.
+// */
+// virtual Point<dim>
+// get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// /**
+// * 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;
+//
+// 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.
+// */
+// template <int dim, int spacedim>
+// class 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);
+//
+// //Boundary Refinenment 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;
+//
+// /**
+// * Construct a new point on a quad.
+// */
+// virtual Point<spacedim>
+// get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
+//
+// /**
+// * 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 ;
+//
+// /**
+// * 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 ;
+//
+// /**
+// * 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 ;
+//
+// 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
template <int dim, int spacedim> class Boundary;
template <int dim, int spacedim> class StraightBoundary;
+template <int spacedim> class Manifold;
+template <int spacedim> class FlatManifold;
template <int, int, int> class TriaAccessor;
template <int spacedim> class TriaAccessor<0,1,spacedim>;
+
namespace internal
{
namespace Triangulation
*/
static const StraightBoundary<dim,spacedim> straight_boundary;
+ /**
+ * Default manifold object. This is used
+ * for those objects for which no
+ * manifold description has been explicitly
+ * set using set_manifold().
+ */
+ static const FlatManifold<spacedim> flat_manifold;
+
+
/**
* Declare some symbolic names
* for mesh smoothing
void set_boundary (const types::manifold_id number,
const Boundary<dim,spacedim> &boundary_object);
- /**
+
+ void set_manifold (const types::manifold_id number,
+ const Manifold<spacedim> &manifold_object);
+
+
+
+/**
* Reset those parts of the boundary with
* the given number to use a straight
* boundary approximation. This is the
*/
void set_boundary (const types::manifold_id number);
+ void set_manifold (const types::manifold_id number);
+
/**
* Return a constant reference to
* a boundary object used for
*/
const Boundary<dim,spacedim> &get_boundary (const types::manifold_id number) const;
+
+ const Manifold<spacedim> &get_manifold (const types::manifold_id number) const;
+
/**
* Returns a vector containing
* all boundary indicators
*/
std::map<types::manifold_id, SmartPointer<const Boundary<dim, spacedim> , Triangulation<dim, spacedim> > > boundary;
+
+ /**
+ * Collection of manifold
+ * objects. We store only
+ * objects, which are not
+ * of type StraightBoundary.
+ */
+ std::map<types::manifold_id, SmartPointer<const Manifold<spacedim> , Triangulation<dim, spacedim> > > manifold;
+
+
/**
* Flag indicating whether
* anisotropic refinement took
template <typename Accessor> class TriaActiveIterator;
template <int dim, int spacedim> class Boundary;
+template <int spacedim> class Manifold;
namespace internal
* for the boundary object.
*/
const Boundary<dim,spacedim> &get_boundary () const;
+
+ /**
+ * Return a constant reference to a
+ * manifold object used for this
+ * object. This function is a shortcut to
+ * retrieving the manifold indicator
+ * using manifold_indicator() and then
+ * asking the
+ * Triangulation::get_manifold() function
+ * for the manifold object.
+ */
+ const Manifold<spacedim> &get_manifold () const;
/**
* @}
}
+
+template <int structdim, int dim, int spacedim>
+const Manifold<spacedim> &
+TriaAccessor<structdim, dim, spacedim>::get_manifold () const
+{
+
+ Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed());
+ // Get the default (manifold_id)
+ types::manifold_id mi = this->objects().manifold_id[this->present_index];
+
+ return this->tria->get_manifold(mi);
+}
+
+
template <int structdim, int dim, int spacedim>
types::manifold_id
TriaAccessor<structdim, dim, spacedim>::manifold_id () const
grid_reordering.cc
grid_tools.cc
intergrid_map.cc
+ manifold.cc
+ manifold_lib.cc
persistent_tria.cc
tria_accessor.cc
tria_boundary.cc
grid_refinement.inst.in
grid_tools.inst.in
intergrid_map.inst.in
+ manifold.inst.in
tria_accessor.inst.in
tria_boundary.inst.in
tria_boundary_lib.inst.in
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id: tria_boundary.cc 30130 2013-07-23 13:01:18Z heltai $
+//
+// Copyright (C) 1998 - 2013 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/point.h>
+#include <deal.II/grid/manifold.h>
+#include <cmath>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+/* -------------------------- Manifold --------------------- */
+
+
+template <int spacedim>
+Manifold<spacedim>::~Manifold ()
+{}
+
+
+
+
+template <int spacedim>
+Point<spacedim>
+Manifold<spacedim>::project_to_manifold (const Point<spacedim> &candidate) const
+{
+ Assert (false, ExcPureFunctionCalled());
+ return Point<spacedim>();
+}
+
+
+/* -------------------------- FlatManifold --------------------- */
+
+
+template <int spacedim>
+FlatManifold<spacedim>::FlatManifold ()
+{}
+
+
+template <int spacedim>
+Point<spacedim>
+FlatManifold<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()));
+
+ Point<spacedim> p;
+ for(unsigned int i=0; i<surrounding_points.size(); ++i)
+ p += surrounding_points[i]*weights[i];
+
+ return p;
+}
+
+// explicit instantiations
+#include "manifold.inst"
+
+DEAL_II_NAMESPACE_CLOSE
+
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id: tria_boundary.inst.in 30130 2013-07-23 13:01:18Z heltai $
+//
+// Copyright (C) 1998 - 2013 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_space_dimension : SPACE_DIMENSIONS)
+ {
+ template class Manifold<deal_II_space_dimension>;
+ template class FlatManifold<deal_II_space_dimension>;
+ }
+
+
+
#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/manifold.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/magic_numbers.h>
#include <deal.II/fe/mapping_q1.h>
DEAL_II_NAMESPACE_OPEN
+/** We create an anonymous namespace to make a helper function which
+ would return default points and weights contained in lines, quads
+ and hexes which are used to query the manifold objects.
+ */
+namespace {
+
+ void get_default_points_and_weights(std::vector<Point<3> > &sp,
+ std::vector<double> &wp,
+ const CellAccessor<3,3>& obj)
+ {
+ const int dim = 3;
+
+ int np = GeometryInfo<dim>::vertices_per_cell+
+ GeometryInfo<dim>::lines_per_cell+
+ GeometryInfo<dim>::faces_per_cell;
+ sp.resize(np);
+ wp.resize(np);
+ unsigned int j=0;
+ for(unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i, ++j)
+ {
+ sp[j] = obj.vertex(i);
+ wp[j] = 1.0/128.0;
+ }
+ for(unsigned int i=0; i<GeometryInfo<dim>::lines_per_cell; ++i, ++j)
+ {
+ sp[j] = obj.line(i)->child(0)->vertex(1);
+ wp[j] = 7.0/192.0;
+ }
+ for(unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i, ++j)
+ {
+ sp[j] = obj.face(i)->isotropic_child(0)->vertex(3);
+ wp[j] = 1.0/12.0;
+ }
+ }
+
+
+ template <typename OBJECT, int spacedim>
+ void get_default_points_and_weights(std::vector<Point<spacedim> > &sp,
+ std::vector<double> &wp,
+ const OBJECT& obj, bool with_laplace = false)
+ {
+ const int dim = OBJECT::AccessorType::structure_dimension;
+ AssertDimension(spacedim, OBJECT::AccessorType::space_dimension);
+ switch(dim)
+ {
+ case 1:
+ sp.resize(2);
+ wp.resize(2);
+ sp[0] = obj->vertex(0); wp[0] = .5;
+ sp[1] = obj->vertex(1); wp[0] = .5;
+ break;
+ case 2:
+ sp.resize(8);
+ wp.resize(8);
+ sp[0] = obj->vertex(0);
+ sp[1] = obj->vertex(1);
+ sp[2] = obj->vertex(2);
+ sp[3] = obj->vertex(3);
+
+ sp[4] = obj->line(0)->has_children() ?
+ obj->line(0)->child(0)->vertex(1) :
+ obj->line(0)->center();
+ sp[5] = obj->line(1)->has_children() ?
+ obj->line(1)->child(0)->vertex(1) :
+ obj->line(1)->center();
+ sp[6] = obj->line(2)->has_children() ?
+ obj->line(2)->child(0)->vertex(1) :
+ obj->line(2)->center();
+ sp[7] = obj->line(3)->has_children() ?
+ obj->line(3)->child(0)->vertex(1) :
+ obj->line(3)->center();
+ if(with_laplace)
+ {
+ std::fill(wp.begin(), wp.begin()+4, 1.0/16.0);
+ std::fill(wp.begin()+4, wp.end(), 3.0/16.0);
+ }
+ else
+ std::fill(wp.begin(), wp.end(), 1.0/8.0);
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ break;
+ }
+ }
+}
+
+
template<int dim>
CellData<dim>::CellData()
{
// not the case, then
// we need to ask a
// boundary object
- if (dim == spacedim)
- {
+ {
// in a first
// step, compute
// the location
// average of the
// 8 surrounding
// vertices
- Point<spacedim> new_point;
- for (unsigned int i=0; i<8; ++i)
- new_point += triangulation.vertices[new_vertices[i]];
- new_point /= 8.0;
-
- triangulation.vertices[next_unused_vertex] = new_point;
+ std::vector<Point<spacedim> > ps(8);
+ std::vector<double> ws(8, 1.0/8);
+ for (unsigned int i=0; i<8; ++i)
+ ps[i] = triangulation.vertices[new_vertices[i]];
+ triangulation.vertices[next_unused_vertex] =
+ cell->get_manifold().get_new_point(ps, ws);
// if the user_flag is set, i.e. if the
// cell is at the boundary, use a
boundary_face=GeometryInfo<dim>::faces_per_cell+1;
}
- if (boundary_face<GeometryInfo<dim>::faces_per_cell)
- // reset the cell's middle vertex
- // to the middle of the straight
- // connection between the new
- // points on this face and on the
- // opposite face
- triangulation.vertices[next_unused_vertex]
- = 0.5*(cell->face(boundary_face)
- ->child(0)->vertex(1)+
- cell->face(GeometryInfo<dim>
+ if (boundary_face<GeometryInfo<dim>::faces_per_cell)
+ {
+ // reset the cell's middle vertex
+ // to the middle of the straight
+ // connection between the new
+ // points on this face and on the
+ // opposite face
+ std::vector<Point<spacedim> > sp(2);
+ std::vector<double> ws(2, .5);
+ sp[0] = cell->face(boundary_face)
+ ->child(0)->vertex(1);
+ sp[1] = cell->face(GeometryInfo<dim>
::opposite_face[boundary_face])
- ->child(0)->vertex(1));
- }
- }
- else
- {
- // if this quad
- // lives in a
- // higher
- // dimensional
- // space then we
- // don't need to
- // worry if it is
- // at the
- // boundary of
- // the manifold
- // -- we always
- // have to use
- // the boundary
- // object anyway;
- // so ignore
- // whether the
- // user flag is
- // set or not
- cell->clear_user_flag();
-
-
- // new vertex is
- // placed on the
- // surface according
- // to the information
- // stored in the
- // boundary class
- triangulation.vertices[next_unused_vertex] =
- cell->get_boundary().get_new_point_on_quad (cell);
- }
+ ->child(0)->vertex(1);
+
+ triangulation.vertices[next_unused_vertex]
+ = cell->get_manifold().get_new_point(ps, ws);
+ }
+ }
+ }
}
const bool /*check_for_distorted_cells*/)
{
const unsigned int dim = 1;
+
+ std::vector<Point<spacedim> > sp(2);
+ std::vector<double> wp(2);
// check whether a new level is
// needed we have to check for this
// midpoint; otherwise we
// have to ask the manifold
// object
- if (dim == spacedim)
- triangulation.vertices[next_unused_vertex] =
- (cell->vertex(0) + cell->vertex(1)) / 2;
- else
- triangulation.vertices[next_unused_vertex] =
- cell->get_boundary().get_new_point_on_line(cell);
+ get_default_points_and_weights(sp, wp, cell);
+
+ triangulation.vertices[next_unused_vertex] =
+ cell->get_manifold().get_new_point(sp, wp);
triangulation.vertices_used[next_unused_vertex] = true;
{
const unsigned int dim = 2;
+ // Helper vectors of points
+ // and weights, to compute
+ // the locations of the newly
+ // created vertices (through
+ // the manifold objects)
+ std::vector<Point<spacedim> > sp;
+ std::vector<double> wp;
+
+
// check whether a new level is
// needed we have to check for this
// on the highest level only (on
// Ask the manifold
// where to put the
// new line
- triangulation.vertices[next_unused_vertex]
- = line->get_boundary().get_new_point_on_line(line);
+ get_default_points_and_weights(sp, wp, line);
+ triangulation.vertices[next_unused_vertex] =
+ line->get_manifold().get_new_point(sp, wp);
// if (spacedim == dim)
{
const unsigned int dim = 3;
+ std::vector<Point<spacedim> > sp;
+ std::vector<double> wp;
+
// this function probably
// also works for spacedim>3
// but it isn't tested. it
ExcTooFewVerticesAllocated());
triangulation.vertices_used[next_unused_vertex] = true;
- triangulation.vertices[next_unused_vertex]
- = line->get_boundary().get_new_point_on_line(line);
+ get_default_points_and_weights(sp, wp, line);
+ triangulation.vertices[next_unused_vertex] =
+ line->get_manifold().get_new_point(sp, wp);
+
// if (line->at_boundary())
// triangulation.vertices[next_unused_vertex]
// the domain, so we need not
// care about boundary quads
// here
- triangulation.vertices[next_unused_vertex]
- = middle_line->get_boundary().get_new_point_on_line(middle_line);
-// = (middle_line->vertex(0) + middle_line->vertex(1)) / 2;
+ get_default_points_and_weights(sp, wp, middle_line);
+ triangulation.vertices[next_unused_vertex] =
+ middle_line->get_manifold().get_new_point(sp, wp);
+
triangulation.vertices_used[next_unused_vertex] = true;
// now search a slot for the two
// isotropic refinement
Assert(quad_ref_case==RefinementCase<dim-1>::no_refinement, ExcInternalError());
- // set the middle vertex
- // appropriately
- if (quad->at_boundary() || quad->manifold_id()!= numbers::invalid_manifold_id)
- triangulation.vertices[next_unused_vertex]
- = quad->get_boundary().get_new_point_on_quad(quad);
-// = triangulation.get_boundary(quad->boundary_indicator()).get_new_point_on_quad (quad);
- else
- // it might be that the
- // quad itself is not
- // at the boundary, but
- // that one of its lines
- // actually is. in this
- // case, the newly
- // created vertices at
- // the centers of the
- // lines are not
- // necessarily the mean
- // values of the
- // adjacent vertices,
- // so do not compute
- // the new vertex as
- // the mean value of
- // the 4 vertices of
- // the face, but rather
- // as a weighted mean
- // value of the 8
- // vertices which we
- // already have (the
- // four old ones, and
- // the four ones
- // inserted as middle
- // points for the four
- // lines). summing up
- // some more points is
- // generally cheaper
- // than first asking
- // whether one of the
- // lines is at the
- // boundary
- //
- // note that the exact
- // weights are chosen
- // such as to minimize
- // the distortion of
- // the four new quads
- // from the optimal
- // shape; their
- // derivation and
- // values is copied
- // over from the
- // @p{MappingQ::set_laplace_on_vector}
- // function
- triangulation.vertices[next_unused_vertex]
- = (quad->vertex(0) + quad->vertex(1) +
- quad->vertex(2) + quad->vertex(3) +
- 3*(quad->line(0)->child(0)->vertex(1) +
- quad->line(1)->child(0)->vertex(1) +
- quad->line(2)->child(0)->vertex(1) +
- quad->line(3)->child(0)->vertex(1)) ) / 16;
+ // // set the middle vertex
+ // // appropriately
+ // get_default_points_and_weights(sp, wp, middle_line);
+ // triangulation.vertices[next_unused_vertex] =
+ // middle_line->get_manifold().get_new_point(sp, wp);
+
+ if (quad->at_boundary() || quad->manifold_id()!= numbers::invalid_manifold_id)
+ {
+ get_default_points_and_weights(sp, wp, quad);
+ triangulation.vertices[next_unused_vertex]
+ = quad->get_manifold().get_new_point(sp, wp);
+ }
+ else
+ {
+
+ // it might be that the quad itself is not at
+ // the boundary, but that one of its lines
+ // actually is. in this case, the newly created
+ // vertices at the centers of the lines are not
+ // necessarily the mean values of the adjacent
+ // vertices, so do not compute the new vertex as
+ // the mean value of the 4 vertices of the face,
+ // but rather as a weighted mean value of the 8
+ // vertices which we already have (the four old
+ // ones, and the four ones inserted as middle
+ // points for the four lines). summing up some
+ // more points is generally cheaper than first
+ // asking whether one of the lines is at the
+ // boundary
+ //
+ // note that the exact weights are chosen such
+ // as to minimize the distortion of the four new
+ // quads from the optimal shape; their
+ // derivation and values is copied over from the
+ // @p{MappingQ::set_laplace_on_vector} function
+ get_default_points_and_weights(sp, wp, quad, true);
+ triangulation.vertices[next_unused_vertex]
+ = quad->get_manifold().get_new_point(sp, wp);
+
+ // triangulation.vertices[next_unused_vertex]
+ // = (quad->vertex(0) + quad->vertex(1) +
+ // quad->vertex(2) + quad->vertex(3) +
+ // 3*(quad->line(0)->child(0)->vertex(1) +
+ // quad->line(1)->child(0)->vertex(1) +
+ // quad->line(2)->child(0)->vertex(1) +
+ // quad->line(3)->child(0)->vertex(1)) ) / 16;
+ }
triangulation.vertices_used[next_unused_vertex] = true;
// vertex at the center of
// the quad is weighted (see
// above)
- triangulation.vertices[next_unused_vertex] = Point<dim>();
- // first add corners of hex
- for (unsigned int vertex=0;
- vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
- triangulation.vertices[next_unused_vertex] += hex->vertex(vertex) / 128;
- // now add center of lines
- for (unsigned int line=0;
- line<GeometryInfo<dim>::lines_per_cell; ++line)
- triangulation.vertices[next_unused_vertex] += hex->line(line)->child(0)->vertex(1) *
- 7./192.;
- // finally add centers of
- // faces. note that vertex 3
- // of child 0 is an invariant
- // with respect to the face
- // orientation, flip and
- // rotation
- for (unsigned int face=0;
- face<GeometryInfo<dim>::faces_per_cell; ++face)
- triangulation.vertices[next_unused_vertex] += hex->face(face)->isotropic_child(0)->vertex(3) *
- 1./12.;
+ get_default_points_and_weights(sp, wp, hex);
+ triangulation.vertices[next_unused_vertex] =
+ hex->get_manifold().get_new_point(sp, wp);
+
+ // // first add corners of hex
+ // for (unsigned int vertex=0;
+ // vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
+ // triangulation.vertices[next_unused_vertex] += hex->vertex(vertex) / 128;
+ // // now add center of lines
+ // for (unsigned int line=0;
+ // line<GeometryInfo<dim>::lines_per_cell; ++line)
+ // triangulation.vertices[next_unused_vertex] += hex->line(line)->child(0)->vertex(1) *
+ // 7./192.;
+ // // finally add centers of
+ // // faces. note that vertex 3
+ // // of child 0 is an invariant
+ // // with respect to the face
+ // // orientation, flip and
+ // // rotation
+ // for (unsigned int face=0;
+ // face<GeometryInfo<dim>::faces_per_cell; ++face)
+ // triangulation.vertices[next_unused_vertex] += hex->face(face)->isotropic_child(0)->vertex(3) *
+ // 1./12.;
// set the data of the
// six lines. first
Triangulation<dim, spacedim>::straight_boundary = StraightBoundary<dim,spacedim>();
+template <int dim, int spacedim>
+const FlatManifold<spacedim>
+Triangulation<dim, spacedim>::flat_manifold = FlatManifold<spacedim>();
+
template <int dim, int spacedim>
const unsigned int
+template <int dim, int spacedim>
+void
+Triangulation<dim, spacedim>::set_manifold (const types::manifold_id m_number,
+ const Manifold<spacedim> &manifold_object)
+{
+ Assert(m_number < numbers::invalid_manifold_id,
+ ExcIndexRange(m_number,0,numbers::invalid_manifold_id));
+
+ manifold[m_number] = &manifold_object;
+}
+
+
template <int dim, int spacedim>
void
Triangulation<dim, spacedim>::set_boundary (const types::manifold_id m_number)
}
+template <int dim, int spacedim>
+void
+Triangulation<dim, spacedim>::set_manifold (const types::manifold_id m_number)
+{
+ Assert(m_number < numbers::invalid_manifold_id,
+ ExcIndexRange(m_number,0,numbers::invalid_manifold_id));
+
+ //delete the entry located at number.
+ manifold.erase(m_number);
+}
+
+
template <int dim, int spacedim>
const Boundary<dim, spacedim> &
}
+template <int dim, int spacedim>
+const Manifold<spacedim> &
+Triangulation<dim, spacedim>::get_manifold (const types::manifold_id m_number) const
+{
+ Assert(m_number < numbers::invalid_manifold_id,
+ ExcIndexRange(m_number,0,numbers::invalid_manifold_id));
+
+ //look, if there is a boundary stored at
+ //boundary_id number.
+ typename std::map<types::manifold_id, SmartPointer<const Manifold<spacedim>, Triangulation<dim, spacedim> > >::const_iterator it
+ = manifold.find(m_number);
+
+ if (it != manifold.end())
+ {
+ //if we have found an entry, return it
+ return *(it->second);
+ }
+ else
+ {
+ //if we have not found an entry
+ //connected with number, we return
+ //flat_manifold
+ return flat_manifold;
+ }
+}
+
+
+
+
template <int dim, int spacedim>
std::vector<types::boundary_id>
Triangulation<dim, spacedim>::get_boundary_indicators () const
template class Triangulation<deal_II_dimension, deal_II_space_dimension>;
#endif
#if deal_II_dimension == deal_II_space_dimension
- template class CellData<deal_II_dimension>;
+ template struct CellData<deal_II_dimension>;
#endif
}