From: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date: Fri, 19 Aug 2016 23:16:36 +0000 (-0600)
Subject: Reduce use of Quadrature objects in Manifolds.
X-Git-Tag: v8.5.0-rc1~708^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3645938c6efe68585b3b32cfe4be6b73484824de;p=dealii.git

Reduce use of Quadrature objects in Manifolds.
---

diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h
index 5f3b7a8ee8..52646e054d 100644
--- a/include/deal.II/grid/manifold.h
+++ b/include/deal.II/grid/manifold.h
@@ -79,7 +79,54 @@ namespace Manifolds
   template <typename MeshIteratorType>
   Quadrature<MeshIteratorType::AccessorType::space_dimension>
   get_default_quadrature(const MeshIteratorType &iterator,
-                         const bool              with_laplace = false);
+                         const bool              with_laplace = false) DEAL_II_DEPRECATED;
+
+  /**
+   * Given a general mesh iterator, construct vectors of quadrature points and
+   * weights that contain the following points:
+   * - If the iterator points to a line, then the quadrature points
+   *   are the two vertices of the line. This results in a point vector
+   *   with two points.
+   * - If the iterator points to a quad, then the quadrature points
+   *   are the vertices and line mid-points. This results in a point vector
+   *   with eight (4+4) points.
+   * - If the iterator points to a hex, then the quadrature points
+   *   are the vertices, the line mid-points, and the face mid-points.
+   *   This results in a points vector with 26 (8+12+6) points.
+   *
+   * The quadrature weights for these points are either chosen identically
+   * and equal to one over the number of quadrature points (if @p with_laplace
+   * is @p false), or in a way that gives points closer to the cell center
+   * (measured on the reference cell) a higher weight. These weights correspond
+   * to solving a Laplace equation and evaluating the solution at the quadrature
+   * points (if @p with_laplace is @p true).
+   *
+   * The function is primarily used to construct the input argument
+   * for the Manifold::get_new_point() function, which computes a new
+   * point on a manifold based on a weighted average of "surrounding"
+   * points represented by the quadrature points and weights stored in the
+   * returned pair of vectors. This function creates such an object based on
+   * the points that "surround" a cell, face, or edge, and weights
+   * are chosen in a way appropriate for computing the new "mid-point"
+   * of the object pointed to. An example of where this is necessary
+   * is for mesh refinement, where (using the 2d situation as an example)
+   * we need to first create new edge mid-points, and then a new cell-point.
+   *
+   * @param[in] iterator A mesh iterator that points to either a line, quad,
+   *   or hex.
+   * @param[in] with_laplace Whether or not to compute the quadrature weights
+   *   by solving a Laplace equation, as discussed above.
+   * @tparam MeshIteratorType An iterator type that corresponds to either
+   *   Triangulation::cell_iterator (or variants such as
+   *   Triangulation::active_cell_iterator or DoFHandler::cell_iterator) or
+   *   that is the result of statements such as
+   *   <code>cell-@>face(f)</code> or <code>cell-@>line(l)</code>.
+   */
+  template <typename MeshIteratorType>
+  std::pair<std::vector<Point<MeshIteratorType::AccessorType::space_dimension> >,
+      std::vector<double> >
+      get_default_points_and_weights(const MeshIteratorType &iterator,
+                                     const bool              with_laplace = false);
 }
 
 
@@ -318,7 +365,28 @@ public:
    */
   virtual
   Point<spacedim>
-  get_new_point (const Quadrature<spacedim> &quad) const;
+  get_new_point (const Quadrature<spacedim> &quad) const DEAL_II_DEPRECATED;
+
+  /**
+   * Return the point which shall become the new vertex surrounded by the
+   * given points @p surrounding_points. @p weights contains appropriate
+   * weights for the surrounding points according to which the manifold
+   * determines the new point's position.
+   *
+   * In its default implementation it uses a pair-wise reduction of
+   * the points by calling the function get_intermediate_point() on the first
+   * two points, then on the resulting point and the next, until all points in
+   * the vector have been taken into account. User classes can get away by
+   * simply implementing the get_intermediate_point() function. Notice that
+   * by default the get_intermediate_point() function calls the
+   * project_to_manifold() function with the convex combination of its
+   * arguments. For simple situations you may get away by implementing
+   * only the project_to_manifold() function.
+   */
+  virtual
+  Point<spacedim>
+  get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+                 const std::vector<double>           &weights) const;
 
   /**
    * Given a point which lies close to the given manifold, it modifies it and
@@ -589,7 +657,33 @@ public:
    */
   virtual
   Point<spacedim>
