]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Break the Manifold interface for performance.
authorDavid Wells <wellsd2@rpi.edu>
Fri, 21 Jul 2017 13:31:43 +0000 (09:31 -0400)
committerDavid Wells <wellsd2@rpi.edu>
Mon, 11 Sep 2017 18:17:09 +0000 (14:17 -0400)
This commit changes the interfaces of

Manifolds::get_default_points_and_weights
Manifold::project_to_manifold
Manifold::get_new_point
Manifold::add_new_points

to use ArrayView instead of std::vector. In addition, the interface of
add_new_points has been changed to populate the ArrayView argument
instead of appending to the end of the array.

This breaks the public interface of Manifold for the sake of improving
performance by about 30%: profiling indicates that, when creating a grid
with a Manifold, we spend about 30% of our time purely calling new and
delete since we must create and destroy so many std::vectors.

27 files changed:
include/deal.II/grid/manifold.h
include/deal.II/grid/manifold_lib.h
include/deal.II/opencascade/boundary_lib.h
source/fe/mapping_manifold.cc
source/fe/mapping_q_generic.cc
source/grid/manifold.cc
source/grid/manifold_lib.cc
source/grid/tria.cc
source/grid/tria_accessor.cc
source/opencascade/boundary_lib.cc
tests/manifold/chart_manifold_03.cc
tests/manifold/chart_manifold_03_embedded.cc
tests/manifold/chart_manifold_09.cc
tests/manifold/composition_manifold_01.cc
tests/manifold/composition_manifold_02.cc
tests/manifold/composition_manifold_03.cc
tests/manifold/composition_manifold_04.cc
tests/manifold/flat_manifold_03.cc
tests/manifold/flat_manifold_07.cc
tests/manifold/function_manifold_03.cc
tests/manifold/polar_manifold_08.cc
tests/manifold/projection_manifold_01.cc
tests/manifold/spherical_manifold_01.cc
tests/manifold/tensor_product_manifold_01.cc
tests/manifold/transfinite_manifold_01.cc
tests/manifold/transfinite_manifold_02.cc
tests/manifold/transfinite_manifold_04.cc

index 5dc6d96f646c7552c6af9a66614147596efa0c7d..aa41005992b338495d21f7f95afbc0fa39d7ffa9 100644 (file)
@@ -20,6 +20,7 @@
 /*----------------------------   manifold.h     ---------------------------*/
 
 #include <deal.II/base/config.h>
+#include <deal.II/base/array_view.h>
 #include <deal.II/base/subscriptor.h>
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/thread_management.h>
@@ -38,6 +39,25 @@ template <int, typename> class Table;
  */
 namespace Manifolds
 {
+  /**
+   * A <code>constexpr</code> helper function that returns the number of
+   * default points for the structure type pointed to by the given
+   * <code>MeshIteratorType</code>. See the documentation of
+   * Manifolds::get_default_points_and_weights() for more information.
+   */
+  template <typename MeshIteratorType>
+  inline
+  constexpr std::size_t n_default_points_per_cell()
+  {
+    // Note that in C++11 a constexpr function can only have a return
+    // statement, so we cannot alias the structure dimension
+    return GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::vertices_per_cell
+           + GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::lines_per_cell
+           + GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::quads_per_cell
+           + GeometryInfo<MeshIteratorType::AccessorType::structure_dimension>::hexes_per_cell
+           - 1; // don't count the cell itself, just the bounding objects
+  }
+
   /**
    * Given a general mesh iterator, construct a quadrature object that
    * contains the following points:
@@ -85,7 +105,7 @@ namespace Manifolds
                          const bool              with_laplace = false) DEAL_II_DEPRECATED;
 
   /**
-   * Given a general mesh iterator, construct vectors of quadrature points and
+   * Given a general mesh iterator, construct arrays 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
@@ -126,13 +146,13 @@ namespace Manifolds
    *   <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> >
+  std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
+      n_default_points_per_cell<MeshIteratorType>()>,
+      std::array<double, n_default_points_per_cell<MeshIteratorType>()> >
       get_default_points_and_weights(const MeshIteratorType &iterator,
                                      const bool              with_laplace = false);
 }
 
-
 /**
  * Manifolds are used to describe the geometry of boundaries of domains as
  * well as the geometry of the interior. Manifold objects are therefore
@@ -373,20 +393,18 @@ public:
    */
   virtual
   Point<spacedim>
-  get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
-                 const std::vector<double>           &weights) const;
+  get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+                 const ArrayView<const double>          &weights) const;
+
 
   /**
-   * Compute a new set of points that interpolate between the given points
-   * @p surrounding_points. @p weights is a table with as many columns as
-   * @p surrounding_points.size(). The number of rows in @p weights determines
-   * how many new points will be computed and appended to the last input
-   * argument @p new_points. After exit of this function, the size of
-   * @p new_points equals the size at entry plus the number of rows in
-   * @p weights.
+   * Compute a new set of points that interpolate between the given points @p
+   * surrounding_points. @p weights is a table with as many columns as @p
+   * surrounding_points.size(). The number of rows in @p weights must match
+   * the length of @p new_points.
    *
    * In its default implementation, this function simply calls get_new_point()
-   * on each row of @p weights and appends those points to the output vector
+   * on each row of @p weights and writes those points into the output array
    * @p new_points. However, this function is more efficient if multiple new
    * points need to be generated like in MappingQGeneric and the manifold does
    * expensive transformations between a chart space and the physical space,
@@ -396,14 +414,14 @@ public:
    * by implementing only the get_new_point() function.
    *
    * The implementation does not allow for @p surrounding_points and
-   * @p new_points to point to the same vector, so make sure to pass different
+   * @p new_points to point to the same array, so make sure to pass different
    * objects into the function.
    */
   virtual
   void
-  add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
-                  const Table<2,double>               &weights,
-                  std::vector<Point<spacedim> >       &new_points) const;
+  add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+                  const Table<2,double>                  &weights,
+                  ArrayView<Point<spacedim>>              new_points) const;
 
   /**
    * Given a point which lies close to the given manifold, it modifies it and
@@ -419,8 +437,8 @@ public:
    * the default behavior should work out of the box.
    */
   virtual
-  Point<spacedim> project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
-                                       const Point<spacedim> &candidate) const;
+  Point<spacedim> project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+                                       const Point<spacedim>                  &candidate) const;
 
   /**
    * Backward compatibility interface.  Return the point which shall become
@@ -704,17 +722,15 @@ public:
    */
   virtual
   Point<spacedim>
-  get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
-                const std::vector<double>           &weights) const;
+  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
+                const ArrayView<const double>          &weights) const override;
+
 
   /**
-   * Compute a new set of points that interpolate between the given points
-   * @p surrounding_points. @p weights is a table with as many columns as
-   * @p surrounding_points.size(). The number of rows in @p weights determines
-   * how many new points will be computed and appended to the last input
-   * argument @p new_points. After exit of this function, the size of
-   * @p new_points equals the size at entry plus the number of rows in
-   * @p weights.
+   * Compute a new set of points that interpolate between the given points @p
+   * surrounding_points. @p weights is a table with as many columns as @p
+   * surrounding_points.size(). The number of rows in @p weights must match
+   * the length of @p new_points.
    *
    * For this particular implementation, the interpolation of the
    * @p surrounding_points according to the @p weights is simply performed in
@@ -722,9 +738,9 @@ public:
    */
   virtual
   void
-  add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
-                  const Table<2,double>               &weights,
-                  std::vector<Point<spacedim> >       &new_points) const;
+  add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+                  const Table<2,double>                  &weights,
+                  ArrayView<Point<spacedim>>              new_points) const;
 
   /**
    * Project to FlatManifold. This is the identity function for flat,
@@ -735,8 +751,8 @@ public:
    */
   virtual
   Point<spacedim>
-  project_to_manifold (const std::vector<Point<spacedim> > &points,
-                       const Point<spacedim> &candidate) const;
+  project_to_manifold (const ArrayView<const Point<spacedim>> &points,
+                       const Point<spacedim>                  &candidate) const;
 
   /**
    * Return a vector that, at $\mathbf x_1$, is tangential to
@@ -922,17 +938,14 @@ public:
    */
   virtual
   Point<spacedim>
