Expose mesh loop bug for anisotropic grids.

diff --git a/doc/news/changes/minor/20201103LucaHeltai b/doc/news/changes/minor/20201103LucaHeltai
new file mode 100644 (file)
index 0000000..05e46af
--- /dev/null
+++ b/doc/news/changes/minor/20201103LucaHeltai
@@ -0,0 +1,3 @@
+Fixed: WorkStream::mesh_loop() should now work on anisotropic grids.
+<br>
+(Luca Heltai, 2020/11/03)

diff --git a/include/deal.II/meshworker/mesh_loop.h b/include/deal.II/meshworker/mesh_loop.h
index 3c9d143a7e6f7f40fb1baa6378dba8d64518b2d3..69fd610285779cfb4a1da3fee013298313c5c4b9 100644 (file)
--- 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 (file)
index 0000000..168d553
--- /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 <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/meshworker/mesh_loop.h>
+
+#include <fstream>
+
+#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 <int dim>
+void
+test()
+{
+  deallog << "Testing dim = " << dim << std::endl
+          << "===============" << std::endl
+          << std::endl;
+  Triangulation<dim> tria;
+  create_grid(tria);
+
+  struct ScratchData
+  {
+    unsigned int foo;
+  };
+  struct CopyData
+  {
+    unsigned int foo;
+  };
+
+  ScratchData scratch;
+  CopyData    copy;
+
+  FE_Q<dim>       fe(1);
+  DoFHandler<dim> dh(tria);
+  dh.distribute_dofs(fe);
+
+  auto workerFoo = [](const typename DoFHandler<dim>::active_cell_iterator &,
+                      ScratchData &,
+                      CopyData &) {};
+  auto copierFoo = [](const CopyData &) {};
+  auto faceFoo = [](const typename DoFHandler<dim>::active_cell_iterator &cell,
+                    const unsigned int                                    f,
+                    const unsigned int                                    sf,
+                    const typename DoFHandler<dim>::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<dim>::active_cell_iterator &,
+                        const unsigned int,
+                        ScratchData &,
+                        CopyData &) {};
+  std::function<void(const typename DoFHandler<dim>::active_cell_iterator &,
+                     ScratchData &,
+                     CopyData &)>
+    emptyCellWorker(workerFoo);
+  std::function<void(const typename DoFHandler<dim>::active_cell_iterator &,
+                     const unsigned int,
+                     const unsigned int,
+                     const typename DoFHandler<dim>::active_cell_iterator &,
+                     const unsigned int,
+                     const unsigned int,
+                     ScratchData &,
+                     CopyData &)>
+    emptyFaceWorker(faceFoo);
+  std::function<void(const typename DoFHandler<dim>::active_cell_iterator &,
+                     const unsigned int,
+                     ScratchData &,
+                     CopyData &)>
+                                        emptyBoundaryWorker(boundaryFoo);
+  std::function<void(const CopyData &)> 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 (file)
index 0000000..24a0ae0
--- /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::