Add topological neighbor method with return of neighbor side

Author: ttruster

include/geom/elem.h

@@ -335,6 +335,18 @@ class Elem : public ReferenceCountedObject<Elem>,
                                MeshBase & mesh,
                                const PointLocatorBase & point_locator,
                                const PeriodicBoundaries * pb);
+                               
+  /**
+   * \returns A pointer to the \f$ i^{th} \f$ neighbor of this element
+   * for interior elements.  If an element is on a periodic
+   * boundary, it will return a corresponding element on the opposite
+   * side, along with the side neigh_side.
+   */
+  const Elem * topological_neighbor_side(const unsigned int i,
+                                         const MeshBase & mesh,
+                                         const PointLocatorBase & point_locator,
+                                         const PeriodicBoundaries * pb,
+                                         unsigned int * neigh_side) const;
 
   /**
    * \returns \p true if the element \p elem in question is a neighbor or

src/geom/elem.C

@@ -1153,6 +1153,45 @@ const Elem * Elem::topological_neighbor (const unsigned int i,
 }
 
 
+const Elem *
+Elem::topological_neighbor_side(const unsigned int i,
+                                const MeshBase & mesh,
+                                const PointLocatorBase & point_locator,
+                                const PeriodicBoundaries * pb,
+                                unsigned int * neigh_side) const
+{
+  libmesh_assert_less (i, this->n_neighbors());
+
+  const Elem * neighbor_i = this->neighbor_ptr(i);
+  if (neighbor_i != nullptr)
+    return neighbor_i;
+
+  if (pb)
+    {
+      // Since the neighbor is nullptr it must be on a boundary. We need
+      // see if this is a periodic boundary in which case it will have a
+      // topological neighbor
+      std::vector<boundary_id_type> bc_ids;
+      mesh.get_boundary_info().boundary_ids(this, cast_int<unsigned short>(i), bc_ids);
+      for (const auto & id : bc_ids)
+        if (pb->boundary(id))
+          {
+            neighbor_i = pb->neighbor(id, point_locator, this, i, neigh_side);
+
+            // Since coarse elements do not have more refined
+            // neighbors we need to make sure that we don't return one
+            // of these types of neighbors.
+            if (neighbor_i)
+              while (level() < neighbor_i->level())
+                neighbor_i = neighbor_i->parent();
+            return neighbor_i;
+          }
+    }
+
+  return nullptr;
+}
+
+
 bool Elem::has_topological_neighbor (const Elem * elem,
                                      const MeshBase & mesh,
                                      const PointLocatorBase & point_locator,