Relax some find_neighbors() assertions that can be overzealous for 1D Elems

src/mesh/unstructured_mesh.C

@@ -1047,17 +1047,34 @@ void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
                       // Let's make sure that "had children made remote"
                       // situation is actually the case
                       libmesh_assert(neigh->has_children());
+
+                      // This assert is overzealous in the case where
+                      // N 1D elements (N>2, even) meet in a common
+                      // node, but an odd number are flagged for
+                      // refinement. In that case, our parent may have
+                      // a coarse level neighbor on this side while we
+                      // do not, because there was no element for us to
+                      // pair up with. Libmesh currently assumes one
+                      // neighbor per side, and "A neighbor-of B" implies
+                      // "B neighbor-of A" so we can't currently handle
+                      // this case correctly. If we are in that case,
+                      // we therefore accept that the neighbor information
+                      // is incomplete/wrong rather than asserting for
+                      // an unrelated reason.
                       bool neigh_has_remote_children = false;
                       for (auto & child : neigh->child_ref_range())
                         if (&child == remote_elem)
                           neigh_has_remote_children = true;
-                      libmesh_assert(neigh_has_remote_children);
+                      libmesh_assert(neigh_has_remote_children ||
+                                     neigh->dim() == 1);
 
-                      // And let's double-check that we don't have
-                      // a remote_elem neighboring an active local element
+                      // And let's double-check that we don't have a
+                      // remote_elem neighboring an active local element
+                      // unless this is a 1D Elem, for the reasons
+                      // discussed above.
                       if (current_elem->active())
-                        libmesh_assert_not_equal_to (current_elem->processor_id(),
-                                                     this->processor_id());
+                        libmesh_assert (current_elem->dim() == 1 ||
+                                        current_elem->processor_id() != this->processor_id());
 #endif // DEBUG
                       neigh = const_cast<RemoteElem *>(remote_elem);
                     }

tests/Makefile.am

@@ -160,7 +160,8 @@ data = matrices/geom_1_extraction_op.m \
        meshes/mesh_with_low_order_edge_blocks.e \
        meshes/shark_tooth_tri6.xda.gz \
        meshes/tetgen_one_tet10.ele \
-       meshes/tetgen_one_tet10.node
+       meshes/tetgen_one_tet10.node \
+       meshes/find_neighbors_junction_1D_only.e
 
 unit_tests_data = $(data)
 

tests/Makefile.in

@@ -2242,7 +2242,8 @@ data = matrices/geom_1_extraction_op.m \
        meshes/mesh_with_low_order_edge_blocks.e \
        meshes/shark_tooth_tri6.xda.gz \
        meshes/tetgen_one_tet10.ele \
-       meshes/tetgen_one_tet10.node
+       meshes/tetgen_one_tet10.node \
+       meshes/find_neighbors_junction_1D_only.e
 
 unit_tests_data = $(data)
 @LIBMESH_DBG_MODE_TRUE@@LIBMESH_ENABLE_CPPUNIT_TRUE@unit_tests_dbg_SOURCES = $(unit_tests_sources)

tests/mesh/nodal_neighbors.C

@@ -4,6 +4,7 @@
 #include <libmesh/mesh_tools.h>
 #include <libmesh/replicated_mesh.h>
 #include <libmesh/elem.h>
+#include <libmesh/mesh_refinement.h>
 
 #include "test_comm.h"
 #include "libmesh_cppunit.h"
@@ -37,6 +38,7 @@ public:
   CPPUNIT_TEST( testEdge3 );
   CPPUNIT_TEST( testEdge4 );
   CPPUNIT_TEST( testOrientation );
+  CPPUNIT_TEST( testMultiNeighborRefinement );
 
   CPPUNIT_TEST_SUITE_END();
 
@@ -242,6 +244,32 @@ public:
         //   }
       }
   }
+
+  void testMultiNeighborRefinement()
+  {
+    LOG_UNIT_TEST;
+
+    // The input mesh has a Node where 6 1D elements all meet.  Then,
+    // we flag three of those elements for refinement, refine the
+    // Mesh, and finally test that we can call MeshBase::prepare_for_use().
+    // Prior to the changes relaxing some assertions on 1D elements in
+    // #3853 (https://github.com/libMesh/libmesh/pull/3853) this test
+    // failed with a segfault.
+    ReplicatedMesh mesh(*TestCommWorld);
+    mesh.read("meshes/find_neighbors_junction_1D_only.e");
+
+    // Mark elements in subdomain 4 for refinement now:
+    for (const auto & elem : mesh.element_ptr_range())
+      if (elem->subdomain_id() == 4)
+        elem->set_refinement_flag(Elem::REFINE);
+
+    // Do the refinement
+    MeshRefinement mesh_refinement(mesh);
+    mesh_refinement.refine_elements();
+
+    // Trigger a new find_neighbors() call
+    mesh.prepare_for_use();
+  }
 };