]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fixed all boundary classes. Removed all references to ancient Boundary class.
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 9 Dec 2013 18:49:53 +0000 (18:49 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 9 Dec 2013 18:49:53 +0000 (18:49 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@31959 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-18/step-18.cc
deal.II/examples/step-47/step-47.cc
deal.II/include/deal.II/grid/manifold.h
deal.II/include/deal.II/grid/tria_boundary_lib.h
deal.II/source/grid/manifold.inst.in
deal.II/source/grid/tria_boundary_lib.cc

index 283419e941c6940fb816f541ee4c9c7953d0326a..f587d740be24c21ced5e2358150636dafbfeeb5c 100644 (file)
@@ -40,7 +40,7 @@
 #include <deal.II/grid/grid_refinement.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/grid/grid_tools.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_accessor.h>
@@ -845,12 +845,12 @@ namespace Step18
     // z-axis. Note that the boundary objects need to live as long as the
     // triangulation does; we can achieve this by making the objects static,
     // which means that they live as long as the program runs:
-    static const CylinderManifold<dim> inner_cylinder (inner_radius, 2);
-    static const CylinderManifold<dim> outer_cylinder (outer_radius, 2);
+    static const CylinderBoundary<dim> inner_cylinder (inner_radius, 2);
+    static const CylinderBoundary<dim> outer_cylinder (outer_radius, 2);
     // We then attach these two objects to the triangulation, and make them
     // correspond to boundary indicators 2 and 3:
-    triangulation.set_manifold (2, inner_cylinder);
-    triangulation.set_manifold (3, outer_cylinder);
+    triangulation.set_boundary (2, inner_cylinder);
+    triangulation.set_boundary (3, outer_cylinder);
 
     // There's one more thing we have to take care of (we should have done so
     // above already, but for didactic reasons it was more appropriate to
@@ -1718,17 +1718,6 @@ namespace Step18
   // to the vertex we consider at present, as they may or may not be ordered
   // in the same order.
   //
-  // This inconvenience could be avoided if finite elements have support
-  // points on vertices (which the one here has; for the concept of support
-  // points, see @ref GlossSupport "support points"). For such a case, one
-  // could construct a custom quadrature rule using
-  // FiniteElement::get_unit_support_points(). The first
-  // <code>GeometryInfo@<dim@>::vertices_per_cell*fe.dofs_per_vertex</code>
-  // quadrature points will then correspond to the vertices of the cell and
-  // are ordered consistent with <code>cell-@>vertex(i)</code>, taking into
-  // account that support points for vector elements will be duplicated
-  // <code>fe.dofs_per_vertex</code> times.
-  //
   // Another point worth explaining about this short function is the way in
   // which the triangulation class exports information about its vertices:
   // through the <code>Triangulation::n_vertices</code> function, it
index 4dfcf9e10d8976b679adbe7b995eaaf025c6a296..87707ba90871e81158ff1da556afe967529f558b 100644 (file)
@@ -33,7 +33,7 @@
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/dofs/dof_tools.h>
 #include <deal.II/fe/fe_q.h>
@@ -1065,8 +1065,8 @@ namespace Step47
             GridGenerator::hyper_ball (triangulation);
             //GridGenerator::hyper_cube (triangulation, -1, 1);
 
-            static const PolarManifold<dim> boundary;
-            triangulation.set_manifold (0, boundary);
+            static const HyperBallBoundary<dim> boundary;
+            triangulation.set_boundary (0, boundary);
 
             triangulation.refine_global (2);
           }
index 36e73a0515ffd8835e36cde8a091edac6f9a8e70..553d9bb0ebe072f54926bf34068867c0ba88612f 100644 (file)
@@ -280,7 +280,7 @@ public:
 
   private:
     FlatManifold<dim> sub_manifold;
-    std::vector<Point<dim> > chart_points;
+    mutable std::vector<Point<dim> > chart_points;
 };
 
 
index 54ad9af979a156e1e9f0f7467548d36b0922c522..a2a6d0ad68617c0c44bf504efc92327e0ab3776a 100644 (file)
@@ -48,7 +48,7 @@ DEAL_II_NAMESPACE_OPEN
  * @author Guido Kanschat, 2001, Wolfgang Bangerth, 2007
  */
 template <int dim, int spacedim = dim>
-class CylinderBoundary : public StraightBoundary<dim,spacedim>
+class CylinderBoundary : public FlatManifold<spacedim>
 {
 public:
   /**
@@ -60,38 +60,16 @@ public:
                     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.
+   * 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 (const std::vector<Point<spacedim> > &surrounding_points,
-//              const std::vector<double> &weights) const;
-
+    
   /**
    * Since this class is derived from FlatManifold, we need to
    * overload only the project_to_manifold function in order for the
@@ -99,66 +77,21 @@ public:
    */
     virtual Point<spacedim>
     project_to_manifold (const Point<spacedim> candidate) 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.
+   * Compute the normal to the surface. 
    */
-  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;
-
+    virtual Point<spacedim>
+    normal_vector (const Point<spacedim> point) 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.
+   * Exception. Thrown by the @p get_radius if the @p
+   * compute_radius_automatically, see below, flag is set true.
    */
   DeclException0 (ExcRadiusNotSet);
 
@@ -182,24 +115,8 @@ protected:
 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.
+   * Given a number for the axis, return a vector that denotes this
+   * direction.
    */
   static Point<spacedim> get_axis_vector (const unsigned int axis);
 };
@@ -240,21 +157,15 @@ private:
  * @author Markus B&uuml;rg, 2009
  */
 template <int dim>
-class ConeBoundary : public StraightBoundary<dim>
+class ConeBoundary : public FlatManifold<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>.
+   * 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,
@@ -262,22 +173,11 @@ public:
                 const Point<dim> x_1);
 
   /**
-   * Return the radius of the
-   * (truncated) cone at given point
+   * 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;
-
+    
   /**
    * Since this class is derived from FlatManifold, we need to
    * overload only the project_to_manifold function in order for the
@@ -288,58 +188,14 @@ public:
 
 
   /**
-   * 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.
+   * Compute the normals to the boundary at the given point
    *
-   * Refer to the general
-   * documentation of this class
-   * and the documentation of the
-   * base class.
+   * 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;
+  Point<dim> 
+  normal_vector (const Point<dim> point) const;
 
 protected:
   /**
@@ -363,24 +219,6 @@ protected:
    * 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;
 };
 
 
@@ -394,16 +232,12 @@ private:
  *   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
+ *   @author Wolfgang Bangerth, 1998, Ralf Hartmann, 2001, Luca Heltai, 2013
  */
 template <int dim, int spacedim=dim>
-class HyperBallBoundary : public StraightBoundary<dim,spacedim>
+class HyperBallBoundary : public FlatManifold<spacedim>
 {
 public:
   /**
@@ -411,17 +245,7 @@ public:
    */
   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
@@ -431,72 +255,17 @@ public:
     get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
                  const std::vector<double> &weights) 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.
+   * Implementation of the function declared in the base class.
    *
-   * Refer to the general
-   * documentation of this class
-   * and the documentation of 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;
+  Point<spacedim>
+  normal_vector (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.
@@ -532,40 +301,15 @@ protected:
   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
+   * 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;
 };
 
 
@@ -592,74 +336,26 @@ public:
   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.
+   * Refer to the general documentation of this class and the
+   * documentation of the base class.
    */
     virtual
     Point<dim>
     get_new_point(const std::vector<Point<dim> > &surrounding_points,
                  const std::vector<double> &weights) 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.
+   * Implementation of the function declared in the base class.
    *
-   * Refer to the general
-   * documentation of this class
-   * and the documentation of the
-   * base class.
+   * 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;
+  virtual
+  Point<dim>
+  normal_vector (const Point<dim> p) const;
 };
 
 
@@ -726,58 +422,13 @@ public:
                           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.
+   * Compute the normal to the surface. 
    */
-  virtual void
-  get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
-                           typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const;
-
+    virtual Point<dim>
+    normal_vector (const Point<dim> point) const;
+    
 private:
   /**
    * Inner and outer radii of the shell.
@@ -796,7 +447,7 @@ private:
  * <tt>dim=2</tt>,<tt>spacedim=3</tt>, that is, just the surface.
  */
 template <int dim, int spacedim>
