]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Avoid more memory allocation
authorRene Gassmoeller <rene.gassmoeller@mailbox.org>
Tue, 5 Dec 2017 00:52:50 +0000 (17:52 -0700)
committerRene Gassmoeller <rene.gassmoeller@mailbox.org>
Wed, 13 Dec 2017 15:51:06 +0000 (08:51 -0700)
include/deal.II/base/array_view.h
include/deal.II/grid/manifold_lib.h
source/grid/manifold_lib.cc

index b3f9f6c02f469e8621c27773594c6fb660fe40e6..013e86a034e368c5fb16f0557fd451a617e8ed2f 100644 (file)
@@ -614,6 +614,54 @@ make_array_view (Table<2,ElementType>                           &table,
 
 
 
+/**
+ * Create a view to an entire Table<2> object. This is equivalent to
+ * initializing an ArrayView object with a pointer to the first element of the
+ * given table, and the number of table entries as the length of the view.
+ *
+ * This function is used for non-@p const references to objects of Table type.
+ * Such objects contain elements that can be written to. Consequently, the
+ * return type of this function is a view to a set of writable objects.
+ *
+ * @param[in] table The Table for which we want to have an array view object.
+ * The array view corresponds to the <em>entire</em> table.
+ *
+ * @relates ArrayView
+ */
+template <typename ElementType>
+inline
+ArrayView<ElementType>
+make_array_view (Table<2,ElementType> &table)
+{
+  return ArrayView<ElementType> (&table[0][0], table.n_elements());
+}
+
+
+
+/**
+ * Create a view to an entire Table<2> object. This is equivalent to
+ * initializing an ArrayView object with a pointer to the first element of the
+ * given table, and the number of table entries as the length of the view.
+ *
+ * This function is used for @p const references to objects of Table type
+ * because they contain immutable elements. Consequently, the return type of
+ * this function is a view to a set of @p const objects.
+ *
+ * @param[in] table The Table for which we want to have an array view object.
+ * The array view corresponds to the <em>entire</em> table.
+ *
+ * @relates ArrayView
+ */
+template <typename ElementType>
+inline
+ArrayView<const ElementType>
+make_array_view (const Table<2,ElementType> &table)
+{
+  return ArrayView<const ElementType> (&table[0][0], table.n_elements());
+}
+
+
+
 /**
  * Create a view to an entire row of a Table<2> object. This is equivalent to
  * initializing an ArrayView object with a pointer to the first element of the
index 2d25d7abfb7136e3005a9e6141d9cfbe3de405c2..89583483701411ed8de420ddf3ee43fc4d957f7a 100644 (file)
@@ -301,6 +301,25 @@ private:
                  const ArrayView<const double> &weights,
                  const Point<spacedim> &candidate_point) const;
 
+  /**
+   * Compute a new set of points that interpolate between the given points @p
+   * surrounding_points. @p weights is an array view with as many entries as @p
+   * surrounding_points.size() times @p new_points.size().
+   *
+   * This function is optimized to perform on a collection
+   * of new points, by collecting operations that are not dependent on the
+   * weights outside of the loop over all new points.
+   *
+   * The implementation does not allow for @p surrounding_points and
+   * @p new_points to point to the same array, so make sure to pass different
+   * objects into the function.
+   */
+  virtual
+  void
+  get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+                  const ArrayView<const double>          &weights,
+                  ArrayView<Point<spacedim>>              new_points) const;
+
   /**
    * A manifold description to be used for get_new_point in 2D.
    **/
index c061714c0693d074ca41587106472a18c08c5114..3bd2c3c25ad9b8c59d75d73904184e0414e3fabc 100644 (file)
@@ -376,15 +376,53 @@ void
 SphericalManifold<dim, spacedim>::
 get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
                 const Table<2,double>                  &weights,
