From 039baa8e514a73147704aa1f3929afc486c0a6b3 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Tue, 3 Nov 2020 10:33:55 +0100 Subject: [PATCH] Expose mesh loop bug for anisotropic grids. --- doc/news/changes/minor/20201103LucaHeltai | 3 + include/deal.II/meshworker/mesh_loop.h | 12 +- tests/meshworker/mesh_loop_anistropic_01.cc | 184 ++++++++++++++++++ .../meshworker/mesh_loop_anistropic_01.output | 61 ++++++ 4 files changed, 255 insertions(+), 5 deletions(-) create mode 100644 doc/news/changes/minor/20201103LucaHeltai create mode 100644 tests/meshworker/mesh_loop_anistropic_01.cc create mode 100644 tests/meshworker/mesh_loop_anistropic_01.output diff --git a/doc/news/changes/minor/20201103LucaHeltai b/doc/news/changes/minor/20201103LucaHeltai new file mode 100644 index 0000000000..05e46affe7 --- /dev/null +++ b/doc/news/changes/minor/20201103LucaHeltai @@ -0,0 +1,3 @@ +Fixed: WorkStream::mesh_loop() should now work on anisotropic grids. +
+(Luca Heltai, 2020/11/03) diff --git a/include/deal.II/meshworker/mesh_loop.h b/include/deal.II/meshworker/mesh_loop.h index 3c9d143a7e..69fd610285 100644 --- a/include/deal.II/meshworker/mesh_loop.h +++ b/include/deal.II/meshworker/mesh_loop.h @@ -394,12 +394,13 @@ namespace MeshWorker cell->has_periodic_neighbor(face_no); if (dim > 1 && ((!periodic_neighbor && - cell->neighbor_is_coarser(face_no)) || + cell->neighbor_is_coarser(face_no) && + neighbor->is_active()) || (periodic_neighbor && - cell->periodic_neighbor_is_coarser(face_no)))) + cell->periodic_neighbor_is_coarser(face_no) && + neighbor->is_active()))) { Assert(cell->is_active(), ExcInternalError()); - Assert(neighbor->is_active(), ExcInternalError()); // skip if only one processor needs to assemble the face // to a ghost cell and the fine cell is not ours. @@ -481,8 +482,9 @@ namespace MeshWorker neighbor->has_children()) continue; - // Now neighbor is on same level, double-check this: - Assert(cell->level() == neighbor->level(), + // Now neighbor is on the same refinement level. + // Double check. + Assert(!cell->neighbor_is_coarser(face_no), ExcInternalError()); // If we own both cells only do faces from one side (unless diff --git a/tests/meshworker/mesh_loop_anistropic_01.cc b/tests/meshworker/mesh_loop_anistropic_01.cc new file mode 100644 index 0000000000..168d55322b --- /dev/null +++ b/tests/meshworker/mesh_loop_anistropic_01.cc @@ -0,0 +1,184 @@ +// --------------------------------------------------------------------- +// +// 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 mesh_loop with anisotropic grids in 2 and 3 dimensions + +#include + +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + + + +// Create 2d grid +void create_grid(Triangulation<2> &tria) +{ + GridGenerator::hyper_cube(tria); + + tria.begin_active()->set_refine_flag(RefinementCase<2>::cut_x); + tria.execute_coarsening_and_refinement(); + tria.begin_active()->set_refine_flag(RefinementCase<2>::cut_x); + tria.execute_coarsening_and_refinement(); + + { + auto it = tria.begin(); + std::advance(it, 4); + it->set_refine_flag(RefinementCase<2>::cut_x); + tria.execute_coarsening_and_refinement(); + } + + { + auto it = tria.begin(); + std::advance(it, 6); + it->set_refine_flag(RefinementCase<2>::cut_y); + tria.execute_coarsening_and_refinement(); + } + + { + auto it = tria.begin(); + std::advance(it, 8); + it->set_refine_flag(RefinementCase<2>::cut_x); + tria.execute_coarsening_and_refinement(); + } + + { + auto it = tria.begin(); + std::advance(it, 2); + it->set_refine_flag(RefinementCase<2>::cut_y); + tria.execute_coarsening_and_refinement(); + } +} + + +// Create 3d grid +void create_grid(Triangulation<3> &tria) +{ + GridGenerator::hyper_cube<3>(tria); + tria.begin_active()->set_refine_flag(dealii::RefinementCase<3>::cut_x); + tria.execute_coarsening_and_refinement(); + tria.begin_active()->set_refine_flag(dealii::RefinementCase<3>::cut_y); + tria.execute_coarsening_and_refinement(); + { + auto it = tria.begin(); + std::advance(it, 2); + it->set_refine_flag(dealii::RefinementCase<3>::cut_z); + tria.execute_coarsening_and_refinement(); + } +} + + + +template +void +test() +{ + deallog << "Testing dim = " << dim << std::endl + << "===============" << std::endl + << std::endl; + Triangulation tria; + create_grid(tria); + + struct ScratchData + { + unsigned int foo; + }; + struct CopyData + { + unsigned int foo; + }; + + ScratchData scratch; + CopyData copy; + + FE_Q fe(1); + DoFHandler dh(tria); + dh.distribute_dofs(fe); + + auto workerFoo = [](const typename DoFHandler::active_cell_iterator &, + ScratchData &, + CopyData &) {}; + auto copierFoo = [](const CopyData &) {}; + auto faceFoo = [](const typename DoFHandler::active_cell_iterator &cell, + const unsigned int f, + const unsigned int sf, + const typename DoFHandler::active_cell_iterator &ncell, + const unsigned int nf, + const unsigned int nsf, + ScratchData &, + CopyData &) { + deallog << "Cell (+): " << cell << " -- " << f << ", " << (int)sf + << std::endl + << "Cell (-): " << ncell << " -- " << nf << ", " << (int)nsf + << std::endl + << std::endl; + }; + auto boundaryFoo = [](const typename DoFHandler::active_cell_iterator &, + const unsigned int, + ScratchData &, + CopyData &) {}; + std::function::active_cell_iterator &, + ScratchData &, + CopyData &)> + emptyCellWorker(workerFoo); + std::function::active_cell_iterator &, + const unsigned int, + const unsigned int, + const typename DoFHandler::active_cell_iterator &, + const unsigned int, + const unsigned int, + ScratchData &, + CopyData &)> + emptyFaceWorker(faceFoo); + std::function::active_cell_iterator &, + const unsigned int, + ScratchData &, + CopyData &)> + emptyBoundaryWorker(boundaryFoo); + std::function emptyCopier(copierFoo); + + MeshWorker::AssembleFlags flags; + flags = MeshWorker::AssembleFlags::assemble_own_cells | + MeshWorker::AssembleFlags::assemble_own_interior_faces_once | + MeshWorker::AssembleFlags::assemble_boundary_faces; + + MeshWorker::mesh_loop(dh.begin_active(), + dh.end(), + emptyCellWorker, + emptyCopier, + scratch, + copy, + flags, + emptyBoundaryWorker, + emptyFaceWorker); +} + + + +int +main() +{ + initlog(); + MultithreadInfo::set_thread_limit(1); + test<2>(); + test<3>(); +} diff --git a/tests/meshworker/mesh_loop_anistropic_01.output b/tests/meshworker/mesh_loop_anistropic_01.output new file mode 100644 index 0000000000..24a0ae0a8f --- /dev/null +++ b/tests/meshworker/mesh_loop_anistropic_01.output @@ -0,0 +1,61 @@ + +DEAL::Testing dim = 2 +DEAL::=============== +DEAL:: +DEAL::Cell (+): 2.0 -- 1, -1 +DEAL::Cell (-): 3.0 -- 0, -1 +DEAL:: +DEAL::Cell (+): 2.2 -- 0, -1 +DEAL::Cell (-): 4.0 -- 1, -1 +DEAL:: +DEAL::Cell (+): 2.2 -- 3, -1 +DEAL::Cell (-): 2.3 -- 2, -1 +DEAL:: +DEAL::Cell (+): 2.3 -- 0, -1 +DEAL::Cell (-): 5.1 -- 1, -1 +DEAL:: +DEAL::Cell (+): 4.0 -- 0, -1 +DEAL::Cell (-): 3.0 -- 1, 0 +DEAL:: +DEAL::Cell (+): 5.0 -- 0, -1 +DEAL::Cell (-): 3.0 -- 1, 1 +DEAL:: +DEAL::Cell (+): 5.0 -- 1, -1 +DEAL::Cell (-): 5.1 -- 0, -1 +DEAL:: +DEAL::Cell (+): 5.0 -- 2, -1 +DEAL::Cell (-): 4.0 -- 3, 0 +DEAL:: +DEAL::Cell (+): 5.1 -- 2, -1 +DEAL::Cell (-): 4.0 -- 3, 1 +DEAL:: +DEAL::Testing dim = 3 +DEAL::=============== +DEAL:: +DEAL::Cell (+): 2.0 -- 3, -1 +DEAL::Cell (-): 2.1 -- 2, -1 +DEAL:: +DEAL::Cell (+): 2.2 -- 0, -1 +DEAL::Cell (-): 2.0 -- 1, 0 +DEAL:: +DEAL::Cell (+): 2.2 -- 3, -1 +DEAL::Cell (-): 2.3 -- 2, -1 +DEAL:: +DEAL::Cell (+): 2.2 -- 5, -1 +DEAL::Cell (-): 2.4 -- 4, -1 +DEAL:: +DEAL::Cell (+): 2.3 -- 0, -1 +DEAL::Cell (-): 2.1 -- 1, 0 +DEAL:: +DEAL::Cell (+): 2.3 -- 5, -1 +DEAL::Cell (-): 2.5 -- 4, -1 +DEAL:: +DEAL::Cell (+): 2.4 -- 0, -1 +DEAL::Cell (-): 2.0 -- 1, 1 +DEAL:: +DEAL::Cell (+): 2.4 -- 3, -1 +DEAL::Cell (-): 2.5 -- 2, -1 +DEAL:: +DEAL::Cell (+): 2.5 -- 0, -1 +DEAL::Cell (-): 2.1 -- 1, 1 +DEAL:: -- 2.39.5