-class TorusBoundary : public Boundary<dim,spacedim>
+class TorusBoundary : public ManifoldChart<dim,spacedim>
 {
 public:
   /**
@@ -805,93 +456,42 @@ public:
    */
   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;
-
-  /**
-   * Refer to the general documentation of
-   * this class and the documentation of the
-   * base class.
-   */
-//   virtual Point<spacedim>
-//   get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
-//              const std::vector<double> &weights) 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.
+   * Pull back the given point in spacedim to the euclidean dim
+   * dimensional space.  Get the surface coordinates of the Torus,
+   * i.e., from <tt>(x,y,z)</tt> to <tt>(theta,phi)</tt>.
+   *
+   * Refer to the general documentation of this class for more
+   * information.
    */
-  virtual void  get_intermediate_points_on_quad (
-    const typename Triangulation< dim, spacedim >::quad_iterator &quad,
-    std::vector< Point< spacedim > >         &points ) const ;
+    virtual Point<dim>
+    pull_back(const Point<spacedim> space_point) const;
+    
+    /**
+     * Given a point in the dim dimensianal Euclidean space, this
+     * method returns a point on the manifold embedded in the spacedim
+     * Euclidean space.  Get the cartesian coordinates of the Torus,
+     * i.e., from <tt>(theta,phi)</tt> to <tt>(x,y,z)</tt>.
+     */
+    virtual Point<spacedim>
+    push_forward(const Point<dim> chart_point) const;
 
   /**
-   * Get the normal from cartesian
-   * coordinates. This normal does not have
-   * unit length.
+   * 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 ;
+  Point<spacedim>  normal_vector(const Point<spacedim> p) 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.
    */
@@ -903,59 +503,6 @@ private:
 
 /* -------------- 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 2e8b1f8bfdcc9e52623a125482d7a250e53354f7..7f99d819b10630fea3c572fa60d7d0d09e25e21c 100644 (file)
@@ -20,6 +20,14 @@ for (deal_II_space_dimension :  SPACE_DIMENSIONS)
   {    
     template class Manifold<deal_II_space_dimension>;
     template class FlatManifold<deal_II_space_dimension>;
+    template class ManifoldChart<deal_II_space_dimension,deal_II_space_dimension>;
+#if deal_II_space_dimension > 1
+    template class ManifoldChart<1,deal_II_space_dimension>;
+#endif
+#if deal_II_space_dimension > 2
+    template class ManifoldChart<2,deal_II_space_dimension>;
+#endif
+
   }
 
 
index 9bd26fb208539c7d6c346baae9a5b2541bbf329d..02975a31625dbb0539c7ce722ee5b1ce6f78ee53 100644 (file)
@@ -59,33 +59,33 @@ CylinderBoundary<dim,spacedim>::get_axis_vector (const unsigned int axis)
 
 
 
-template <int dim, int spacedim>
-Point<spacedim>
-CylinderBoundary<dim,spacedim>::
-get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const
-{
-  // compute a proposed new point
-  const Point<spacedim> middle = StraightBoundary<dim,spacedim>::get_new_point_on_line (line);
-
-  // we then have to project this
-  // point out to the given radius
-  // from the axis. to this end, we
-  // have to take into account the
-  // offset point_on_axis and the
-  // direction of the axis
-  const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
-                                           ((middle-point_on_axis) * direction) * direction;
-  // scale it to the desired length
-  // and put everything back
-  // together, unless we have a point
-  // on the axis
-  if (vector_from_axis.norm() <= 1e-10 * middle.norm())
-    return middle;
-  else
-    return (vector_from_axis / vector_from_axis.norm() * radius +
-            ((middle-point_on_axis) * direction) * direction +
-            point_on_axis);
-}
+// template <int dim, int spacedim>
+// Point<spacedim>
+// CylinderBoundary<dim,spacedim>::
+// get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const
+// {
+//   // compute a proposed new point
+//   const Point<spacedim> middle = StraightBoundary<dim,spacedim>::get_new_point_on_line (line);
+
+//   // we then have to project this
+//   // point out to the given radius
+//   // from the axis. to this end, we
+//   // have to take into account the
+//   // offset point_on_axis and the
+//   // direction of the axis
+//   const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
+//                                            ((middle-point_on_axis) * direction) * direction;
+//   // scale it to the desired length
+//   // and put everything back
+//   // together, unless we have a point
+//   // on the axis
+//   if (vector_from_axis.norm() <= 1e-10 * middle.norm())
+//     return middle;
+//   else
+//     return (vector_from_axis / vector_from_axis.norm() * radius +
+//             ((middle-point_on_axis) * direction) * direction +
+//             point_on_axis);
+// }
 
 template <int dim, int spacedim>
 Point<spacedim>
@@ -114,173 +114,168 @@ project_to_manifold (const Point<spacedim> middle) const
 
 
 
-template<>
-Point<3>
-CylinderBoundary<3>::
-get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
-{
-  const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
-
-  // same algorithm as above
-  const unsigned int spacedim = 3;
-
-  const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
-                                           ((middle-point_on_axis) * direction) * direction;
-  if (vector_from_axis.norm() <= 1e-10 * middle.norm())
-    return middle;
-  else
-    return (vector_from_axis / vector_from_axis.norm() * radius +
-            ((middle-point_on_axis) * direction) * direction +
-            point_on_axis);
-}
-
-template<>
-Point<3>
-CylinderBoundary<2,3>::
-get_new_point_on_quad (const Triangulation<2,3>::quad_iterator &quad) const
-{
-  const Point<3> middle = StraightBoundary<2,3>::get_new_point_on_quad (quad);
-
-  // same algorithm as above
-  const unsigned int spacedim = 3;
-  const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
-                                           ((middle-point_on_axis) * direction) * direction;
-  if (vector_from_axis.norm() <= 1e-10 * middle.norm())
-    return middle;
-  else
-    return (vector_from_axis / vector_from_axis.norm() * radius +
-            ((middle-point_on_axis) * direction) * direction +
-            point_on_axis);
-}
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-CylinderBoundary<dim,spacedim>::
-get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &) const
-{
-  Assert (false, ExcImpossibleInDim(dim));
-  return Point<spacedim>();
-}
-
-
-
-template <int dim, int spacedim>
-void
-CylinderBoundary<dim,spacedim>::get_intermediate_points_on_line (
-  const typename Triangulation<dim,spacedim>::line_iterator &line,
-  std::vector<Point<spacedim> > &points) const
-{
-  if (points.size()==1)
-    points[0]=get_new_point_on_line(line);
-  else
-    get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points);
-}
-
+// template<>
+// Point<3>
+// CylinderBoundary<3>::
+// get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
+// {
+//   const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
+
+//   // same algorithm as above
+//   const unsigned int spacedim = 3;
+
+//   const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
+//                                            ((middle-point_on_axis) * direction) * direction;
+//   if (vector_from_axis.norm() <= 1e-10 * middle.norm())
+//     return middle;
+//   else
+//     return (vector_from_axis / vector_from_axis.norm() * radius +
+//             ((middle-point_on_axis) * direction) * direction +
+//             point_on_axis);
+// }
 
-template <int dim, int spacedim>
-void
-CylinderBoundary<dim,spacedim>::get_intermediate_points_between_points (
-  const Point<spacedim> &v0,
-  const Point<spacedim> &v1,
-  std::vector<Point<spacedim> > &points) const
-{
-  const unsigned int n=points.size();
-  Assert(n>0, ExcInternalError());
+// template<>
+// Point<3>
+// CylinderBoundary<2,3>::
+// get_new_point_on_quad (const Triangulation<2,3>::quad_iterator &quad) const
+// {
+//   const Point<3> middle = StraightBoundary<2,3>::get_new_point_on_quad (quad);
+
+//   // same algorithm as above
+//   const unsigned int spacedim = 3;
+//   const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
+//                                            ((middle-point_on_axis) * direction) * direction;
+//   if (vector_from_axis.norm() <= 1e-10 * middle.norm())
+//     return middle;
+//   else
+//     return (vector_from_axis / vector_from_axis.norm() * radius +
+//             ((middle-point_on_axis) * direction) * direction +
+//             point_on_axis);
+// }
 
-  // Do a simple linear interpolation followed by projection, using the same
-  // algorithm as above
-  const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
 
-  for (unsigned int i=0; i<n; ++i)
-    {
-      const double x = line_points[i+1][0];
-      const Point<spacedim> middle = (1-x)*v0 + x*v1;
+// template <int dim, int spacedim>
+// Point<spacedim>
+// CylinderBoundary<dim,spacedim>::
+// get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &) const
+// {
+//   Assert (false, ExcImpossibleInDim(dim));
+//   return Point<spacedim>();
+// }
 
-      const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
-                                               ((middle-point_on_axis) * direction) * direction;
-      if (vector_from_axis.norm() <= 1e-10 * middle.norm())
-        points[i] = middle;
-      else
-        points[i] = (vector_from_axis / vector_from_axis.norm() * radius +
-                     ((middle-point_on_axis) * direction) * direction +
-                     point_on_axis);
-    }
-}
 
 
+// template <int dim, int spacedim>
+// void
+// CylinderBoundary<dim,spacedim>::get_intermediate_points_on_line (
+//   const typename Triangulation<dim,spacedim>::line_iterator &line,
+//   std::vector<Point<spacedim> > &points) const
+// {
+//   if (points.size()==1)
+//     points[0]=get_new_point_on_line(line);
+//   else
+//     get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points);
+// }
 
-template <>
-void
-CylinderBoundary<3>::get_intermediate_points_on_quad (
-  const Triangulation<3>::quad_iterator &quad,
-  std::vector<Point<3> > &points) const
-{
-  if (points.size()==1)
-    points[0]=get_new_point_on_quad(quad);
-  else
-    {
-      unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
-      Assert(points.size()==m*m, ExcInternalError());
 
-      std::vector<Point<3> > lp0(m);
-      std::vector<Point<3> > lp1(m);
+// template <int dim, int spacedim>
+// void
+// CylinderBoundary<dim,spacedim>::get_intermediate_points_between_points (
+//   const Point<spacedim> &v0,
+//   const Point<spacedim> &v1,
+//   std::vector<Point<spacedim> > &points) const
+// {
+//   const unsigned int n=points.size();
+//   Assert(n>0, ExcInternalError());
+
+//   // Do a simple linear interpolation followed by projection, using the same
+//   // algorithm as above
+//   const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
+
+//   for (unsigned int i=0; i<n; ++i)
+//     {
+//       const double x = line_points[i+1][0];
+//       const Point<spacedim> middle = (1-x)*v0 + x*v1;
+
+//       const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
+//                                                ((middle-point_on_axis) * direction) * direction;
+//       if (vector_from_axis.norm() <= 1e-10 * middle.norm())
+//         points[i] = middle;
+//       else
+//         points[i] = (vector_from_axis / vector_from_axis.norm() * radius +
+//                      ((middle-point_on_axis) * direction) * direction +
+//                      point_on_axis);
+//     }
+// }
 
-      get_intermediate_points_on_line(quad->line(0), lp0);
-      get_intermediate_points_on_line(quad->line(1), lp1);
 
-      std::vector<Point<3> > lps(m);
-      for (unsigned int i=0; i<m; ++i)
-        {
-          get_intermediate_points_between_points(lp0[i], lp1[i], lps);
 
-          for (unsigned int j=0; j<m; ++j)
-            points[i*m+j]=lps[j];
-        }
-    }
-}
+// template <>
+// void
+// CylinderBoundary<3>::get_intermediate_points_on_quad (
+//   const Triangulation<3>::quad_iterator &quad,
+//   std::vector<Point<3> > &points) const
+// {
+//   if (points.size()==1)
+//     points[0]=get_new_point_on_quad(quad);
+//   else
+//     {
+//       unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
+//       Assert(points.size()==m*m, ExcInternalError());
+
+//       std::vector<Point<3> > lp0(m);
+//       std::vector<Point<3> > lp1(m);
+
+//       get_intermediate_points_on_line(quad->line(0), lp0);
+//       get_intermediate_points_on_line(quad->line(1), lp1);
+
+//       std::vector<Point<3> > lps(m);
+//       for (unsigned int i=0; i<m; ++i)
+//         {
+//           get_intermediate_points_between_points(lp0[i], lp1[i], lps);
+
+//           for (unsigned int j=0; j<m; ++j)
+//             points[i*m+j]=lps[j];
+//         }
+//     }
+// }
 
 
 
-template <int dim, int spacedim>
-void
-CylinderBoundary<dim,spacedim>::get_intermediate_points_on_quad (
-  const typename Triangulation<dim,spacedim>::quad_iterator &,
-  std::vector<Point<spacedim> > &) const
-{
-  Assert (false, ExcImpossibleInDim(dim));
-}
+// template <int dim, int spacedim>
+// void
+// CylinderBoundary<dim,spacedim>::get_intermediate_points_on_quad (
+//   const typename Triangulation<dim,spacedim>::quad_iterator &,
+//   std::vector<Point<spacedim> > &) const
+// {
+//   Assert (false, ExcImpossibleInDim(dim));
+// }
 
 
 
 
-template <>
-void
-CylinderBoundary<1>::
-get_normals_at_vertices (const Triangulation<1>::face_iterator &,
-                         Boundary<1,1>::FaceVertexNormals &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
+// template <>
+// void
+// CylinderBoundary<1>::
+// get_normals_at_vertices (const Triangulation<1>::face_iterator &,
+//                          Boundary<1,1>::FaceVertexNormals &) const
+// {
+//   Assert (false, ExcImpossibleInDim(1));
+// }
 
 
 
 
 template <int dim, int spacedim>
-void
+Point<spacedim>
 CylinderBoundary<dim,spacedim>::
-get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
-                         typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const
+normal_vector (const Point<spacedim> vertex) const
 {
-  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
-    {
-      const Point<spacedim> vertex = face->vertex(v);
-
-      const Point<spacedim> vector_from_axis = (vertex-point_on_axis) -
-                                               ((vertex-point_on_axis) * direction) * direction;
+  const Point<spacedim> vector_from_axis = (vertex-point_on_axis) -
+                                          ((vertex-point_on_axis) * direction) * direction;
 
-      face_vertex_normals[v] = (vector_from_axis / vector_from_axis.norm());
-    }
+  Point<spacedim> normal = (vector_from_axis / vector_from_axis.norm());
+  return normal;
 }
 
 
@@ -321,51 +316,51 @@ double ConeBoundary<dim>::get_radius (Point<dim> x) const
 
 
 
-template<int dim>
-void
-ConeBoundary<dim>::
-get_intermediate_points_between_points (const Point<dim> &p0,
-                                        const Point<dim> &p1,
-                                        std::vector<Point<dim> > &points) const
-{
-  const unsigned int n = points.size ();
-  const Point<dim> axis = x_1 - x_0;
-
-  Assert (n > 0, ExcInternalError ());
-
-  const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
-
-  for (unsigned int i=0; i<n; ++i)
-    {
-      const double x = line_points[i+1][0];
-
-      // Compute the current point.
-      const Point<dim> x_i = (1-x)*p0 + x*p1;
-      // To project this point on the boundary of the cone we first compute
-      // the orthogonal projection of this point onto the axis of the cone.
-      const double c = (x_i - x_0) * axis / axis.square ();
-      const Point<dim> x_ip = x_0 + c * axis;
-      // Compute the projection of the middle point on the boundary of the
-      // cone.
-      points[i] = x_ip + get_radius (x_ip) *  (x_i - x_ip) / (x_i - x_ip).norm ();
-    }
-}
+// template<int dim>
+// void
+// ConeBoundary<dim>::
+// get_intermediate_points_between_points (const Point<dim> &p0,
+//                                         const Point<dim> &p1,
+//                                         std::vector<Point<dim> > &points) const
+// {
+//   const unsigned int n = points.size ();
+//   const Point<dim> axis = x_1 - x_0;
+
+//   Assert (n > 0, ExcInternalError ());
+
+//   const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
+
+//   for (unsigned int i=0; i<n; ++i)
+//     {
+//       const double x = line_points[i+1][0];
+
+//       // Compute the current point.
+//       const Point<dim> x_i = (1-x)*p0 + x*p1;
+//       // To project this point on the boundary of the cone we first compute
+//       // the orthogonal projection of this point onto the axis of the cone.
+//       const double c = (x_i - x_0) * axis / axis.square ();
+//       const Point<dim> x_ip = x_0 + c * axis;
+//       // Compute the projection of the middle point on the boundary of the
+//       // cone.
+//       points[i] = x_ip + get_radius (x_ip) *  (x_i - x_ip) / (x_i - x_ip).norm ();
+//     }
+// }
 
-template<int dim>
-Point<dim>
-ConeBoundary<dim>::
-get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
-{
-  const Point<dim> axis = x_1 - x_0;
-  // Compute the middle point of the line.
-  const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
-  // To project it on the boundary of the cone we first compute the orthogonal
-  // projection of the middle point onto the axis of the cone.
-  const double c = (middle - x_0) * axis / axis.square ();
-  const Point<dim> middle_p = x_0 + c * axis;
-  // Compute the projection of the middle point on the boundary of the cone.
-  return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
-}
+// template<int dim>
+// Point<dim>
+// ConeBoundary<dim>::
+// get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+// {
+//   const Point<dim> axis = x_1 - x_0;
+//   // Compute the middle point of the line.
+//   const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
+//   // To project it on the boundary of the cone we first compute the orthogonal
+//   // projection of the middle point onto the axis of the cone.
+//   const double c = (middle - x_0) * axis / axis.square ();
+//   const Point<dim> middle_p = x_0 + c * axis;
+//   // Compute the projection of the middle point on the boundary of the cone.
+//   return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
+// }
 
 
 template<int dim>
@@ -384,133 +379,130 @@ project_to_manifold (const Point<dim> middle) const
 
 
 
-template <>
-Point<3>
-ConeBoundary<3>::
-get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
-{
-  const int dim = 3;
-
-  const Point<dim> axis = x_1 - x_0;
-  // Compute the middle point of the quad.
-  const Point<dim> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
-  // Same algorithm as above: To project it on the boundary of the cone we
-  // first compute the orthogonal projection of the middle point onto the axis
-  // of the cone.
-  const double c = (middle - x_0) * axis / axis.square ();
-  const Point<dim> middle_p = x_0 + c * axis;
-  // Compute the projection of the middle point on the boundary of the cone.
-  return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
-}
+// template <>
+// Point<3>
+// ConeBoundary<3>::
+// get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
+// {
+//   const int dim = 3;
+
+//   const Point<dim> axis = x_1 - x_0;
+//   // Compute the middle point of the quad.
+//   const Point<dim> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
+//   // Same algorithm as above: To project it on the boundary of the cone we
+//   // first compute the orthogonal projection of the middle point onto the axis
+//   // of the cone.
+//   const double c = (middle - x_0) * axis / axis.square ();
+//   const Point<dim> middle_p = x_0 + c * axis;
+//   // Compute the projection of the middle point on the boundary of the cone.
+//   return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
+// }
 
 
 
-template<int dim>
-Point<dim>
-ConeBoundary<dim>::
-get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &) const
-{
-  Assert (false, ExcImpossibleInDim (dim));
+// template<int dim>
+// Point<dim>
+// ConeBoundary<dim>::
+// get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &) const
+// {
+//   Assert (false, ExcImpossibleInDim (dim));
 
-  return Point<dim>();
-}
+//   return Point<dim>();
+// }
 
 
 
-template<int dim>
-void
-ConeBoundary<dim>::
-get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
-                                 std::vector<Point<dim> > &points) const
-{
-  if (points.size () == 1)
-    points[0] = get_new_point_on_line (line);
-  else
-    get_intermediate_points_between_points (line->vertex (0), line->vertex (1), points);
-}
+// template<int dim>
+// void
+// ConeBoundary<dim>::
+// get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
+//                                  std::vector<Point<dim> > &points) const
+// {
+//   if (points.size () == 1)
+//     points[0] = get_new_point_on_line (line);
+//   else
+//     get_intermediate_points_between_points (line->vertex (0), line->vertex (1), points);
+// }
 
 
 
 
-template<>
-void
-ConeBoundary<3>::
-get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad,
-                                 std::vector<Point<3> > &points) const
-{
-  if (points.size () == 1)
-    points[0] = get_new_point_on_quad (quad);
-  else
-    {
-      unsigned int n = static_cast<unsigned int> (std::sqrt (static_cast<double> (points.size ())));
+// template<>
+// void
+// ConeBoundary<3>::
+// get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad,
+//                                  std::vector<Point<3> > &points) const
+// {
+//   if (points.size () == 1)
+//     points[0] = get_new_point_on_quad (quad);
+//   else
+//     {
+//       unsigned int n = static_cast<unsigned int> (std::sqrt (static_cast<double> (points.size ())));
 
-      Assert (points.size () == n * n, ExcInternalError ());
+//       Assert (points.size () == n * n, ExcInternalError ());
 
-      std::vector<Point<3> > points_line_0 (n);
-      std::vector<Point<3> > points_line_1 (n);
+//       std::vector<Point<3> > points_line_0 (n);
+//       std::vector<Point<3> > points_line_1 (n);
 
-      get_intermediate_points_on_line (quad->line (0), points_line_0);
-      get_intermediate_points_on_line (quad->line (1), points_line_1);
+//       get_intermediate_points_on_line (quad->line (0), points_line_0);
+//       get_intermediate_points_on_line (quad->line (1), points_line_1);
 
-      std::vector<Point<3> > points_line_segment (n);
+//       std::vector<Point<3> > points_line_segment (n);
 
-      for (unsigned int i = 0; i < n; ++i)
-        {
-          get_intermediate_points_between_points (points_line_0[i],
-                                                  points_line_1[i],
-                                                  points_line_segment);
+//       for (unsigned int i = 0; i < n; ++i)
+//         {
+//           get_intermediate_points_between_points (points_line_0[i],
+//                                                   points_line_1[i],
+//                                                   points_line_segment);
 
-          for (unsigned int j = 0; j < n; ++j)
-            points[i * n + j] = points_line_segment[j];
-        }
-    }
-}
+//           for (unsigned int j = 0; j < n; ++j)
+//             points[i * n + j] = points_line_segment[j];
+//         }
+//     }
+// }
 
 
 
-template <int dim>
-void
-ConeBoundary<dim>::
-get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &,
-                                 std::vector<Point<dim> > &) const
-{
-  Assert (false, ExcImpossibleInDim (dim));
-}
+// template <int dim>
+// void
+// ConeBoundary<dim>::
+// get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &,
+//                                  std::vector<Point<dim> > &) const
+// {
+//   Assert (false, ExcImpossibleInDim (dim));
+// }
 
 
 
 
-template<>
-void
-ConeBoundary<1>::
-get_normals_at_vertices (const Triangulation<1>::face_iterator &,
-                         Boundary<1,1>::FaceVertexNormals &) const
-{
-  Assert (false, ExcImpossibleInDim (1));
-}
+// template<>
+// void
+// ConeBoundary<1>::
+// get_normals_at_vertices (const Triangulation<1>::face_iterator &,
+//                          Boundary<1,1>::FaceVertexNormals &) const
+// {
+//   Assert (false, ExcImpossibleInDim (1));
+// }
 
 
 
 template<int dim>
-void
+Point<dim>
 ConeBoundary<dim>::
-get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
-                         typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
+normal_vector (const Point<dim> vertex) const 
 {
   const Point<dim> axis = x_1 - x_0;
 
-  for (unsigned int vertex = 0; vertex < GeometryInfo<dim>::vertices_per_cell; ++vertex)
-    {
-      // Compute the orthogonal projection of the vertex onto the axis of the
-      // cone.
-      const double c = (face->vertex (vertex) - x_0) * axis / axis.square ();
-      const Point<dim> vertex_p = x_0 + c * axis;
+  // Compute the orthogonal projection of the vertex onto the axis of
+  // the cone.
+  const double c = (vertex - x_0) * axis / axis.square ();
+  const Point<dim> vertex_p = x_0 + c * axis;
       // Then compute the vector pointing from the point <tt>vertex_p</tt> on
       // the axis to the vertex.
-      const Point<dim> axis_to_vertex = face->vertex (vertex) - vertex_p;
+  const Point<dim> axis_to_vertex = vertex - vertex_p;
 
-      face_vertex_normals[vertex] = axis_to_vertex / axis_to_vertex.norm ();
-    }
+  Point<dim> normal = axis_to_vertex / axis_to_vertex.norm ();
+  return normal;
 }
 
 
@@ -544,27 +536,27 @@ HyperBallBoundary<dim,spacedim>::HyperBallBoundary (const Point<spacedim> p,
 // }
 
 
-template <int dim, int spacedim>
-Point<spacedim>
-HyperBallBoundary<dim,spacedim>::get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const
-{
-  Point<spacedim> middle = StraightBoundary<dim,spacedim>::get_new_point_on_line (line);
-
-  middle -= center;
-
-  double r=0;
-  if (compute_radius_automatically)
-    {
-      const Point<spacedim> vertex_relative = line->vertex(0) - center;
-      r = std::sqrt(vertex_relative.square());
-    }
-  else
-    r=radius;
-  // project to boundary
-  middle *= r / std::sqrt(middle.square());
-  middle += center;
-  return middle;
-}
+// template <int dim, int spacedim>
+// Point<spacedim>
+// HyperBallBoundary<dim,spacedim>::get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const
+// {
+//   Point<spacedim> middle = StraightBoundary<dim,spacedim>::get_new_point_on_line (line);
+
+//   middle -= center;
+
+//   double r=0;
+//   if (compute_radius_automatically)
+//     {
+//       const Point<spacedim> vertex_relative = line->vertex(0) - center;
+//       r = std::sqrt(vertex_relative.square());
+//     }
+//   else
+//     r=radius;
+//   // project to boundary
+//   middle *= r / std::sqrt(middle.square());
+//   middle += center;
+//   return middle;
+// }
 
 
 template <int dim, int spacedim>
@@ -590,257 +582,253 @@ HyperBallBoundary<dim,spacedim>::get_new_point(const std::vector<Point<spacedim>
 
 
 
-template <>
-Point<1>
-HyperBallBoundary<1,1>::
-get_new_point_on_quad (const Triangulation<1,1>::quad_iterator &) const
-{
-  Assert (false, ExcInternalError());
-  return Point<1>();
-}
-
-
-template <>
-Point<2>
-HyperBallBoundary<1,2>::
-get_new_point_on_quad (const Triangulation<1,2>::quad_iterator &) const
-{
-  Assert (false, ExcInternalError());
-  return Point<2>();
-}
-
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-HyperBallBoundary<dim,spacedim>::
-get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const
-{
-  Point<spacedim> middle = StraightBoundary<dim,spacedim>::get_new_point_on_quad (quad);
-
-  middle -= center;
-
-  double r=0;
-  if (compute_radius_automatically)
-    {
-      const Point<spacedim> vertex_relative = quad->vertex(0) - center;
-      r = std::sqrt(vertex_relative.square());
-    }
-  else
-    r=radius;
-  // project to boundary
-  middle *= r / std::sqrt(middle.square());
-
-  middle += center;
-  return middle;
-}
-
-
-
-template <>
-void
-HyperBallBoundary<1>::get_intermediate_points_on_line (
-  const Triangulation<1>::line_iterator &,
-  std::vector<Point<1> > &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
-
-
-
-template <int dim, int spacedim>
-void
-HyperBallBoundary<dim,spacedim>::get_intermediate_points_on_line (
-  const typename Triangulation<dim,spacedim>::line_iterator &line,
-  std::vector<Point<spacedim> > &points) const
-{
-  if (points.size()==1)
-    points[0]=get_new_point_on_line(line);
-  else
-    get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points);
-}
-
-
-
-template <int dim, int spacedim>
-void
-HyperBallBoundary<dim,spacedim>::get_intermediate_points_between_points (
-  const Point<spacedim> &p0, const Point<spacedim> &p1,
-  std::vector<Point<spacedim> > &points) const
-{
-  const unsigned int n=points.size();
-  Assert(n>0, ExcInternalError());
-
-  const Point<spacedim> v0=p0-center,
-                        v1=p1-center;
-  const double length=std::sqrt((v1-v0).square());
-
-  double eps=1e-12;
-  double r=0;
-  if (compute_radius_automatically)
-    {
-      const Point<spacedim> vertex_relative = p0 - center;
-      r = std::sqrt(vertex_relative.square());
-    }
-  else
-    r=radius;
-
-
-  const double r2=r*r;
-  Assert(std::fabs(v0.square()-r2)<eps*r2, ExcInternalError());
-  Assert(std::fabs(v1.square()-r2)<eps*r2, ExcInternalError());
+// template <>
+// Point<1>
+// HyperBallBoundary<1,1>::
+// get_new_point_on_quad (const Triangulation<1,1>::quad_iterator &) const
+// {
+//   Assert (false, ExcInternalError());
+//   return Point<1>();
+// }
 
-  const double alpha=std::acos((v0*v1)/std::sqrt(v0.square()*v1.square()));
-  const Point<spacedim> pm=0.5*(v0+v1);
 
-  const double h=std::sqrt(pm.square());
+// template <>
+// Point<2>
+// HyperBallBoundary<1,2>::
+// get_new_point_on_quad (const Triangulation<1,2>::quad_iterator &) const
+// {
+//   Assert (false, ExcInternalError());
+//   return Point<2>();
+// }
 
-  // n even:  m=n/2,
-  // n odd:   m=(n-1)/2
-  const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
-  const unsigned int m=n/2;
-  for (unsigned int i=0; i<m ; ++i)
-    {
-      const double beta = alpha * (line_points[i+1][0]-0.5);
-      const double d=h*std::tan(beta);
-      points[i]=pm+d/length*(v1-v0);
-      points[n-1-i]=pm-d/length*(v1-v0);
-    }
 
-  if ((n+1)%2==0)
-    // if the number of parts is even insert the midpoint
-    points[(n-1)/2]=pm;
 
+// template <int dim, int spacedim>
+// Point<spacedim>
+// HyperBallBoundary<dim,spacedim>::
+// get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const
+// {
+//   Point<spacedim> middle = StraightBoundary<dim,spacedim>::get_new_point_on_quad (quad);
+
+//   middle -= center;
+
+//   double r=0;
+//   if (compute_radius_automatically)
+//     {
+//       const Point<spacedim> vertex_relative = quad->vertex(0) - center;
+//       r = std::sqrt(vertex_relative.square());
+//     }
+//   else
+//     r=radius;
+//   // project to boundary
+//   middle *= r / std::sqrt(middle.square());
+
+//   middle += center;
+//   return middle;
+// }
 
-  // project the points from the straight line to the HyperBallBoundary
-  for (unsigned int i=0; i<n; ++i)
-    {
-      points[i] *= r / std::sqrt(points[i].square());
-      points[i] += center;
-    }
-}
 
 
+// template <>
+// void
+// HyperBallBoundary<1>::get_intermediate_points_on_line (
+//   const Triangulation<1>::line_iterator &,
+//   std::vector<Point<1> > &) const
+// {
+//   Assert (false, ExcImpossibleInDim(1));
+// }
 
-template <>
-void
-HyperBallBoundary<3>::get_intermediate_points_on_quad (
-  const Triangulation<3>::quad_iterator &quad,
-  std::vector<Point<3> > &points) const
-{
-  if (points.size()==1)
-    points[0]=get_new_point_on_quad(quad);
-  else
-    {
-      unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
-      Assert(points.size()==m*m, ExcInternalError());
 
-      std::vector<Point<3> > lp0(m);
-      std::vector<Point<3> > lp1(m);
 
-      get_intermediate_points_on_line(quad->line(0), lp0);
-      get_intermediate_points_on_line(quad->line(1), lp1);
+// template <int dim, int spacedim>
+// void
+// HyperBallBoundary<dim,spacedim>::get_intermediate_points_on_line (
+//   const typename Triangulation<dim,spacedim>::line_iterator &line,
+//   std::vector<Point<spacedim> > &points) const
+// {
+//   if (points.size()==1)
+//     points[0]=get_new_point_on_line(line);
+//   else
+//     get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points);
+// }
 
-      std::vector<Point<3> > lps(m);
-      for (unsigned int i=0; i<m; ++i)
-        {
-          get_intermediate_points_between_points(lp0[i], lp1[i], lps);
 
-          for (unsigned int j=0; j<m; ++j)
-            points[i*m+j]=lps[j];
-        }
-    }
-}
 
+// template <int dim, int spacedim>
+// void
+// HyperBallBoundary<dim,spacedim>::get_intermediate_points_between_points (
+//   const Point<spacedim> &p0, const Point<spacedim> &p1,
+//   std::vector<Point<spacedim> > &points) const
+// {
+//   const unsigned int n=points.size();
+//   Assert(n>0, ExcInternalError());
+
+//   const Point<spacedim> v0=p0-center,
+//                         v1=p1-center;
+//   const double length=std::sqrt((v1-v0).square());
+
+//   double eps=1e-12;
+//   double r=0;
+//   if (compute_radius_automatically)
+//     {
+//       const Point<spacedim> vertex_relative = p0 - center;
+//       r = std::sqrt(vertex_relative.square());
+//     }
+//   else
+//     r=radius;
+
+
+//   const double r2=r*r;
+//   Assert(std::fabs(v0.square()-r2)<eps*r2, ExcInternalError());
+//   Assert(std::fabs(v1.square()-r2)<eps*r2, ExcInternalError());
+
+//   const double alpha=std::acos((v0*v1)/std::sqrt(v0.square()*v1.square()));
+//   const Point<spacedim> pm=0.5*(v0+v1);
+
+//   const double h=std::sqrt(pm.square());
+
+//   // n even:  m=n/2,
+//   // n odd:   m=(n-1)/2
+//   const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
+//   const unsigned int m=n/2;
+//   for (unsigned int i=0; i<m ; ++i)
+//     {
+//       const double beta = alpha * (line_points[i+1][0]-0.5);
+//       const double d=h*std::tan(beta);
+//       points[i]=pm+d/length*(v1-v0);
+//       points[n-1-i]=pm-d/length*(v1-v0);
+//     }
+
+//   if ((n+1)%2==0)
+//     // if the number of parts is even insert the midpoint
+//     points[(n-1)/2]=pm;
+
+
+//   // project the points from the straight line to the HyperBallBoundary
+//   for (unsigned int i=0; i<n; ++i)
+//     {
+//       points[i] *= r / std::sqrt(points[i].square());
+//       points[i] += center;
+//     }
+// }
 
 
-template <>
-void
-HyperBallBoundary<2,3>::get_intermediate_points_on_quad (
-  const Triangulation<2,3>::quad_iterator &quad,
-  std::vector<Point<3> > &points) const
-{
-  if (points.size()==1)
-    points[0]=get_new_point_on_quad(quad);
-  else
-    {
-      unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
-      Assert(points.size()==m*m, ExcInternalError());
 
-      std::vector<Point<3> > lp0(m);
-      std::vector<Point<3> > lp1(m);
+// template <>
+// void
+// HyperBallBoundary<3>::get_intermediate_points_on_quad (
+//   const Triangulation<3>::quad_iterator &quad,
+//   std::vector<Point<3> > &points) const
+// {
+//   if (points.size()==1)
+//     points[0]=get_new_point_on_quad(quad);
+//   else
+//     {
+//       unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
+//       Assert(points.size()==m*m, ExcInternalError());
+
+//       std::vector<Point<3> > lp0(m);
+//       std::vector<Point<3> > lp1(m);
+
+//       get_intermediate_points_on_line(quad->line(0), lp0);
+//       get_intermediate_points_on_line(quad->line(1), lp1);
+
+//       std::vector<Point<3> > lps(m);
+//       for (unsigned int i=0; i<m; ++i)
+//         {
+//           get_intermediate_points_between_points(lp0[i], lp1[i], lps);
+
+//           for (unsigned int j=0; j<m; ++j)
+//             points[i*m+j]=lps[j];
+//         }
+//     }
+// }
 
-      get_intermediate_points_on_line(quad->line(0), lp0);
-      get_intermediate_points_on_line(quad->line(1), lp1);
 
-      std::vector<Point<3> > lps(m);
-      for (unsigned int i=0; i<m; ++i)
-        {
-          get_intermediate_points_between_points(lp0[i], lp1[i], lps);
 
-          for (unsigned int j=0; j<m; ++j)
-            points[i*m+j]=lps[j];
-        }
-    }
-}
+// template <>
+// void
+// HyperBallBoundary<2,3>::get_intermediate_points_on_quad (
+//   const Triangulation<2,3>::quad_iterator &quad,
+//   std::vector<Point<3> > &points) const
+// {
+//   if (points.size()==1)
+//     points[0]=get_new_point_on_quad(quad);
+//   else
+//     {
+//       unsigned int m=static_cast<unsigned int> (std::sqrt(static_cast<double>(points.size())));
+//       Assert(points.size()==m*m, ExcInternalError());
+
+//       std::vector<Point<3> > lp0(m);
+//       std::vector<Point<3> > lp1(m);
+
+//       get_intermediate_points_on_line(quad->line(0), lp0);
+//       get_intermediate_points_on_line(quad->line(1), lp1);
+
+//       std::vector<Point<3> > lps(m);
+//       for (unsigned int i=0; i<m; ++i)
+//         {
+//           get_intermediate_points_between_points(lp0[i], lp1[i], lps);
+
+//           for (unsigned int j=0; j<m; ++j)
+//             points[i*m+j]=lps[j];
+//         }
+//     }
+// }
 
 
 
-template <int dim, int spacedim>
-void
-HyperBallBoundary<dim,spacedim>::get_intermediate_points_on_quad (
-  const typename Triangulation<dim,spacedim>::quad_iterator &,
-  std::vector<Point<spacedim> > &) const
-{
-  Assert(false, ExcImpossibleInDim(dim));
-}
+// template <int dim, int spacedim>
+// void
+// HyperBallBoundary<dim,spacedim>::get_intermediate_points_on_quad (
+//   const typename Triangulation<dim,spacedim>::quad_iterator &,
+//   std::vector<Point<spacedim> > &) const
+// {
+//   Assert(false, ExcImpossibleInDim(dim));
+// }
 
 
 
-template <int dim, int spacedim>
-Tensor<1,spacedim>
-HyperBallBoundary<dim,spacedim>::
-normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &,
-               const Point<spacedim> &p) const
-{
-  const Tensor<1,spacedim> unnormalized_normal = p-center;
-  return unnormalized_normal/unnormalized_normal.norm();
-}
+// template <int dim, int spacedim>
+// Tensor<1,spacedim>
+// HyperBallBoundary<dim,spacedim>::
+// normal_vector (const typename Triangulation<dim,spacedim>::face_iterator &,
+//                const Point<spacedim> &p) const
+// {
+//   const Tensor<1,spacedim> unnormalized_normal = p-center;
+//   return unnormalized_normal/unnormalized_normal.norm();
+// }
 
 
 
-template <>
-void
-HyperBallBoundary<1>::
-get_normals_at_vertices (const Triangulation<1>::face_iterator &,
-                         Boundary<1,1>::FaceVertexNormals &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
+// template <>
+// void
+// HyperBallBoundary<1>::
+// get_normals_at_vertices (const Triangulation<1>::face_iterator &,
+//                          Boundary<1,1>::FaceVertexNormals &) const
+// {
+//   Assert (false, ExcImpossibleInDim(1));
+// }
 
-template <>
-void
-HyperBallBoundary<1,2>::
-get_normals_at_vertices (const Triangulation<1,2>::face_iterator &,
-                         Boundary<1,2>::FaceVertexNormals &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
+// template <>
+// void
+// HyperBallBoundary<1,2>::
+// get_normals_at_vertices (const Triangulation<1,2>::face_iterator &,
+//                          Boundary<1,2>::FaceVertexNormals &) const
+// {
+//   Assert (false, ExcImpossibleInDim(1));
+// }
 
 
 
 template <int dim, int spacedim>
-void
+Point<spacedim>
 HyperBallBoundary<dim,spacedim>::
-get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
-                         typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const
+normal_vector (const Point<spacedim> vertex) const
 {
-  for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_face; ++vertex)
-    face_vertex_normals[vertex] = face->vertex(vertex)-center;
+  return vertex-center;
 }
 
-
-
 template <int dim, int spacedim>
 Point<spacedim>
 HyperBallBoundary<dim,spacedim>::get_center () const
@@ -869,13 +857,13 @@ HalfHyperBallBoundary<dim>::HalfHyperBallBoundary (const Point<dim> center,
 {}
 
 
-template <int dim>
-Point<dim>
-HalfHyperBallBoundary<dim>::get_new_point_on_line(const typename Triangulation<dim>::line_iterator &line) const
-{
-  Assert(false, ExcInternalError());
-  return Point<dim>();
-}
+// template <int dim>
+// Point<dim>
+// HalfHyperBallBoundary<dim>::get_new_point_on_line(const typename Triangulation<dim>::line_iterator &line) const
+// {
+//   Assert(false, ExcInternalError());
+//   return Point<dim>();
+// }
 
   
 template <int dim>
@@ -883,117 +871,157 @@ Point<dim>
 HalfHyperBallBoundary<dim>::get_new_point(const std::vector<Point<dim> > &surrounding_points,
                                          const std::vector<double> &weights) const
 {
-  // check whether center of object is at x==x_center, since then it belongs
-  // to the plane part of the boundary. however, this is not the case if it is
-  // at the outer perimeter
+  // check whether every point is at distance R
+  bool all_at_distance_R = true;
+  for(unsigned int i=0; i<surrounding_points.size(); ++i)
+    if( (surrounding_points[i].distance(this->center)-this->radius) >1e-5)
+      {
+       all_at_distance_R = false;
+       break;
+      }
+  
   const Point<dim> object_center = FlatManifold<dim>::get_new_point(surrounding_points, weights);
-  if (object_center.distance(this->center) < 1e-5)
+  if (std::abs(object_center(0)-this->center(0)) < 1e-5  && !all_at_distance_R)
     return object_center;
   else
-    return HyperBallBoundary<dim>::get_new_point (surrounding_points,weights);
+    if(all_at_distance_R)
+      return HyperBallBoundary<dim>::get_new_point (surrounding_points,weights);
+    else
+      {
+       Assert(false, ExcInternalError());
+       return Point<dim>();
+      }
 }
 
 
 
-template <>
-Point<1>
-HalfHyperBallBoundary<1>::
-get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
-{
-  Assert (false, ExcInternalError());
-  return Point<1>();
-}
+// template <>
+// Point<1>
+// HalfHyperBallBoundary<1>::
+// get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
+// {
+//   Assert (false, ExcInternalError());
+//   return Point<1>();
+// }
 
 
 
-template <int dim>
-Point<dim>
-HalfHyperBallBoundary<dim>::
-get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
-{
-  const Point<dim> quad_center = quad->center();
-  if (quad_center(0) == this->center(0))
-    return quad_center;
-  else
-    return HyperBallBoundary<dim>::get_new_point_on_quad (quad);
-}
+// template <int dim>
+// Point<dim>
+// HalfHyperBallBoundary<dim>::
+// get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
+// {
+//   const Point<dim> quad_center = quad->center();
+//   if (quad_center(0) == this->center(0))
+//     return quad_center;
+//   else
+//     return HyperBallBoundary<dim>::get_new_point_on_quad (quad);
+// }
 
 
 
-template <int dim>
-void
-HalfHyperBallBoundary<dim>::
-get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
-                                 std::vector<Point<dim> > &points) const
-{
-  // check whether center of object is at x==0, since then it belongs to the
-  // plane part of the boundary
-  const Point<dim> line_center = line->center();
-  if (line_center(0) == this->center(0))
-    return StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
-  else
-    return HyperBallBoundary<dim>::get_intermediate_points_on_line (line, points);
-}
+// template <int dim>
+// void
+// HalfHyperBallBoundary<dim>::
+// get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
+//                                  std::vector<Point<dim> > &points) const
+// {
+//   // check whether center of object is at x==0, since then it belongs to the
+//   // plane part of the boundary
+//   const Point<dim> line_center = line->center();
+//   if (line_center(0) == this->center(0))
+//     return StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
+//   else
+//     return HyperBallBoundary<dim>::get_intermediate_points_on_line (line, points);
+// }
 
 
 
-template <int dim>
-void
-HalfHyperBallBoundary<dim>::
-get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
-                                 std::vector<Point<dim> > &points) const
-{
-  if (points.size()==1)
-    points[0]=get_new_point_on_quad(quad);
-  else
-    {
-      // check whether center of object is at x==0, since then it belongs to
-      // the plane part of the boundary
-      const Point<dim> quad_center = quad->center();
-      if (quad_center(0) == this->center(0))
-        StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
-      else
-        HyperBallBoundary<dim>::get_intermediate_points_on_quad (quad, points);
-    }
-}
+// template <int dim>
+// void
+// HalfHyperBallBoundary<dim>::
+// get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
+//                                  std::vector<Point<dim> > &points) const
+// {
+//   if (points.size()==1)
+//     points[0]=get_new_point_on_quad(quad);
+//   else
+//     {
+//       // check whether center of object is at x==0, since then it belongs to
+//       // the plane part of the boundary
+//       const Point<dim> quad_center = quad->center();
+//       if (quad_center(0) == this->center(0))
+//         StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
+//       else
+//         HyperBallBoundary<dim>::get_intermediate_points_on_quad (quad, points);
+//     }
+// }
 
 
 
-template <>
-void
-HalfHyperBallBoundary<1>::
-get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &,
-                                 std::vector<Point<1> > &) const
-{
-  Assert (false, ExcInternalError());
-}
+// template <>
+// void
+// HalfHyperBallBoundary<1>::
+// get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &,
+//                                  std::vector<Point<1> > &) const
+// {
+//   Assert (false, ExcInternalError());
+// }
 
 
 
-template <>
-void
-HalfHyperBallBoundary<1>::
-get_normals_at_vertices (const Triangulation<1>::face_iterator &,
-                         Boundary<1,1>::FaceVertexNormals &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
+// template <>
+// void
+// HalfHyperBallBoundary<1>::
+// get_normals_at_vertices (const Triangulation<1>::face_iterator &,
+//                          Boundary<1,1>::FaceVertexNormals &) const
+// {
+//   Assert (false, ExcImpossibleInDim(1));
+// }
 
 
 
+// template <int dim>
+// void
+// HalfHyperBallBoundary<dim>::
+// get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+//                          typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
+// {
+//   // check whether center of object is at x==0, since then it belongs to the
+//   // plane part of the boundary
+//   const Point<dim> quad_center = face->center();
+//   if (quad_center(0) == this->center(0))
+//     StraightBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
+//   else
+//     HyperBallBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
+// }
+
+
 template <int dim>
-void
+Point<dim>
 HalfHyperBallBoundary<dim>::
-get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
-                         typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
+normal_vector (const Point<dim> vertex) const
 {
-  // check whether center of object is at x==0, since then it belongs to the
-  // plane part of the boundary
-  const Point<dim> quad_center = face->center();
-  if (quad_center(0) == this->center(0))
-    StraightBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
+  // Check whether the point is inside the flat part, or outside. If
+  // it is both at distance R from the center and with zero first
+  // component, throw an exception, since we are on a edge, and normal
+  // is not well defined.
+  bool at_distance_R = (std::abs(vertex.distance(this->center)-this->radius)<1e-10);
+  
+  if (abs(vertex(0)) < 1e-10 && !at_distance_R)
+    {
+      Point<dim> normal;
+      normal[0] = -1;
+      return normal;
+    }
+  else if(at_distance_R)
+    return vertex-this->center;
   else
-    HyperBallBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
+    {
+      Assert(false,
+            ExcMessage("Normals at edge is not well defined. Use two distinct manifolds!"));
+      return Point<dim>();
+    }
 }
 
 
@@ -1033,241 +1061,257 @@ HalfHyperShellBoundary<dim>::HalfHyperShellBoundary (const Point<dim> &center,
 
 
 
-template <int dim>
-Point<dim>
-HalfHyperShellBoundary<dim>::
-get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
-{
-  switch (dim)
-    {
-      // in 2d, first check whether the two end points of the line are on the
-      // axis of symmetry. if so, then return the mid point
-    case 2:
-    {
-      if ((line->vertex(0)(0) == this->center(0))
-          &&
-          (line->vertex(1)(0) == this->center(0)))
-        return (line->vertex(0) + line->vertex(1))/2;
-      else
-        // otherwise we are on the outer or inner part of the shell. proceed
-        // as in the base class
-        return HyperShellBoundary<dim>::get_new_point_on_line (line);
-    }
-
-    // in 3d, a line is a straight line if it is on the symmetry plane and if
-    // not both of its end points are on either the inner or outer sphere
-    case 3:
-    {
-
-      if (((line->vertex(0)(0) == this->center(0))
-           &&
-           (line->vertex(1)(0) == this->center(0)))
-          &&
-          !(((std::fabs (line->vertex(0).distance (this->center)
-                         - inner_radius) < 1e-12 * outer_radius)
-             &&
-             (std::fabs (line->vertex(1).distance (this->center)
-                         - inner_radius) < 1e-12 * outer_radius))
-            ||
-            ((std::fabs (line->vertex(0).distance (this->center)
-                         - outer_radius) < 1e-12 * outer_radius)
-             &&
-             (std::fabs (line->vertex(1).distance (this->center)
-                         - outer_radius) < 1e-12 * outer_radius))))
-        return (line->vertex(0) + line->vertex(1))/2;
-      else
-        // otherwise we are on the outer or inner part of the shell. proceed
-        // as in the base class
-        return HyperShellBoundary<dim>::get_new_point_on_line (line);
-    }
-
-    default:
-      Assert (false, ExcNotImplemented());
-    }
-
-  return Point<dim>();
-}
-
-
-
-template <>
-Point<1>
-HalfHyperShellBoundary<1>::
-get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
-{
-  Assert (false, ExcInternalError());
-  return Point<1>();
-}
-
+// template <int dim>
+// Point<dim>
+// HalfHyperShellBoundary<dim>::
+// get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const
+// {
+//   switch (dim)
+//     {
+//       // in 2d, first check whether the two end points of the line are on the
+//       // axis of symmetry. if so, then return the mid point
+//     case 2:
+//     {
+//       if ((line->vertex(0)(0) == this->center(0))
+//           &&
+//           (line->vertex(1)(0) == this->center(0)))
+//         return (line->vertex(0) + line->vertex(1))/2;
+//       else
+//         // otherwise we are on the outer or inner part of the shell. proceed
+//         // as in the base class
+//         return HyperShellBoundary<dim>::get_new_point_on_line (line);
+//     }
+
+//     // in 3d, a line is a straight line if it is on the symmetry plane and if
+//     // not both of its end points are on either the inner or outer sphere
+//     case 3:
+//     {
+
+//       if (((line->vertex(0)(0) == this->center(0))
+//            &&
+//            (line->vertex(1)(0) == this->center(0)))
+//           &&
+//           !(((std::fabs (line->vertex(0).distance (this->center)
+//                          - inner_radius) < 1e-12 * outer_radius)
+//              &&
+//              (std::fabs (line->vertex(1).distance (this->center)
+//                          - inner_radius) < 1e-12 * outer_radius))
+//             ||
+//             ((std::fabs (line->vertex(0).distance (this->center)
+//                          - outer_radius) < 1e-12 * outer_radius)
+//              &&
+//              (std::fabs (line->vertex(1).distance (this->center)
+//                          - outer_radius) < 1e-12 * outer_radius))))
+//         return (line->vertex(0) + line->vertex(1))/2;
+//       else
+//         // otherwise we are on the outer or inner part of the shell. proceed
+//         // as in the base class
+//         return HyperShellBoundary<dim>::get_new_point_on_line (line);
+//     }
+
+//     default:
+//       Assert (false, ExcNotImplemented());
+//     }
+
+//   return Point<dim>();
+// }
 
 
 
-template <int dim>
-Point<dim>
-HalfHyperShellBoundary<dim>::
-get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
-{
-  // if this quad is on the symmetry plane, take the center point and project
-  // it outward to the same radius as the centers of the two radial lines
-  if ((quad->vertex(0)(0) == this->center(0)) &&
-      (quad->vertex(1)(0) == this->center(0)) &&
-      (quad->vertex(2)(0) == this->center(0)) &&
-      (quad->vertex(3)(0) == this->center(0)))
-    {
-      const Point<dim> quad_center = (quad->vertex(0) + quad->vertex(1) +
-                                      quad->vertex(2) + quad->vertex(3)   )/4;
-      const Point<dim> quad_center_offset = quad_center - this->center;
-
-
-      if (std::fabs (quad->line(0)->center().distance(this->center) -
-                     quad->line(1)->center().distance(this->center))
-          < 1e-12 * outer_radius)
-        {
-          // lines 0 and 1 are radial
-          const double needed_radius
-            = quad->line(0)->center().distance(this->center);
-
-          return (this->center +
-                  quad_center_offset/quad_center_offset.norm() * needed_radius);
-        }
-      else if (std::fabs (quad->line(2)->center().distance(this->center) -
-                          quad->line(3)->center().distance(this->center))
-               < 1e-12 * outer_radius)
-        {
-          // lines 2 and 3 are radial
-          const double needed_radius
-            = quad->line(2)->center().distance(this->center);
-
-          return (this->center +
-                  quad_center_offset/quad_center_offset.norm() * needed_radius);
-        }
-      else
-        Assert (false, ExcInternalError());
-    }
+// template <>
+// Point<1>
+// HalfHyperShellBoundary<1>::
+// get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const
+// {
+//   Assert (false, ExcInternalError());
+//   return Point<1>();
+// }
 
-  // otherwise we are on the outer or inner part of the shell. proceed as in
-  // the base class
-  return HyperShellBoundary<dim>::get_new_point_on_quad (quad);
-}
 
 
 
-template <int dim>
-void
-HalfHyperShellBoundary<dim>::
-get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
-                                 std::vector<Point<dim> > &points) const
-{
-  switch (dim)
-    {
-      // in 2d, first check whether the two end points of the line are on the
-      // axis of symmetry. if so, then return the mid point
-    case 2:
-    {
-      if ((line->vertex(0)(0) == this->center(0))
-          &&
-          (line->vertex(1)(0) == this->center(0)))
-        StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
-      else
-        // otherwise we are on the outer or inner part of the shell. proceed
-        // as in the base class
-        HyperShellBoundary<dim>::get_intermediate_points_on_line (line, points);
-      break;
-    }
+// template <int dim>
+// Point<dim>
+// HalfHyperShellBoundary<dim>::
+// get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
+// {
+//   // if this quad is on the symmetry plane, take the center point and project
+//   // it outward to the same radius as the centers of the two radial lines
+//   if ((quad->vertex(0)(0) == this->center(0)) &&
+//       (quad->vertex(1)(0) == this->center(0)) &&
+//       (quad->vertex(2)(0) == this->center(0)) &&
+//       (quad->vertex(3)(0) == this->center(0)))
+//     {
+//       const Point<dim> quad_center = (quad->vertex(0) + quad->vertex(1) +
+//                                       quad->vertex(2) + quad->vertex(3)   )/4;
+//       const Point<dim> quad_center_offset = quad_center - this->center;
+
+
+//       if (std::fabs (quad->line(0)->center().distance(this->center) -
+//                      quad->line(1)->center().distance(this->center))
+//           < 1e-12 * outer_radius)
+//         {
+//           // lines 0 and 1 are radial
+//           const double needed_radius
+//             = quad->line(0)->center().distance(this->center);
+
+//           return (this->center +
+//                   quad_center_offset/quad_center_offset.norm() * needed_radius);
+//         }
+//       else if (std::fabs (quad->line(2)->center().distance(this->center) -
+//                           quad->line(3)->center().distance(this->center))
+//                < 1e-12 * outer_radius)
+//         {
+//           // lines 2 and 3 are radial
+//           const double needed_radius
+//             = quad->line(2)->center().distance(this->center);
+
+//           return (this->center +
+//                   quad_center_offset/quad_center_offset.norm() * needed_radius);
+//         }
+//       else
+//         Assert (false, ExcInternalError());
+//     }
+
+//   // otherwise we are on the outer or inner part of the shell. proceed as in
+//   // the base class
+//   return HyperShellBoundary<dim>::get_new_point_on_quad (quad);
+// }
 
-    // in 3d, a line is a straight line if it is on the symmetry plane and if
-    // not both of its end points are on either the inner or outer sphere
-    case 3:
-    {
-      if (((line->vertex(0)(0) == this->center(0))
-           &&
-           (line->vertex(1)(0) == this->center(0)))
-          &&
-          !(((std::fabs (line->vertex(0).distance (this->center)
-                         - inner_radius) < 1e-12 * outer_radius)
-             &&
-             (std::fabs (line->vertex(1).distance (this->center)
-                         - inner_radius) < 1e-12 * outer_radius))
-            ||
-            ((std::fabs (line->vertex(0).distance (this->center)
-                         - outer_radius) < 1e-12 * outer_radius)
-             &&
-             (std::fabs (line->vertex(1).distance (this->center)
-                         - outer_radius) < 1e-12 * outer_radius))))
-        StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
-      else
-        // otherwise we are on the outer or inner part of the shell. proceed
-        // as in the base class
-        HyperShellBoundary<dim>::get_intermediate_points_on_line (line, points);
 
-      break;
-    }
 
-    default:
-      Assert (false, ExcNotImplemented());
-    }
-}
+// template <int dim>
+// void
+// HalfHyperShellBoundary<dim>::
+// get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
+//                                  std::vector<Point<dim> > &points) const
+// {
+//   switch (dim)
+//     {
+//       // in 2d, first check whether the two end points of the line are on the
+//       // axis of symmetry. if so, then return the mid point
+//     case 2:
+//     {
+//       if ((line->vertex(0)(0) == this->center(0))
+//           &&
+//           (line->vertex(1)(0) == this->center(0)))
+//         StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
+//       else
+//         // otherwise we are on the outer or inner part of the shell. proceed
+//         // as in the base class
+//         HyperShellBoundary<dim>::get_intermediate_points_on_line (line, points);
+//       break;
+//     }
+
+//     // in 3d, a line is a straight line if it is on the symmetry plane and if
+//     // not both of its end points are on either the inner or outer sphere
+//     case 3:
+//     {
+//       if (((line->vertex(0)(0) == this->center(0))
+//            &&
+//            (line->vertex(1)(0) == this->center(0)))
+//           &&
+//           !(((std::fabs (line->vertex(0).distance (this->center)
+//                          - inner_radius) < 1e-12 * outer_radius)
+//              &&
+//              (std::fabs (line->vertex(1).distance (this->center)
+//                          - inner_radius) < 1e-12 * outer_radius))
+//             ||
+//             ((std::fabs (line->vertex(0).distance (this->center)
+//                          - outer_radius) < 1e-12 * outer_radius)
+//              &&
+//              (std::fabs (line->vertex(1).distance (this->center)
+//                          - outer_radius) < 1e-12 * outer_radius))))
+//         StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
+//       else
+//         // otherwise we are on the outer or inner part of the shell. proceed
+//         // as in the base class
+//         HyperShellBoundary<dim>::get_intermediate_points_on_line (line, points);
+
+//       break;
+//     }
+
+//     default:
+//       Assert (false, ExcNotImplemented());
+//     }
+// }
 
 
 
-template <int dim>
-void
-HalfHyperShellBoundary<dim>::
-get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
-                                 std::vector<Point<dim> > &points) const
-{
-  Assert (dim < 3, ExcNotImplemented());
+// template <int dim>
+// void
+// HalfHyperShellBoundary<dim>::
+// get_intermediate_points_on_quad (const typename Triangulation<dim>::quad_iterator &quad,
+//                                  std::vector<Point<dim> > &points) const
+// {
+//   Assert (dim < 3, ExcNotImplemented());
+
+//   // check whether center of object is at x==0, since then it belongs to the
+//   // plane part of the boundary
+//   const Point<dim> quad_center = quad->center();
+//   if (quad_center(0) == this->center(0))
+//     StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
+//   else
+//     HyperShellBoundary<dim>::get_intermediate_points_on_quad (quad, points);
+// }
 
-  // check whether center of object is at x==0, since then it belongs to the
-  // plane part of the boundary
-  const Point<dim> quad_center = quad->center();
-  if (quad_center(0) == this->center(0))
-    StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
-  else
-    HyperShellBoundary<dim>::get_intermediate_points_on_quad (quad, points);
-}
 
 
+// template <>
+// void
+// HalfHyperShellBoundary<1>::
+// get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &,
+//                                  std::vector<Point<1> > &) const
+// {
+//   Assert (false, ExcInternalError());
+// }
 
-template <>
-void
-HalfHyperShellBoundary<1>::
-get_intermediate_points_on_quad (const Triangulation<1>::quad_iterator &,
-                                 std::vector<Point<1> > &) const
-{
-  Assert (false, ExcInternalError());
-}
 
 
+// template <>
+// void
+// HalfHyperShellBoundary<1>::
+// get_normals_at_vertices (const Triangulation<1>::face_iterator &,
+//                          Boundary<1,1>::FaceVertexNormals &) const
+// {
+//   Assert (false, ExcImpossibleInDim(1));
+// }
 
-template <>
-void
-HalfHyperShellBoundary<1>::
-get_normals_at_vertices (const Triangulation<1>::face_iterator &,
-                         Boundary<1,1>::FaceVertexNormals &) const
-{
-  Assert (false, ExcImpossibleInDim(1));
-}
 
 
 
 
+// template <int dim>
+// void
+// HalfHyperShellBoundary<dim>::
+// get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
+//                          typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
+// {
+//   if (face->center()(0) == this->center(0))
+//     StraightBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
+//   else
+//     HyperShellBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
+// }
 
 template <int dim>
-void
+Point<dim>
 HalfHyperShellBoundary<dim>::
-get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
-                         typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
+normal_vector (const Point<dim> vertex) const
 {
-  if (face->center()(0) == this->center(0))
-    StraightBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
+  if (vertex(0) == this->center(0)) 
+    {
+      Point<dim> n;
+      n[0] = -1;
+      return n;
+    }
   else
-    HyperShellBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
+    return HyperShellBoundary<dim>::normal_vector (vertex);
 }
 
 
 
 
+
 template <int dim, int spacedim>
 TorusBoundary<dim,spacedim>::TorusBoundary (const double R__,
                                             const double r__)
@@ -1316,7 +1360,7 @@ TorusBoundary<dim,spacedim>::get_correct_angle(const double angle,
 
 template <>
 Point<3>
-TorusBoundary<2,3>::get_real_coord (const Point<2> &surfP) const
+TorusBoundary<2,3>::push_forward (const Point<2> surfP) const
 {
   const double theta=surfP(0);
   const double phi=surfP(1);
@@ -1330,7 +1374,7 @@ TorusBoundary<2,3>::get_real_coord (const Point<2> &surfP) const
 
 template <>
 Point<2>
-TorusBoundary<2,3>::get_surf_coord(const Point<3> &p) const
+TorusBoundary<2,3>::pull_back(const Point<3> p) const
 {
   const double phi=std::asin(std::abs(p(1))/r);
   const double Rr_2=p(0)*p(0)+p(2)*p(2);
@@ -1356,63 +1400,63 @@ TorusBoundary<2,3>::get_surf_coord(const Point<3> &p) const
 
 
 
-template <>
-Point<3>
-TorusBoundary<2,3>::get_new_point_on_line (const Triangulation<2,3>::line_iterator &line) const
-{
-  //Just get the average
-  Point<2>  p0=get_surf_coord(line->vertex(0));
-  Point<2>  p1=get_surf_coord(line->vertex(1));
+// template <>
+// Point<3>
+// TorusBoundary<2,3>::get_new_point_on_line (const Triangulation<2,3>::line_iterator &line) const
+// {
+//   //Just get the average
+//   Point<2>  p0=get_surf_coord(line->vertex(0));
+//   Point<2>  p1=get_surf_coord(line->vertex(1));
 
-  Point<2>  middle(0,0);
+//   Point<2>  middle(0,0);
 
-  //Take care for periodic conditions, For instance phi0= 0, phi1= 3/2*Pi
-  //middle has to be 7/4*Pi not 3/4*Pi. This also works for -Pi/2 + Pi, middle
-  //is 5/4*Pi
-  for (unsigned int i=0; i<2; i++)
-    if (std::abs(p0(i)-p1(i))> numbers::PI)
-      middle(i)=2*numbers::PI;
+//   //Take care for periodic conditions, For instance phi0= 0, phi1= 3/2*Pi
+//   //middle has to be 7/4*Pi not 3/4*Pi. This also works for -Pi/2 + Pi, middle
+//   //is 5/4*Pi
+//   for (unsigned int i=0; i<2; i++)
+//     if (std::abs(p0(i)-p1(i))> numbers::PI)
+//       middle(i)=2*numbers::PI;
 
-  middle+=  p0 + p1;
-  middle*=0.5;
+//   middle+=  p0 + p1;
+//   middle*=0.5;
 
-  Point<3> midReal=get_real_coord(middle);
-  return midReal;
-}
+//   Point<3> midReal=get_real_coord(middle);
+//   return midReal;
+// }
 
 
 
-template <>
-Point<3>
-TorusBoundary<2,3>::get_new_point_on_quad (const Triangulation<2,3>::quad_iterator &quad) const
-{
-  //Just get the average
-  Point<2> p[4];
+// template <>
+// Point<3>
+// TorusBoundary<2,3>::get_new_point_on_quad (const Triangulation<2,3>::quad_iterator &quad) const
+// {
+//   //Just get the average
+//   Point<2> p[4];
 
-  for (unsigned int i=0; i<4; i++)
-    p[i]=get_surf_coord(quad->vertex(i));
+//   for (unsigned int i=0; i<4; i++)
+//     p[i]=get_surf_coord(quad->vertex(i));
 
-  Point<2>  middle(0,0);
+//   Point<2>  middle(0,0);
 
-  //Take care for periodic conditions, see get_new_point_on_line() above
-  //For instance phi0= 0, phi1= 3/2*Pi  middle has to be 7/4*Pi not 3/4*Pi
-  //This also works for -Pi/2 + Pi + Pi- Pi/2, middle is 5/4*Pi
-  for (unsigned int i=0; i<2; i++)
-    for (unsigned int j=1; j<4; j++)
-      {
-        if (std::abs(p[0](i)-p[j](i))> numbers::PI)
-          {
-            middle(i)+=2*numbers::PI;
-          }
-      }
+//   //Take care for periodic conditions, see get_new_point_on_line() above
+//   //For instance phi0= 0, phi1= 3/2*Pi  middle has to be 7/4*Pi not 3/4*Pi
+//   //This also works for -Pi/2 + Pi + Pi- Pi/2, middle is 5/4*Pi
+//   for (unsigned int i=0; i<2; i++)
+//     for (unsigned int j=1; j<4; j++)
+//       {
+//         if (std::abs(p[0](i)-p[j](i))> numbers::PI)
+//           {
+//             middle(i)+=2*numbers::PI;
+//           }
+//       }
 
-  for (unsigned int i=0; i<4; i++)
-    middle+=p[i];
+//   for (unsigned int i=0; i<4; i++)
+//     middle+=p[i];
 
-  middle*= 0.25;
+//   middle*= 0.25;
 
-  return get_real_coord(middle);
-}
+//   return get_real_coord(middle);
+// }
 
 
 
@@ -1440,136 +1484,136 @@ TorusBoundary<2,3>:: get_surf_norm_from_sp(const Point<2> &surfP) const
 //Normal field without unit length
 template <>
 Point<3>
-TorusBoundary<2,3>::get_surf_norm(const Point<3> &p) const
+TorusBoundary<2,3>::normal_vector(const Point<3> p) const
 {
 
-  Point<2> surfP=get_surf_coord(p);
+  Point<2> surfP=pull_back(p);
   return get_surf_norm_from_sp(surfP);
 
 }
 
 
 
-template<>
-void
-TorusBoundary<2,3>::
-get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator   &line,
-                                 std::vector< Point< 3 > > &points) const
-{
-  //Almost the same implementation as StraightBoundary<2,3>
-  unsigned int npoints=points.size();
-  if (npoints==0) return;
-
-  Point<2> p[2];
-
-  for (unsigned int i=0; i<2; i++)
-    p[i]=get_surf_coord(line->vertex(i));
-
-  unsigned int offset[2];
-  offset[0]=0;
-  offset[1]=0;
-
-  //Take care for periodic conditions & negative angles, see
-  //get_new_point_on_line() above. Because we dont have a symmetric
-  //interpolation (just the middle) we need to add 2*Pi to each almost zero
-  //and negative angles.
-  for (unsigned int i=0; i<2; i++)
-    for (unsigned int j=1; j<2; j++)
-      {
-        if (std::abs(p[0](i)-p[j](i))> numbers::PI)
-          {
-            offset[i]++;
-            break;
-          }
-      }
-
-  for (unsigned int i=0; i<2; i++)
-    for (unsigned int j=0; j<2; j++)
-      if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
-        p[j](i)+=2*numbers::PI*offset[i];
-
-
-  Point<2>  target;
-  const std::vector<Point<1> > &line_points = this->get_line_support_points(npoints);
-  for (unsigned int i=0; i<npoints; i++)
-    {
-      const double x = line_points[i+1][0];
-      target=  (1-x)*p[0] + x*p[1];
-      points[i]=get_real_coord(target);
-    }
-}
-
-
+// template<>
+// void
+// TorusBoundary<2,3>::
+// get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator   &line,
+//                                  std::vector< Point< 3 > > &points) const
+// {
+//   //Almost the same implementation as StraightBoundary<2,3>
+//   unsigned int npoints=points.size();
+//   if (npoints==0) return;
+
+//   Point<2> p[2];
+
+//   for (unsigned int i=0; i<2; i++)
+//     p[i]=get_surf_coord(line->vertex(i));
+
+//   unsigned int offset[2];
+//   offset[0]=0;
+//   offset[1]=0;
+
+//   //Take care for periodic conditions & negative angles, see
+//   //get_new_point_on_line() above. Because we dont have a symmetric
+//   //interpolation (just the middle) we need to add 2*Pi to each almost zero
+//   //and negative angles.
+//   for (unsigned int i=0; i<2; i++)
+//     for (unsigned int j=1; j<2; j++)
+//       {
+//         if (std::abs(p[0](i)-p[j](i))> numbers::PI)
+//           {
+//             offset[i]++;
+//             break;
+//           }
+//       }
+
+//   for (unsigned int i=0; i<2; i++)
+//     for (unsigned int j=0; j<2; j++)
+//       if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
+//         p[j](i)+=2*numbers::PI*offset[i];
+
+
+//   Point<2>  target;
+//   const std::vector<Point<1> > &line_points = this->get_line_support_points(npoints);
+//   for (unsigned int i=0; i<npoints; i++)
+//     {
+//       const double x = line_points[i+1][0];
+//       target=  (1-x)*p[0] + x*p[1];
+//       points[i]=get_real_coord(target);
+//     }
+// }
 
-template<>
-void
-TorusBoundary<2,3>::
-get_intermediate_points_on_quad (const Triangulation< 2, 3 >::quad_iterator &quad,
-                                 std::vector< Point< 3 > > &points )const
-{
-  //Almost the same implementation as  StraightBoundary<2,3>
-  const unsigned int n=points.size(),
-                     m=static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
-  // is n a square number
-  Assert(m*m==n, ExcInternalError());
-
-  Point<2>  p[4];
-
-  for (unsigned int i=0; i<4; i++)
-    p[i]=get_surf_coord(quad->vertex(i));
-
-  Point<2>  target;
-  unsigned int offset[2];
-  offset[0]=0;
-  offset[1]=0;
-
-  //Take care for periodic conditions & negative angles, see
-  //get_new_point_on_line() above.  Because we dont have a symmetric
-  //interpolation (just the middle) we need to add 2*Pi to each almost zero
-  //and negative angles.
-  for (unsigned int i=0; i<2; i++)
-    for (unsigned int j=1; j<4; j++)
-      {
-        if (std::abs(p[0](i)-p[j](i))> numbers::PI)
-          {
-            offset[i]++;
-            break;
-          }
-      }
 
-  for (unsigned int i=0; i<2; i++)
-    for (unsigned int j=0; j<4; j++)
-      if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
-        p[j](i)+=2*numbers::PI*offset[i];
 
-  const std::vector<Point<1> > &line_points = this->get_line_support_points(m);
-  for (unsigned int i=0; i<m; ++i)
-    {
-      const double y=line_points[i+1][0];
-      for (unsigned int j=0; j<m; ++j)
-        {
-          const double x=line_points[j+1][0];
-          target=((1-x) * p[0] +
-                  x     * p[1]) * (1-y) +
-                 ((1-x) * p[2] +
-                  x     * p[3]) * y;
-
-          points[i*m+j]=get_real_coord(target);
-        }
-    }
-}
+// template<>
+// void
+// TorusBoundary<2,3>::
+// get_intermediate_points_on_quad (const Triangulation< 2, 3 >::quad_iterator &quad,
+//                                  std::vector< Point< 3 > > &points )const
+// {
+//   //Almost the same implementation as  StraightBoundary<2,3>
+//   const unsigned int n=points.size(),
+//                      m=static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
+//   // is n a square number
+//   Assert(m*m==n, ExcInternalError());
+
+//   Point<2>  p[4];
+
+//   for (unsigned int i=0; i<4; i++)
+//     p[i]=get_surf_coord(quad->vertex(i));
+
+//   Point<2>  target;
+//   unsigned int offset[2];
+//   offset[0]=0;
+//   offset[1]=0;
+
+//   //Take care for periodic conditions & negative angles, see
+//   //get_new_point_on_line() above.  Because we dont have a symmetric
+//   //interpolation (just the middle) we need to add 2*Pi to each almost zero
+//   //and negative angles.
+//   for (unsigned int i=0; i<2; i++)
+//     for (unsigned int j=1; j<4; j++)
+//       {
+//         if (std::abs(p[0](i)-p[j](i))> numbers::PI)
+//           {
+//             offset[i]++;
+//             break;
+//           }
+//       }
+
+//   for (unsigned int i=0; i<2; i++)
+//     for (unsigned int j=0; j<4; j++)
+//       if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
+//         p[j](i)+=2*numbers::PI*offset[i];
+
+//   const std::vector<Point<1> > &line_points = this->get_line_support_points(m);
+//   for (unsigned int i=0; i<m; ++i)
+//     {
+//       const double y=line_points[i+1][0];
+//       for (unsigned int j=0; j<m; ++j)
+//         {
+//           const double x=line_points[j+1][0];
+//           target=((1-x) * p[0] +
+//                   x     * p[1]) * (1-y) +
+//                  ((1-x) * p[2] +
+//                   x     * p[3]) * y;
+
+//           points[i*m+j]=get_real_coord(target);
+//         }
+//     }
+// }
 
 
 
-template<>
-void
-TorusBoundary<2,3>::
-get_normals_at_vertices (const Triangulation<2,3 >::face_iterator &face,
-                         Boundary<2,3>::FaceVertexNormals &face_vertex_normals) const
-{
-  for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_face; i++)
-    face_vertex_normals[i]=get_surf_norm(face->vertex(i));
-}
+// template<>
+// void
+// TorusBoundary<2,3>::
+// get_normals_at_vertices (const Triangulation<2,3 >::face_iterator &face,
+//                          Boundary<2,3>::FaceVertexNormals &face_vertex_normals) const
+// {
+//   for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_face; i++)
+//     face_vertex_normals[i]=get_surf_norm(face->vertex(i));
+// }
 
 
 

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.