-  get_new_point(const Quadrature<spacedim> &quad) const;
+  get_new_point(const Quadrature<spacedim> &quad) const DEAL_II_DEPRECATED;
+
+  /**
+   * Let the new point be the average sum of surrounding vertices.
+   *
+   * This particular implementation constructs the weighted average of the
+   * surrounding points, and then calls internally the function
+   * project_to_manifold(). The reason why we do it this way, is to allow lazy
+   * programmers to implement only the project_to_manifold() function for their
+   * own Manifold classes which are small (or trivial) perturbations of a flat
+   * manifold. This is the case whenever the coarse mesh is a decent
+   * approximation of the manifold geometry. In this case, the middle point of
+   * a cell is close to true middle point of the manifold, and a projection
+   * may suffice.
+   *
+   * For most simple geometries, it is possible to get reasonable results by
+   * deriving your own Manifold class from FlatManifold, and write a new
+   * interface only for the project_to_manifold function. You will have good
+   * approximations also with large deformations, as long as in the coarsest
+   * mesh size you are trying to refine, the middle point is not too far from
+   * the manifold mid point, i.e., as long as the coarse mesh size is small
+   * enough.
+   */
+  virtual
+  Point<spacedim>
+  get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
+                const std::vector<double>           &weights) const;
 
 
   /**
@@ -784,7 +878,16 @@ public:
    */
   virtual
   Point<spacedim>
