Overwrite LocationToLevelSet values with unassigned in MeshClassifier::reclassify()

Fixes a bug where the faces can get the wrong LocationToLevelSet values
if reclassify() is called multiple times and the level set function has
been updated in between.

diff --git a/source/non_matching/mesh_classifier.cc b/source/non_matching/mesh_classifier.cc
index 3d919a69d5bd8163601b96321b4259f8ceef77c5..c0c92fa64d252a58102d432ffef02745a55d1514 100644 (file)
--- a/source/non_matching/mesh_classifier.cc
+++ b/source/non_matching/mesh_classifier.cc
@@ -349,9 +349,9 @@ namespace NonMatching
   MeshClassifier<dim>::reclassify()
   {
     initialize();
-    cell_locations.resize(triangulation->n_active_cells(),
+    cell_locations.assign(triangulation->n_active_cells(),
                           LocationToLevelSet::unassigned);
-    face_locations.resize(triangulation->n_raw_faces(),
+    face_locations.assign(triangulation->n_raw_faces(),
                           LocationToLevelSet::unassigned);
 
     // Loop over all cells and determine the location of all non artificial

diff --git a/tests/non_matching/mesh_classifier.cc b/tests/non_matching/mesh_classifier.cc
index 5188877a6b2384c15db7e1c03aefa59dfa67d75b..3330a08ce2686251086ed6b7cbfaaca003c35195 100644 (file)
--- a/tests/non_matching/mesh_classifier.cc
+++ b/tests/non_matching/mesh_classifier.cc
@@ -315,6 +315,48 @@ test_lagrange_coefficents_positive()
 
 
 
+// Check that the values of LocationToLevelSet for the cells and faces get
+// updated correctly when calling reclassify() multiple times.
+//
+// First, make the level set function all negative, call reclassify(), and check
+// that the values of LocationToLevelSet for all cells and faces equals
+// LocationToLevelSet::inside. Then, change the level set function to all
+// positive, call reclassify() again, and check that all values have been
+// changed to LocationToLevelSet::outside.
+template <int dim>
+void
+test_reclassify_called_multiple_times()
+{
+  deallog << "test_reclassify_called_multiple_times" << std::endl;
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation);
+
+  const FE_Q<dim> element(1);
+
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(element);
+
+  Vector<double>                   level_set(element.dofs_per_cell);
+  NonMatching::MeshClassifier<dim> classifier(dof_handler, level_set);
+
+  const typename Triangulation<dim>::active_cell_iterator cell =
+    triangulation.begin_active();
+
+  deallog << "Level set negative" << std::endl;
+  level_set = -1;
+  classifier.reclassify();
+  print_cell_and_face_locations(classifier, cell);
+
+  deallog << "Level set positive" << std::endl;
+  level_set = 1;
+  classifier.reclassify();
+  print_cell_and_face_locations(classifier, cell);
+
+  deallog << std::endl;
+}
+
+
+
 template <int dim>
 void
 run_test()
@@ -327,6 +369,8 @@ run_test()
   // This test doesn't make sense in 1D.
   if (dim != 1)
     test_lagrange_coefficents_positive<dim>();
+
+  test_reclassify_called_multiple_times<dim>();
 }
 
 

diff --git a/tests/non_matching/mesh_classifier.output b/tests/non_matching/mesh_classifier.output
index a29c6fa6bec4fe4f91c72248d40196942c0d72bd..338309e69ea16b400631c75f94d7aff4339e1a70 100644 (file)
--- a/tests/non_matching/mesh_classifier.output
+++ b/tests/non_matching/mesh_classifier.output
@@ -36,6 +36,16 @@ DEAL::cell intersected
 DEAL::face 0 inside
 DEAL::face 1 outside
 DEAL::
+DEAL::test_reclassify_called_multiple_times
+DEAL::Level set negative
+DEAL::cell inside
+DEAL::face 0 inside
+DEAL::face 1 inside
+DEAL::Level set positive
+DEAL::cell outside
+DEAL::face 0 outside
+DEAL::face 1 outside
+DEAL::
 DEAL::dim = 2
 DEAL::test_negative_function
 DEAL::
@@ -92,6 +102,20 @@ DEAL::face 1 outside
 DEAL::face 2 intersected
 DEAL::face 3 outside
 DEAL::
+DEAL::test_reclassify_called_multiple_times
+DEAL::Level set negative
+DEAL::cell inside
+DEAL::face 0 inside
+DEAL::face 1 inside
+DEAL::face 2 inside
+DEAL::face 3 inside
+DEAL::Level set positive
+DEAL::cell outside
+DEAL::face 0 outside
+DEAL::face 1 outside
+DEAL::face 2 outside
+DEAL::face 3 outside
+DEAL::
 DEAL::dim = 3
 DEAL::test_negative_function
 DEAL::
@@ -162,3 +186,21 @@ DEAL::face 3 outside
 DEAL::face 4 intersected
 DEAL::face 5 outside
 DEAL::
+DEAL::test_reclassify_called_multiple_times
+DEAL::Level set negative
+DEAL::cell inside
+DEAL::face 0 inside
+DEAL::face 1 inside
+DEAL::face 2 inside
+DEAL::face 3 inside
+DEAL::face 4 inside
+DEAL::face 5 inside
+DEAL::Level set positive
+DEAL::cell outside
+DEAL::face 0 outside
+DEAL::face 1 outside
+DEAL::face 2 outside
+DEAL::face 3 outside
+DEAL::face 4 outside
+DEAL::face 5 outside
+DEAL::