Multiple points in Mapping::transform_real_to_unit_cell

diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h
index 8b302d8ba74d26b64a78e5fb78c4bd39931d576c..ac4f6652b0a3d96f10ee66fbf2a3c9126fb5e8eb 100644 (file)
--- a/include/deal.II/fe/mapping.h
+++ b/include/deal.II/fe/mapping.h
@@ -452,6 +452,23 @@ public:
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const Point<spacedim> &                                     p) const = 0;
 
+  /**
+   * Map multiple points from the real point locations to points in reference
+   * locations. The functionality is essentially the same as looping over all
+   * points and calling the Mapping::transform_real_to_unit_cell() function
+   * for each point individually, but it can be much faster for certain
+   * mappings that implement a more specialized version such as
+   * MappingQGeneric. The only difference in behavior is that this function
+   * will never throw an ExcTransformationFailed() exception. If the
+   * transformation fails for `real_points[i]`, the returned `unit_points[i]`
+   * contains std::numeric_limits<double>::infinity() as the first entry.
+   */
+  virtual void
+  transform_points_real_to_unit_cell(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const ArrayView<const Point<spacedim>> &                    real_points,
+    ArrayView<Point<dim>> &unit_points) const;
+
   /**
    * Transform the point @p p on the real @p cell to the corresponding point
    * on the reference cell, and then project this point to a (dim-1)-dimensional

diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h
index 6dec09464394092024a8069027f55a78f4682896..501ef16289846dce5a5c59d97eac094dba4f1dac 100644 (file)
--- a/include/deal.II/fe/mapping_q_generic.h
+++ b/include/deal.II/fe/mapping_q_generic.h
@@ -181,6 +181,13 @@ public:
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const Point<spacedim> &p) const override;
 
+  // for documentation, see the Mapping base class
+  virtual void
+  transform_points_real_to_unit_cell(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const ArrayView<const Point<spacedim>> &                    real_points,
+    ArrayView<Point<dim>> &unit_points) const override;
+
   /**
    * @}
    */

diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc
index 5433bf7471a06a27d7ff513c44f251dcc9d70574..a70503a51ba4be233ce429766feedf6a3b7e2141 100644 (file)
--- a/source/fe/mapping.cc
+++ b/source/fe/mapping.cc
@@ -76,6 +76,30 @@ Mapping<dim, spacedim>::get_bounding_box(
 
 
 
+template <int dim, int spacedim>
+void
+Mapping<dim, spacedim>::transform_points_real_to_unit_cell(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const ArrayView<const Point<spacedim>> &                    real_points,
+  ArrayView<Point<dim>> &                                     unit_points) const
+{
+  AssertDimension(real_points.size(), unit_points.size());
+  for (unsigned int i = 0; i < real_points.size(); ++i)
+    {
+      try
+        {
+          unit_points[i] = transform_real_to_unit_cell(cell, real_points[i]);
+        }
+      catch (typename Mapping<dim>::ExcTransformationFailed &)
+        {
+          unit_points[i]    = Point<dim>();
+          unit_points[i][0] = std::numeric_limits<double>::infinity();
+        }
+    }
+}
+
+
+
 template <int dim, int spacedim>
 Point<dim - 1>
 Mapping<dim, spacedim>::project_real_point_to_unit_point_on_face(

diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
index 9e4eab6b044c45ec3c6f257441078a61ffdf0b2c..ab1058c97e1eb61423edb5525277ce54f5a7a922 100644 (file)
--- a/source/fe/mapping_q_generic.cc
+++ b/source/fe/mapping_q_generic.cc
@@ -29,6 +29,7 @@
 #include <deal.II/fe/mapping_q1.h>
 #include <deal.II/fe/mapping_q_generic.h>
 
+#include <deal.II/grid/grid_tools.h>
 #include <deal.II/grid/manifold_lib.h>
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/tria_iterator.h>
@@ -893,17 +894,15 @@ namespace internal
       /**
        * Implementation of transform_real_to_unit_cell for dim==spacedim
        */
-      template <int dim>
+      template <int dim, int spacedim>
       Point<dim>
       do_transform_real_to_unit_cell_internal(
-        const Point<dim> &                                  p,
+        const Point<spacedim> &                             p,
         const Point<dim> &                                  initial_p_unit,
-        const std::vector<Point<dim>> &                     points,
+        const std::vector<Point<spacedim>> &                points,
         const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
         const std::vector<unsigned int> &                   renumber)
       {
-        const unsigned int spacedim = dim;
-
         AssertDimension(points.size(),
                         Utilities::pow(polynomials_1d.size(), dim));
 
@@ -917,8 +916,8 @@ namespace internal
         // The shape values and derivatives of the mapping at this point are
         // previously computed.
 
-        Point<dim>                            p_unit = initial_p_unit;
-        std::pair<Point<dim>, Tensor<2, dim>> p_real =
+        Point<dim>                                      p_unit = initial_p_unit;
+        std::pair<Point<spacedim>, Tensor<2, spacedim>> p_real =
           compute_mapped_location_of_point(points,
                                            polynomials_1d,
                                            renumber,
@@ -967,6 +966,9 @@ namespace internal
         const double       eps                    = 1.e-11;
         const unsigned int newton_iteration_limit = 20;
 
+        Point<dim> invalid_point;
+        invalid_point[0] = std::numeric_limits<double>::infinity();
+
         unsigned int newton_iteration = 0;
         double       last_f_weighted_norm_square;
         do
@@ -979,9 +981,9 @@ namespace internal
             const Tensor<2, spacedim> &df = p_real.second;
 
             // Solve  [f'(x)]d=f(x)
-            AssertThrow(
-              determinant(df) > 0,
-              (typename Mapping<dim, spacedim>::ExcTransformationFailed()));
+            if (determinant(df) <= 0)
+              return invalid_point;
+
             const Tensor<2, spacedim> df_inverse = invert(df);
             const Tensor<1, spacedim> delta      = df_inverse * f;
 
@@ -1035,18 +1037,13 @@ namespace internal
                 else if (step_length > 0.05)
                   step_length /= 2;
                 else
-                  AssertThrow(
-                    false,
-                    (typename Mapping<dim,
-                                      spacedim>::ExcTransformationFailed()));
+                  return invalid_point;
               }
             while (true);
 
             ++newton_iteration;
             if (newton_iteration > newton_iteration_limit)
-              AssertThrow(
-                false,
-                (typename Mapping<dim, spacedim>::ExcTransformationFailed()));
+              return invalid_point;
             last_f_weighted_norm_square = (df_inverse * f).norm_square();
           }
         while (last_f_weighted_norm_square > eps * eps);
@@ -2295,7 +2292,58 @@ MappingQGeneric<dim, spacedim>::transform_real_to_unit_cell(
 
   // perform the Newton iteration and return the result. note that this
   // statement may throw an exception, which we simply pass up to the caller
-  return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
+  const Point<dim> p_unit =
+    this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
+  if (p_unit[0] == std::numeric_limits<double>::infinity())
+    AssertThrow(false,
+                (typename Mapping<dim, spacedim>::ExcTransformationFailed()));
+  return p_unit;
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingQGeneric<dim, spacedim>::transform_points_real_to_unit_cell(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const ArrayView<const Point<spacedim>> &                    real_points,
+  ArrayView<Point<dim>> &                                     unit_points) const
+{
+  AssertDimension(real_points.size(), unit_points.size());
+  const std::vector<Point<spacedim>> support_points =
+    this->compute_mapping_support_points(cell);
+
+  // From the chosen (high-order) support points, now only pick the first
+  // 2^dim points and construct an affine approximation from those.
+  const std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>>
+    affine_factors = GridTools::affine_cell_approximation<dim>(
+      ArrayView<const Point<spacedim>>(support_points.data(),
+                                       GeometryInfo<dim>::vertices_per_cell));
+  const DerivativeForm<1, spacedim, dim> A_inv =
+    affine_factors.first.covariant_form().transpose();
+
+  for (unsigned int i = 0; i < real_points.size(); ++i)
+    {
+      try
+        {
+          // Compute an initial guess by inverting the affine approximation
+          // A * x_unit + b = x_real
+          Point<dim> initial_guess(
+            apply_transformation(A_inv,
+                                 real_points[i] - affine_factors.second));
+          unit_points[i] = internal::MappingQGenericImplementation::
+            do_transform_real_to_unit_cell_internal<dim, spacedim>(
+              real_points[i],
+              GeometryInfo<dim>::project_to_unit_cell(initial_guess),
+              support_points,
+              polynomials_1d,
+              renumber_lexicographic_to_hierarchic);
+        }
+      catch (typename Mapping<dim>::ExcTransformationFailed &)
+        {
+          unit_points[i][0] = std::numeric_limits<double>::infinity();
+        }
+    }
 }