-  get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
-                const std::vector<double>           &weights) const;
+  get_new_point(const ArrayView<const Point<spacedim>> &surrounding_points,
+                const ArrayView<const double>          &weights) const override;
 
   /**
-   * Compute a new set of points that interpolate between the given points
-   * @p surrounding_points. @p weights is a table with as many columns as
-   * @p surrounding_points.size(). The number of rows in @p weights determines
-   * how many new points will be computed and appended to the last input
-   * argument @p new_points. After exit of this function, the size of
-   * @p new_points equals the size at entry plus the number of rows in
-   * @p weights.
+   * Compute a new set of points that interpolate between the given points @p
+   * surrounding_points. @p weights is a table with as many columns as @p
+   * surrounding_points.size(). The number of rows in @p weights must match
+   * the length of @p new_points.
    *
    * The implementation of this function first transforms the
    * @p surrounding_points to the chart space by calling pull_back(). Then, new
@@ -951,10 +964,9 @@ public:
    */
   virtual
   void
-  add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
-                  const Table<2,double>               &weights,
-                  std::vector<Point<spacedim> >       &new_points) const;
-
+  add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+                  const Table<2,double>                  &weights,
+                  ArrayView<Point<spacedim>>              new_points) const;
   /**
    * Pull back the given point in spacedim to the Euclidean chartdim
    * dimensional space.
@@ -1076,8 +1088,6 @@ private:
 };
 
 
-
-
 /* -------------- declaration of explicit specializations ------------- */
 
 #ifndef DOXYGEN
@@ -1130,24 +1140,30 @@ namespace Manifolds
   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);
+    const auto points_and_weights = get_default_points_and_weights(iterator, with_laplace);
+    static const int spacedim = MeshIteratorType::AccessorType::space_dimension;
+    return Quadrature<spacedim>
+           (std::vector<Point<spacedim>>(points_and_weights.first.begin(),
+                                         points_and_weights.first.end()),
+            std::vector<double>(points_and_weights.second.begin(),
+                                points_and_weights.second.end()));
   }
 
+
+
   template <typename MeshIteratorType>
-  std::pair<std::vector<Point<MeshIteratorType::AccessorType::space_dimension> >,
-      std::vector<double> >
+  std::pair<std::array<Point<MeshIteratorType::AccessorType::space_dimension>,
+      n_default_points_per_cell<MeshIteratorType>()>,
+      std::array<double, n_default_points_per_cell<MeshIteratorType>()> >
       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;
+    const int spacedim = MeshIteratorType::AccessorType::space_dimension;
+    constexpr std::size_t points_per_cell = n_default_points_per_cell<MeshIteratorType>();
 
-    std::pair<std::vector<Point<spacedim> >,
-        std::vector<double> > points_weights;
+    std::pair<std::array<Point<spacedim>, points_per_cell>, std::array<double, points_per_cell> >
+    points_weights;
 
 
     // note that the exact weights are chosen such as to minimize the
@@ -1157,16 +1173,16 @@ namespace Manifolds
     switch (dim)
       {
       case 1:
-        points_weights.first.resize(2);
-        points_weights.second.resize(2);
+        Assert(points_weights.first.size() == 2, ExcInternalError());
+        Assert(points_weights.second.size() == 2, ExcInternalError());
         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:
-        points_weights.first.resize(8);
-        points_weights.second.resize(8);
+        Assert(points_weights.first.size() == 8, ExcInternalError());
+        Assert(points_weights.second.size() == 8, ExcInternalError());
 
         for (unsigned int i=0; i<4; ++i)
           {
@@ -1192,9 +1208,10 @@ namespace Manifolds
           GeometryInfo<dim>::vertices_per_cell+
           GeometryInfo<dim>::lines_per_cell+
           GeometryInfo<dim>::faces_per_cell;
-        points_weights.first.resize(np);
-        points_weights.second.resize(np);
-        std::vector<Point<3> > *sp3 = reinterpret_cast<std::vector<Point<3> > *>(&points_weights.first);
+        Assert(points_weights.first.size() == np, ExcInternalError());
+        Assert(points_weights.second.size() == np, ExcInternalError());
+        auto *sp3 = reinterpret_cast<std::array<Point<3>, n_default_points_per_cell<decltype(hex)>()> *>
+                    (&points_weights.first);
 
         unsigned int j=0;
 
index 80a7b05bbf5ff206f88cb165c0d68a3a3757e443..62c2200c5464cdbcb513f42a23b0938379a22dbb 100644 (file)
@@ -22,6 +22,8 @@
 #include <deal.II/base/function.h>
 #include <deal.II/base/function_parser.h>
 
+#include <boost/container/small_vector.hpp>
+
 DEAL_II_NAMESPACE_OPEN
 
 /**
@@ -233,8 +235,8 @@ public:
    */
   virtual
   Point<spacedim>
-  get_new_point (const std::vector<Point<spacedim> > &vertices,
-                 const std::vector<double> &weights) const;
+  get_new_point (const ArrayView<const Point<spacedim>> &vertices,
+                 const ArrayView<const double>          &weights) const override;
 
   /**
    * The center of the spherical coordinate system.
@@ -303,12 +305,12 @@ public:
    * 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 override;
+  get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+                 const ArrayView<const double>          &weights) const override;
 
 protected:
   /**
-   * A vector orthogonal to direcetion.
+   * A vector orthogonal to the normal direction.
    */
   const Tensor<1,spacedim> normal_direction;
 
@@ -617,11 +619,10 @@ private:
  * the case for PolarManifold but not for Spherical manifold, so be careful
  * when using the latter. In case the quality of the manifold is not good
  * enough, upon mesh refinement it may happen that the transformation to a
- * chart inside the get_new_point() or add_new_points() methods produces points
- * that are outside the unit cell. Then this class throws an exception of
- * type Manifold@<dim,spacedim@>::ExcTransformationFailed. In that case,
- * the mesh should be refined before attaching this class, as done in the
- * following example:
+ * chart inside the get_new_point() or add_new_points() methods produces
+ * points that are outside the unit cell. Then this class throws an exception
+ * of type Mapping::ExcTransformationFailed. In that case, the mesh should be
+ * refined before attaching this class, as done in the following example:
  *
  * @code
  * SphericalManifold<dim> spherical_manifold;
@@ -696,17 +697,14 @@ public:
    */
   virtual
   Point<spacedim>
-  get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
-                 const std::vector<double>           &weights) const;
+  get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+                 const ArrayView<const double>          &weights) const override;
 
   /**
-   * Compute a new set of points that interpolate between the given points
-   * @p surrounding_points. @p weights is a table with as many columns as
-   * @p surrounding_points.size(). The number of rows in @p weights determines
-   * how many new points will be computed and appended to the last input
-   * argument @p new_points. After exit of this function, the size of
-   * @p new_points equals the size at entry plus the number of rows in
-   * @p weights.
+   * Compute a new set of points that interpolate between the given points @p
+   * surrounding_points. @p weights is a table with as many columns as @p
+   * surrounding_points.size(). The number of columns in @p weights must match
+   * the length of @p new_points.
    *
    * The implementation in this class overrides the method in the base class
    * and computes the new point by a transfinite interpolation. The first step
@@ -723,9 +721,9 @@ public:
    */
   virtual
   void
-  add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
-                  const Table<2,double>               &weights,
-                  std::vector<Point<spacedim> >       &new_points) const;
+  add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+                  const Table<2,double>                  &weights,
+                  ArrayView<Point<spacedim>>              new_points) const;
 
 private:
   /**
@@ -739,16 +737,18 @@ private:
    * 10 entries of a the indices <tt>cell->index()</tt>.
    */
   std::array<unsigned int, 10>
-  get_possible_cells_around_points(const std::vector<Point<spacedim> > &surrounding_points) const;
+  get_possible_cells_around_points(const ArrayView<const Point<spacedim>> &surrounding_points) const;
 
   /**
-   * Finalizes the identification of the correct chart and returns the location
-   * of the surrounding points on the chart. This method internally calls
-   * @p get_possible_cells_around_points().
+   * Finalizes the identification of the correct chart and populates @p
+   * chart_points with the pullbacks of the surrounding points. This method
+   * internally calls @p get_possible_cells_around_points().
+   *
+   * Returns an iterator to the cell on which the chart is defined.
    */