-  get_new_point(const Quadrature<spacedim> &quad) const;
+  get_new_point(const Quadrature<spacedim> &quad) const DEAL_II_DEPRECATED;
+
+  /**
+   * Refer to the general documentation of this class and the documentation of
+   * the base class for more information.
+   */
+  virtual
+  Point<spacedim>
+  get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
+                const std::vector<double>           &weights) const;
 
   /**
    * Pull back the given point in spacedim to the Euclidean chartdim
@@ -953,12 +1056,25 @@ namespace Manifolds
   Quadrature<MeshIteratorType::AccessorType::space_dimension>
   get_default_quadrature(const MeshIteratorType &iterator,
                          const bool              with_laplace)
+  {
+    const std::pair<std::vector<Point<MeshIteratorType::AccessorType::space_dimension> >,
+          std::vector<double> > points_and_weights = get_default_points_and_weights(iterator,
+                                                     with_laplace);
+    return Quadrature<MeshIteratorType::AccessorType::space_dimension>(points_and_weights.first,
+           points_and_weights.second);
+  }
+
+  template <typename MeshIteratorType>
+  std::pair<std::vector<Point<MeshIteratorType::AccessorType::space_dimension> >,
+      std::vector<double> >
+      get_default_points_and_weights(const MeshIteratorType &iterator,
+                                     const bool              with_laplace)
   {
     const int spacedim = MeshIteratorType::AccessorType::space_dimension;
     const int dim = MeshIteratorType::AccessorType::structure_dimension;
 
-    std::vector<Point<spacedim> > sp;
-    std::vector<double> wp;
+    std::pair<std::vector<Point<spacedim> >,
+        std::vector<double> > points_weights;
 
 
     // note that the exact weights are chosen such as to minimize the
@@ -968,32 +1084,32 @@ namespace Manifolds
     switch (dim)
       {
       case 1:
-        sp.resize(2);
-        wp.resize(2);
-        sp[0] = iterator->vertex(0);
-        wp[0] = .5;
-        sp[1] = iterator->vertex(1);
-        wp[1] = .5;
+        points_weights.first.resize(2);
+        points_weights.second.resize(2);
+        points_weights.first[0] = iterator->vertex(0);
+        points_weights.second[0] = .5;
+        points_weights.first[1] = iterator->vertex(1);
+        points_weights.second[1] = .5;
         break;
       case 2:
-        sp.resize(8);
-        wp.resize(8);
+        points_weights.first.resize(8);
+        points_weights.second.resize(8);
 
         for (unsigned int i=0; i<4; ++i)
           {
-            sp[i] = iterator->vertex(i);
-            sp[4+i] = ( iterator->line(i)->has_children() ?
-                        iterator->line(i)->child(0)->vertex(1) :
-                        iterator->line(i)->get_manifold().get_new_point_on_line(iterator->line(i)) );
+            points_weights.first[i] = iterator->vertex(i);
+            points_weights.first[4+i] = ( iterator->line(i)->has_children() ?
+                                          iterator->line(i)->child(0)->vertex(1) :
+                                          iterator->line(i)->get_manifold().get_new_point_on_line(iterator->line(i)) );
           }
 
         if (with_laplace)
           {
-            std::fill(wp.begin(), wp.begin()+4, 1.0/16.0);
-            std::fill(wp.begin()+4, wp.end(), 3.0/16.0);
+            std::fill(points_weights.second.begin(), points_weights.second.begin()+4, 1.0/16.0);
+            std::fill(points_weights.second.begin()+4, points_weights.second.end(), 3.0/16.0);
           }
         else
-          std::fill(wp.begin(), wp.end(), 1.0/8.0);
+          std::fill(points_weights.second.begin(), points_weights.second.end(), 1.0/8.0);
         break;
       case 3:
       {
@@ -1003,9 +1119,9 @@ namespace Manifolds
           GeometryInfo<dim>::vertices_per_cell+
           GeometryInfo<dim>::lines_per_cell+
           GeometryInfo<dim>::faces_per_cell;
-        sp.resize(np);
-        wp.resize(np);
-        std::vector<Point<3> > *sp3 = reinterpret_cast<std::vector<Point<3> > *>(&sp);
+        points_weights.first.resize(np);
+        points_weights.second.resize(np);
+        std::vector<Point<3> > *sp3 = reinterpret_cast<std::vector<Point<3> > *>(&points_weights.first);
 
         unsigned int j=0;
 
@@ -1016,33 +1132,33 @@ namespace Manifolds
         for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i, ++j)
           {
             (*sp3)[j] = hex->vertex(i);
-            wp[j] = 1.0/128.0;
+            points_weights.second[j] = 1.0/128.0;
           }
         for (unsigned int i=0; i<GeometryInfo<dim>::lines_per_cell; ++i, ++j)
           {
             (*sp3)[j] = (hex->line(i)->has_children() ?
                          hex->line(i)->child(0)->vertex(1) :
                          hex->line(i)->get_manifold().get_new_point_on_line(hex->line(i)));
-            wp[j] = 7.0/192.0;
+            points_weights.second[j] = 7.0/192.0;
           }
         for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i, ++j)
           {
             (*sp3)[j] = (hex->quad(i)->has_children() ?
                          hex->quad(i)->isotropic_child(0)->vertex(3) :
                          hex->quad(i)->get_manifold().get_new_point_on_quad(hex->quad(i)));
-            wp[j] = 1.0/12.0;
+            points_weights.second[j] = 1.0/12.0;
           }
         // Overwrite the weights with 1/np if we don't want to use
         // laplace vectors.
         if (with_laplace == false)
-          std::fill(wp.begin(), wp.end(), 1.0/np);
+          std::fill(points_weights.second.begin(), points_weights.second.end(), 1.0/np);
       }
       break;
       default:
         Assert(false, ExcInternalError());
         break;
       }
-    return Quadrature<spacedim>(sp,wp);
+    return points_weights;
   }
 }
 
diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h
index 40a383aea4..f9583582bd 100644
--- a/include/deal.II/grid/manifold_lib.h
+++ b/include/deal.II/grid/manifold_lib.h
@@ -289,7 +289,15 @@ public:
    * the base class for a detailed description of what this function does.
    */
   virtual Point<spacedim>
