From a797dda7cf3588f1a3024c4e8d30cc3600a9698c Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 10 Aug 2017 16:12:29 -0600 Subject: [PATCH] First version of mesh_loop. --- include/deal.II/meshworker/mesh_loop.h | 232 +++++++++++++++++++++++++ tests/meshworker/mesh_loop_02.cc | 130 ++++++++++++++ tests/meshworker/mesh_loop_02.output | 162 +++++++++++++++++ 3 files changed, 524 insertions(+) create mode 100644 include/deal.II/meshworker/mesh_loop.h create mode 100644 tests/meshworker/mesh_loop_02.cc create mode 100644 tests/meshworker/mesh_loop_02.output diff --git a/include/deal.II/meshworker/mesh_loop.h b/include/deal.II/meshworker/mesh_loop.h new file mode 100644 index 0000000000..db52f446c0 --- /dev/null +++ b/include/deal.II/meshworker/mesh_loop.h @@ -0,0 +1,232 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +#ifndef dealii__mesh_loop_h +#define dealii__mesh_loop_h + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +template class TriaActiveIterator; + +namespace MeshWorker +{ + /** + * A low level work function of this namespace. + * + * @ingroup MeshWorker + * @author Luca Heltai, 2017 + */ + template + void mesh_loop(const Iterator & begin, + const typename identity::type &end, + +// const std::function &cell_worker, +// const std::function ©er, + CellWorker cell_worker, + Copyer copyer, + + ScratchData &scratch_data, + CopyData ©_data, + + const AssembleFlags flags = assemble_own_cells, + +// const std::function &boundary_worker= +// std::function(), + +// const std::function &face_worker= +// std::function(), + + BoundaryWorker boundary_worker = BoundaryWorker(), + FaceWorker face_worker = FaceWorker(), + + const unsigned int queue_length = 2*MultithreadInfo::n_threads(), + const unsigned int chunk_size = 8) + { + auto cell_action = [&] (const Iterator &cell, ScratchData &scratch, CopyData ©) + { + const bool ignore_subdomain = (cell->get_triangulation().locally_owned_subdomain() + == numbers::invalid_subdomain_id); + + types::subdomain_id csid = /*(cell->is_level_cell()) + ? cell->level_subdomain_id() + : */cell->subdomain_id(); + + const bool own_cell = ignore_subdomain || (csid == cell->get_triangulation().locally_owned_subdomain()); + + if ((!ignore_subdomain) && (csid == numbers::artificial_subdomain_id)) + return; + if( (flags & (assemble_cells_first)) && + ( ((flags & (assemble_own_cells)) && own_cell) + || ( (flags & assemble_ghost_cells) && !own_cell) ) ) + cell_worker(cell, scratch, copy); + + if (flags & assemble_own_faces) + for (unsigned int face_no=0; face_no < GeometryInfo::faces_per_cell; ++face_no) + { + typename Iterator::AccessorType::Container::face_iterator face = cell->face(face_no); + if (cell->at_boundary(face_no) && !cell->has_periodic_neighbor(face_no)) + { + // only integrate boundary faces of own cells + if ( (flags & assemble_boundary_faces) && own_cell) + { + boundary_worker(cell, face_no, scratch, copy); + } + } + else if (flags & assemble_own_interior_faces) + { + // Interior face + TriaIterator neighbor = cell->neighbor_or_periodic_neighbor(face_no); + + types::subdomain_id neighbid = numbers::artificial_subdomain_id; +// if (neighbor->is_level_cell()) +// neighbid = neighbor->level_subdomain_id(); + //subdomain id is only valid for active cells + /*else */if (neighbor->active()) + neighbid = neighbor->subdomain_id(); + + const bool own_neighbor = ignore_subdomain || + (neighbid == cell->get_triangulation().locally_owned_subdomain()); + + // skip all faces between two ghost cells + if (!own_cell && !own_neighbor) + continue; + + // skip if the user doesn't want faces between own cells + if (own_cell && own_neighbor && !(flags & assemble_own_interior_faces)) + continue; + + // skip face to ghost + if (own_cell != own_neighbor && !(flags & assemble_ghost_faces)) + continue; + + // Deal with refinement edges from the refined side. Assuming one-irregular + // meshes, this situation should only occur if both cells are active. + const bool periodic_neighbor = cell->has_periodic_neighbor(face_no); + + if ((!periodic_neighbor && cell->neighbor_is_coarser(face_no)) + || (periodic_neighbor && cell->periodic_neighbor_is_coarser(face_no))) + { + Assert(!cell->has_children(), ExcInternalError()); + Assert(!neighbor->has_children(), ExcInternalError()); + + // skip if only one processor needs to assemble the face + // to a ghost cell and the fine cell is not ours. + if (!own_cell && (flags & assemble_ghost_faces_once)) + continue; + + const std::pair neighbor_face_no + = periodic_neighbor? + cell->periodic_neighbor_of_coarser_periodic_neighbor(face_no): + cell->neighbor_of_coarser_neighbor(face_no); + const typename Iterator::AccessorType::Container::face_iterator nface + = neighbor->face(neighbor_face_no.first); + + face_worker(cell, face_no, numbers::invalid_unsigned_int, + neighbor, neighbor_face_no.first, neighbor_face_no.second, + scratch, copy); + + if(flags & assemble_own_interior_faces_both) { + face_worker(neighbor, neighbor_face_no.first, neighbor_face_no.second, + cell, face_no, numbers::invalid_unsigned_int, + scratch, copy); + } + } + else + { + // If iterator is active and neighbor is refined, skip + // internal face. + if (internal::is_active_iterator(cell) && neighbor->has_children()) + { + continue; + } + + // Now neighbor is on same level, double-check this: + Assert(cell->level()==neighbor->level(), ExcInternalError()); + + // If we own both cells only do faces from one side (unless + // AssembleFlags says otherwise). Here, we rely on cell comparison + // that will look at cell->index(). + if (own_cell && own_neighbor + && (flags & assemble_own_interior_faces_once) + && (neighbor < cell)) + continue; + + // independent of loop_control.faces_to_ghost, + // we only look at faces to ghost on the same level once + // (only where own_cell=true and own_neighbor=false) + if (!own_cell) + continue; + + // now only one processor assembles faces_to_ghost. We let the + // processor with the smaller (level-)subdomain id assemble the + // face. + if (own_cell && !own_neighbor + && (flags & assemble_ghost_faces_once) + && (neighbid < csid)) + continue; + + const unsigned int neighbor_face_no = periodic_neighbor? + cell->periodic_neighbor_face_no(face_no): + cell->neighbor_face_no(face_no); + Assert (periodic_neighbor || neighbor->face(neighbor_face_no) == face, ExcInternalError()); + + face_worker(cell, face_no, numbers::invalid_unsigned_int, + neighbor, neighbor_face_no, numbers::invalid_unsigned_int, + scratch, copy); + + } + } + } // faces + + // Execute this, if faces + // have to be handled first + if ((flags & assemble_own_cells) && !(flags & assemble_cells_first) && + ( ((flags & assemble_own_cells) && own_cell) || ((flags & assemble_ghost_cells) && !own_cell))) + cell_worker(cell, scratch, copy); + }; + + + // Loop over all cells + WorkStream::run(begin, end, + cell_action, copyer, + scratch_data, copy_data, + queue_length, chunk_size); + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/tests/meshworker/mesh_loop_02.cc b/tests/meshworker/mesh_loop_02.cc new file mode 100644 index 0000000000..985be0a1a7 --- /dev/null +++ b/tests/meshworker/mesh_loop_02.cc @@ -0,0 +1,130 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2007 - 2016 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// test assemble_flags.h + +#include "../tests.h" +#include +#include +#include +#include +#include + +#include + +struct ScratchData {}; +struct CopyData {}; + +using namespace MeshWorker; + +template +void test() { + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(1); + tria.begin_active()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + + ScratchData scratch; + CopyData copy; + + auto cell=tria.begin_active(); + auto endc=tria.end(); + + typedef decltype(cell) Iterator; + + + + auto cell_worker = [] (const Iterator &cell, ScratchData &s, CopyData &c) { + deallog << "Cell worker on : " << cell << std::endl; + }; + + auto boundary_worker = [] (const Iterator &cell, const unsigned int &f, ScratchData &, CopyData &) { + deallog << "Boundary worker on : " << cell << ", Face : "<< f << std::endl; + }; + + 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) { + deallog << "Face worker on : " << cell << ", Neighbor cell : " << ncell + << ", Face : "<< f << ", Neighbor Face : " << nf + << ", Subface: " << sf + << ", Neighbor Subface: " << nsf << std::endl; + }; + + auto copyer = [](const CopyData &) { + deallog << "Copyer" << std::endl; + }; + + deallog << "CELLS ONLY" << std::endl << std::endl; + + mesh_loop(cell, endc, cell_worker, copyer, scratch, copy, + assemble_own_cells, + boundary_worker, face_worker); + + + deallog << "BOUNDARY ONLY" << std::endl << std::endl; + + mesh_loop(cell, endc, cell_worker, copyer, scratch, copy, + assemble_boundary_faces, + boundary_worker, face_worker); + + deallog << "CELLS AND BOUNDARY" << std::endl << std::endl; + + mesh_loop(cell, endc, cell_worker, copyer, scratch, copy, + assemble_own_cells | assemble_boundary_faces, + boundary_worker, face_worker); + + deallog << "CELLS FIRST AND BOUNDARY" << std::endl << std::endl; + + mesh_loop(cell, endc, cell_worker, copyer, scratch, copy, + assemble_own_cells | assemble_cells_first | assemble_boundary_faces, + boundary_worker, face_worker); + + + deallog << "ONLY FACES ONCE" << std::endl << std::endl; + + mesh_loop(cell, endc, cell_worker, copyer, scratch, copy, + assemble_own_interior_faces_once, + boundary_worker, face_worker); + + + deallog << "ONLY FACES BOTH" << std::endl << std::endl; + + mesh_loop(cell, endc, cell_worker, copyer, scratch, copy, + assemble_own_interior_faces_both, + boundary_worker, face_worker); + + + + + deallog << "CELLS FIRST AND BOUNDARY" << std::endl << std::endl; + + mesh_loop(cell, endc, cell_worker, copyer, scratch, copy, + assemble_own_cells | assemble_cells_first | assemble_boundary_faces, + boundary_worker, face_worker); + +} + + +int main() +{ + initlog(); + + test<2,2>(); +} diff --git a/tests/meshworker/mesh_loop_02.output b/tests/meshworker/mesh_loop_02.output new file mode 100644 index 0000000000..f63bd6b747 --- /dev/null +++ b/tests/meshworker/mesh_loop_02.output @@ -0,0 +1,162 @@ + +DEAL::CELLS ONLY +DEAL:: +DEAL::Cell worker on : 1.1 +DEAL::Cell worker on : 1.2 +DEAL::Cell worker on : 1.3 +DEAL::Cell worker on : 2.0 +DEAL::Cell worker on : 2.1 +DEAL::Cell worker on : 2.2 +DEAL::Cell worker on : 2.3 +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::BOUNDARY ONLY +DEAL:: +DEAL::Boundary worker on : 1.1, Face : 1 +DEAL::Boundary worker on : 1.1, Face : 2 +DEAL::Boundary worker on : 1.2, Face : 0 +DEAL::Boundary worker on : 1.2, Face : 3 +DEAL::Boundary worker on : 1.3, Face : 1 +DEAL::Boundary worker on : 1.3, Face : 3 +DEAL::Boundary worker on : 2.0, Face : 0 +DEAL::Boundary worker on : 2.0, Face : 2 +DEAL::Boundary worker on : 2.1, Face : 2 +DEAL::Boundary worker on : 2.2, Face : 0 +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::CELLS AND BOUNDARY +DEAL:: +DEAL::Boundary worker on : 1.1, Face : 1 +DEAL::Boundary worker on : 1.1, Face : 2 +DEAL::Cell worker on : 1.1 +DEAL::Boundary worker on : 1.2, Face : 0 +DEAL::Boundary worker on : 1.2, Face : 3 +DEAL::Cell worker on : 1.2 +DEAL::Boundary worker on : 1.3, Face : 1 +DEAL::Boundary worker on : 1.3, Face : 3 +DEAL::Cell worker on : 1.3 +DEAL::Boundary worker on : 2.0, Face : 0 +DEAL::Boundary worker on : 2.0, Face : 2 +DEAL::Cell worker on : 2.0 +DEAL::Boundary worker on : 2.1, Face : 2 +DEAL::Cell worker on : 2.1 +DEAL::Boundary worker on : 2.2, Face : 0 +DEAL::Cell worker on : 2.2 +DEAL::Cell worker on : 2.3 +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::CELLS FIRST AND BOUNDARY +DEAL:: +DEAL::Cell worker on : 1.1 +DEAL::Boundary worker on : 1.1, Face : 1 +DEAL::Boundary worker on : 1.1, Face : 2 +DEAL::Cell worker on : 1.2 +DEAL::Boundary worker on : 1.2, Face : 0 +DEAL::Boundary worker on : 1.2, Face : 3 +DEAL::Cell worker on : 1.3 +DEAL::Boundary worker on : 1.3, Face : 1 +DEAL::Boundary worker on : 1.3, Face : 3 +DEAL::Cell worker on : 2.0 +DEAL::Boundary worker on : 2.0, Face : 0 +DEAL::Boundary worker on : 2.0, Face : 2 +DEAL::Cell worker on : 2.1 +DEAL::Boundary worker on : 2.1, Face : 2 +DEAL::Cell worker on : 2.2 +DEAL::Boundary worker on : 2.2, Face : 0 +DEAL::Cell worker on : 2.3 +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::ONLY FACES ONCE +DEAL:: +DEAL::Face worker on : 1.1, Neighbor cell : 1.3, Face : 3, Neighbor Face : 2, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 1.2, Neighbor cell : 1.3, Face : 1, Neighbor Face : 0, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.0, Neighbor cell : 2.1, Face : 1, Neighbor Face : 0, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.0, Neighbor cell : 2.2, Face : 3, Neighbor Face : 2, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.1, Neighbor cell : 1.1, Face : 1, Neighbor Face : 0, Subface: 4294967295, Neighbor Subface: 0 +DEAL::Face worker on : 2.1, Neighbor cell : 2.3, Face : 3, Neighbor Face : 2, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.2, Neighbor cell : 2.3, Face : 1, Neighbor Face : 0, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.2, Neighbor cell : 1.2, Face : 3, Neighbor Face : 2, Subface: 4294967295, Neighbor Subface: 0 +DEAL::Face worker on : 2.3, Neighbor cell : 1.1, Face : 1, Neighbor Face : 0, Subface: 4294967295, Neighbor Subface: 1 +DEAL::Face worker on : 2.3, Neighbor cell : 1.2, Face : 3, Neighbor Face : 2, Subface: 4294967295, Neighbor Subface: 1 +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::ONLY FACES BOTH +DEAL:: +DEAL::Face worker on : 1.1, Neighbor cell : 1.3, Face : 3, Neighbor Face : 2, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 1.2, Neighbor cell : 1.3, Face : 1, Neighbor Face : 0, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 1.3, Neighbor cell : 1.2, Face : 0, Neighbor Face : 1, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 1.3, Neighbor cell : 1.1, Face : 2, Neighbor Face : 3, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.0, Neighbor cell : 2.1, Face : 1, Neighbor Face : 0, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.0, Neighbor cell : 2.2, Face : 3, Neighbor Face : 2, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.1, Neighbor cell : 2.0, Face : 0, Neighbor Face : 1, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.1, Neighbor cell : 1.1, Face : 1, Neighbor Face : 0, Subface: 4294967295, Neighbor Subface: 0 +DEAL::Face worker on : 1.1, Neighbor cell : 2.1, Face : 0, Neighbor Face : 1, Subface: 0, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.1, Neighbor cell : 2.3, Face : 3, Neighbor Face : 2, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.2, Neighbor cell : 2.3, Face : 1, Neighbor Face : 0, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.2, Neighbor cell : 2.0, Face : 2, Neighbor Face : 3, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.2, Neighbor cell : 1.2, Face : 3, Neighbor Face : 2, Subface: 4294967295, Neighbor Subface: 0 +DEAL::Face worker on : 1.2, Neighbor cell : 2.2, Face : 2, Neighbor Face : 3, Subface: 0, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.3, Neighbor cell : 2.2, Face : 0, Neighbor Face : 1, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.3, Neighbor cell : 1.1, Face : 1, Neighbor Face : 0, Subface: 4294967295, Neighbor Subface: 1 +DEAL::Face worker on : 1.1, Neighbor cell : 2.3, Face : 0, Neighbor Face : 1, Subface: 1, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.3, Neighbor cell : 2.1, Face : 2, Neighbor Face : 3, Subface: 4294967295, Neighbor Subface: 4294967295 +DEAL::Face worker on : 2.3, Neighbor cell : 1.2, Face : 3, Neighbor Face : 2, Subface: 4294967295, Neighbor Subface: 1 +DEAL::Face worker on : 1.2, Neighbor cell : 2.3, Face : 2, Neighbor Face : 3, Subface: 1, Neighbor Subface: 4294967295 +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::CELLS FIRST AND BOUNDARY +DEAL:: +DEAL::Cell worker on : 1.1 +DEAL::Boundary worker on : 1.1, Face : 1 +DEAL::Boundary worker on : 1.1, Face : 2 +DEAL::Cell worker on : 1.2 +DEAL::Boundary worker on : 1.2, Face : 0 +DEAL::Boundary worker on : 1.2, Face : 3 +DEAL::Cell worker on : 1.3 +DEAL::Boundary worker on : 1.3, Face : 1 +DEAL::Boundary worker on : 1.3, Face : 3 +DEAL::Cell worker on : 2.0 +DEAL::Boundary worker on : 2.0, Face : 0 +DEAL::Boundary worker on : 2.0, Face : 2 +DEAL::Cell worker on : 2.1 +DEAL::Boundary worker on : 2.1, Face : 2 +DEAL::Cell worker on : 2.2 +DEAL::Boundary worker on : 2.2, Face : 0 +DEAL::Cell worker on : 2.3 +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer +DEAL::Copyer -- 2.39.5