From c0e80eeb8fbb57cf99e9c7ad4055df04a1cebc60 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 6 Apr 2020 07:25:08 +0200 Subject: [PATCH] Break up first part of MappingQ setup to make it parallelizable --- .../matrix_free/mapping_info.templates.h | 272 +++++++++++------- 1 file changed, 166 insertions(+), 106 deletions(-) diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 00b2f3f19f..b5460636c5 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -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 + void + mapping_q_query_fe_values( + const unsigned int begin_cell, + const unsigned int end_cell, + const MappingQGeneric & mapping_q, + const dealii::Triangulation & tria, + const std::vector> &cell_array, + const double jacobian_size, + std::vector &preliminary_cell_type, + AlignedVector & plain_quadrature_points, + AlignedVector, dim + 1>> + &jacobians_on_stencil) + { + const unsigned int mapping_degree = mapping_q.get_degree(); + FE_Nothing dummy_fe; + QGaussLobatto quadrature(mapping_degree + 1); + const unsigned int n_mapping_points = + Utilities::pow(mapping_degree + 1, dim); + + FEValues 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::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, 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 + std::vector + mapping_q_find_compression( + const double jacobian_size, + const AlignedVector, dim + 1>> + & jacobians_on_stencil, + const unsigned int n_mapping_points, + const AlignedVector &plain_quadrature_points, + std::vector & preliminary_cell_type) + { + std::vector 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 comparator(1e4 * jacobian_size); + std::map, dim + 1>, + unsigned int, + FPArrayComparator> + 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 cell_data_index(cell_array.size()); + std::vector cell_data_index; std::vector preliminary_cell_type(cell_array.size()); { - FE_Nothing dummy_fe; - QGaussLobatto quadrature(mapping_degree + 1); - - FEValues 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 comparator(1e4 * jacobian_size); - std::map, dim + 1>, - unsigned int, - FPArrayComparator> - compressed_jacobians(comparator); - - unsigned int n_data_buckets = 0; - for (unsigned int cell = 0; cell < cell_array.size(); ++cell) - { - typename dealii::Triangulation::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, 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, 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 -- 2.39.5