-  std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
-      std::vector<Point<dim> > >
-      compute_chart_points(const std::vector<Point<spacedim> > &surrounding_points) const;
+  typename Triangulation<dim,spacedim>::cell_iterator
+  compute_chart_points(const ArrayView<const Point<spacedim>> &surrounding_points,
+                       ArrayView<Point<dim>>                   chart_points) const;
 
   /**
    * Pull back operation into the unit coordinates on the given coarse cell.
index 7055db2880332d7631e9ef76ef11587ab3ae5f1c..6f576b4c43d6c0e9efc5ef78cc357e68f1b36e8e 100644 (file)
@@ -87,8 +87,8 @@ namespace OpenCASCADE
      * algorithms.
      */
     virtual Point<spacedim>
-    project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
-                         const Point<spacedim> &candidate) const;
+    project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+                         const Point<spacedim>                  &candidate) const;
 
 
   private:
@@ -149,8 +149,8 @@ namespace OpenCASCADE
      * projection algorithms.
      */
     virtual Point<spacedim>
-    project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
-                         const Point<spacedim> &candidate) const;
+    project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+                         const Point<spacedim>                  &candidate) const;
 
   private:
     /**
@@ -236,8 +236,8 @@ namespace OpenCASCADE
      * exception is thrown.
      */
     virtual Point<spacedim>
-    project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
-                         const Point<spacedim> &candidate) const;
+    project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+                         const Point<spacedim>                  &candidate) const;
 
   private:
     /**
index 2b98672b4f7266ded10b8bcee2ca0178a22e9924..f8279275172e89de018df05365cf47247fd1dd78 100644 (file)
@@ -213,14 +213,18 @@ MappingManifold<dim,spacedim>::
 transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                              const Point<dim> &p) const
 {
-  std::vector<Point<spacedim> > vertices;
-  std::vector<double> weights;
+  std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertices;
+  std::array<double, GeometryInfo<dim>::vertices_per_cell> weights;
+
   for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
     {
-      vertices.push_back(cell->vertex(v));
-      weights.push_back(GeometryInfo<dim>::d_linear_shape_function(p,v));
+      vertices[v] = cell->vertex(v);
+      weights[v] = GeometryInfo<dim>::d_linear_shape_function(p, v);
     }
-  return cell->get_manifold().get_new_point(vertices, weights);
+  return cell->get_manifold().get_new_point(make_array_view(vertices.begin(),
+                                                            vertices.end()),
+                                            make_array_view(weights.begin(),
+                                                            weights.end()));
 }
 
 
@@ -375,9 +379,9 @@ namespace internal
           {
             for (unsigned int point=0; point<quadrature_points.size(); ++point)
               {
-                quadrature_points[point] = data.manifold->
-                                           get_new_point(data.vertices,
-                                                         data.cell_manifold_quadrature_weights[point+data_set]);
+                quadrature_points[point] = data.manifold->get_new_point
+                                           (make_array_view(data.vertices),
+                                            make_array_view(data.cell_manifold_quadrature_weights[point+data_set]));
               }
           }
       }
@@ -413,9 +417,9 @@ namespace internal
                 const Point<dim> &p = data.quad.point(point+data_set);
 
                 // And get its image on the manifold:
-                const Point<spacedim> P = data.manifold->
-                                          get_new_point(data.vertices,
-                                                        data.cell_manifold_quadrature_weights[point+data_set]);
+                const Point<spacedim> P = data.manifold->get_new_point
+                                          (make_array_view(data.vertices),
+                                           make_array_view(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
@@ -443,8 +447,8 @@ namespace internal
                       data.vertex_weights[j] = GeometryInfo<dim>::d_linear_shape_function(np, j);
 
                     const Point<spacedim> NP=
-                      data.manifold->get_new_point(data.vertices,
-                                                   data.vertex_weights);
+                      data.manifold->get_new_point(make_array_view(data.vertices),
+                                                   make_array_view(data.vertex_weights));
 
                     const Tensor<1,spacedim> T = data.manifold->get_tangent_vector(P, NP);
 
index 72ab60eaa91466cf3de172d381fbbc3b04c1944f..4a3987ffb4b492ad0f1b8a523c102cd188a143c9 100644 (file)
@@ -3810,11 +3810,21 @@ add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterat
             }
           else
             {
-              tmp_points.resize(2);
-              tmp_points[0] = cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0));
-              tmp_points[1] = cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 1));
-              manifold.add_new_points(tmp_points,
-                                      support_point_weights_perimeter_to_interior[0], a);
+              const std::array<Point<spacedim>, 2> vertices
+              {
+                {
+                  cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 0)),
+                  cell->vertex(GeometryInfo<dim>::line_to_cell_vertices(line_no, 1))
+                }
+              };
+
+              const std::size_t n_rows = support_point_weights_perimeter_to_interior[0].size(0);
+              a.resize(a.size() + n_rows);
+              auto a_view = make_array_view(a.end() - n_rows, a.end());
+              manifold.add_new_points(make_array_view(vertices.begin(),
+                                                      vertices.end()),
+                                      support_point_weights_perimeter_to_interior[0],
+                                      a_view);
             }
         }
     }
@@ -3890,7 +3900,9 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell,
           // extract the points surrounding a quad from the points
           // already computed. First get the 4 vertices and then the points on
           // the four lines
-          tmp_points.resize(4 + 4*(polynomial_degree-1));
+          boost::container::small_vector<Point<3>, 200>
+          tmp_points(GeometryInfo<2>::vertices_per_cell
+                     + GeometryInfo<2>::lines_per_cell*(polynomial_degree-1));
           for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
             tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no,v)];
           if (polynomial_degree > 1)
@@ -3900,9 +3912,14 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell,
                   a[GeometryInfo<3>::vertices_per_cell +
                     (polynomial_degree-1)*
                     GeometryInfo<3>::face_to_cell_lines(face_no,line) + i];
-          face->get_manifold().add_new_points (tmp_points,
+
+          const std::size_t n_rows = support_point_weights_perimeter_to_interior[1].size(0);
+          a.resize(a.size() + n_rows);
+          auto a_view = make_array_view(a.end() - n_rows, a.end());
+          face->get_manifold().add_new_points (make_array_view(tmp_points.begin(),
+                                                               tmp_points.end()),
                                                support_point_weights_perimeter_to_interior[1],
-                                               a);
+                                               a_view);
         }
     }
 }
@@ -3924,9 +3941,10 @@ add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell,
     }
   else
     {
-      std::vector<Point<3> > vertices;
+      std::array<Point<3>, GeometryInfo<2>::vertices_per_cell> vertices;
       for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
-        vertices.push_back(cell->vertex(i));
+        vertices[i] = cell->vertex(i);
+
       Table<2,double> weights(Utilities::fixed_power<2>(polynomial_degree-1),
                               GeometryInfo<2>::vertices_per_cell);
       for (unsigned int q=0, q2=0; q2<polynomial_degree-1; ++q2)
@@ -3938,7 +3956,13 @@ add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell,
               weights(q,i) = GeometryInfo<2>::d_linear_shape_function(point, i);
           }
       // TODO: use all surrounding points once Boundary path is removed
-      cell->get_manifold().add_new_points(vertices, weights, a);
+      const std::size_t n_rows = weights.size(0);
+      a.resize(a.size() + n_rows);
+      auto a_view = make_array_view(a.end() - n_rows, a.end());
+      cell->get_manifold().add_new_points(make_array_view(vertices.begin(),
+                                                          vertices.end()),
+                                          weights,
+                                          a_view);
     }
 }
 
@@ -3989,8 +4013,13 @@ compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_
 
       if (all_manifold_ids_are_equal)
         {
-          std::vector<Point<spacedim> > vertices(a);
-          cell->get_manifold().add_new_points(vertices, support_point_weights_cell, a);
+          const std::size_t n_rows = support_point_weights_cell.size(0);
+          a.resize(a.size() + n_rows);
+          auto a_view = make_array_view(a.end() - n_rows, a.end());
+          cell->get_manifold().add_new_points(make_array_view(a.begin(),
+                                                              a.end() - n_rows),
+                                              support_point_weights_cell,
+                                              a_view);
         }
       else
         switch (dim)
@@ -4008,10 +4037,13 @@ compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_
               add_quad_support_points(cell, a);
             else
               {
-                std::vector<Point<spacedim> > tmp_points(a);
-                cell->get_manifold().add_new_points(tmp_points,
+                const std::size_t n_rows = support_point_weights_perimeter_to_interior[1].size(0);
+                a.resize(a.size() + n_rows);
+                auto a_view = make_array_view(a.end() - n_rows, a.end());
+                cell->get_manifold().add_new_points(make_array_view(a.begin(),
+                                                                    a.end() - n_rows),
                                                     support_point_weights_perimeter_to_interior[1],
-                                                    a);
+                                                    a_view);
               }
             break;
 
@@ -4022,10 +4054,13 @@ compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_
 
             // then compute the interior points
             {
-              std::vector<Point<spacedim> > tmp_points(a);
-              cell->get_manifold().add_new_points(tmp_points,
+              const std::size_t n_rows = support_point_weights_perimeter_to_interior[2].size(0);
+              a.resize(a.size() + n_rows);
+              auto a_view = make_array_view(a.end() - n_rows, a.end());
+              cell->get_manifold().add_new_points(make_array_view(a.begin(),
+                                                                  a.end() - n_rows),
                                                   support_point_weights_perimeter_to_interior[2],
-                                                  a);
+                                                  a_view);
             }
             break;
 
index c2d2238f515f18dd7a5a779d74292c8c3eb7c504..1968c74bbcc6a9de83bb6d2fb55ba599eebda8f3 100644 (file)
@@ -29,27 +29,7 @@ DEAL_II_NAMESPACE_OPEN
 
 using namespace Manifolds;
 
-// This structure is used as comparison function for std::sort when sorting
-// points according to their weight.
-struct CompareWeights
-{
-public:
-  CompareWeights(const std::vector<double> &weights)
-    :
-    compare_weights(weights)
-  {}
-
-  bool operator() (unsigned int a, unsigned int b) const
-  {
-    return compare_weights[a] < compare_weights[b];
-  }
-
-private:
-  const std::vector<double> &compare_weights;
-};
-
 /* -------------------------- Manifold --------------------- */