-                ArrayView<Point<spacedim>>              new_points) const
+                ArrayView<Point<spacedim>>             new_points) const
 {
   AssertDimension(new_points.size(), weights.size(0));
   AssertDimension(surrounding_points.size(), weights.size(1));
 
+  get_new_points(surrounding_points,
+                 make_array_view(weights),
+                 new_points);
+
+  return;
+}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+SphericalManifold<dim,spacedim>::
+get_new_point (const ArrayView<const Point<spacedim>> &vertices,
+               const ArrayView<const double>          &weights) const
+{
+  // To avoid duplicating all of the logic in get_new_points, simply call it
+  // for one position.
+  Point<spacedim> new_point;
+  get_new_points(vertices,
+                 weights,
+                 make_array_view(&new_point,&new_point+1));
+
+  return new_point;
+}
+
+
+
+template <int dim, int spacedim>
+void
+SphericalManifold<dim,spacedim>::
+get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
+                const ArrayView<const double>          &weights,
+                ArrayView<Point<spacedim>>              new_points) const
+{
+  AssertDimension(weights.size(), new_points.size() * surrounding_points.size());
+  const unsigned int weight_rows = new_points.size();
+  const unsigned int weight_columns = surrounding_points.size();
+
   if (surrounding_points.size() == 2)
     {
-      for (unsigned int row=0; row<weights.size(0); ++row)
-        new_points[row] = get_intermediate_point(surrounding_points[0], surrounding_points[1], weights[row][1]);
+      for (unsigned int row=0; row<weight_rows; ++row)
+        new_points[row] = get_intermediate_point(surrounding_points[0], surrounding_points[1], weights[row * weight_columns + 1]);
       return;
     }
 
@@ -418,11 +456,11 @@ get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
   boost::container::small_vector<bool,100> accurate_point_was_found(new_points.size(),false);
   const ArrayView<const Tensor<1,spacedim>> array_directions = make_array_view(directions.begin(),directions.end());
   const ArrayView<const double> array_distances = make_array_view(distances.begin(),distances.end());
-  for (unsigned int row=0; row<weights.size(0); ++row)
+  for (unsigned int row=0; row<weight_rows; ++row)
     {
       new_candidates[row] = guess_new_point(array_directions,
                                             array_distances,
-                                            make_array_view(weights,row));
+                                            ArrayView<const double>(&weights[row*weight_columns],weight_columns));
 
       // If the candidate is the center, mark it as found to avoid entering
       // the Newton iteration in step 2, which would crash.
@@ -437,7 +475,7 @@ get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
       // after we verified that the candidate is not the center.
       if (spacedim<3)
         {
-          new_points[row] = polar_manifold.get_new_point(surrounding_points, make_array_view(weights,row));
+          new_points[row] = polar_manifold.get_new_point(surrounding_points, ArrayView<const double>(&weights[row*weight_columns],weight_columns));
           accurate_point_was_found[row] = true;
           continue;
         }
@@ -450,7 +488,7 @@ get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
   // doesn't already succeed.
   if (max_distance < 2e-2 || spacedim<3)
     {
-      for (unsigned int row=0; row<weights.size(0); ++row)
+      for (unsigned int row=0; row<weight_rows; ++row)
         new_points[row] = center + new_candidates[row].first * new_candidates[row].second;
 
       return;
@@ -461,7 +499,7 @@ get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
 
   // Search for duplicate directions and merge them to minimize the cost of
   // the get_new_point function call below.
-  boost::container::small_vector<double, 1000>             merged_weights(weights.size(0)*weights.size(1));
+  boost::container::small_vector<double, 1000>             merged_weights(weights.size());
   boost::container::small_vector<Tensor<1, spacedim>, 100> merged_directions(surrounding_points.size(),Point<spacedim>());
   boost::container::small_vector<double, 100>              merged_distances(surrounding_points.size(),0.0);
 
@@ -478,8 +516,8 @@ get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
           if (!found_duplicate && squared_distance < 1e-28)
             {
               found_duplicate = true;
-              for (unsigned int row = 0; row < weights.size(0); ++row)
-                merged_weights[row*weights.size(1) + j] += weights[row][i];
+              for (unsigned int row = 0; row < weight_rows; ++row)
+                merged_weights[row*weight_columns + j] += weights[row*weight_columns + i];
             }
         }
 
@@ -487,8 +525,8 @@ get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
         {
           merged_directions[n_unique_directions] = directions[i];
           merged_distances[n_unique_directions] = distances[i];
-          for (unsigned int row = 0; row < weights.size(0); ++row)
-            merged_weights[row*weights.size(1) + n_unique_directions] = weights[row][i];
+          for (unsigned int row = 0; row < weight_rows; ++row)
+            merged_weights[row*weight_columns + n_unique_directions] = weights[row*weight_columns + i];
 
           ++n_unique_directions;
         }
@@ -498,7 +536,7 @@ get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
   // the get_new_point function call below.
   boost::container::small_vector<unsigned int, 100> merged_weights_index(new_points.size(),numbers::invalid_unsigned_int);
   unsigned int unique_weight_rows = 0;
-  for (unsigned int row = 0; row < new_points.size(); ++row)
+  for (unsigned int row = 0; row < weight_rows; ++row)
     {
       bool found_identical_row = false;
       for (unsigned int existing_row = 0; existing_row < row; ++existing_row)
@@ -506,7 +544,7 @@ get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
           bool identical_weights = true;
 
           for (unsigned int weight_index = 0; weight_index < n_unique_directions; ++weight_index)
-            if (merged_weights[row*weights.size(1) + weight_index] != merged_weights[existing_row*weights.size(1) + weight_index])
+            if (std::abs(merged_weights[row*weight_columns + weight_index] - merged_weights[existing_row*weight_columns + weight_index]) > tolerance)
               {
                 identical_weights = false;
                 break;
@@ -527,12 +565,12 @@ get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
   const ArrayView<const double> array_merged_distances = make_array_view(merged_distances.begin(),
                                                          merged_distances.begin()+n_unique_directions);
 
-  for (unsigned int row=0; row<weights.size(0); ++row)
+  for (unsigned int row=0; row<weight_rows; ++row)
     if (!accurate_point_was_found[row])
       {
         if (merged_weights_index[row] == numbers::invalid_unsigned_int)
           {
-            const ArrayView<const double> array_merged_weights (&merged_weights[row*weights.size(1)],n_unique_directions);
+            const ArrayView<const double> array_merged_weights (&merged_weights[row*weight_columns],n_unique_directions);
             new_candidates[row].second = get_new_point(array_merged_directions,
                                                        array_merged_distances,
                                                        array_merged_weights,
@@ -547,24 +585,6 @@ get_new_points (const ArrayView<const Point<spacedim>> &surrounding_points,
 
 
 
-template <int dim, int spacedim>
-Point<spacedim>
-SphericalManifold<dim,spacedim>::
-get_new_point (const ArrayView<const Point<spacedim>> &vertices,
-               const ArrayView<const double>          &weights) const
-{
-  // To avoid duplicating all of the logic in get_new_points, simply call it
-  // for one position.
-  const Table<2,double> table_weights(1,weights.size(),weights.begin());
-  Point<spacedim> new_point;
-  ArrayView<Point<spacedim>> new_points(&new_point,1);
-
-  get_new_points(vertices,table_weights,new_points);
-  return new_points[0];
-}
-
-
-
 template <int dim, int spacedim>
 std::pair<double, Tensor<1,spacedim> >
 SphericalManifold<dim,spacedim>::

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.