Reverted back the way HyperBallBoundary was creating intermediate points.

git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32261 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 971558bcef9b0612d1950b79a19ec876bd167e7d..ee1a6319b4a324c8ef26e97e132fb51dfd100b34 100644 (file)
--- a/deal.II/include/deal.II/grid/tria_boundary_lib.h
+++ b/deal.II/include/deal.II/grid/tria_boundary_lib.h
@@ -312,6 +312,19 @@ protected:
    * called.
    */
   bool compute_radius_automatically;
+
+
+private:
+
+  /**                                                                           
+   * Called by @p get_new_point for special cases where surrounding
+   * points are two. Returns the point which is at xi between
+   * p0 and p1, where xi is between 0 and 1. 
+   */
+  Point<spacedim> get_intermediate_point (const Point<spacedim> &p0, 
+                                         const Point<spacedim> &p1,
+                                         const double xi) const;
+
 };
 
 

diff --git a/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/source/grid/tria_boundary_lib.cc
index 42c35d8b0c4517a14e5590ae01956da2245ca966..e0f2017682d4adca4d7d78b62f2ab4d94ea103be 100644 (file)
--- a/deal.II/source/grid/tria_boundary_lib.cc
+++ b/deal.II/source/grid/tria_boundary_lib.cc
@@ -183,21 +183,87 @@ 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);
-  middle -= center;
+  // Special cases where the points are vertices in one or two dimensions
+  switch(surrounding_points.size()) {
+  case 2:
+    return get_intermediate_point(surrounding_points[0], surrounding_points[1], weights[1]);
+    break;
+  case 4:
+    {
+      Point<spacedim> v[2];
+      // Coordinates of the point in the reference element
+      double xi1 = weights[1]+weights[3]; // = x
+      double xi2 = weights[2]+weights[3]; // = y
+      
+      v[0] =  get_intermediate_point(surrounding_points[0], 
+                                    surrounding_points[1],
+                                    xi1);
+      v[1] =  get_intermediate_point(surrounding_points[2], 
+                                    surrounding_points[3],
+                                    xi1);
+      return get_intermediate_point(v[0], v[1], xi2);
+    }
+    break;
+  default: 
+    {
+      Point<spacedim> middle = FlatManifold<spacedim>::get_new_point(surrounding_points, weights);
+      middle -= center;
+      
+      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;
+    }
+    break;
+  }
+}
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+HyperBallBoundary<dim,spacedim>::get_intermediate_point(const Point<spacedim> &p0,
+                                                       const Point<spacedim> &p1,
+                                                       const double xi) const {
   
+  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 = surrounding_points[0] - center;
+      const Point<spacedim> vertex_relative = p0 - center;
       r = std::sqrt(vertex_relative.square());
     }
   else
     r=radius;
-  // project to boundary
-  middle *= r / std::sqrt(middle.square());
-  middle += center;
-  return middle;
+
+  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());
+
+  const double beta = alpha * (xi-0.5);
+  const double d=h*std::tan(beta);
+  Point<spacedim> p=pm+d/length*(v1-v0);
+  
+  p *= r / std::sqrt(p.square());
+  p += center;
+
+  return p;
 }