-
 template <int dim, int spacedim>
 Manifold<dim, spacedim>::~Manifold ()
 {}
@@ -59,7 +39,7 @@ Manifold<dim, spacedim>::~Manifold ()
 template <int dim, int spacedim>
 Point<spacedim>
 Manifold<dim, spacedim>::
-project_to_manifold (const std::vector<Point<spacedim> > &,
+project_to_manifold (const ArrayView<const Point<spacedim>> &,
                      const Point<spacedim> &) const
 {
   Assert (false, ExcPureFunctionCalled());
@@ -75,10 +55,9 @@ get_intermediate_point (const Point<spacedim> &p1,
                         const Point<spacedim> &p2,
                         const double w) const
 {
-  std::vector<Point<spacedim> > vertices;
-  vertices.push_back(p1);
-  vertices.push_back(p2);
-  return project_to_manifold(vertices, w * p2 + (1-w)*p1 );
+  const std::array<Point<spacedim>, 2> vertices {{p1, p2}};
+  return project_to_manifold(make_array_view(vertices.begin(), vertices.end()),
+                             w * p2 + (1-w)*p1);
 }
 
 
@@ -86,8 +65,8 @@ get_intermediate_point (const Point<spacedim> &p1,
 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
+get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+               const ArrayView<const double>          &weights) const
 {
   const double tol = 1e-10;
   const unsigned int n_points = surrounding_points.size();
@@ -106,7 +85,11 @@ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
   // associative (as for the SphericalManifold).
   boost::container::small_vector<unsigned int, 100> permutation(n_points);
   std::iota(permutation.begin(), permutation.end(), 0u);
-  std::sort(permutation.begin(), permutation.end(), CompareWeights(weights));
+  std::sort(permutation.begin(), permutation.end(),
+            [&weights](const std::size_t a, const std::size_t b)
+  {
+    return weights[a] < weights[b];
+  });
 
   // Now loop over points in the order of their associated weight
   Point<spacedim> p = surrounding_points[permutation[0]];
@@ -133,22 +116,24 @@ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
 template <int dim, int spacedim>
 void
 Manifold<dim, spacedim>::
-add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
-                const Table<2,double>               &weights,
-                std::vector<Point<spacedim> >       &new_points) const
+add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+                const Table<2,double>                  &weights,
+                ArrayView<Point<spacedim>>              new_points) const
 {
   AssertDimension(surrounding_points.size(), weights.size(1));
-  Assert(&surrounding_points != &new_points,
-         ExcMessage("surrounding_points and new_points cannot be the same "
-                    "array"));
 
-  const unsigned int n_points = surrounding_points.size();
-  std::vector<double> local_weights(n_points);
+  const std::size_t n_points = surrounding_points.size();
+  // TODO find a better dimension-dependent estimate for the size of this
+  // vector
+  boost::container::small_vector<double, 20> local_weights(n_points);
   for (unsigned int row=0; row<weights.size(0); ++row)
     {
       for (unsigned int i=0; i<n_points; ++i)
         local_weights[i] = weights(row,i);
-      new_points.push_back(get_new_point(surrounding_points, local_weights));
+      new_points[row] = get_new_point(make_array_view(surrounding_points.begin(),
+                                                      surrounding_points.end()),
+                                      make_array_view(local_weights.begin(),
+                                                      local_weights.end()));
     }
 }
 
@@ -353,8 +338,11 @@ Point<spacedim>
 Manifold<dim, spacedim>::
 get_new_point_on_line (const typename Triangulation<dim, spacedim>::line_iterator &line) const
 {
-  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);
+  const auto points_weights = get_default_points_and_weights(line);
+  return get_new_point (make_array_view(points_weights.first.begin(),
+                                        points_weights.first.end()),
+                        make_array_view(points_weights.second.begin(),
+                                        points_weights.second.end()));
 }
 
 
@@ -364,8 +352,11 @@ Point<spacedim>
 Manifold<dim, spacedim>::
 get_new_point_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
 {
-  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);
+  const auto points_weights = get_default_points_and_weights(quad);
+  return get_new_point (make_array_view(points_weights.first.begin(),
+                                        points_weights.first.end()),
+                        make_array_view(points_weights.second.begin(),
+                                        points_weights.second.end()));
 }
 
 
@@ -492,8 +483,11 @@ Point<3>
 Manifold<3,3>::
 get_new_point_on_hex (const Triangulation<3, 3>::hex_iterator &hex) const
 {
-  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);
+  const auto points_weights = get_default_points_and_weights(hex, true);
+  return get_new_point (make_array_view(points_weights.first.begin(),
+                                        points_weights.first.end()),
+                        make_array_view(points_weights.second.begin(),
+                                        points_weights.second.end()));
 }
 
 
@@ -505,15 +499,12 @@ Manifold<dim,spacedim>::get_tangent_vector(const Point<spacedim> &x1,
 {
   const double epsilon = 1e-8;
 
-  std::vector<Point<spacedim> > q;
-  q.push_back(x1);
-  q.push_back(x2);
-
-  std::vector<double> w;
-  w.push_back(epsilon);
-  w.push_back(1.0-epsilon);
-
-  const Tensor<1,spacedim> neighbor_point = get_new_point (q, w);
+  const std::array<Point<spacedim>, 2> points {{x1, x2}};
+  const std::array<double, 2> weights {{epsilon, 1.0 - epsilon}};
+  const Point<spacedim> neighbor_point = get_new_point (make_array_view(points.begin(),
+                                                        points.end()),
+                                                        make_array_view(weights.begin(),
+                                                            weights.end()));
   return (neighbor_point-x1)/epsilon;
 }
 
@@ -533,8 +524,8 @@ FlatManifold<dim,spacedim>::FlatManifold (const Tensor<1,spacedim> &periodicity,
 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
+get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+               const ArrayView<const 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!"));
@@ -588,15 +579,15 @@ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
 template <int dim, int spacedim>
 void
 FlatManifold<dim, spacedim>::
-add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
-                const Table<2,double>               &weights,
-                std::vector<Point<spacedim> >       &new_points) const
+add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+                const Table<2,double>                  &weights,
+                ArrayView<Point<spacedim>>              new_points) const
 {
   AssertDimension(surrounding_points.size(), weights.size(1));
   if (weights.size(0) == 0)
     return;
 
-  const unsigned int n_points = surrounding_points.size();
+  const std::size_t n_points = surrounding_points.size();
 
   Tensor<1,spacedim> minP = periodicity;
   for (unsigned int d=0; d<spacedim; ++d)
@@ -612,7 +603,8 @@ add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
   // check whether periodicity shifts some of the points. Only do this if
   // necessary to avoid memory allocation
   const Point<spacedim> *surrounding_points_start = &surrounding_points[0];
-  std::vector<Point<spacedim> > modified_points;
+
+  boost::container::small_vector<Point<spacedim>, 200> modified_points;
   bool adjust_periodicity = false;
   for (unsigned int d=0; d<spacedim; ++d)
     if (periodicity[d] > 0)
@@ -624,7 +616,9 @@ add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
           }
   if (adjust_periodicity == true)
     {
-      modified_points = surrounding_points;
+      modified_points.resize(surrounding_points.size());
+      std::copy(surrounding_points.begin(), surrounding_points.end(),
+                modified_points.begin());
       for (unsigned int d=0; d<spacedim; ++d)
         if (periodicity[d] > 0)
           for (unsigned int i=0; i<n_points; ++i)
@@ -648,16 +642,20 @@ add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
           if (new_point[d] < 0)
             new_point[d] += periodicity[d];
 
-      new_points.push_back(project_to_manifold(surrounding_points, new_point));
+      // TODO should this use surrounding_points_start or surrounding_points?
+      // The older version used surrounding_points
+      new_points[row] = project_to_manifold(make_array_view(surrounding_points.begin(),
+                                                            surrounding_points.end()),
+                                            new_point);
     }
 }
 
 
