Add new function to query multiple points to Manifold

diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h
index 8e97d71509655dcb353e8fccf8e963c0fb19bfa9..424b331bf44c26f4b5ae9edaee46225514df5040 100644 (file)
--- a/include/deal.II/grid/manifold.h
+++ b/include/deal.II/grid/manifold.h
@@ -29,6 +29,9 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+// forward declaration
+template <int, typename> class Table;
+
 /**
  * We collect here some helper functions used in the Manifold<dim,spacedim>
  * classes.
@@ -396,6 +399,33 @@ public:
   get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
                  const std::vector<double>           &weights) const;
 
+  /**
+   * Compute a new set of points around 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.
+   *
+   * 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.
+   *
+   * The implementation does not allow for @p surrounding_points and
+   * @p new_points to point to the same vector, so make sure 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;
+
   /**
    * Given a point which lies close to the given manifold, it modifies it and
    * projects it to manifold itself.
@@ -723,6 +753,23 @@ public:
   get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
                 const std::vector<double>           &weights) const;
 
+  /**
+   * Compute a new set of points around 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.
+   *
+   * For this particular implementation, an interpolation of the
+   * @p surrounding_points according to the @p weights is performed.
+   */
+  virtual
+  void
+  add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
+                  const Table<2,double>               &weights,
+                  std::vector<Point<spacedim> >       &new_points) const;
 
   /**
    * Project to FlatManifold. This is the identity function for flat,
@@ -934,6 +981,33 @@ public:
   get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
                 const std::vector<double>           &weights) const;
 
+  /**
+   * Compute a new set of points around 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.
+   *
+   * The implementation of this function first transforms the
+   * @p surrounding_points to the chart 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.
+   */
+  virtual
+  void
+  add_new_points (const std::vector<Point<spacedim> > &surrounding_points,
+                  const Table<2,double>               &weights,
+                  std::vector<Point<spacedim> >       &new_points) const;
+
   /**
    * Pull back the given point in spacedim to the Euclidean chartdim
    * dimensional space.

diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc
index f36b8ca3eaf2d25db7327a793e13383cbc824723..93a949cd00521b11272bbd153fcf32d8f33ed04e 100644 (file)
--- a/source/grid/manifold.cc
+++ b/source/grid/manifold.cc
@@ -14,6 +14,7 @@
 // ---------------------------------------------------------------------
 
 #include <deal.II/base/tensor.h>
+#include <deal.II/base/table.h>
 #include <deal.II/grid/manifold.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_iterator.h>
@@ -50,6 +51,8 @@ template <int dim, int spacedim>
 Manifold<dim, spacedim>::~Manifold ()
 {}
 
+
+
 template <int dim, int spacedim>
 Point<spacedim>
 Manifold<dim, spacedim>::
@@ -60,6 +63,8 @@ project_to_manifold (const std::vector<Point<spacedim> > &,
   return Point<spacedim>();
 }
 
+
+
 template <int dim, int spacedim>
 Point<spacedim>
 Manifold<dim, spacedim>::
@@ -73,6 +78,8 @@ get_intermediate_point (const Point<spacedim> &p1,
   return project_to_manifold(vertices, w * p2 + (1-w)*p1 );
 }
 
+
+
 template <int dim, int spacedim>
 Point<spacedim>
 Manifold<dim, spacedim>::
@@ -81,6 +88,8 @@ get_new_point (const Quadrature<spacedim> &quad) const
   return get_new_point(quad.get_points(),quad.get_weights());
 }
 
+
+
 template <int dim, int spacedim>
 Point<spacedim>
 Manifold<dim, spacedim>::
@@ -131,6 +140,29 @@ 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
+{
+  AssertDimension(surrounding_points.size(), weights.size(1));
+  Assert(&surrounding_points != &new_points,
+         ExcMessage("surrounding_points and new_points cannot be the same "
+                    "array"));
+
+  std::vector<double> local_weights(surrounding_points.size());
+  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];
+      new_points.push_back(get_new_point(surrounding_points, local_weights));
+    }
+}
+
+
+
 template <>
 Tensor<1,2>
 Manifold<2, 2>::
@@ -244,6 +276,7 @@ normal_vector (const Triangulation<3, 3>::face_iterator &face,
 }
 
 
+
 template <int dim, int spacedim>
 Tensor<1,spacedim>
 Manifold<dim, spacedim>::
@@ -275,6 +308,7 @@ get_normals_at_vertices (const Triangulation<2, 2>::face_iterator &face,
 }
 
 
+
 template <>
 void
 Manifold<3, 3>::
@@ -307,6 +341,7 @@ get_normals_at_vertices (const Triangulation<3, 3>::face_iterator &face,
 }
 
 
+
 template <int dim, int spacedim>
 void
 Manifold<dim, spacedim>::
@@ -322,7 +357,6 @@ get_normals_at_vertices (const typename Triangulation<dim, spacedim>::face_itera
 
 
 
-
 template <int dim, int spacedim>
 Point<spacedim>
 Manifold<dim, spacedim>::
@@ -344,6 +378,7 @@ get_new_point_on_quad (const typename Triangulation<dim, spacedim>::quad_iterato
 }
 
 
+
 template <int dim, int spacedim>
 Point<spacedim>
 Manifold<dim,spacedim>::
@@ -383,6 +418,7 @@ get_new_point_on_cell (const typename Triangulation<dim,spacedim>::cell_iterator
 }
 
 
+
 template <>
 Point<1>
 Manifold<1,1>::
@@ -393,6 +429,7 @@ get_new_point_on_face (const Triangulation<1,1>::face_iterator &) const
 }
 
 
+
 template <>
 Point<2>
 Manifold<1,2>::
@@ -414,6 +451,7 @@ get_new_point_on_face (const Triangulation<1,3>::face_iterator &) const
 }
 
 
+
 template <>
 Point<1>
 Manifold<1,1>::
@@ -499,6 +537,8 @@ FlatManifold<dim,spacedim>::FlatManifold (const Tensor<1,spacedim> &periodicity,
   tolerance(tolerance)
 {}
 
+
+
 template <int dim, int spacedim>
 Point<spacedim>
 FlatManifold<dim, spacedim>::
@@ -518,11 +558,69 @@ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
   Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0)-1.0) < 1e-10,
          ExcMessage("The weights for the individual points should sum to 1!"));
 
-  Tensor<1,spacedim> minP = periodicity;
+  Point<spacedim> p;
 
+  // if there is no periodicity, use a shortcut
+  if (periodicity == Tensor<1,spacedim>())
+    {
+      for (unsigned int i=0; i<surrounding_points.size(); ++i)
+        p += surrounding_points[i] * weights[i];
+    }
+  else
+    {
+      Tensor<1,spacedim> minP = periodicity;
+
+      for (unsigned int d=0; d<spacedim; ++d)
+        if (periodicity[d] > 0)
+          for (unsigned int i=0; i<surrounding_points.size(); ++i)
+            {
+              minP[d] = std::min(minP[d], surrounding_points[i][d]);
+              Assert( (surrounding_points[i][d] < periodicity[d]+tolerance*periodicity[d]) ||
+                      (surrounding_points[i][d] >= -tolerance*periodicity[d]),
+                      ExcPeriodicBox(d, surrounding_points[i], periodicity[i]));
+            }
+
+      // compute the weighted average point, possibly taking into account periodicity
+      for (unsigned int i=0; i<surrounding_points.size(); ++i)
+        {
+          Point<spacedim> dp;
+          for (unsigned int d=0; d<spacedim; ++d)
+            if (periodicity[d] > 0)
+              dp[d] = ( (surrounding_points[i][d]-minP[d]) > periodicity[d]/2.0 ?
+                        -periodicity[d] : 0.0 );
+
+          p += (surrounding_points[i]+dp)*weights[i];
+        }
+
+      // if necessary, also adjust the weighted point by the periodicity
+      for (unsigned int d=0; d<spacedim; ++d)
+        if (periodicity[d] > 0)
+          if (p[d] < 0)
+            p[d] += periodicity[d];
+    }
+
+  return project_to_manifold(surrounding_points, p);
+}
+
+
+
+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
+{
+  AssertDimension(surrounding_points.size(), weights.size(1));
+  if (weights.size(0) == 0)
+    return;
+
+  const unsigned int n_points = surrounding_points.size();
+
+  Tensor<1,spacedim> minP = periodicity;
   for (unsigned int d=0; d<spacedim; ++d)
     if (periodicity[d] > 0)
-      for (unsigned int i=0; i<surrounding_points.size(); ++i)
+      for (unsigned int i=0; i<n_points; ++i)
         {
           minP[d] = std::min(minP[d], surrounding_points[i][d]);
           Assert( (surrounding_points[i][d] < periodicity[d]+tolerance*periodicity[d]) ||
@@ -530,26 +628,47 @@ get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
                   ExcPeriodicBox(d, surrounding_points[i], periodicity[i]));
         }
 
-  // compute the weighted average point, possibly taking into account periodicity
-  Point<spacedim> p;
-  for (unsigned int i=0; i<surrounding_points.size(); ++i)
+  // 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;
+  bool adjust_periodicity = false;
+  for (unsigned int d=0; d<spacedim; ++d)
+    if (periodicity[d] > 0)
+      for (unsigned int i=0; i<n_points; ++i)
+        if ((surrounding_points[i][d]-minP[d]) > periodicity[d]/2.0)
+          {
+            adjust_periodicity = true;
+            break;
+          }
+  if (adjust_periodicity == true)
     {
-      Point<spacedim> dp;
+      modified_points = surrounding_points;
       for (unsigned int d=0; d<spacedim; ++d)
         if (periodicity[d] > 0)
-          dp[d] = ( (surrounding_points[i][d]-minP[d]) > periodicity[d]/2.0 ?
-                    -periodicity[d] : 0.0 );
-
-      p += (surrounding_points[i]+dp)*weights[i];
+          for (unsigned int i=0; i<n_points; ++i)
+            if ((surrounding_points[i][d]-minP[d]) > periodicity[d]/2.0)
+              modified_points[i][d] -= periodicity[d];
+      surrounding_points_start = &modified_points[0];
     }
 
-  // if necessary, also adjust the weighted point by the periodicity
-  for (unsigned int d=0; d<spacedim; ++d)
-    if (periodicity[d] > 0)
-      if (p[d] < 0)
-        p[d] += periodicity[d];
+  // Now perform the interpolation
+  for (unsigned int row=0; row<weights.size(0); ++row)
+    {
+      Assert(std::abs(std::accumulate(&weights(row,0), &weights(row,0)+n_points, 0.0)-1.0) < 1e-10,
+             ExcMessage("The weights for the individual points should sum to 1!"));
+      Point<spacedim> new_point;
+      for (unsigned int p=0; p<n_points; ++p)
+        new_point += surrounding_points_start[p] * weights(row,p);
 
-  return project_to_manifold(surrounding_points, p);
+      // if necessary, also adjust the weighted point by the periodicity
+      for (unsigned int d=0; d<spacedim; ++d)
+        if (periodicity[d] > 0)
+          if (new_point[d] < 0)
+            new_point[d] += periodicity[d];
+
+      new_points.push_back(project_to_manifold(surrounding_points, new_point));
+    }
 }
 
 
@@ -642,6 +761,36 @@ 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
+{
+  Assert(weights.size(0) > 0, ExcEmptyObject());
+  AssertDimension(surrounding_points.size(), weights.size(1));
+
+  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]);
+
+  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));
+    }
+}
+
+
+
 template <int dim, int spacedim, int chartdim>
 DerivativeForm<1,chartdim,spacedim>
 ChartManifold<dim,spacedim,chartdim>::
@@ -695,4 +844,3 @@ ChartManifold<dim, spacedim, chartdim>::get_periodicity() const
 #include "manifold.inst"
 
 DEAL_II_NAMESPACE_CLOSE
-