From: Jean-Paul Pelteret Date: Thu, 2 May 2019 17:32:41 +0000 (+0200) Subject: Implement iterator range support in MeshWorker::mesh_loop() X-Git-Tag: v9.1.0-rc1~126^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6ea60d738f6f51a0c140de099615ae1dd1945183;p=dealii.git Implement iterator range support in MeshWorker::mesh_loop() --- diff --git a/doc/news/changes/minor/20190502Jean-PaulPelteret-2 b/doc/news/changes/minor/20190502Jean-PaulPelteret-2 new file mode 100644 index 0000000000..ae68edd751 --- /dev/null +++ b/doc/news/changes/minor/20190502Jean-PaulPelteret-2 @@ -0,0 +1,5 @@ +Improved: The `Meshworker::mesh_loop()` function is now capable of working with an +IteratorRange, and also supports iterator ranges constructed from +FilteredIterators. +
+(Jean-Paul Pelteret, 2019/05/02) diff --git a/include/deal.II/meshworker/mesh_loop.h b/include/deal.II/meshworker/mesh_loop.h index 9aa0ff4369..d1c7ba3169 100644 --- a/include/deal.II/meshworker/mesh_loop.h +++ b/include/deal.II/meshworker/mesh_loop.h @@ -32,6 +32,7 @@ #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -40,6 +41,82 @@ class TriaActiveIterator; namespace MeshWorker { + namespace internal + { + /** + * A helper class to provide a type definition for the underlying cell + * iterator type. + */ + template + struct CellIteratorBaseType + { + /** + * Type definition for the cell iterator type. + */ + using type = CellIteratorType; + }; + + /** + * A helper class to provide a type definition for the underlying cell + * iterator type. + * + * This specialization is for IteratorRange, which may have either a + * TriaActiveIterator or a FilteredIterator as its base type. + */ + template + struct CellIteratorBaseType< + typename IteratorRange::IteratorOverIterators> + { + /** + * Type definition for the cell iterator type. + */ + using type = typename IteratorRange< + CellIteratorType>::IteratorOverIterators::BaseIterator; + + static_assert(std::is_same::value, + "Expected another type as the base iterator."); + }; + + /** + * A helper class to provide a type definition for the underlying cell + * iterator type. + * + * This specialization is for FilteredIterator, which may have either a + * TriaActiveIterator as its base type, or may be nested with another + * FilteredIterator as the type to iterate over. + */ + template + struct CellIteratorBaseType> + { + /** + * Type definition for the cell iterator type. + */ + // Since we can have nested filtered iterators, we recursively + // remove the template layers to retrieve the underlying iterator type. + using type = typename CellIteratorBaseType::type; + }; + + /** + * A helper class to provide a type definition for the underlying cell + * iterator type. + * + * This specialization is for an IteratorRange that iterates over + * FilteredIterators. + */ + template + struct CellIteratorBaseType>::IteratorOverIterators> + { + /** + * Type definition for the cell iterator type. + */ + // Since we can have nested filtered iterators, we recursively + // remove the template layers to retrieve the underlying iterator type. + using type = + typename CellIteratorBaseType>::type; + }; + } // namespace internal + /** * This function extends the WorkStream concept to work on meshes * (cells and/or faces) and handles the complicated logic for @@ -115,50 +192,55 @@ namespace MeshWorker * @ingroup MeshWorker * @author Luca Heltai and Timo Heister, 2017 */ - template + template ::type> void - mesh_loop(const CellIteratorType & begin, - const typename identity::type &end, + mesh_loop( + const CellIteratorType & begin, + const typename identity::type &end, - const typename identity>::type - &cell_worker, - const typename identity>::type - &copier, - - const ScratchData &sample_scratch_data, - const CopyData & sample_copy_data, - - const AssembleFlags flags = assemble_own_cells, - - const typename identity>::type - &boundary_worker = std::function(), - - const typename identity>::type - &face_worker = std::function(), + const typename identity>::type + &cell_worker, + const typename identity>::type + &copier, - const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), - const unsigned int chunk_size = 8) + const ScratchData &sample_scratch_data, + const CopyData & sample_copy_data, + + const AssembleFlags flags = assemble_own_cells, + + const typename identity>::type + &boundary_worker = std::function(), + + const typename identity>::type + &face_worker = std::function(), + + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) { Assert( (!cell_worker) == !(flags & work_on_cells), @@ -193,9 +275,9 @@ namespace MeshWorker ExcMessage( "If you specify a boundary_worker, assemble_boundary_faces needs to be set.")); - auto cell_action = [&](const CellIteratorType &cell, - ScratchData & scratch, - CopyData & copy) { + auto cell_action = [&](const CellIteratorBaseType &cell, + ScratchData & scratch, + CopyData & copy) { // First reset the CopyData class to the empty copy_data given by the // user. copy = sample_copy_data; @@ -224,8 +306,8 @@ namespace MeshWorker if (flags & (work_on_faces | work_on_boundary)) for (unsigned int face_no = 0; - face_no < GeometryInfo::faces_per_cell; + face_no < GeometryInfo::faces_per_cell; ++face_no) { if (cell->at_boundary(face_no) && @@ -238,8 +320,8 @@ namespace MeshWorker else { // interior face, potentially assemble - TriaIterator neighbor = - cell->neighbor_or_periodic_neighbor(face_no); + TriaIterator + neighbor = cell->neighbor_or_periodic_neighbor(face_no); types::subdomain_id neighbor_subdomain_id = numbers::artificial_subdomain_id; @@ -325,7 +407,7 @@ namespace MeshWorker { // If iterator is active and neighbor is refined, skip // internal face. - if (internal::is_active_iterator(cell) && + if (dealii::internal::is_active_iterator(cell) && neighbor->has_children()) continue; @@ -393,6 +475,75 @@ namespace MeshWorker chunk_size); } + /** + * Same as the function above, but for iterator ranges (and, therefore, + * filtered iterators). + */ + template ::type> + void + mesh_loop( + IteratorRange iterator_range, + const typename identity>::type + &cell_worker, + const typename identity>::type + &copier, + + const ScratchData &sample_scratch_data, + const CopyData & sample_copy_data, + + const AssembleFlags flags = assemble_own_cells, + + const typename identity>::type + &boundary_worker = std::function(), + + const typename identity>::type + &face_worker = std::function(), + + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + // Call the function above + mesh_loop::IteratorOverIterators, + ScratchData, + CopyData, + CellIteratorBaseType>(iterator_range.begin(), + iterator_range.end(), + cell_worker, + copier, + sample_scratch_data, + sample_copy_data, + flags, + boundary_worker, + face_worker, + queue_length, + chunk_size); + } + /** * This is a variant of the mesh_loop() function, that can be used for worker * and copier functions that are member functions of a class. @@ -539,6 +690,60 @@ namespace MeshWorker queue_length, chunk_size); } + + /** + * Same as the function above, but for iterator ranges (and, therefore, + * filtered iterators). + */ + template ::type> + void + mesh_loop(IteratorRange iterator_range, + MainClass & main_class, + void (MainClass::*cell_worker)(const CellIteratorBaseType &, + ScratchData &, + CopyData &), + void (MainClass::*copier)(const CopyData &), + const ScratchData & sample_scratch_data, + const CopyData & sample_copy_data, + const AssembleFlags flags = assemble_own_cells, + void (MainClass::*boundary_worker)(const CellIteratorBaseType &, + const unsigned int, + ScratchData &, + CopyData &) = nullptr, + void (MainClass::*face_worker)(const CellIteratorBaseType &, + const unsigned int, + const unsigned int, + const CellIteratorBaseType &, + const unsigned int, + const unsigned int, + ScratchData &, + CopyData &) = nullptr, + const unsigned int queue_length = 2 * MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + // Call the function above + mesh_loop::IteratorOverIterators, + ScratchData, + CopyData, + MainClass, + CellIteratorBaseType>(iterator_range.begin(), + iterator_range.end(), + main_class, + cell_worker, + copier, + sample_scratch_data, + sample_copy_data, + flags, + boundary_worker, + face_worker, + queue_length, + chunk_size); + } } // namespace MeshWorker DEAL_II_NAMESPACE_CLOSE diff --git a/tests/meshworker/mesh_loop_04.cc b/tests/meshworker/mesh_loop_04.cc new file mode 100644 index 0000000000..7decee699b --- /dev/null +++ b/tests/meshworker/mesh_loop_04.cc @@ -0,0 +1,140 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2019 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 mesh_loop for iterator ranges +// This test is based off meshworker/scratch_data_01.cc, but for iterator ranges + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include "../tests.h" + +using namespace MeshWorker; + +template +void +test() +{ + Triangulation tria; + FE_Q fe(1); + DoFHandler dh(tria); + GridGenerator::hyper_cube(tria); + tria.refine_global(1); + tria.begin_active()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + dh.distribute_dofs(fe); + + QGauss quad(3); + QGauss face_quad(3); + + UpdateFlags cell_flags = update_JxW_values; + UpdateFlags face_flags = update_JxW_values; + + using ScratchData = MeshWorker::ScratchData; + using CopyData = MeshWorker::CopyData<0, 1, 0>; + + ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags); + CopyData copy(3); // store measures + + using Iterator = decltype(dh.begin_active()); + + auto measures = std::make_tuple(0.0, 0.0, 0.0); + + + auto cell_worker = [](const Iterator &cell, ScratchData &s, CopyData &c) { + const auto &fev = s.reinit(cell); + const auto &JxW = s.get_JxW_values(); + c = 0; + for (auto w : JxW) + c.vectors[0][0] += w; + }; + + auto boundary_worker = [](const Iterator & cell, + const unsigned int &f, + ScratchData & s, + CopyData & c) { + const auto &fev = s.reinit(cell, f); + const auto &JxW = s.get_JxW_values(); + for (auto w : JxW) + c.vectors[0][1] += w; + }; + + auto face_worker = [](const Iterator & cell, + const unsigned int &f, + const unsigned int &sf, + const Iterator & ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData & s, + CopyData & c) { + const auto &fev = s.reinit(cell, f, sf); + const auto &nfev = s.reinit_neighbor(ncell, nf, nsf); + + const auto &JxW = s.get_JxW_values(); + const auto &nJxW = s.get_neighbor_JxW_values(); + for (auto w : JxW) + c.vectors[0][2] += w; + }; + + auto copier = [&measures](const CopyData &c) { + std::get<0>(measures) += c.vectors[0][0]; + std::get<1>(measures) += c.vectors[0][1]; + std::get<2>(measures) += c.vectors[0][2]; + }; + + mesh_loop(dh.active_cell_iterators(), + cell_worker, + copier, + scratch, + copy, + assemble_own_cells | assemble_boundary_faces | + assemble_own_interior_faces_once, + boundary_worker, + face_worker); + + deallog << "Testing <" << dim << "," << spacedim << ">" << std::endl; + + deallog << "Volume: " << std::get<0>(measures) << std::endl; + + deallog << "Boundary surface: " << std::get<1>(measures) << std::endl; + + deallog << "Interior surface: " << std::get<2>(measures) << std::endl; +} + + + +int +main() +{ + initlog(); + + test<2, 2>(); + test<2, 3>(); + test<3, 3>(); +} diff --git a/tests/meshworker/mesh_loop_04.output b/tests/meshworker/mesh_loop_04.output new file mode 100644 index 0000000000..b9c1ecd3d3 --- /dev/null +++ b/tests/meshworker/mesh_loop_04.output @@ -0,0 +1,13 @@ + +DEAL::Testing <2,2> +DEAL::Volume: 1.00000 +DEAL::Boundary surface: 4.00000 +DEAL::Interior surface: 3.00000 +DEAL::Testing <2,3> +DEAL::Volume: 1.00000 +DEAL::Boundary surface: 4.00000 +DEAL::Interior surface: 3.00000 +DEAL::Testing <3,3> +DEAL::Volume: 1.00000 +DEAL::Boundary surface: 6.00000 +DEAL::Interior surface: 3.75000 diff --git a/tests/meshworker/mesh_loop_05.cc b/tests/meshworker/mesh_loop_05.cc new file mode 100644 index 0000000000..d251046114 --- /dev/null +++ b/tests/meshworker/mesh_loop_05.cc @@ -0,0 +1,157 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2019 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 mesh_loop for filtered iterator ranges +// This test is based off base/work_stream_06.cc, but for parallel distributed +// Triangulations. + +#include + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include "../tests.h" + +using namespace MeshWorker; + +template +void +test() +{ + parallel::distributed::Triangulation tria( + MPI_COMM_WORLD, + Triangulation::limit_level_difference_at_vertices); + FE_Q fe(1); + DoFHandler dh(tria); + GridGenerator::hyper_cube(tria); + tria.refine_global(1); + tria.begin_active()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + dh.distribute_dofs(fe); + + QGauss quad(3); + QGauss face_quad(3); + + UpdateFlags cell_flags = update_JxW_values; + UpdateFlags face_flags = update_JxW_values; + + using ScratchData = MeshWorker::ScratchData; + using CopyData = MeshWorker::CopyData<0, 1, 0>; + + ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags); + CopyData copy(3); // store measures + + using Iterator = decltype(dh.begin_active()); + + auto measures = std::make_tuple(0.0, 0.0, 0.0); + + + auto cell_worker = [](const Iterator &cell, ScratchData &s, CopyData &c) { + const auto &fev = s.reinit(cell); + const auto &JxW = s.get_JxW_values(); + c = 0; + for (auto w : JxW) + c.vectors[0][0] += w; + }; + + auto boundary_worker = [](const Iterator & cell, + const unsigned int &f, + ScratchData & s, + CopyData & c) { + const auto &fev = s.reinit(cell, f); + const auto &JxW = s.get_JxW_values(); + for (auto w : JxW) + c.vectors[0][1] += w; + }; + + auto face_worker = [](const Iterator & cell, + const unsigned int &f, + const unsigned int &sf, + const Iterator & ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData & s, + CopyData & c) { + const auto &fev = s.reinit(cell, f, sf); + const auto &nfev = s.reinit_neighbor(ncell, nf, nsf); + + const auto &JxW = s.get_JxW_values(); + const auto &nJxW = s.get_neighbor_JxW_values(); + for (auto w : JxW) + c.vectors[0][2] += w; + }; + + auto copier = [&measures](const CopyData &c) { + std::get<0>(measures) += c.vectors[0][0]; + std::get<1>(measures) += c.vectors[0][1]; + std::get<2>(measures) += c.vectors[0][2]; + }; + + const auto filtered_iterator_range = + filter_iterators(dh.active_cell_iterators(), + IteratorFilters::LocallyOwnedCell()); + + mesh_loop(filtered_iterator_range, + cell_worker, + copier, + scratch, + copy, + assemble_own_cells | assemble_boundary_faces | + assemble_own_interior_faces_once, + boundary_worker, + face_worker); + + std::get<0>(measures) = + Utilities::MPI::sum(std::get<0>(measures), MPI_COMM_WORLD); + std::get<1>(measures) = + Utilities::MPI::sum(std::get<1>(measures), MPI_COMM_WORLD); + std::get<2>(measures) = + Utilities::MPI::sum(std::get<2>(measures), MPI_COMM_WORLD); + + deallog << "Testing <" << dim << "," << spacedim << ">" << std::endl; + + deallog << "Volume: " << std::get<0>(measures) << std::endl; + + deallog << "Boundary surface: " << std::get<1>(measures) << std::endl; + + deallog << "Interior surface: " << std::get<2>(measures) << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + mpi_initlog(); + + test<2, 2>(); + test<2, 3>(); + test<3, 3>(); +} diff --git a/tests/meshworker/mesh_loop_05.with_p4est=true.mpirun=2.output b/tests/meshworker/mesh_loop_05.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..7c2f39fa9c --- /dev/null +++ b/tests/meshworker/mesh_loop_05.with_p4est=true.mpirun=2.output @@ -0,0 +1,13 @@ + +DEAL::Testing <2,2> +DEAL::Volume: 1.00000 +DEAL::Boundary surface: 4.00000 +DEAL::Interior surface: 2.00000 +DEAL::Testing <2,3> +DEAL::Volume: 1.00000 +DEAL::Boundary surface: 4.00000 +DEAL::Interior surface: 2.00000 +DEAL::Testing <3,3> +DEAL::Volume: 1.00000 +DEAL::Boundary surface: 6.00000 +DEAL::Interior surface: 3.00000 diff --git a/tests/meshworker/mesh_loop_06.cc b/tests/meshworker/mesh_loop_06.cc new file mode 100644 index 0000000000..3fb235e5a3 --- /dev/null +++ b/tests/meshworker/mesh_loop_06.cc @@ -0,0 +1,157 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2019 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 mesh_loop for filtered iterator ranges +// This test is based off base/work_stream_06.cc, but for filtered iterator +// ranges. In particular, we test the material ID predicate. + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include "../tests.h" + +using namespace MeshWorker; + +template +void +test(const types::material_id test_id) +{ + Assert(test_id == 0 || test_id == 1, + ExcMessage("Only material IDs 0 and 1 are valid.")); + + Triangulation tria; + FE_Q fe(1); + DoFHandler dh(tria); + GridGenerator::hyper_cube(tria); + tria.refine_global(1); + tria.begin_active()->set_material_id(1); + tria.begin_active()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + dh.distribute_dofs(fe); + + QGauss quad(3); + QGauss face_quad(3); + + UpdateFlags cell_flags = update_JxW_values; + UpdateFlags face_flags = update_JxW_values; + + using ScratchData = MeshWorker::ScratchData; + using CopyData = MeshWorker::CopyData<0, 1, 0>; + + ScratchData scratch(fe, quad, cell_flags, face_quad, face_flags); + CopyData copy(3); // store measures + + using Iterator = decltype(dh.begin_active()); + + auto measures = std::make_tuple(0.0, 0.0, 0.0); + + + auto cell_worker = [](const Iterator &cell, ScratchData &s, CopyData &c) { + const auto &fev = s.reinit(cell); + const auto &JxW = s.get_JxW_values(); + c = 0; + for (auto w : JxW) + c.vectors[0][0] += w; + }; + + auto boundary_worker = [](const Iterator & cell, + const unsigned int &f, + ScratchData & s, + CopyData & c) { + const auto &fev = s.reinit(cell, f); + const auto &JxW = s.get_JxW_values(); + for (auto w : JxW) + c.vectors[0][1] += w; + }; + + auto face_worker = [](const Iterator & cell, + const unsigned int &f, + const unsigned int &sf, + const Iterator & ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData & s, + CopyData & c) { + const auto &fev = s.reinit(cell, f, sf); + const auto &nfev = s.reinit_neighbor(ncell, nf, nsf); + + const auto &JxW = s.get_JxW_values(); + const auto &nJxW = s.get_neighbor_JxW_values(); + for (auto w : JxW) + c.vectors[0][2] += w; + }; + + auto copier = [&measures](const CopyData &c) { + std::get<0>(measures) += c.vectors[0][0]; + std::get<1>(measures) += c.vectors[0][1]; + std::get<2>(measures) += c.vectors[0][2]; + }; + + const auto filtered_iterator_range = + filter_iterators(dh.active_cell_iterators(), + IteratorFilters::MaterialIdEqualTo(test_id)); + + mesh_loop(filtered_iterator_range, + cell_worker, + copier, + scratch, + copy, + assemble_own_cells | assemble_boundary_faces | + assemble_own_interior_faces_once, + boundary_worker, + face_worker); + + deallog << "Testing <" << dim << "," << spacedim << ">" << std::endl; + + deallog << "Volume: " << std::get<0>(measures) << std::endl; + + deallog << "Boundary surface: " << std::get<1>(measures) << std::endl; + + deallog << "Interior surface: " << std::get<2>(measures) << std::endl; +} + + + +int +main() +{ + initlog(); + + for (auto mat_id : {0, 1}) + { + deallog.push("Material ID " + Utilities::to_string(mat_id)); + { + test<2, 2>(mat_id); + test<2, 3>(mat_id); + test<3, 3>(mat_id); + } + deallog.pop(); + } +} diff --git a/tests/meshworker/mesh_loop_06.output b/tests/meshworker/mesh_loop_06.output new file mode 100644 index 0000000000..784fe6e462 --- /dev/null +++ b/tests/meshworker/mesh_loop_06.output @@ -0,0 +1,25 @@ + +DEAL:Material ID 0::Testing <2,2> +DEAL:Material ID 0::Volume: 0.750000 +DEAL:Material ID 0::Boundary surface: 3.00000 +DEAL:Material ID 0::Interior surface: 1.00000 +DEAL:Material ID 0::Testing <2,3> +DEAL:Material ID 0::Volume: 0.750000 +DEAL:Material ID 0::Boundary surface: 3.00000 +DEAL:Material ID 0::Interior surface: 1.00000 +DEAL:Material ID 0::Testing <3,3> +DEAL:Material ID 0::Volume: 0.875000 +DEAL:Material ID 0::Boundary surface: 5.25000 +DEAL:Material ID 0::Interior surface: 2.25000 +DEAL:Material ID 1::Testing <2,2> +DEAL:Material ID 1::Volume: 0.250000 +DEAL:Material ID 1::Boundary surface: 1.00000 +DEAL:Material ID 1::Interior surface: 2.00000 +DEAL:Material ID 1::Testing <2,3> +DEAL:Material ID 1::Volume: 0.250000 +DEAL:Material ID 1::Boundary surface: 1.00000 +DEAL:Material ID 1::Interior surface: 2.00000 +DEAL:Material ID 1::Testing <3,3> +DEAL:Material ID 1::Volume: 0.125000 +DEAL:Material ID 1::Boundary surface: 0.750000 +DEAL:Material ID 1::Interior surface: 1.50000