]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Break up first part of MappingQ setup to make it parallelizable
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 6 Apr 2020 05:25:08 +0000 (07:25 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 9 May 2020 10:34:51 +0000 (12:34 +0200)
include/deal.II/matrix_free/mapping_info.templates.h

index 00b2f3f19f1617434802da4bca6445183e06604c..b5460636c5e3f8fc1e06ac06d790846a580688f5 100644 (file)
@@ -1099,6 +1099,154 @@ namespace internal
 
 
 
+      /**
+       * This invokes the FEValues part of the initialization of MappingQ,
+       * storing the resulting quadrature points and an initial representation
+       * of Jacobians in two arrays.
+       */
+      template <int dim>
+      void
+      mapping_q_query_fe_values(
+        const unsigned int                                        begin_cell,
+        const unsigned int                                        end_cell,
+        const MappingQGeneric<dim> &                              mapping_q,
+        const dealii::Triangulation<dim> &                        tria,
+        const std::vector<std::pair<unsigned int, unsigned int>> &cell_array,
+        const double                                              jacobian_size,
+        std::vector<GeometryType> &preliminary_cell_type,
+        AlignedVector<double> &    plain_quadrature_points,
+        AlignedVector<std::array<Tensor<2, dim>, dim + 1>>
+          &jacobians_on_stencil)
+      {
+        const unsigned int mapping_degree = mapping_q.get_degree();
+        FE_Nothing<dim>    dummy_fe;
+        QGaussLobatto<dim> quadrature(mapping_degree + 1);
+        const unsigned int n_mapping_points =
+          Utilities::pow(mapping_degree + 1, dim);
+
+        FEValues<dim> fe_values(mapping_q,
+                                dummy_fe,
+                                quadrature,
+                                update_quadrature_points | update_jacobians);
+
+        for (unsigned int cell = begin_cell; cell < end_cell; ++cell)
+          {
+            typename dealii::Triangulation<dim>::cell_iterator cell_it(
+              &tria, cell_array[cell].first, cell_array[cell].second);
+            fe_values.reinit(cell_it);
+            for (unsigned int d = 0; d < dim; ++d)
+              for (unsigned int q = 0; q < n_mapping_points; ++q)
+                plain_quadrature_points[(cell * dim + d) * n_mapping_points +
+                                        q] = fe_values.quadrature_point(q)[d];
+
+            // store the first, second, n-th and n^2-th one along a
+            // stencil-like pattern
+            std::array<Tensor<2, dim, double>, dim + 1> &my_jacobians =
+              jacobians_on_stencil[cell];
+            my_jacobians[0] = Tensor<2, dim, double>(fe_values.jacobian(0));
+            for (unsigned int d = 0, skip = 1; d < dim;
+                 ++d, skip *= (mapping_degree + 1))
+              my_jacobians[1 + d] =
+                Tensor<2, dim, double>(fe_values.jacobian(skip));
+
+            // check whether cell is Cartesian/affine/general
+            GeometryType type = cartesian;
+            for (unsigned int d = 0; d < dim; ++d)
+              for (unsigned int e = 0; e < dim; ++e)
+                if (d != e)
+                  if (std::abs(my_jacobians[0][d][e]) > 1e-12 * jacobian_size)
+                    type = affine;
+
+            for (unsigned int q = 1; q < n_mapping_points; ++q)
+              for (unsigned int d = 0; d < dim; ++d)
+                for (unsigned int e = 0; e < dim; ++e)
+                  if (std::abs(fe_values.jacobian(q)[d][e] -
+                               fe_values.jacobian(0)[d][e]) >
+                      1e-12 * jacobian_size)
+                    {
+                      type = general;
+                      goto endloop;
+                    }
+          endloop:
+            preliminary_cell_type[cell] = type;
+          }
+      }
+
+
+
+      template <int dim>
+      std::vector<unsigned int>
+      mapping_q_find_compression(
+        const double jacobian_size,
+        const AlignedVector<std::array<Tensor<2, dim>, dim + 1>>
+          &                          jacobians_on_stencil,
+        const unsigned int           n_mapping_points,
+        const AlignedVector<double> &plain_quadrature_points,
+        std::vector<GeometryType> &  preliminary_cell_type)
+      {
+        std::vector<unsigned int> cell_data_index(jacobians_on_stencil.size());
+
+        // we include a map to store some compressed information about the
+        // Jacobians which we collect by a stencil-like pattern around the
+        // first quadrature point on the cell - we use a relatively coarse
+        // tolerance to account for some inaccuracies in the manifold
+        // evaluation
+        const FPArrayComparator<double> comparator(1e4 * jacobian_size);
+        std::map<std::array<Tensor<2, dim>, dim + 1>,
+                 unsigned int,
+                 FPArrayComparator<double>>
+          compressed_jacobians(comparator);
+
+        unsigned int n_data_buckets = 0;
+        for (unsigned int cell = 0; cell < jacobians_on_stencil.size(); ++cell)
+          {
+            // check in the map for the index of this cell
+            auto inserted = compressed_jacobians.insert(
+              std::make_pair(jacobians_on_stencil[cell], cell));
+            bool add_this_cell = inserted.second;
+            if (inserted.second == false)
+              {
+                // check if the found duplicate really is a translation and
+                // the similarity identified by the map is not by accident
+                double        max_distance = 0;
+                const double *ptr_origin =
+                  plain_quadrature_points.data() +
+                  inserted.first->second * dim * n_mapping_points;
+                const double *ptr_mine = plain_quadrature_points.data() +
+                                         cell * dim * n_mapping_points;
+                for (unsigned int d = 0; d < dim; ++d)
+                  {
+                    const double translate_d =
+                      ptr_origin[d * n_mapping_points] -
+                      ptr_mine[d * n_mapping_points];
+                    for (unsigned int q = 1; q < n_mapping_points; ++q)
+                      max_distance =
+                        std::max(std::abs(ptr_origin[d * n_mapping_points + q] -
+                                          ptr_mine[d * n_mapping_points + q] -
+                                          translate_d),
+                                 max_distance);
+                  }
+
+                // this is not a duplicate, must add it again
+                if (max_distance > 1e-10 * jacobian_size)
+                  add_this_cell = true;
+              }
+            if (add_this_cell)
+              cell_data_index[cell] = n_data_buckets++;
+            else
+              {
+                cell_data_index[cell] = cell_data_index[inserted.first->second];
+                // make sure that the cell type is the same as in the original
+                // field, despite possible small differences due to roundoff
+                // and the tolerances we use
+                preliminary_cell_type[cell] =
+                  preliminary_cell_type[inserted.first->second];
+              }
+          }
+        return cell_data_index;
+      }
+
+
       /**
        * This evaluates the mapping information on a range of cells calling
        * into the tensor product interpolators of the matrix-free framework,
@@ -2366,114 +2514,26 @@ namespace internal
 
       const double jacobian_size = ExtractCellHelper::get_jacobian_size(tria);
 
-      std::vector<unsigned int> cell_data_index(cell_array.size());
+      std::vector<unsigned int> cell_data_index;
       std::vector<GeometryType> preliminary_cell_type(cell_array.size());
       {
-        FE_Nothing<dim>    dummy_fe;
-        QGaussLobatto<dim> quadrature(mapping_degree + 1);
-
-        FEValues<dim> fe_values(*mapping_q,
-                                dummy_fe,
-                                quadrature,
-                                update_quadrature_points | update_jacobians);
-
-        // we include a map to store some compressed information about the
-        // Jacobians which we collect by a stencil-like pattern around the
-        // first quadrature point on the cell - we use a relatively coarse
-        // tolerance to account for some inaccuracies in the manifold
-        // evaluation
-        const FPArrayComparator<double> comparator(1e4 * jacobian_size);
-        std::map<std::array<Tensor<2, dim>, dim + 1>,
-                 unsigned int,
-                 FPArrayComparator<double>>
-          compressed_jacobians(comparator);
-
-        unsigned int n_data_buckets = 0;
-        for (unsigned int cell = 0; cell < cell_array.size(); ++cell)
-          {
-            typename dealii::Triangulation<dim>::cell_iterator cell_it(
-              &tria, cell_array[cell].first, cell_array[cell].second);
-            fe_values.reinit(cell_it);
-            for (unsigned int d = 0; d < dim; ++d)
-              for (unsigned int q = 0; q < n_mapping_points; ++q)
-                plain_quadrature_points[(cell * dim + d) * n_mapping_points +
-                                        q] = fe_values.quadrature_point(q)[d];
-
-            // store the first, second, n-th and n^2-th one along a
-            // stencil-like pattern
-            std::array<Tensor<2, dim, double>, dim + 1> jacobians_on_stencil;
-            jacobians_on_stencil[0] =
-              Tensor<2, dim, double>(fe_values.jacobian(0));
-            for (unsigned int d = 0, skip = 1; d < dim;
-                 ++d, skip *= (mapping_degree + 1))
-              jacobians_on_stencil[1 + d] =
-                Tensor<2, dim, double>(fe_values.jacobian(skip));
-
-            // check in the map for the index of this cell
-            auto inserted = compressed_jacobians.insert(
-              std::make_pair(jacobians_on_stencil, cell));
-            bool add_this_cell = inserted.second;
-            if (inserted.second == false)
-              {
-                // check if the found duplicate really is a translation and
-                // the similarity identified by the map is not by accident
-                double        max_distance = 0;
-                const double *ptr_origin =
-                  plain_quadrature_points.data() +
-                  inserted.first->second * dim * n_mapping_points;
-                const double *ptr_mine = plain_quadrature_points.data() +
-                                         cell * dim * n_mapping_points;
-                for (unsigned int d = 0; d < dim; ++d)
-                  {
-                    const double translate_d =
-                      ptr_origin[d * n_mapping_points] -
-                      ptr_mine[d * n_mapping_points];
-                    for (unsigned int q = 1; q < n_mapping_points; ++q)
-                      max_distance =
-                        std::max(std::abs(ptr_origin[d * n_mapping_points + q] -
-                                          ptr_mine[d * n_mapping_points + q] -
-                                          translate_d),
-                                 max_distance);
-                  }
-
-                // this is not a duplicate, must add it again
-                if (max_distance > 1e-10 * jacobian_size)
-                  add_this_cell = true;
-              }
-
-            if (add_this_cell == true)
-              {
-                // check whether cell is Cartesian/affine/general
-                GeometryType type = cartesian;
-                for (unsigned int d = 0; d < dim; ++d)
-                  for (unsigned int e = 0; e < dim; ++e)
-                    if (d != e)
-                      if (std::abs(inserted.first->first[0][d][e]) >
-                          1e-12 * jacobian_size)
-                        type = affine;
-
-                for (unsigned int q = 1; q < n_mapping_points; ++q)
-                  for (unsigned int d = 0; d < dim; ++d)
-                    for (unsigned int e = 0; e < dim; ++e)
-                      if (std::abs(fe_values.jacobian(q)[d][e] -
-                                   fe_values.jacobian(0)[d][e]) >
-                          1e-12 * jacobian_size)
-                        {
-                          type = general;
-                          goto endloop;
-                        }
-              endloop:
-                cell_data_index[cell]       = n_data_buckets;
-                preliminary_cell_type[cell] = type;
-                ++n_data_buckets;
-              }
-            else
-              {
-                cell_data_index[cell] = cell_data_index[inserted.first->second];
-                preliminary_cell_type[cell] =
-                  preliminary_cell_type[inserted.first->second];
-              }
-          }
+        AlignedVector<std::array<Tensor<2, dim>, dim + 1>> jacobians_on_stencil(
+          cell_array.size());
+        ExtractCellHelper::mapping_q_query_fe_values(0,
+                                                     cell_array.size(),
+                                                     *mapping_q,
+                                                     tria,
+                                                     cell_array,
+                                                     jacobian_size,
+                                                     preliminary_cell_type,
+                                                     plain_quadrature_points,
+                                                     jacobians_on_stencil);
+        cell_data_index =
+          ExtractCellHelper::mapping_q_find_compression(jacobian_size,
+                                                        jacobians_on_stencil,
+                                                        n_mapping_points,
+                                                        plain_quadrature_points,
+                                                        preliminary_cell_type);
       }
 
       // step 2: compute the appropriate evaluation matrices for cells and

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.