-  get_new_point(const Quadrature<spacedim> &quad) const;
+  get_new_point(const Quadrature<spacedim> &quad) const DEAL_II_DEPRECATED;
+
+  /**
+   * Compute new points on the CylindricalManifold. See the documentation of
+   * the base class for a detailed description of what this function does.
+   */
+  virtual Point<spacedim>
+  get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
+                const std::vector<double>           &weights) const;
 
 protected:
   /**
diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc
index 481155db61..e0b7e4e1ec 100644
--- a/source/fe/mapping_manifold.cc
+++ b/source/fe/mapping_manifold.cc
@@ -219,7 +219,7 @@ transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_it
       vertices.push_back(cell->vertex(v));
       weights.push_back(GeometryInfo<dim>::d_linear_shape_function(p,v));
     }
-  return cell->get_manifold().get_new_point(Quadrature<spacedim>(vertices, weights));
+  return cell->get_manifold().get_new_point(vertices, weights);
 }
 
 
@@ -403,8 +403,8 @@ namespace internal
           for (unsigned int point=0; point<quadrature_points.size(); ++point)
             {
               quadrature_points[point] = data.manifold->
-                                         get_new_point(Quadrature<spacedim>(data.vertices,
-                                                       data.cell_manifold_quadrature_weights[point+data_set]));
+                                         get_new_point(data.vertices,
+                                                       data.cell_manifold_quadrature_weights[point+data_set]);
             }
         }
     }
@@ -440,8 +440,8 @@ namespace internal
 
               // And get its image on the manifold:
               const Point<spacedim> P = data.manifold->
-                                        get_new_point(Quadrature<spacedim>(data.vertices,
-                                                      data.cell_manifold_quadrature_weights[point+data_set]));
+                                        get_new_point(data.vertices,
+                                                      data.cell_manifold_quadrature_weights[point+data_set]);
 
               // To compute the Jacobian, we choose dim points aligned
               // with the dim reference axes, which are still in the
@@ -469,8 +469,8 @@ namespace internal
                     data.vertex_weights[j] = GeometryInfo<dim>::d_linear_shape_function(np, j);
 
                   const Point<spacedim> NP=
-                    data.manifold->get_new_point(Quadrature<spacedim>(data.vertices,
-                                                                      data.vertex_weights));
+                    data.manifold->get_new_point(data.vertices,
+                                                 data.vertex_weights);
 
                   const Tensor<1,spacedim> T = data.manifold->get_tangent_vector(P, NP);
 
diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
index 317dc0b5c1..5de96793cc 100644
--- a/source/fe/mapping_q_generic.cc
+++ b/source/fe/mapping_q_generic.cc
@@ -3484,8 +3484,7 @@ namespace
             const double x = line_support_points.point(i+1)[0];
             w[1] = x;
             w[0] = (1-x);
-            Quadrature<spacedim> quadrature(surrounding_points, w);
-            points[i] = manifold.get_new_point(quadrature);
+            points[i] = manifold.get_new_point(surrounding_points, w);
           }
         break;
       }
@@ -3511,8 +3510,7 @@ namespace
                 for (unsigned int l=0; l<4; ++l)
                   w[l] = GeometryInfo<2>::d_linear_shape_function(p, l);
 
-                Quadrature<spacedim> quadrature(surrounding_points, w);
-                points[c]=manifold.get_new_point(quadrature);
+                points[c]=manifold.get_new_point(surrounding_points, w);
               }
           }
         break;
@@ -3542,8 +3540,7 @@ namespace
                     for (unsigned int l=0; l<8; ++l)
                       w[l] = GeometryInfo<3>::d_linear_shape_function(p, l);
 
-                    Quadrature<spacedim> quadrature(surrounding_points, w);
-                    points[c]=manifold.get_new_point(quadrature);
+                    points[c]=manifold.get_new_point(surrounding_points, w);
                   }
               }
           }
diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc
index b318665216..d4f9a9fe9f 100644
--- a/source/grid/manifold.cc
+++ b/source/grid/manifold.cc
@@ -59,29 +59,67 @@ template <int dim, int spacedim>
 Point<spacedim>
 Manifold<dim, spacedim>::
 get_new_point (const Quadrature<spacedim> &quad) const
+{
+  return get_new_point(quad.get_points(),quad.get_weights());
+}
+
+template <int dim, int spacedim>
+Point<spacedim>
+Manifold<dim, spacedim>::
+get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+               const std::vector<double>           &weights) const
 {
   const double tol = 1e-10;
+  const unsigned int n_points = surrounding_points.size();
 
-  Assert(quad.size() > 0,
-         ExcMessage("Quadrature should have at least one point."));
+  Assert(n_points > 0,
+         ExcMessage("There should be at least one point."));
 
-  Assert(std::abs(std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.0)-1.0) < tol,
+  Assert(n_points == weights.size(),
+         ExcMessage("There should be as many surrounding points as weights given."));
+
+  Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0)-1.0) < tol,
          ExcMessage("The weights for the individual points should sum to 1!"));
 
-  QSorted<spacedim> sorted_quad(quad);
-  Point<spacedim> p = sorted_quad.point(0);
-  double w = sorted_quad.weight(0);
+  // The number of points is small (at most 26 in 3D), therefore try to minize
+  // overhead and use a very simple search of O(N^2) to find the order in which
+  // to loop over points.
+  std::vector<unsigned int> permutation(n_points);
+  std::vector<bool> found(n_points,false);
+
+  unsigned int min_index = numbers::invalid_unsigned_int;
+  for (unsigned int i=0; i < n_points; ++i)
+    {
+      double min_weight = std::numeric_limits<double>::max();
+
+      for (unsigned int j=0; j < n_points; ++j)
+        {
+          if ((!found[j]) && (weights[j] < min_weight))
+            {
+              min_weight = weights[j];
+              min_index = j;
+            }
+        }
+      permutation[i] = min_index;
+      found[min_index] = true;
+    }
+
+  // Now loop over points in the order of their weights. This is done to
+  // produce unique points even if get_intermediate_points is not
+  // associative (as for the SphericalManifold).
+  Point<spacedim> p = surrounding_points[permutation[0]];
+  double w = weights[permutation[0]];
 
-  for (unsigned int i=1; i<sorted_quad.size(); ++i)
+  for (unsigned int i=1; i<n_points; ++i)
     {
       double weight = 0.0;
-      if ( (sorted_quad.weight(i) + w) < tol )
+      if ( (weights[permutation[i]] + w) < tol )
         weight = 0.0;
       else
-        weight =  w/(sorted_quad.weight(i) + w);
+        weight =  w/(weights[permutation[i]] + w);
 
-      p = get_intermediate_point(p, sorted_quad.point(i),1.0 - weight );
-      w += sorted_quad.weight(i);
+      p = get_intermediate_point(p, surrounding_points[permutation[i]],1.0 - weight );
+      w += weights[permutation[i]];
     }
 
   return p;
@@ -274,7 +312,8 @@ Point<spacedim>
 Manifold<dim, spacedim>::
 get_new_point_on_line (const typename Triangulation<dim, spacedim>::line_iterator &line) const
 {
-  return get_new_point (get_default_quadrature(line));
+  const std::pair<std::vector<Point<spacedim> >, std::vector<double> > points_weights(get_default_points_and_weights(line));
+  return get_new_point (points_weights.first,points_weights.second);
 }
 
 
@@ -284,7 +323,8 @@ Point<spacedim>
 Manifold<dim, spacedim>::
 get_new_point_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
 {
-  return get_new_point (get_default_quadrature(quad));
+  const std::pair<std::vector<Point<spacedim> >, std::vector<double> > points_weights(get_default_points_and_weights(quad));
+  return get_new_point (points_weights.first,points_weights.second);
 }
 
 
@@ -407,7 +447,8 @@ Point<3>
 Manifold<3,3>::
 get_new_point_on_hex (const Triangulation<3, 3>::hex_iterator &hex) const
 {
-  return get_new_point(get_default_quadrature(hex, true));
+  const std::pair<std::vector<Point<3> >, std::vector<double> > points_weights(get_default_points_and_weights(hex,true));
+  return get_new_point (points_weights.first,points_weights.second);
 }
 
 
@@ -427,7 +468,7 @@ Manifold<dim,spacedim>::get_tangent_vector(const Point<spacedim> &x1,
   w.push_back(epsilon);
   w.push_back(1.0-epsilon);
 
-  const Tensor<1,spacedim> neighbor_point = get_new_point (Quadrature<spacedim>(q, w));
+  const Tensor<1,spacedim> neighbor_point = get_new_point (q, w);
   return (neighbor_point-x1)/epsilon;
 }
 
@@ -447,11 +488,19 @@ Point<spacedim>
 FlatManifold<dim, spacedim>::
 get_new_point (const Quadrature<spacedim> &quad) const
 {
-  Assert(std::abs(std::accumulate(quad.get_weights().begin(), quad.get_weights().end(), 0.0)-1.0) < 1e-10,
-         ExcMessage("The weights for the individual points should sum to 1!"));
+  return get_new_point(quad.get_points(),quad.get_weights());
+}
 
-  const std::vector<Point<spacedim> > &surrounding_points = quad.get_points();
-  const std::vector<double> &weights = quad.get_weights();
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+FlatManifold<dim, spacedim>::
+get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+               const std::vector<double>           &weights) const
+{
+  Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0)-1.0) < 1e-10,
+         ExcMessage("The weights for the individual points should sum to 1!"));
 
   Tensor<1,spacedim> minP = periodicity;
 
@@ -554,15 +603,23 @@ Point<spacedim>
 ChartManifold<dim,spacedim,chartdim>::
 get_new_point (const Quadrature<spacedim> &quad) const
 {
-  const std::vector<Point<spacedim> > &surrounding_points = quad.get_points();
-  const std::vector<double> &weights = quad.get_weights();
+  return get_new_point(quad.get_points(),quad.get_weights());
+}
+
+
+
+template <int dim, int spacedim, int chartdim>
+Point<spacedim>
+ChartManifold<dim,spacedim,chartdim>::
+get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+               const std::vector<double>           &weights) const
+{
   std::vector<Point<chartdim> > chart_points(surrounding_points.size());
 
   for (unsigned int i=0; i<surrounding_points.size(); ++i)
     chart_points[i] = pull_back(surrounding_points[i]);
 
-  const Quadrature<chartdim> chart_quad(chart_points, weights);
-  const Point<chartdim> p_chart = sub_manifold.get_new_point(chart_quad);
+  const Point<chartdim> p_chart = sub_manifold.get_new_point(chart_points,weights);
 
   return push_forward(p_chart);
 }
diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index 0b0fd94ae0..e278616f97 100644
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -308,11 +308,20 @@ Point<spacedim>
 CylindricalManifold<dim,spacedim>::
 get_new_point (const Quadrature<spacedim> &quad) const
 {
-  const std::vector<Point<spacedim> > &surrounding_points = quad.get_points();
-  const std::vector<double> &weights = quad.get_weights();
+  return get_new_point(quad.get_points(),quad.get_weights());
+}
+
 
+
+template <int dim, int spacedim>
+Point<spacedim>
+CylindricalManifold<dim,spacedim>::
+get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
+               const std::vector<double>           &weights) const
+{
   // compute a proposed new point
-  Point<spacedim> middle = flat_manifold.get_new_point(quad);
+  const Point<spacedim> middle = flat_manifold.get_new_point(surrounding_points,
+                                                             weights);
 
   double radius = 0;
   Tensor<1,spacedim> on_plane;
diff --git a/source/grid/tria.cc b/source/grid/tria.cc
index a70616c3c6..0dc044a918 100644
--- a/source/grid/tria.cc
+++ b/source/grid/tria.cc
@@ -4098,9 +4098,8 @@ namespace internal
                         ps[1] = cell->face(GeometryInfo<dim>
                                            ::opposite_face[boundary_face])
                                 ->child(0)->vertex(1);
-                        Quadrature<spacedim> qs(ps,ws);
                         triangulation.vertices[next_unused_vertex]
-                          = cell->get_manifold().get_new_point(qs);
+                          = cell->get_manifold().get_new_point(ps,ws);
                       }
                   }
               }
diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc
index 594fede818..4f3b53a30c 100644
--- a/source/grid/tria_accessor.cc
+++ b/source/grid/tria_accessor.cc
@@ -1050,8 +1050,10 @@ namespace
     else
       {
         TriaRawIterator<TriaAccessor<structdim, dim, spacedim> > it(obj);
-        Quadrature<spacedim> quadrature = Manifolds::get_default_quadrature(it, use_laplace);
-        return obj.get_manifold().get_new_point(quadrature);
+        const std::pair<std::vector<Point<spacedim> >,
+              std::vector<double> > points_and_weights = Manifolds::get_default_points_and_weights(it, use_laplace);
+        return obj.get_manifold().get_new_point(points_and_weights.first,
+                                                points_and_weights.second);
       }
   }
 }
@@ -1216,8 +1218,7 @@ TriaAccessor<structdim, dim, spacedim>::intermediate_point (const Point<structdi
       w[i] = fe.shape_value(i, coordinates);
     }
 
-  Quadrature<spacedim> quadrature(p, w);
-  return this->get_manifold().get_new_point(quadrature);
+  return this->get_manifold().get_new_point(p, w);
 }
 
 
diff --git a/tests/manifold/chart_manifold_03.cc b/tests/manifold/chart_manifold_03.cc
index 3e992b01c6..31c35e91b6 100644
--- a/tests/manifold/chart_manifold_03.cc
+++ b/tests/manifold/chart_manifold_03.cc
@@ -119,8 +119,7 @@ void test(unsigned int ref=1)
 
   for (unsigned int i=0; i<ps.size(); ++i)
     {
-      quad = Quadrature<spacedim>(ps[i],ws);
-      middle = manifold.get_new_point(quad);
+      middle = manifold.get_new_point(ps[i],ws);
       deallog << "P0: " << ps[i][0] << " , P1: " << ps[i][1] << " , Middle: " << middle << std::endl;
     }
 
diff --git a/tests/manifold/chart_manifold_03_embedded.cc b/tests/manifold/chart_manifold_03_embedded.cc
index bc59316482..a1c4bd2b9d 100644
--- a/tests/manifold/chart_manifold_03_embedded.cc
+++ b/tests/manifold/chart_manifold_03_embedded.cc
@@ -127,8 +127,7 @@ void test(unsigned int ref=1)
 
   for (unsigned int i=0; i<ps.size(); ++i)
     {
-      quad = Quadrature<spacedim>(ps[i],ws);
-      middle = manifold.get_new_point(quad);
+      middle = manifold.get_new_point(ps[i],ws);
       deallog << "P0: " << ps[i][0] << " , P1: " << ps[i][1] << " , Middle: " << middle << std::endl;
     }
 
diff --git a/tests/manifold/chart_manifold_09.cc b/tests/manifold/chart_manifold_09.cc
index 66fdebd283..c994caf9cf 100644
--- a/tests/manifold/chart_manifold_09.cc
+++ b/tests/manifold/chart_manifold_09.cc
@@ -60,7 +60,7 @@ main()
           weights[0] = (double)i/((double)n_intermediates-1);
           weights[1] = 1.0-weights[0];
 
-          deallog << manifold.get_new_point(Quadrature<2>(points, weights)) << std::endl;
+          deallog << manifold.get_new_point(points, weights) << std::endl;
         }
     }
 
diff --git a/tests/manifold/composition_manifold_01.cc b/tests/manifold/composition_manifold_01.cc
index b6bfc236e6..4351b5cb75 100644
--- a/tests/manifold/composition_manifold_01.cc
+++ b/tests/manifold/composition_manifold_01.cc
@@ -59,7 +59,7 @@ int main ()
       w[0] = 1.0-(double)i/((double)n_intermediates);
       w[1] = 1.0 - w[0];
 
-      Point<spacedim> ip = manifold.get_new_point(Quadrature<spacedim>(sp, w));
+      Point<spacedim> ip = manifold.get_new_point(sp, w);
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
 
diff --git a/tests/manifold/composition_manifold_02.cc b/tests/manifold/composition_manifold_02.cc
index 6b91128bab..98c40d82fa 100644
--- a/tests/manifold/composition_manifold_02.cc
+++ b/tests/manifold/composition_manifold_02.cc
@@ -70,7 +70,7 @@ int main ()
       w[0] = 1.0-(double)i/((double)n_intermediates);
       w[1] = 1.0 - w[0];
 
-      Point<spacedim> ip = manifold.get_new_point(Quadrature<spacedim>(sp, w));
+      Point<spacedim> ip = manifold.get_new_point(sp, w);
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
 
diff --git a/tests/manifold/composition_manifold_03.cc b/tests/manifold/composition_manifold_03.cc
index bfbd11361d..2766d5e863 100644
--- a/tests/manifold/composition_manifold_03.cc
+++ b/tests/manifold/composition_manifold_03.cc
@@ -74,7 +74,7 @@ int main ()
       w[0] = 1.0-(double)i/((double)n_intermediates);
       w[1] = 1.0 - w[0];
 
-      Point<spacedim> ip = manifold.get_new_point(Quadrature<spacedim>(sp, w));
+      Point<spacedim> ip = manifold.get_new_point(sp, w);
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
 
diff --git a/tests/manifold/composition_manifold_04.cc b/tests/manifold/composition_manifold_04.cc
index a442601691..6ee8081701 100644
--- a/tests/manifold/composition_manifold_04.cc
+++ b/tests/manifold/composition_manifold_04.cc
@@ -75,7 +75,7 @@ int main ()
       w[0] = 1.0-(double)i/((double)n_intermediates);
       w[1] = 1.0 - w[0];
 
-      Point<spacedim> ip = manifold.get_new_point(Quadrature<spacedim>(sp, w));
+      Point<spacedim> ip = manifold.get_new_point(sp, w);
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
 
@@ -95,7 +95,7 @@ int main ()
       w[1] = 1.0 - w[0];
 
       Point<spacedim> ip = manifold.
-                           pull_back(manifold.get_new_point(Quadrature<spacedim>(sp, w)));
+                           pull_back(manifold.get_new_point(sp, w));
 
       ip[0] = w[1];
 
diff --git a/tests/manifold/flat_manifold_03.cc b/tests/manifold/flat_manifold_03.cc
index 946142ac56..52a4f0612c 100644
--- a/tests/manifold/flat_manifold_03.cc
+++ b/tests/manifold/flat_manifold_03.cc
@@ -78,8 +78,7 @@ void test(unsigned int ref=1)
 
   for (unsigned int i=0; i<ps.size(); ++i)
     {
-      quad = Quadrature<spacedim>(ps[i],ws);
-      middle = manifold.get_new_point(quad);
+      middle = manifold.get_new_point(ps[i],ws);
       deallog << "P0: " << ps[i][0] << " , P1: " << ps[i][1] << " , Middle: " << middle << std::endl;
     }
 
diff --git a/tests/manifold/flat_manifold_07.cc b/tests/manifold/flat_manifold_07.cc
index 2b7e0ad662..623eae5acf 100644
--- a/tests/manifold/flat_manifold_07.cc
+++ b/tests/manifold/flat_manifold_07.cc
@@ -60,7 +60,7 @@ main()
           weights[0] = (double)i/((double)n_intermediates-1);
           weights[1] = 1.0-weights[0];
 
-          deallog << manifold.get_new_point(Quadrature<2>(points, weights)) << std::endl;
+          deallog << manifold.get_new_point(points, weights) << std::endl;
         }
     }
 
diff --git a/tests/manifold/function_manifold_03.cc b/tests/manifold/function_manifold_03.cc
index a698f7ca1b..136c8d767f 100644
--- a/tests/manifold/function_manifold_03.cc
+++ b/tests/manifold/function_manifold_03.cc
@@ -69,7 +69,7 @@ void test(unsigned int ref=1)
       w[0] = 1.0-(double)i/((double)n_intermediates);
       w[1] = 1.0 - w[0];
 
-      Point<spacedim> ip = manifold.get_new_point(Quadrature<spacedim>(p, w));
+      Point<spacedim> ip = manifold.get_new_point(p, w);
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, p[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, p[1]);
 
diff --git a/tests/manifold/polar_manifold_08.cc b/tests/manifold/polar_manifold_08.cc
index 4638a037df..84c64c718e 100644
--- a/tests/manifold/polar_manifold_08.cc
+++ b/tests/manifold/polar_manifold_08.cc
@@ -66,7 +66,7 @@ main()
           weights[0] = (double)i/((double)n_intermediates-1);
           weights[1] = 1.0-weights[0];
 
-          deallog << manifold.get_new_point(Quadrature<2>(points, weights)) << std::endl;
+          deallog << manifold.get_new_point(points, weights) << std::endl;
         }
     }
 
diff --git a/tests/manifold/spherical_manifold_01.cc b/tests/manifold/spherical_manifold_01.cc
index 1cef855d29..ecba0318cc 100644
--- a/tests/manifold/spherical_manifold_01.cc
+++ b/tests/manifold/spherical_manifold_01.cc
@@ -158,13 +158,9 @@ main()
     weights[1] = 1.0/3.0;
     weights[2] = 1.0/3.0;
 
-    Quadrature<3> quad1(points1, weights);
-    Quadrature<3> quad2(points2, weights);
-    Quadrature<3> quad3(points3, weights);
-
-    Point<3> Q = manifold.get_new_point(quad1);
-    Point<3> S = manifold.get_new_point(quad2);
-    Point<3> T = manifold.get_new_point(quad3);
+    Point<3> Q = manifold.get_new_point(points1, weights);
+    Point<3> S = manifold.get_new_point(points2, weights);
+    Point<3> T = manifold.get_new_point(points3, weights);
 
     Point<3> P5(0.707107, 0.707107, 0.0);
     Point<3> P4(0.0, 0.0, 1.0);
diff --git a/tests/manifold/tensor_product_manifold_01.cc b/tests/manifold/tensor_product_manifold_01.cc
index a1ad09f494..bde990c3c7 100644
--- a/tests/manifold/tensor_product_manifold_01.cc
+++ b/tests/manifold/tensor_product_manifold_01.cc
@@ -63,7 +63,7 @@ void test1()
       w[0] = 1.0-(double)i/((double)n_intermediates);
       w[1] = 1.0 - w[0];
 
-      Point<spacedim> ip = manifold.get_new_point(Quadrature<spacedim>(sp, w));
+      Point<spacedim> ip = manifold.get_new_point(sp, w);
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);