-
 template <int dim, int spacedim>
 Point<spacedim>
-FlatManifold<dim, spacedim>::project_to_manifold (const std::vector<Point<spacedim> > &/*vertices*/,
-                                                  const Point<spacedim> &candidate) const
+FlatManifold<dim, spacedim>::project_to_manifold
+(const ArrayView<const Point<spacedim>> &/*vertices*/,
+ const Point<spacedim>                  &candidate) const
 {
   return candidate;
 }
@@ -717,15 +715,19 @@ ChartManifold<dim,spacedim,chartdim>::ChartManifold (const Tensor<1,chartdim> &p
 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
+get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+               const ArrayView<const double>          &weights) const
 {
-  std::vector<Point<chartdim> > chart_points(surrounding_points.size());
+  const std::size_t n_points = surrounding_points.size();
+
+  boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
 
-  for (unsigned int i=0; i<surrounding_points.size(); ++i)
+  for (unsigned int i=0; i<n_points; ++i)
     chart_points[i] = pull_back(surrounding_points[i]);
 
-  const Point<chartdim> p_chart = sub_manifold.get_new_point(chart_points,weights);
+  const Point<chartdim> p_chart = sub_manifold.get_new_point
+                                  (make_array_view(chart_points.begin(), chart_points.end()),
+                                   weights);
 
   return push_forward(p_chart);
 }
@@ -735,25 +737,28 @@ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
 template <int dim, int spacedim, int chartdim>
 void
 ChartManifold<dim,spacedim,chartdim>::
-add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
-                const Table<2,double>               &weights,
-                std::vector<Point<spacedim> >       &new_points) const
+add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+                const Table<2,double>                  &weights,
+                ArrayView<Point<spacedim>>              new_points) const
 {
   Assert(weights.size(0) > 0, ExcEmptyObject());
   AssertDimension(surrounding_points.size(), weights.size(1));
 
-  const unsigned int n_points = surrounding_points.size();
+  const std::size_t n_points = surrounding_points.size();
 
-  std::vector<Point<chartdim> > chart_points(n_points);
-  for (unsigned int i=0; i<n_points; ++i)
+  boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
+  for (std::size_t i=0; i<n_points; ++i)
     chart_points[i] = pull_back(surrounding_points[i]);
 
-  std::vector<Point<chartdim> > new_points_on_chart;
-  new_points_on_chart.reserve(weights.size(0));
-  sub_manifold.add_new_points(chart_points, weights, new_points_on_chart);
+  boost::container::small_vector<Point<chartdim>, 200> new_points_on_chart(weights.size(0));
+  sub_manifold.add_new_points(make_array_view(chart_points.begin(),
+                                              chart_points.end()),
+                              weights,
+                              make_array_view(new_points_on_chart.begin(),
+                                              new_points_on_chart.end()));
 
-  for (unsigned int row=0; row<weights.size(0); ++row)
-    new_points.push_back(push_forward(new_points_on_chart[row]));
+  for (std::size_t row=0; row<weights.size(0); ++row)
+    new_points[row] = push_forward(new_points_on_chart[row]);
 }
 
 
index 6b1ee081458f75986a59c95336b00f734f7d766d..fd8c68c321e2c6c836981c1e5e56124711b1304a 100644 (file)
@@ -276,8 +276,8 @@ get_tangent_vector (const Point<spacedim> &p1,
 template <int dim, int spacedim>
 Point<spacedim>
 SphericalManifold<dim,spacedim>::
-get_new_point (const std::vector<Point<spacedim> > &vertices,
-               const std::vector<double> &weights) const
+get_new_point (const ArrayView<const Point<spacedim>> &vertices,
+               const ArrayView<const double>          &weights) const
 {
   const unsigned int n_points = vertices.size();
 
@@ -368,8 +368,8 @@ CylindricalManifold<dim, spacedim>::CylindricalManifold(const Point<spacedim> &d
 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
+get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+               const ArrayView<const double>          &weights) const
 {
   // First check if the average in space lies on the axis.
   Point<spacedim> middle;
@@ -536,11 +536,12 @@ DerivativeForm<1,chartdim, spacedim>
 FunctionManifold<dim,spacedim,chartdim>::push_forward_gradient(const Point<chartdim> &chart_point) const
 {
   DerivativeForm<1, chartdim, spacedim> DF;
-  std::vector<Tensor<1, chartdim> > gradients(spacedim);
-  push_forward_function->vector_gradient(chart_point, gradients);
   for (unsigned int i=0; i<spacedim; ++i)
-    for (unsigned int j=0; j<chartdim; ++j)
-      DF[i][j] = gradients[i][j];
+    {
+      const auto gradient = push_forward_function->gradient(chart_point, i);
+      for (unsigned int j=0; j<chartdim; ++j)
+        DF[i][j] = gradient[j];
+    }
   return DF;
 }
 
@@ -724,8 +725,12 @@ namespace
 
         // add the contribution from the lines around the cell (first line in
         // formula)
-        std::vector<double> weights(GeometryInfo<2>::vertices_per_face);
-        std::vector<Point<spacedim> > points(GeometryInfo<2>::vertices_per_face);
+        std::array<double, GeometryInfo<2>::vertices_per_face> weights;
+        std::array<Point<spacedim>, GeometryInfo<2>::vertices_per_face> points;
+        // note that the views are immutable, but the arrays are not
+        const auto weights_view = make_array_view(weights.begin(), weights.end());
+        const auto points_view = make_array_view(points.begin(), points.end());
+
         for (unsigned int line=0; line<GeometryInfo<2>::lines_per_cell; ++line)
           {
             const double my_weight = line%2 ? chart_point[line/2] : 1-chart_point[line/2];
@@ -748,7 +753,8 @@ namespace
                 weights[0] = 1. - line_point;
                 weights[1] = line_point;
                 new_point += my_weight *
-                             cell.line(line)->get_manifold().get_new_point(points, weights);
+                             cell.line(line)->get_manifold().get_new_point(points_view,
+                                                                           weights_view);
               }
           }
 
@@ -790,8 +796,12 @@ namespace
           weights_lines[line] = 0;
 
         // start with the contributions of the faces
-        std::vector<double> weights;
-        std::vector<Point<spacedim> > points;
+        std::array<double, GeometryInfo<2>::vertices_per_cell> weights;
+        std::array<Point<spacedim>, GeometryInfo<2>::vertices_per_cell> points;
+        // note that the views are immutable, but the arrays are not
+        const auto weights_view = make_array_view(weights.begin(), weights.end());
+        const auto points_view = make_array_view(points.begin(), points.end());
+
         for (unsigned int face=0; face<GeometryInfo<3>::faces_per_cell; ++face)
           {
             Point<2> quad_point(chart_point[(face/2+1)%3], chart_point[(face/2+2)%3]);
@@ -817,15 +827,14 @@ namespace
               }
             else
               {
-                points.resize(GeometryInfo<2>::vertices_per_cell);
-                weights.resize(GeometryInfo<2>::vertices_per_cell);
                 for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
                   {
                     points[v] = cell.vertex(GeometryInfo<3>::face_to_cell_vertices(face,v));
                     weights[v] = GeometryInfo<2>::d_linear_shape_function(quad_point, v);
                   }
                 new_point += my_weight *
-                             cell.face(face)->get_manifold().get_new_point(points, weights);
+                             cell.face(face)->get_manifold().get_new_point(points_view,
+                                                                           weights_view);
               }
           }
 
@@ -860,14 +869,13 @@ namespace
               }
             else
               {
-                points.resize(2);
-                weights.resize(2);
                 points[0] = cell.vertex(GeometryInfo<3>::line_to_cell_vertices(line,0));
                 points[1] = cell.vertex(GeometryInfo<3>::line_to_cell_vertices(line,1));
                 weights[0] = 1. - line_point;
                 weights[1] = line_point;
                 new_point -= my_weight *
-                             cell.line(line)->get_manifold().get_new_point(points, weights);
+                             cell.line(line)->get_manifold().get_new_point(points_view,
+                                                                           weights_view);
               }
           }
 
