Allow 1D elements to break symmetry of the "is neighbor of" relationship

src/geom/elem.C

@@ -1219,26 +1219,37 @@ void Elem::libmesh_assert_valid_neighbors() const
           else
             {
               unsigned int rev = neigh->which_neighbor_am_i(elem);
-              libmesh_assert_less (rev, neigh->n_neighbors());
 
-              if (this->subactive() && !neigh->subactive())
+              // We relax the valid neighbor requirements for 1D elements
+              // since multiple 1D elements may be connected to the same
+              // node, and in this case, they do not satisfy the
+              // "A is a neighbor of B" -> "B is a neighbor of A"
+              // assumption.
+              libmesh_assert (rev < neigh->n_neighbors() ||
+                              this->dim() == 1);
+
+              if (rev < neigh->n_neighbors())
                 {
-                  while (neigh->neighbor_ptr(rev) != elem)
+                  if (this->subactive() && !neigh->subactive())
                     {
-                      libmesh_assert(elem->parent());
-                      elem = elem->parent();
+                      while (neigh->neighbor_ptr(rev) != elem)
+                        {
+                          libmesh_assert(elem->parent());
+                          elem = elem->parent();
+                        }
                     }
-                }
-              else
-                {
-                  const Elem * nn = neigh->neighbor_ptr(rev);
-                  libmesh_assert(nn);
+                  else
+                    {
+                      const Elem * nn = neigh->neighbor_ptr(rev);
+                      libmesh_assert(nn);
 
-                  for (; elem != nn; elem = elem->parent())
-                    libmesh_assert(elem);
+                      for (; elem != nn; elem = elem->parent())
+                        libmesh_assert(elem);
+                    }
                 }
             }
         }
+
       // If we don't have a neighbor and we're not subactive, our
       // ancestors shouldn't have any neighbors in this same
       // direction.

src/mesh/unstructured_mesh.C

@@ -882,7 +882,9 @@ void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
   // with identical side keys and then check to see if they
   // are neighbors
   {
-    // data structures -- Use the hash_multimap if available
+    // map from side_key -> (Elem, side_index)
+    // We use a multimap since side keys are not guaranteed to be
+    // unique, then we verify sides actually match using Elem::operator==()
     typedef dof_id_type                     key_type;
     typedef std::pair<Elem *, unsigned char> val_type;
     typedef std::unordered_multimap<key_type, val_type> map_type;
@@ -963,7 +965,6 @@ void UnstructuredMesh::find_neighbors (const bool reset_remote_elements,
                               {
                                 neighbor->set_neighbor(ns,element);
                               }
-                            side_to_elem_map.erase (bounds.first);
 
                             // get out of this nested crap
                             goto next_side;