]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added Manifold interface. The branch compiles, no tests have been performed. Document...
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 Aug 2013 15:35:27 +0000 (15:35 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 Aug 2013 15:35:27 +0000 (15:35 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@30358 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/grid/manifold.h [new file with mode: 0644]
deal.II/include/deal.II/grid/manifold_lib.h [new file with mode: 0644]
deal.II/include/deal.II/grid/tria.h
deal.II/include/deal.II/grid/tria_accessor.h
deal.II/include/deal.II/grid/tria_accessor.templates.h
deal.II/source/grid/CMakeLists.txt
deal.II/source/grid/manifold.cc [new file with mode: 0644]
deal.II/source/grid/manifold.inst.in [new file with mode: 0644]
deal.II/source/grid/tria.cc
deal.II/source/grid/tria.inst.in

diff --git a/deal.II/include/deal.II/grid/manifold.h b/deal.II/include/deal.II/grid/manifold.h
new file mode 100644 (file)
index 0000000..8354314
--- /dev/null
@@ -0,0 +1,172 @@
+// ---------------------------------------------------------------------
+// $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
diff --git a/deal.II/include/deal.II/grid/manifold_lib.h b/deal.II/include/deal.II/grid/manifold_lib.h
new file mode 100644 (file)
index 0000000..7bd1cf5
--- /dev/null
@@ -0,0 +1,900 @@
+// ---------------------------------------------------------------------
+// $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&uuml;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> &center = 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> &center = 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
index adb735e58add8405b3d93bfe92d12f71db2e3f94..db9f3669d64bae3f89e7004110f2ddaf1e02a507 100644 (file)
@@ -42,10 +42,13 @@ DEAL_II_NAMESPACE_OPEN
 
 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
@@ -1242,6 +1245,15 @@ public:
    */
   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
@@ -1808,7 +1820,13 @@ public:
   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
@@ -1823,6 +1841,8 @@ public:
    */
   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
@@ -1836,6 +1856,9 @@ public:
    */
   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
@@ -3579,6 +3602,16 @@ private:
    */
   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
index 984374f98f4242365221c2c58b4064beebbbd0d9..6c2dcc55b9ae2dabeefac6898e5d2f82f56f7925 100644 (file)
@@ -37,6 +37,7 @@ template <typename Accessor> class TriaIterator;
 template <typename Accessor> class TriaActiveIterator;
 
 template <int dim, int spacedim> class Boundary;
+template <int spacedim> class Manifold;
 
 
 namespace internal
@@ -1161,6 +1162,18 @@ public:
    * 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;
 
   /**
    * @}
index 6e153eb07523334d41f57308cb6821e1b7214520..e7d0d573256ed4c4cc35e29ebadcd3b75c4d742a 100644 (file)
@@ -1959,6 +1959,20 @@ TriaAccessor<structdim, dim, spacedim>::get_boundary () 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
index 5e13fa938c63e5535b2620728c1b23bbc27ff991..4968dd0fe6dc4566c11e228ca39582591467d77f 100644 (file)
@@ -28,6 +28,8 @@ SET(_src
   grid_reordering.cc
   grid_tools.cc
   intergrid_map.cc
+  manifold.cc
+  manifold_lib.cc
   persistent_tria.cc
   tria_accessor.cc
   tria_boundary.cc
@@ -45,6 +47,7 @@ SET(_inst
   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
diff --git a/deal.II/source/grid/manifold.cc b/deal.II/source/grid/manifold.cc
new file mode 100644 (file)
index 0000000..5135baf
--- /dev/null
@@ -0,0 +1,72 @@
+// ---------------------------------------------------------------------
+// $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
+
diff --git a/deal.II/source/grid/manifold.inst.in b/deal.II/source/grid/manifold.inst.in
new file mode 100644 (file)
index 0000000..2e8b1f8
--- /dev/null
@@ -0,0 +1,26 @@
+// ---------------------------------------------------------------------
+// $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>;
+  }
+
+
+
index a29e247a5edd9d50a1dfece41d7cce5c549ae8d9..9fc1ae83edc1cf78dfcc783a14445b8a4bc6934c 100644 (file)
@@ -26,6 +26,7 @@
 #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()
 {
@@ -3960,8 +4048,7 @@ namespace internal
             // not the case, then
             // we need to ask a
             // boundary object
-            if (dim == spacedim)
-              {
+           {
                 // in a first
                 // step, compute
                 // the location
@@ -3970,12 +4057,12 @@ namespace internal
                 // 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
@@ -4014,52 +4101,26 @@ namespace internal
                             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);
+                     }
+                 }             
+           }
           }
 
 
@@ -4316,6 +4377,9 @@ namespace internal
                           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
@@ -4444,12 +4508,10 @@ namespace internal
                    // 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;
 
@@ -4605,6 +4667,15 @@ namespace internal
       {
         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
@@ -4862,8 +4933,9 @@ namespace internal
                                                   // 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)
@@ -5043,6 +5115,9 @@ namespace internal
       {
         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
@@ -5423,8 +5498,10 @@ namespace internal
                           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]
@@ -6033,9 +6110,10 @@ namespace internal
                             // 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
@@ -6099,65 +6177,53 @@ namespace internal
                     // 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;
 
@@ -8423,26 +8489,29 @@ namespace internal
                       // 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
@@ -9501,6 +9570,10 @@ const StraightBoundary<dim,spacedim>
 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
@@ -9610,6 +9683,18 @@ 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,
+                                            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)
@@ -9622,6 +9707,18 @@ 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> &
@@ -9650,6 +9747,35 @@ Triangulation<dim, spacedim>::get_boundary (const types::manifold_id m_number) c
 }
 
 
+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
index 66e61c2864ffcd7b54043927fcf70ab2a689782a..c1f3faa1039ea1809923cc3d4cd0744c1b7a01f7 100644 (file)
@@ -22,6 +22,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
     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
   }

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.