From: Peter Munch Date: Tue, 1 Dec 2020 18:54:26 +0000 (+0100) Subject: MatrixFree: split up ranges according to active_fe_indices X-Git-Tag: v9.3.0-rc1~816^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11260%2Fhead;p=dealii.git MatrixFree: split up ranges according to active_fe_indices --- diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 4438f81211..e04a41b761 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -1481,6 +1481,51 @@ public: const unsigned int fe_index, const unsigned int dof_handler_index = 0) const; + /** + * In the hp adaptive case, a subrange of internal faces as computed during + * loop() might contain internal faces with elements of different active + * fe indices. Use this function to compute what the subrange for a given pair + * of active fe indices is. + */ + std::pair + create_inner_face_subrange_hp_by_index( + const std::pair &range, + const unsigned int fe_index_interior, + const unsigned int fe_index_exterior, + const unsigned int dof_handler_index = 0) const; + + /** + * In the hp adaptive case, a subrange of boundary faces as computed during + * loop() might contain boundary faces with elements of different active + * fe indices. Use this function to compute what the subrange for a given + * active fe indices is. + */ + std::pair + create_boundary_face_subrange_hp_by_index( + const std::pair &range, + const unsigned int fe_index, + const unsigned int dof_handler_index = 0) const; + + /** + * In the hp adaptive case, return number of active_fe_indices. + */ + unsigned int + n_active_fe_indices() const; + + /** + * In the hp adaptive case, return the active_fe_index of a cell range. + */ + unsigned int + get_cell_active_fe_index( + const std::pair range) const; + + /** + * In the hp adaptive case, return the active_fe_index of a face range. + */ + unsigned int + get_face_active_fe_index(const std::pair range, + const bool is_interior_face = true) const; + //@} /** @@ -2471,6 +2516,75 @@ MatrixFree::at_irregular_cell( +template +unsigned int +MatrixFree::n_active_fe_indices() const +{ + return shape_info.size(2); +} + + +template +unsigned int +MatrixFree::get_cell_active_fe_index( + const std::pair range) const +{ + const auto &fe_indices = dof_info[0].cell_active_fe_index; + + if (fe_indices.empty() == true) + return 0; + + const auto index = fe_indices[range.first]; + + for (unsigned int i = range.first; i < range.second; ++i) + AssertDimension(index, fe_indices[i]); + + return index; +} + + + +template +unsigned int +MatrixFree::get_face_active_fe_index( + const std::pair range, + const bool is_interior_face) const +{ + const auto &fe_indices = dof_info[0].cell_active_fe_index; + + if (fe_indices.empty() == true) + return 0; + + if (is_interior_face) + { + const unsigned int index = + fe_indices[face_info.faces[range.first].cells_interior[0] / + VectorizedArrayType::size()]; + + for (unsigned int i = range.first; i < range.second; ++i) + AssertDimension(index, + fe_indices[face_info.faces[i].cells_interior[0] / + VectorizedArrayType::size()]); + + return index; + } + else + { + const unsigned int index = + fe_indices[face_info.faces[range.first].cells_exterior[0] / + VectorizedArrayType::size()]; + + for (unsigned int i = range.first; i < range.second; ++i) + AssertDimension(index, + fe_indices[face_info.faces[i].cells_exterior[0] / + VectorizedArrayType::size()]); + + return index; + } +} + + + template inline unsigned int MatrixFree::n_components_filled( @@ -4525,8 +4639,17 @@ namespace internal cell(const std::pair &cell_range) override { if (cell_function != nullptr && cell_range.second > cell_range.first) - (container.* - cell_function)(matrix_free, this->dst, this->src, cell_range); + for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i) + { + const auto cell_subrange = + matrix_free.create_cell_subrange_hp_by_index(cell_range, i); + + if (cell_subrange.second <= cell_subrange.first) + continue; + + (container.* + cell_function)(matrix_free, this->dst, this->src, cell_subrange); + } } // Runs the assembler on interior faces. If no function is given, nothing @@ -4535,8 +4658,19 @@ namespace internal face(const std::pair &face_range) override { if (face_function != nullptr && face_range.second > face_range.first) - (container.* - face_function)(matrix_free, this->dst, this->src, face_range); + for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i) + for (unsigned int j = 0; j < matrix_free.n_active_fe_indices(); ++j) + { + const auto face_subrange = + matrix_free.create_inner_face_subrange_hp_by_index(face_range, + i, + j); + + if (face_subrange.second <= face_subrange.first) + continue; + (container.* + face_function)(matrix_free, this->dst, this->src, face_subrange); + } } // Runs the assembler on boundary faces. If no function is given, nothing @@ -4545,8 +4679,20 @@ namespace internal boundary(const std::pair &face_range) override { if (boundary_function != nullptr && face_range.second > face_range.first) - (container.* - boundary_function)(matrix_free, this->dst, this->src, face_range); + for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i) + { + const auto face_subrange = + matrix_free.create_boundary_face_subrange_hp_by_index(face_range, + i); + + if (face_subrange.second <= face_subrange.first) + continue; + + (container.*boundary_function)(matrix_free, + this->dst, + this->src, + face_subrange); + } } // Starts the communication for the update ghost values operation. We diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index e2c0e006f3..de3b9ef39a 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -85,12 +85,17 @@ std::pair MatrixFree::create_cell_subrange_hp_by_index( const std::pair &range, const unsigned int fe_index, - const unsigned int vector_component) const + const unsigned int dof_handler_index) const { - AssertIndexRange(fe_index, dof_info[vector_component].max_fe_index); + if (dof_info[dof_handler_index].max_fe_index == 0) + return range; + + AssertIndexRange(fe_index, dof_info[dof_handler_index].max_fe_index); const std::vector &fe_indices = - dof_info[vector_component].cell_active_fe_index; - if (fe_indices.empty() == true) + dof_info[dof_handler_index].cell_active_fe_index; + + if (fe_indices.empty() == true || + dof_handlers[dof_handler_index]->get_fe_collection().size() == 1) return range; else { @@ -124,14 +129,145 @@ MatrixFree::create_cell_subrange_hp_by_index( +namespace +{ + class FaceRangeCompartor + { + public: + FaceRangeCompartor(const std::vector &fe_indices, + const bool include) + : fe_indices(fe_indices) + , include(include) + {} + + template + bool + operator()(const internal::MatrixFreeFunctions::FaceToCellTopology< + vectorization_width> &face, + const unsigned int & fe_index) + { + const unsigned int fe_index_face = + fe_indices[face.cells_interior[0] / vectorization_width]; + + return fe_index_face < fe_index || + (include ? (fe_index_face == fe_index) : false); + } + + template + bool + operator()(const internal::MatrixFreeFunctions::FaceToCellTopology< + vectorization_width> & face, + const std::pair &fe_index) + { + const std::pair fe_index_face = { + fe_indices[face.cells_interior[0] / vectorization_width], + fe_indices[face.cells_exterior[0] / vectorization_width]}; + + return include ? (fe_index_face <= fe_index) : (fe_index_face < fe_index); + } + + private: + const std::vector &fe_indices; + const bool include; + }; +} // namespace + + + +template +std::pair +MatrixFree:: + create_inner_face_subrange_hp_by_index( + const std::pair &range, + const unsigned int fe_index_interior, + const unsigned int fe_index_exterior, + const unsigned int dof_handler_index) const +{ + if (dof_info[dof_handler_index].max_fe_index == 0 || + dof_handlers[dof_handler_index]->get_fe_collection().size() == 1) + return range; + + AssertIndexRange(fe_index_interior, dof_info[dof_handler_index].max_fe_index); + AssertIndexRange(fe_index_exterior, dof_info[dof_handler_index].max_fe_index); + const std::vector &fe_indices = + dof_info[dof_handler_index].cell_active_fe_index; + if (fe_indices.empty() == true) + return range; + else + { + std::pair return_range; + return_range.first = + std::lower_bound(face_info.faces.begin() + range.first, + face_info.faces.begin() + range.second, + std::pair{ + fe_index_interior, fe_index_exterior}, + FaceRangeCompartor(fe_indices, false)) - + face_info.faces.begin(); + return_range.second = + std::lower_bound(face_info.faces.begin() + return_range.first, + face_info.faces.begin() + range.second, + std::pair{ + fe_index_interior, fe_index_exterior}, + FaceRangeCompartor(fe_indices, true)) - + face_info.faces.begin(); + Assert(return_range.first >= range.first && + return_range.second <= range.second, + ExcInternalError()); + return return_range; + } +} + + + +template +std::pair +MatrixFree:: + create_boundary_face_subrange_hp_by_index( + const std::pair &range, + const unsigned int fe_index, + const unsigned int dof_handler_index) const +{ + if (dof_info[dof_handler_index].max_fe_index == 0 || + dof_handlers[dof_handler_index]->get_fe_collection().size() == 1) + return range; + + AssertIndexRange(fe_index, dof_info[dof_handler_index].max_fe_index); + const std::vector &fe_indices = + dof_info[dof_handler_index].cell_active_fe_index; + if (fe_indices.empty() == true) + return range; + else + { + std::pair return_range; + return_range.first = + std::lower_bound(face_info.faces.begin() + range.first, + face_info.faces.begin() + range.second, + fe_index, + FaceRangeCompartor(fe_indices, false)) - + face_info.faces.begin(); + return_range.second = + std::lower_bound(face_info.faces.begin() + return_range.first, + face_info.faces.begin() + range.second, + fe_index, + FaceRangeCompartor(fe_indices, true)) - + face_info.faces.begin(); + Assert(return_range.first >= range.first && + return_range.second <= range.second, + ExcInternalError()); + return return_range; + } +} + + + template void MatrixFree::renumber_dofs( std::vector &renumbering, - const unsigned int vector_component) + const unsigned int dof_handler_index) { - AssertIndexRange(vector_component, dof_info.size()); - dof_info[vector_component].compute_dof_renumbering(renumbering); + AssertIndexRange(dof_handler_index, dof_info.size()); + dof_info[dof_handler_index].compute_dof_renumbering(renumbering); } @@ -1546,7 +1682,7 @@ MatrixFree::initialize_indices( #ifdef DEAL_II_WITH_MPI { // non-buffering mode is only supported if the indices of all cells are - // contiguous for all dof_info objects and hp is not enabled. + // contiguous for all dof_info objects. bool is_non_buffering_sm_supported = true; for (const auto &di : dof_info) { diff --git a/tests/simplex/matrix_free_02.cc b/tests/simplex/matrix_free_02.cc index 7cb9ade094..d884c35976 100644 --- a/tests/simplex/matrix_free_02.cc +++ b/tests/simplex/matrix_free_02.cc @@ -47,15 +47,122 @@ #include "../tests.h" -#include "./simplex_grids.h" - using namespace dealii; +namespace dealii +{ + namespace GridGenerator + { + template + void + subdivided_hyper_rectangle_with_simplices_mix( + Triangulation & tria, + const std::vector &repetitions, + const Point & p1, + const Point & p2, + const bool colorize = false) + { + AssertDimension(dim, spacedim); + + AssertThrow(colorize == false, ExcNotImplemented()); + + std::vector> vertices; + std::vector> cells; + + if (dim == 2) + { + // determine cell sizes + const Point dx((p2[0] - p1[0]) / repetitions[0], + (p2[1] - p1[1]) / repetitions[1]); + + // create vertices + for (unsigned int j = 0; j <= repetitions[1]; ++j) + for (unsigned int i = 0; i <= repetitions[0]; ++i) + vertices.push_back( + Point(p1[0] + dx[0] * i, p1[1] + dx[1] * j)); + + // create cells + for (unsigned int j = 0; j < repetitions[1]; ++j) + for (unsigned int i = 0; i < repetitions[0]; ++i) + { + // create reference QUAD cell + std::array quad{{ + (j + 0) * (repetitions[0] + 1) + i + 0, // + (j + 0) * (repetitions[0] + 1) + i + 1, // + (j + 1) * (repetitions[0] + 1) + i + 0, // + (j + 1) * (repetitions[0] + 1) + i + 1 // + }}; // + + if (j < repetitions[1] / 2 && i < repetitions[0] / 2) + { + CellData quad_; + quad_.vertices = {quad[0], quad[1], quad[2], quad[3]}; + cells.push_back(quad_); + + continue; + } + + // TRI cell 0 + { + CellData tri; + tri.vertices = {quad[0], quad[1], quad[2]}; + cells.push_back(tri); + } + + // TRI cell 1 + { + CellData tri; + tri.vertices = {quad[3], quad[2], quad[1]}; + cells.push_back(tri); + } + } + } + else + { + AssertThrow(colorize == false, ExcNotImplemented()); + } + + // actually create triangulation + tria.create_triangulation(vertices, cells, SubCellData()); + } + + + template + void + subdivided_hyper_cube_with_simplices_mix(Triangulation &tria, + const unsigned int repetitions, + const double p1 = 0.0, + const double p2 = 1.0, + const bool colorize = false) + { + if (dim == 2) + { + subdivided_hyper_rectangle_with_simplices_mix( + tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2}, colorize); + } + else if (dim == 3) + { + subdivided_hyper_rectangle_with_simplices_mix( + tria, + {{repetitions, repetitions, repetitions}}, + {p1, p1, p1}, + {p2, p2, p2}, + colorize); + } + else + { + AssertThrow(false, ExcNotImplemented()) + } + } + } // namespace GridGenerator +} // namespace dealii + template class PoissonOperator { public: - using VectorType = LinearAlgebra::distributed::Vector; + using VectorType = LinearAlgebra::distributed::Vector; + using FECellIntegrator = FEEvaluation; PoissonOperator(const MatrixFree &matrix_free, const bool do_helmholtz) @@ -75,24 +182,17 @@ public: const int dummy = 0; matrix_free.template cell_loop( - [&](const auto &data, auto &dst, const auto &, const auto cell_range) { - for (unsigned int i = 0; i < 2; ++i) - { - const auto cell_subrange = - data.create_cell_subrange_hp_by_index(cell_range, i); + [&](const auto &data, auto &dst, const auto &, const auto range) { + const auto i = data.get_cell_active_fe_index(range); + FECellIntegrator phi(matrix_free, 0, 0, 0, i, i); - FEEvaluation phi(matrix_free, 0, 0, 0, i, i); - - for (unsigned int cell = cell_subrange.first; - cell < cell_subrange.second; - ++cell) - { - phi.reinit(cell); - for (unsigned int q = 0; q < phi.n_q_points; ++q) - phi.submit_value(1.0, q); + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_value(1.0, q); - phi.integrate_scatter(true, false, dst); - } + phi.integrate_scatter(true, false, dst); } }, vec, @@ -105,30 +205,24 @@ public: vmult(VectorType &dst, const VectorType &src) const { matrix_free.template cell_loop( - [&](const auto &data, auto &dst, const auto &src, const auto cell_range) { - for (unsigned int i = 0; i < 2; ++i) + [&](const auto &data, auto &dst, const auto &src, const auto range) { + const auto i = data.get_cell_active_fe_index(range); + FECellIntegrator phi(matrix_free, 0, 0, 0, i, i); + + for (unsigned int cell = range.first; cell < range.second; ++cell) { - const auto cell_subrange = - data.create_cell_subrange_hp_by_index(cell_range, i); + phi.reinit(cell); + phi.gather_evaluate(src, do_helmholtz, true); - FEEvaluation phi(matrix_free, 0, 0, 0, i, i); - for (unsigned int cell = cell_subrange.first; - cell < cell_subrange.second; - ++cell) + for (unsigned int q = 0; q < phi.n_q_points; ++q) { - phi.reinit(cell); - phi.gather_evaluate(src, do_helmholtz, true); + if (do_helmholtz) + phi.submit_value(phi.get_value(q), q); - for (unsigned int q = 0; q < phi.n_q_points; ++q) - { - if (do_helmholtz) - phi.submit_value(phi.get_value(q), q); - - phi.submit_gradient(phi.get_gradient(q), q); - } - - phi.integrate_scatter(do_helmholtz, true, dst); + phi.submit_gradient(phi.get_gradient(q), q); } + + phi.integrate_scatter(do_helmholtz, true, dst); } }, dst, diff --git a/tests/simplex/matrix_free_range_iteration_01.cc b/tests/simplex/matrix_free_range_iteration_01.cc new file mode 100644 index 0000000000..54479d8815 --- /dev/null +++ b/tests/simplex/matrix_free_range_iteration_01.cc @@ -0,0 +1,164 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// Test ShapeData for Simplex::FE_P and Simplex::QGauss + +#include + +#include + +#include +#include + +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +test() +{ + const unsigned int degree = 1; + + Triangulation tria; + GridGenerator::subdivided_hyper_cube(tria, 4); + + hp::FECollection fe{FE_Q<2>(degree), FE_Q<2>(degree)}; + MappingQGeneric mapping(1); + QGauss quadrature(degree + 1); + + DoFHandler dof_handler(tria); + + unsigned int counter = 0; + + for (auto &cell : dof_handler) + if (counter++ < tria.n_cells() / 2) + cell.set_active_fe_index(1); + + dof_handler.distribute_dofs(fe); + + + typename MatrixFree::AdditionalData additional_data; + additional_data.mapping_update_flags_boundary_faces = + update_gradients | update_values; + additional_data.mapping_update_flags_inner_faces = + update_gradients | update_values; + additional_data.initialize_mapping = false; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + + AffineConstraints constraints; + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraints, quadrature, additional_data); + + using VectorType = Vector; + + VectorType src, dst; + matrix_free.initialize_dof_vector(src); + matrix_free.initialize_dof_vector(dst); + + std::vector> cells(2); + std::vector>> boundary_faces(2); + + std::vector>>> inner_faces( + 2); + for (unsigned int i = 0; i < 2; ++i) + inner_faces[i].resize(2); + + matrix_free.template loop( + [&](const auto &data, auto &dst, const auto &src, const auto range) { + const auto i = data.get_cell_active_fe_index(range); + + for (unsigned int cell = range.first; cell < range.second; ++cell) + for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell); + ++v) + cells[i].push_back(data.get_cell_iterator(cell, v)->id()); + }, + [&](const auto &data, auto &dst, const auto &src, const auto range) { + const auto i = data.get_face_active_fe_index(range, true); + const auto j = data.get_face_active_fe_index(range, false); + + for (unsigned int face = range.first; face < range.second; ++face) + for (unsigned int v = 0; v < data.n_active_entries_per_face_batch(face); + ++v) + inner_faces[i][j].emplace_back( + data.get_face_iterator(face, v, true).first->id(), + data.get_face_iterator(face, v, false).first->id()); + }, + [&](const auto &data, auto &dst, const auto &src, const auto range) { + const auto i = data.get_face_active_fe_index(range); + + for (unsigned int face = range.first; face < range.second; ++face) + for (unsigned int v = 0; v < data.n_active_entries_per_face_batch(face); + ++v) + boundary_faces[i].emplace_back( + data.get_face_iterator(face, v).first->id(), + data.get_face_iterator(face, v).second); + }, + dst, + src); + + deallog << "CELL categories:" << std::endl; + for (auto &i : cells) + { + std::sort(i.begin(), i.end()); + + for (const auto &j : i) + deallog << j << " "; + deallog << std::endl; + } + + deallog << "BOUNDARY_FACE categories:" << std::endl; + for (auto &i : boundary_faces) + { + std::sort(i.begin(), i.end()); + for (const auto &j : i) + deallog << j.first << "@" << j.second << " "; + deallog << std::endl; + } + + deallog << "INNER_FACE categories:" << std::endl; + for (auto &k : inner_faces) + { + for (auto &i : k) + { + std::sort(i.begin(), i.end()); + for (const auto &j : i) + deallog << j.first << "@" << j.second << " "; + deallog << std::endl; + } + } +} + +int +main() +{ + initlog(); + + test<2>(); +} diff --git a/tests/simplex/matrix_free_range_iteration_01.output b/tests/simplex/matrix_free_range_iteration_01.output new file mode 100644 index 0000000000..3229c0c705 --- /dev/null +++ b/tests/simplex/matrix_free_range_iteration_01.output @@ -0,0 +1,12 @@ + +DEAL::CELL categories: +DEAL::8_0: 9_0: 10_0: 11_0: 12_0: 13_0: 14_0: 15_0: +DEAL::0_0: 1_0: 2_0: 3_0: 4_0: 5_0: 6_0: 7_0: +DEAL::BOUNDARY_FACE categories: +DEAL::8_0:@0 11_0:@1 12_0:@0 12_0:@3 13_0:@3 14_0:@3 15_0:@1 15_0:@3 +DEAL::0_0:@0 0_0:@2 1_0:@2 2_0:@2 3_0:@1 3_0:@2 4_0:@0 7_0:@1 +DEAL::INNER_FACE categories: +DEAL::9_0:@8_0: 10_0:@9_0: 11_0:@10_0: 12_0:@8_0: 13_0:@9_0: 13_0:@12_0: 14_0:@10_0: 14_0:@13_0: 15_0:@11_0: 15_0:@14_0: +DEAL:: +DEAL::4_0:@8_0: 5_0:@9_0: 6_0:@10_0: 7_0:@11_0: +DEAL::1_0:@0_0: 2_0:@1_0: 3_0:@2_0: 4_0:@0_0: 5_0:@1_0: 5_0:@4_0: 6_0:@2_0: 6_0:@5_0: 7_0:@3_0: 7_0:@6_0: