From cccdb93577ccbcc98bf88a0b035daceb88143d63 Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@sissa.it>
Date: Wed, 23 Mar 2016 16:26:07 +0100
Subject: [PATCH] Use Manifold, not boundary, functionality.

---
 source/fe/mapping_c1.cc       | 28 ++++++++++++++++------------
 source/grid/grid_generator.cc |  8 ++++----
 source/grid/tria.cc           | 20 ++++++++++++++------
 3 files changed, 34 insertions(+), 22 deletions(-)

diff --git a/source/fe/mapping_c1.cc b/source/fe/mapping_c1.cc
index 216d11aeb2..f20482269e 100644
--- a/source/fe/mapping_c1.cc
+++ b/source/fe/mapping_c1.cc
@@ -66,7 +66,9 @@ MappingC1<2>::MappingC1Generic::add_line_support_points (const Triangulation<2>:
                                                          std::vector<Point<2> > &a) const
 {
   const unsigned int dim = 2;
-  std::vector<Point<dim> > line_points (2);
+  const std::array<double, 2> interior_gl_points
+  {{0.5 - 0.5*std::sqrt(1.0/5.0), 0.5 + 0.5*std::sqrt(1.0/5.0)}};
+
 
   // loop over each of the lines, and if it is at the boundary, then first get
   // the boundary description and second compute the points on it. if not at
@@ -79,11 +81,10 @@ MappingC1<2>::MappingC1Generic::add_line_support_points (const Triangulation<2>:
         {
           // first get the normal vectors at the two vertices of this line
           // from the boundary description
-          const Boundary<dim> &boundary
-            = line->get_triangulation().get_boundary(line->boundary_id());
+          const Manifold<dim> &manifold = line->get_manifold();
 
-          Boundary<dim>::FaceVertexNormals face_vertex_normals;
-          boundary.get_normals_at_vertices (line, face_vertex_normals);
+          Manifold<dim>::FaceVertexNormals face_vertex_normals;
+          manifold.get_normals_at_vertices (line, face_vertex_normals);
 
           // then transform them into interpolation points for a cubic
           // polynomial
@@ -124,9 +125,8 @@ MappingC1<2>::MappingC1Generic::add_line_support_points (const Triangulation<2>:
                              -face_vertex_normals[1][0] * std::sin(alpha)))
                            -2*c;
 
-          QGaussLobatto<1> quad_points(4);
-          const double t1 = quad_points.point(1)[0];
-          const double t2 = quad_points.point(2)[0];
+          const double t1 = interior_gl_points[0];
+          const double t2 = interior_gl_points[1];
           const double s_t1 = (((-b-c)*t1+b)*t1+c)*t1;
           const double s_t2 = (((-b-c)*t2+b)*t2+c)*t2;
 
@@ -149,11 +149,15 @@ MappingC1<2>::MappingC1Generic::add_line_support_points (const Triangulation<2>:
             };
         }
       else
-        // not at boundary
+        // not at boundary, so just use scaled Gauss-Lobatto points (i.e., use
+        // plain straight lines).
         {
-          static const StraightBoundary<dim> straight_boundary;
-          straight_boundary.get_intermediate_points_on_line (line, line_points);
-          a.insert (a.end(), line_points.begin(), line_points.end());
+          // Note that the zeroth Gauss-Lobatto point is a boundary point, so
+          // we push back mapped versions of the first and second.
+          a.push_back((1.0 - interior_gl_points[0])*line->vertex(0) +
+                      (interior_gl_points[0]*line->vertex(1)));
+          a.push_back((1.0 - interior_gl_points[1])*line->vertex(0) +
+                      (interior_gl_points[1]*line->vertex(1)));
         }
     }
 }
diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index f839683775..9307cf69f5 100644
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -16,10 +16,10 @@
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_reordering.h>
 #include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_accessor.h>
 #include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
 #include <deal.II/grid/intergrid_map.h>
 
 #include <deal.II/distributed/tria.h>
@@ -3484,11 +3484,11 @@ namespace GridGenerator
         // once and then re-arrange all
         // interior nodes so that the mesh is
         // the least distorted
-        HyperShellBoundary<3> boundary (p);
+        SphericalManifold<3> boundary (p);
         Triangulation<3> tmp;
         hyper_shell (tmp, p, inner_radius, outer_radius, 12);
-        tmp.set_boundary(0, boundary);
-        tmp.set_boundary(1, boundary);
+        tmp.set_manifold(0, boundary);
+        tmp.set_manifold(1, boundary);
         tmp.refine_global (1);
 
         // let's determine the distance at
diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index 5bdb643229..d7fdfb3687 100644
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -8932,6 +8932,16 @@ namespace internal
         return true;
       }
     };
+
+
+
+    template <int dim, int spacedim>
+    const Manifold<dim, spacedim> &
+    get_default_flat_manifold()
+    {
+      static const FlatManifold<dim, spacedim> flat_manifold;
+      return flat_manifold;
+    }
   }
 }
 
@@ -9219,12 +9229,10 @@ Triangulation<dim, spacedim>::get_manifold (const types::manifold_id m_number) c
       //if we have found an entry, return it
       return *(it->second);
     }
-  else
-    {
-      //if we have not found an entry connected with number, we return
-      //straight_boundary
-      return straight_boundary;
-    }
+
+  // if we have not found an entry connected with number, we return
+  // the default (flat) manifold
+  return internal::Triangulation::get_default_flat_manifold<dim,spacedim>();
 }
 
 
-- 
2.39.5