@@ -995,7 +1003,7 @@ TransfiniteInterpolationManifold<dim,spacedim>
 template <int dim, int spacedim>
 std::array<unsigned int, 10>
 TransfiniteInterpolationManifold<dim,spacedim>
-::get_possible_cells_around_points(const std::vector<Point<spacedim> > &points) const
+::get_possible_cells_around_points(const ArrayView<const Point<spacedim>> &points) const
 {
   // The methods to identify cells around points in GridTools are all written
   // for the active cells, but we are here looking at some cells at the coarse
@@ -1012,7 +1020,7 @@ TransfiniteInterpolationManifold<dim,spacedim>
   typename Triangulation<dim,spacedim>::cell_iterator
   cell = triangulation->begin(level_coarse),
   endc = triangulation->end(level_coarse);
-  std::vector<std::pair<double, unsigned int> > distances_and_cells;
+  boost::container::small_vector<std::pair<double, unsigned int>, 200> distances_and_cells;
   for ( ; cell != endc; ++cell)
     {
       // only consider cells where the current manifold is attached
@@ -1069,14 +1077,14 @@ TransfiniteInterpolationManifold<dim,spacedim>
 
 
 template <int dim, int spacedim>
-std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
-    std::vector<Point<dim> > >
-    TransfiniteInterpolationManifold<dim, spacedim>
-    ::compute_chart_points (const std::vector<Point<spacedim> > &surrounding_points) const
+typename Triangulation<dim,spacedim>::cell_iterator
+TransfiniteInterpolationManifold<dim, spacedim>
+::compute_chart_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+                        ArrayView<Point<dim>>                   chart_points) const
 {
-  std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
-      std::vector<Point<dim> > > chart_points;
-  chart_points.second.resize(surrounding_points.size());
+  Assert(surrounding_points.size() == chart_points.size(),
+         ExcMessage("The chart points array view must be as large as the "
+                    "surrounding points array view."));
 
   std::array<unsigned int,10> nearby_cells =
     get_possible_cells_around_points(surrounding_points);
@@ -1092,11 +1100,11 @@ std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
       bool inside_unit_cell = true;
       for (unsigned int i=0; i<surrounding_points.size(); ++i)
         {
-          chart_points.second[i] = pull_back(cell, surrounding_points[i]);
+          chart_points[i] = pull_back(cell, surrounding_points[i]);
 
           // Tolerance 1e-6 chosen that the method also works with
           // SphericalManifold
-          if (GeometryInfo<dim>::is_inside_unit_cell(chart_points.second[i],
+          if (GeometryInfo<dim>::is_inside_unit_cell(chart_points[i],
                                                      1e-6) == false)
             {
               inside_unit_cell = false;
@@ -1105,16 +1113,14 @@ std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
         }
       if (inside_unit_cell == true)
         {
-          chart_points.first = cell;
-          return chart_points;
+          return cell;
         }
     }
 
   // a valid inversion should have returned a point above.
   AssertThrow(false,
               (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
-  chart_points.second.clear();
-  return chart_points;
+  return typename Triangulation<dim,spacedim>::cell_iterator();
 }
 
 
@@ -1122,16 +1128,17 @@ std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
 template <int dim, int spacedim>
 Point<spacedim>
 TransfiniteInterpolationManifold<dim, spacedim>
-::get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
-                 const std::vector<double>           &weights) const
+::get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
+                 const ArrayView<const double>          &weights) const
 {
-  const std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
-        std::vector<Point<dim> > > chart_points =
-          compute_chart_points(surrounding_points);
+  boost::container::small_vector<Point<dim>, 100> chart_points(surrounding_points.size());
+  ArrayView<Point<dim>> chart_points_view = make_array_view(chart_points.begin(),
+                                                            chart_points.end());
+  const auto cell = compute_chart_points(surrounding_points, chart_points_view);
 
-  const Point<dim> p_chart = chart_manifold.get_new_point(chart_points.second,weights);
+  const Point<dim> p_chart = chart_manifold.get_new_point (chart_points_view, weights);
 
-  return push_forward(chart_points.first, p_chart);
+  return push_forward(cell, p_chart);
 }
 
 
@@ -1139,23 +1146,26 @@ TransfiniteInterpolationManifold<dim, spacedim>
 template <int dim, int spacedim>
 void
 TransfiniteInterpolationManifold<dim,spacedim>::
-add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
-                const Table<2,double>               &weights,
-                std::vector<Point<spacedim> >       &new_points) const
+add_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+                const Table<2,double>                  &weights,
+                ArrayView<Point<spacedim> >             new_points) const
 {
   Assert(weights.size(0) > 0, ExcEmptyObject());
   AssertDimension(surrounding_points.size(), weights.size(1));
 
-  const std::pair<typename Triangulation<dim,spacedim>::cell_iterator,
-        std::vector<Point<dim> > > chart_points =
-          compute_chart_points(surrounding_points);
+  boost::container::small_vector<Point<dim>, 100> chart_points(surrounding_points.size());
+  ArrayView<Point<dim>> chart_points_view = make_array_view(chart_points.begin(),
+                                                            chart_points.end());
+  const auto cell = compute_chart_points(surrounding_points, chart_points_view);
 
-  std::vector<Point<dim> > new_points_on_chart;
-  new_points_on_chart.reserve(weights.size(0));
-  chart_manifold.add_new_points(chart_points.second, weights, new_points_on_chart);
+  boost::container::small_vector<Point<dim>, 100> new_points_on_chart(weights.size(0));
+  chart_manifold.add_new_points(chart_points_view,
+                                weights,
+                                make_array_view(new_points_on_chart.begin(),
+                                                new_points_on_chart.end()));
 
   for (unsigned int row=0; row<weights.size(0); ++row)
-    new_points.push_back(push_forward(chart_points.first, new_points_on_chart[row]));
+    new_points[row] = push_forward(cell, new_points_on_chart[row]);
 }
 
 
index 6e0c4fabe258f9491a1f5f40cc446f9b66677624..5bdb643229239dfe2030b4a00dd126367a56d0e5 100644 (file)
@@ -4044,15 +4044,20 @@ namespace internal
                       // as returned by the underlying manifold
                       // object.
                       {
-                        std::vector<Point<spacedim> > ps(2);
-                        std::vector<double> ws(2, 0.5);
-                        ps[0] = cell->face(boundary_face)
-                                ->child(0)->vertex(1);
-                        ps[1] = cell->face(GeometryInfo<dim>
-                                           ::opposite_face[boundary_face])
-                                ->child(0)->vertex(1);
+                        const std::array<Point<spacedim>, 2> ps
+                        {
+                          {
+                            cell->face(boundary_face)->child(0)->vertex(1),
+                            cell->face(GeometryInfo<dim>::opposite_face[boundary_face])
+                            ->child(0)->vertex(1)
+                          }
+                        };
+                        const std::array<double, 2> ws {{0.5, 0.5}};
                         triangulation.vertices[next_unused_vertex]
-                          = cell->get_manifold().get_new_point(ps,ws);
+                          = cell->get_manifold().get_new_point(make_array_view(ps.begin(),
+                                                                               ps.end()),
+                                                               make_array_view(ws.begin(),
+                                                                               ws.end()));
                       }
                   }
               }
index 59f3f1fa2355290e7ea9afdc692cd1c1b25fe757..106b9615a9bf825f3f930c6aac083a927d11bd39 100644 (file)
@@ -1051,10 +1051,11 @@ namespace
     else
       {
         TriaRawIterator<TriaAccessor<structdim, dim, spacedim> > it(obj);
-        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);
+        const auto points_and_weights = Manifolds::get_default_points_and_weights(it, use_laplace);
+        return obj.get_manifold().get_new_point(make_array_view(points_and_weights.first.begin(),
+                                                                points_and_weights.first.end()),
+                                                make_array_view(points_and_weights.second.begin(),
+                                                                points_and_weights.second.end()));
       }
   }
 }
@@ -1206,8 +1207,8 @@ Point<spacedim>
 TriaAccessor<structdim, dim, spacedim>::intermediate_point (const Point<structdim> &coordinates) const
 {
   // Surrounding points and weights.
-  std::vector<Point<spacedim> > p(GeometryInfo<structdim>::vertices_per_cell);
-  std::vector<double>   w(GeometryInfo<structdim>::vertices_per_cell);
+  std::array<Point<spacedim>, GeometryInfo<structdim>::vertices_per_cell> p;
+  std::array<double, GeometryInfo<structdim>::vertices_per_cell> w;
 
   for (unsigned int i=0; i<GeometryInfo<structdim>::vertices_per_cell; ++i)
     {
@@ -1215,7 +1216,8 @@ TriaAccessor<structdim, dim, spacedim>::intermediate_point (const Point<structdi
       w[i] = GeometryInfo<structdim>::d_linear_shape_function(coordinates, i);
     }
 
-  return this->get_manifold().get_new_point(p, w);
+  return this->get_manifold().get_new_point(make_array_view(p.begin(), p.end()),
+                                            make_array_view(w.begin(), w.end()));
 }
 
 
index 98284f832632fd4ab6f093917f4e60473fd4eee3..8125c5bfd233cdd6a991e72f2676127e9a495097 100644 (file)
@@ -90,8 +90,8 @@ namespace OpenCASCADE
 
   template <int dim, int spacedim>
   Point<spacedim>  NormalProjectionBoundary<dim,spacedim>::
-  project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
-                       const Point<spacedim> &candidate) const
+  project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+                       const Point<spacedim>                  &candidate) const
   {
     (void)surrounding_points;
 #ifdef DEBUG
@@ -120,8 +120,8 @@ namespace OpenCASCADE
 
   template <int dim, int spacedim>
   Point<spacedim>  DirectionalProjectionBoundary<dim,spacedim>::
-  project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
-                       const Point<spacedim> &candidate) const
+  project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+                       const Point<spacedim>                  &candidate) const
   {
     (void)surrounding_points;
 #ifdef DEBUG
@@ -151,8 +151,8 @@ namespace OpenCASCADE
 
   template <int dim, int spacedim>
   Point<spacedim>  NormalToMeshProjectionBoundary<dim,spacedim>::
-  project_to_manifold (const std::vector<Point<spacedim> > &surrounding_points,
-                       const Point<spacedim> &candidate) const
+  project_to_manifold (const ArrayView<const Point<spacedim>> &surrounding_points,
+                       const Point<spacedim>                  &candidate) const
   {
     TopoDS_Shape out_shape;
     Tensor<1,3> average_normal;
index 5ce140d19dc8e3a5926edd27713d9aa8e00547db..b7ee85c83ff26fd45dbb68bec4f023aa194d0c37 100644 (file)
@@ -117,7 +117,8 @@ void test(unsigned int ref=1)
 
   for (unsigned int i=0; i<ps.size(); ++i)
     {
-      middle = manifold.get_new_point(ps[i],ws);
+      middle = manifold.get_new_point(make_array_view(ps[i]),
+                                      make_array_view(ws));
       deallog << "P0: " << ps[i][0] << " , P1: " << ps[i][1] << " , Middle: " << middle << std::endl;
     }
 
@@ -135,4 +136,3 @@ int main ()
 
   return 0;
 }
-
index 46cee38235746a5579ffa450f3563d6c6e0ccf9d..a244c2c845b5d9a2417ab354acb91baf474e4f0b 100644 (file)
@@ -125,7 +125,8 @@ void test(unsigned int ref=1)
 
   for (unsigned int i=0; i<ps.size(); ++i)
     {
-      middle = manifold.get_new_point(ps[i],ws);
+      middle = manifold.get_new_point(make_array_view(ps[i]),
+                                      make_array_view(ws));
       deallog << "P0: " << ps[i][0] << " , P1: " << ps[i][1] << " , Middle: " << middle << std::endl;
     }
 
@@ -141,4 +142,3 @@ int main ()
 
   return 0;
 }
-
index c994caf9cfc81b0564f4d842644f83fb5429bd04..bb1f74d5e61e450a30c922b002b8cb853a5a5e77 100644 (file)
@@ -60,12 +60,11 @@ main()
           weights[0] = (double)i/((double)n_intermediates-1);
           weights[1] = 1.0-weights[0];
 
-          deallog << manifold.get_new_point(points, weights) << std::endl;
+          deallog << manifold.get_new_point(make_array_view(points),
+                                            make_array_view(weights))
+                  << std::endl;
         }
     }
 
   return 0;
 }
-
-
-
index d3577a2567851e769388a11bb7035389fc3029a0..4c3cfe9066d4a38b31f02bfdd62cf290658d2477 100644 (file)
@@ -57,7 +57,8 @@ int main ()
       w[0] = 1.0-(double)i/((double)n_intermediates);
       w[1] = 1.0 - w[0];
 
-      Point<spacedim> ip = manifold.get_new_point(sp, w);
+      Point<spacedim> ip = manifold.get_new_point(make_array_view(sp),
+                                                  make_array_view(w));
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
 
@@ -70,4 +71,3 @@ int main ()
 
   return 0;
 }
-
index 79dacbf933a269118db49a0502ad6d812f1e0406..2b9660fc32fd497b39d8474c25051830f4f901c0 100644 (file)
@@ -68,7 +68,8 @@ int main ()
       w[0] = 1.0-(double)i/((double)n_intermediates);
       w[1] = 1.0 - w[0];
 
-      Point<spacedim> ip = manifold.get_new_point(sp, w);
+      Point<spacedim> ip = manifold.get_new_point(make_array_view(sp),
+                                                  make_array_view(w));
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
 
@@ -80,4 +81,3 @@ int main ()
 
   return 0;
 }
-
index 986fb87b71fccf6c88b0209781a9833da71bea30..4fc292eaa939d4821d54178cdd78f7d7cdf0195a 100644 (file)
@@ -72,7 +72,8 @@ int main ()
       w[0] = 1.0-(double)i/((double)n_intermediates);
       w[1] = 1.0 - w[0];
 
-      Point<spacedim> ip = manifold.get_new_point(sp, w);
+      Point<spacedim> ip = manifold.get_new_point(make_array_view(sp),
+                                                  make_array_view(w));
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
 
@@ -84,4 +85,3 @@ int main ()
 
   return 0;
 }
-
index c38d7193e94e936caf160a9077c3ef8b1b01bae7..b4b0cd2fc02d1ca9d4c8f7c7937353a515a5f460 100644 (file)
@@ -73,7 +73,8 @@ int main ()
       w[0] = 1.0-(double)i/((double)n_intermediates);
       w[1] = 1.0 - w[0];
 
-      Point<spacedim> ip = manifold.get_new_point(sp, w);
+      Point<spacedim> ip = manifold.get_new_point(make_array_view(sp),
+                                                  make_array_view(w));
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
 
@@ -93,7 +94,8 @@ int main ()
       w[1] = 1.0 - w[0];
 
       Point<spacedim> ip = manifold.
-                           pull_back(manifold.get_new_point(sp, w));
+                           pull_back(manifold.get_new_point(make_array_view(sp),
+                                                            make_array_view(w)));
 
       ip[0] = w[1];
 
