Backward compatibility for Cylinder and Ball. Still missing the Half-Thingy family.

git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@31888 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/include/deal.II/grid/tria_boundary_lib.h b/deal.II/include/deal.II/grid/tria_boundary_lib.h
index 20e2b5416f0679bfefaef8c3582db6717d29586b..54ad9af979a156e1e9f0f7467548d36b0922c522 100644 (file)
--- a/deal.II/include/deal.II/grid/tria_boundary_lib.h
+++ b/deal.II/include/deal.II/grid/tria_boundary_lib.h
@@ -92,6 +92,15 @@ public:
 //   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
+   * boundary to function properly.
+   */
+    virtual Point<spacedim>
+    project_to_manifold (const Point<spacedim> candidate) const;
+    
+    
   /**
    * Refer to the general documentation of
    * this class and the documentation of the
@@ -270,13 +279,12 @@ public:
   get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
 
   /**
-   * Refer to the general documentation of
-   * this class and the documentation of the
-   * base class.
+   * Since this class is derived from FlatManifold, we need to
+   * overload only the project_to_manifold function in order for the
+   * boundary to function properly.
    */
-//   virtual Point<dim>
-//   get_new_point (const std::vector<Point<dim> > &surrounding_points,
-//              const std::vector<double> &weights) const;
+    virtual Point<dim>
+    project_to_manifold (const Point<dim> candidate) const;
 
 
   /**
@@ -413,16 +421,17 @@ public:
   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;
-
+    virtual Point<spacedim>
+    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
@@ -599,15 +608,18 @@ public:
   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
+   * 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;
-
+    virtual
+    Point<dim>
+    get_new_point(const std::vector<Point<dim> > &surrounding_points,
+                 const std::vector<double> &weights) const;
+      
     
   /**
    * Refer to the general
@@ -719,16 +731,7 @@ public:
    */
   virtual Point<dim>
   get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) const;
-  /**
-   * Refer to the general documentation of
-   * this class and the documentation of the
-   * base class.
-   */
-//   virtual Point<dim>
-//   get_new_point (const std::vector<Point<dim> > &surrounding_points,
-//              const std::vector<double> &weights) const;
-
-
+    
   /**
    * Construct a new point on a quad.
    */

diff --git a/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/source/grid/tria_boundary_lib.cc
index dfe6410b08e3399c18a4b8aacb6d18a130420568..9bd26fb208539c7d6c346baae9a5b2541bbf329d 100644 (file)
--- a/deal.II/source/grid/tria_boundary_lib.cc
+++ b/deal.II/source/grid/tria_boundary_lib.cc
@@ -87,6 +87,31 @@ get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator
             point_on_axis);
 }
 
+template <int dim, int spacedim>
+Point<spacedim>
+CylinderBoundary<dim,spacedim>::
+project_to_manifold (const Point<spacedim> middle) const
+{
+  // The temptative point is computed from the get_new_point function
+  // of the FlatManifold class, which in turns calls this class
+  // internally. We have 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<>
@@ -343,6 +368,21 @@ get_new_point_on_line (const typename Triangulation<dim>::line_iterator &line) c
 }
 
 
+template<int dim>
+Point<dim>
+ConeBoundary<dim>::
+project_to_manifold (const Point<dim> middle) const
+{
+  const Point<dim> axis = x_1 - x_0;
+  // 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>
@@ -527,6 +567,28 @@ HyperBallBoundary<dim,spacedim>::get_new_point_on_line (const typename Triangula
 }
 
 
+template <int dim, int spacedim>
+Point<spacedim>
+HyperBallBoundary<dim,spacedim>::get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
+                                              const std::vector<double> &weights) const
+{
+  Point<spacedim> middle = FlatManifold<spacedim>::get_new_point(surrounding_points, weights);
+
+  double r=0;
+  if (compute_radius_automatically)
+    {
+      const Point<spacedim> vertex_relative = surrounding_points[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 <>
 Point<1>
@@ -807,28 +869,28 @@ 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
+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
-  const Point<dim> line_center = line->center();
-  const Point<dim> vertices[2] = { line->vertex(0), line->vertex(1) };
-
-  if ((line_center(0) == this->center(0))
-      &&
-      ((std::fabs(vertices[0].distance(this->center)-this->radius) >
-        1e-5*this->radius)
-       ||
-       (std::fabs(vertices[1].distance(this->center)-this->radius) >
-        1e-5*this->radius)))
-    return line_center;
+  const Point<dim> object_center = FlatManifold<dim>::get_new_point(surrounding_points, weights);
+  if (object_center.distance(this->center) < 1e-5)
+    return object_center;
   else
-    return HyperBallBoundary<dim>::get_new_point_on_line (line);
+    return HyperBallBoundary<dim>::get_new_point (surrounding_points,weights);
 }