]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Performance improvements. Augment documentation.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sun, 19 Mar 2017 09:16:53 +0000 (10:16 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 22 Mar 2017 19:34:41 +0000 (20:34 +0100)
include/deal.II/grid/manifold.h
source/grid/manifold.cc

index 424b331bf44c26f4b5ae9edaee46225514df5040..fd355a490fec91930e6cfcc17a874bf7b5ed2c61 100644 (file)
@@ -400,7 +400,7 @@ public:
                  const std::vector<double>           &weights) const;
 
   /**
-   * Compute a new set of points around the given points
+   * 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
@@ -408,16 +408,18 @@ public:
    * @p new_points equals the size at entry plus the number of rows in
    * @p weights.
    *
-   * In its default implementation, this function simply calls
-   * get_new_point() on each row of @weights and appends those points to the
-   * output vector @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, such as ChartManifold. If efficiency is not important,
-   * you may get away by implementing only the get_new_point() function.
+   * In its default implementation, this function simply calls get_new_point()
+   * on each row of @weights and appends those points to the output vector
+   * @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,
+   * such as ChartManifold. For this function, the surrounding points need to
+   * be transformed back to the chart sparse only once, rather than for every
+   * call to get_new_point(). If efficiency is not important, you may get away
+   * 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 pass different
+   * @p new_points to point to the same vector, so make sure to pass different
    * objects into the function.
    */
   virtual
@@ -754,7 +756,7 @@ public:
                 const std::vector<double>           &weights) const;
 
   /**
-   * Compute a new set of points around the given points
+   * 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
@@ -762,8 +764,9 @@ public:
    * @p new_points equals the size at entry plus the number of rows in
    * @p weights.
    *
-   * For this particular implementation, an interpolation of the
-   * @p surrounding_points according to the @p weights is performed.
+   * For this particular implementation, the interpolation of the
+   * @p surrounding_points according to the @p weights is simply performed in
+   * Cartesian space.
    */
   virtual
   void
@@ -982,7 +985,7 @@ public:
                 const std::vector<double>           &weights) const;
 
   /**
-   * Compute a new set of points around the given points
+   * 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
@@ -991,16 +994,19 @@ public:
    * @p weights.
    *
    * The implementation of this function first transforms the
-   * @p surrounding_points to the chart by calling pull_back(). Then, new
+   * @p surrounding_points to the chart space by calling pull_back(). Then, new
    * points are computed on the chart by usual interpolation according to the
    * given @p weights, which are finally transformed to the image space by
    * push_forward().
    *
    * This implementation can be much more efficient for computing multiple new
    * points from the same surrounding points than separate calls to
-   * get_new_point() in case the pull_back() operation is expensive. Often,
-   * this is indeed the case because pull_back() might involve Newton
-   * iterations or something similar in non-trivial manifolds.
+   * get_new_point() in case the pull_back() operation is expensive. This is
+   * because pull_back() is only called once for the surrounding points and
+   * the interpolation is done for all given weights using this set of
+   * points. Often, pull_back() is also more expensive than push_forward()
+   * because the former might involve some kind of Newton iteration in
+   * non-trivial manifolds.
    */
   virtual
   void
index 93a949cd00521b11272bbd153fcf32d8f33ed04e..eb88fe9e35408bec949261dbbe173c5bbe81a441 100644 (file)
@@ -33,16 +33,16 @@ struct CompareWeights
 public:
   CompareWeights(const std::vector<double> &weights)
     :
-    compare_weights(&weights)
+    compare_weights(weights)
   {}
 
   bool operator() (unsigned int a, unsigned int b) const
   {
-    return (*compare_weights)[a] < (*compare_weights)[b];
+    return compare_weights[a] < compare_weights[b];
   }
 
 private:
-  const std::vector<double> *compare_weights;
+  const std::vector<double> &compare_weights;
 };
 
 /* -------------------------- Manifold --------------------- */
@@ -111,12 +111,22 @@ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
   // First sort 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).
-  std::vector<unsigned int> permutation(n_points);
+  unsigned int permutation_short[30];
+  std::vector<unsigned int> permutation_long;
+  unsigned int *permutation;
+  if (n_points > 30)
+    {
+      permutation_long.resize(n_points);
+      permutation = &permutation_long[0];
+    }
+  else
+    permutation = &permutation_short[0];
+
   for (unsigned int i=0; i<n_points; ++i)
     permutation[i] = i;
 
-  std::sort(permutation.begin(),
-            permutation.end(),
+  std::sort(permutation,
+            permutation + n_points,
             CompareWeights(weights));
 
   // Now loop over points in the order of their associated weight
@@ -131,7 +141,8 @@ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
       else
         weight =  w/(weights[permutation[i]] + w);
 
-      p = get_intermediate_point(p, surrounding_points[permutation[i]],1.0 - weight );
+      if (std::abs(weight) > 1e-14)
+        p = get_intermediate_point(p, surrounding_points[permutation[i]],1.0 - weight );
       w += weights[permutation[i]];
     }
 
@@ -152,11 +163,12 @@ add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
          ExcMessage("surrounding_points and new_points cannot be the same "
                     "array"));
 
-  std::vector<double> local_weights(surrounding_points.size());
+  const unsigned int n_points = surrounding_points.size();
+  std::vector<double> local_weights(n_points);
   for (unsigned int row=0; row<weights.size(0); ++row)
     {
-      for (unsigned int p=0; p<surrounding_points.size(); ++p)
-        local_weights[p] = weights[row][p];
+      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));
     }
 }
@@ -774,19 +786,15 @@ add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
   const unsigned int n_points = surrounding_points.size();
 
   std::vector<Point<chartdim> > chart_points(n_points);
-  std::vector<double> local_weights(n_points);
-
   for (unsigned int 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);
+
   for (unsigned int row=0; row<weights.size(0); ++row)
-    {
-      for (unsigned int p=0; p<n_points; ++p)
-        local_weights[p] = weights[row][p];
-      const Point<chartdim> p_chart = sub_manifold.get_new_point(chart_points,
-                                                                 local_weights);
-      new_points.push_back(push_forward(p_chart));
-    }
+    new_points.push_back(push_forward(new_points_on_chart[row]));
 }
 
 

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.