index bde2e81ed8c6f93e354f9a1364df378b8299e893..58d7058daea12d7809e43c20f4d894a02c1a6bac 100644 (file)
@@ -76,7 +76,8 @@ void test(unsigned int ref=1)
 
   for (unsigned int i=0; i<ps.size(); ++i)
     {
-      middle = manifold.get_new_point(ps[i],ws);
+      middle = manifold.get_new_point(make_array_view(ps[i]),
+                                      make_array_view(ws));
       deallog << "P0: " << ps[i][0] << " , P1: " << ps[i][1] << " , Middle: " << middle << std::endl;
     }
 
@@ -94,4 +95,3 @@ int main ()
 
   return 0;
 }
-
index 623eae5acf4c7fd5a409d3635f6aaf4a3a1168b8..12ddf1ddae66b0b79f1aade3c2289ca7c2360f4e 100644 (file)
@@ -60,12 +60,11 @@ main()
           weights[0] = (double)i/((double)n_intermediates-1);
           weights[1] = 1.0-weights[0];
 
-          deallog << manifold.get_new_point(points, weights) << std::endl;
+          deallog << manifold.get_new_point(make_array_view(points),
+                                            make_array_view(weights))
+                  << std::endl;
         }
     }
 
   return 0;
 }
-
-
-
index a01d825ba156f12c14868be79ac59a3dd69f04a3..3b164fbb8cebdc7f2eca425bfdb82e8f7ec89cfd 100644 (file)
@@ -67,7 +67,8 @@ 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(p, w);
+      Point<spacedim> ip = manifold.get_new_point(make_array_view(p),
+                                                  make_array_view(w));
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, p[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, p[1]);
 
@@ -89,4 +90,3 @@ int main ()
 
   return 0;
 }
-
index 84c64c718e34b74348d0479d9d13e894b09cc183..0ab7a76dc84d3ba663fa85be2f8ff4d3841be88a 100644 (file)
@@ -66,12 +66,11 @@ main()
           weights[0] = (double)i/((double)n_intermediates-1);
           weights[1] = 1.0-weights[0];
 
-          deallog << manifold.get_new_point(points, weights) << std::endl;
+          deallog << manifold.get_new_point(make_array_view(points),
+                                            make_array_view(weights))
+                  << std::endl;
         }
     }
 
   return 0;
 }
-
-
-
index 64d65dcc2a0e7cd8b421a987f18ce2e81bb0b595..b7d53db2f45e23dd23d55170646e18dcdeea53cb 100644 (file)
@@ -31,8 +31,8 @@ class MyManifold : public Manifold<dim, spacedim>
 {
 public:
   Point<spacedim>
-  project_to_manifold (const std::vector<Point<spacedim> > &vertices,
-                       const Point<spacedim> &candidate) const
+  project_to_manifold (const ArrayView<const Point<spacedim> > &vertices,
+                       const Point<spacedim> &candidate) const override
   {
     // Shift the y coordinate to 4*x*(1-x)
     Point<spacedim> p = candidate;
@@ -73,4 +73,3 @@ int main ()
 
   return 0;
 }
-
index ecba0318cc2d89155d83e1a95f5dd51abb0388e9..93bbbbaaebef913638649065cab77e988b0adeb0 100644 (file)
@@ -158,9 +158,12 @@ main()
     weights[1] = 1.0/3.0;
     weights[2] = 1.0/3.0;
 
-    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> Q = manifold.get_new_point(make_array_view(points1),
+                                        make_array_view(weights));
+    Point<3> S = manifold.get_new_point(make_array_view(points2),
+                                        make_array_view(weights));
+    Point<3> T = manifold.get_new_point(make_array_view(points3),
+                                        make_array_view(weights));
 
     Point<3> P5(0.707107, 0.707107, 0.0);
     Point<3> P4(0.0, 0.0, 1.0);
@@ -177,6 +180,3 @@ main()
   // Quadrature (const std::vector< Point< dim > > &points, const std::vector< double > &weights);
   return 0;
 }
-
-
-
index 5301ac42507d0f0f4094828eabb3f32bc4fdd2f6..598f003d241009c36bffe0df8d4f91e42ff17e56 100644 (file)
@@ -61,7 +61,8 @@ void test1()
       w[0] = 1.0-(double)i/((double)n_intermediates);
       w[1] = 1.0 - w[0];
 
-      Point<spacedim> ip = manifold.get_new_point(sp, w);
+      Point<spacedim> ip = manifold.get_new_point(make_array_view(sp),
+                                                  make_array_view(w));
       Tensor<1,spacedim> t1 = manifold.get_tangent_vector(ip, sp[0]);
       Tensor<1,spacedim> t2 = manifold.get_tangent_vector(ip, sp[1]);
 
@@ -83,4 +84,3 @@ int main ()
 
   return 0;
 }
-
index 7c83b2819fbb17ddfecec63831d4ac30b82f9612..1abfbecb701a7a52614eddd06f44dd2c39487f91 100644 (file)
@@ -43,14 +43,18 @@ void do_test(const Triangulation<dim,spacedim> &tria)
             std::vector<double> weights(2);
             weights[0] = 0.1;
             weights[1] = 0.9;
-            Point<spacedim> p = cell->get_manifold().get_new_point(points, weights);
-            Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+            Point<spacedim> p = cell->get_manifold().get_new_point(make_array_view(points),
+                                                                   make_array_view(weights));
+            Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+                                   make_array_view(weights));
             deallog << "Distance between cell manifold and face manifold: "
                     << (pref-p) << std::endl;
             weights[0] = 0.55;
             weights[1] = 0.45;
-            p = cell->get_manifold().get_new_point(points, weights);
-            pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+            p = cell->get_manifold().get_new_point(make_array_view(points),
+                                                   make_array_view(weights));
+            pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+                                                                  make_array_view(weights));
             deallog << "Distance between cell manifold and face manifold: "
                     << (pref-p) << std::endl;
           }
@@ -147,4 +151,3 @@ int main ()
 
   return 0;
 }
-
index b3d47d4c076394da9293d0c622bcdd4324ccbc39..32f2c0dc23143c2e079cb5bf3e35aaaee7dc05b4 100644 (file)
@@ -43,14 +43,18 @@ void do_test(const Triangulation<dim,spacedim> &tria)
             std::vector<double> weights(2);
             weights[0] = 0.1;
             weights[1] = 0.9;
-            Point<spacedim> p = cell->get_manifold().get_new_point(points, weights);
-            Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+            Point<spacedim> p = cell->get_manifold().get_new_point(make_array_view(points),
+                                                                   make_array_view(weights));
+            Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+                                   make_array_view(weights));
             deallog << "Distance between cell manifold and face manifold: "
                     << (pref-p) << std::endl;
             weights[0] = 0.55;
             weights[1] = 0.45;
-            p = cell->get_manifold().get_new_point(points, weights);
-            pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+            p = cell->get_manifold().get_new_point(make_array_view(points),
+                                                   make_array_view(weights));
+            pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+                                                                  make_array_view(weights));
             deallog << "Distance between cell manifold and face manifold: "
                     << (pref-p) << std::endl;
           }
@@ -156,4 +160,3 @@ int main ()
 
   return 0;
 }
-
index c427246c889377652c99bd15fa1268c8b6060dcf..37e635d82fdea2e57edd745cf7e47bd1b0e47963 100644 (file)
@@ -43,14 +43,18 @@ void do_test(const Triangulation<dim,spacedim> &tria)
             std::vector<double> weights(2);
             weights[0] = 0.1;
             weights[1] = 0.9;
-            Point<spacedim> p = cell->get_manifold().get_new_point(points, weights);
-            Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+            Point<spacedim> p = cell->get_manifold().get_new_point(make_array_view(points),
+                                                                   make_array_view(weights));
+            Point<spacedim> pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+                                   make_array_view(weights));
             deallog << "Distance between cell manifold and face manifold: "
                     << (pref-p) << std::endl;
             weights[0] = 0.55;
             weights[1] = 0.45;
-            p = cell->get_manifold().get_new_point(points, weights);
-            pref = cell->face(face)->get_manifold().get_new_point(points, weights);
+            p = cell->get_manifold().get_new_point(make_array_view(points),
+                                                   make_array_view(weights));
+            pref = cell->face(face)->get_manifold().get_new_point(make_array_view(points),
+                                                                  make_array_view(weights));
             deallog << "Distance between cell manifold and face manifold: "
                     << (pref-p) << std::endl;
           }
@@ -118,5 +122,3 @@ int main ()
 
   return 0;
 }
-
-

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.