diff --git a/src/mesh/mesh_tools.C b/src/mesh/mesh_tools.C
index e014c5f7dee..405b5c4506b 100644
--- a/src/mesh/mesh_tools.C
+++ b/src/mesh/mesh_tools.C
@@ -350,6 +350,8 @@ void find_nodal_neighbors_helper(const dof_id_type global_id,
                 }
             }
 
+          const auto elem_order = Elem::type_to_default_order_map[elem->type()];
+
           // Index of the current edge
           unsigned current_edge = 0;
 
@@ -373,18 +375,23 @@ void find_nodal_neighbors_helper(const dof_id_type global_id,
 
                   // Find another node in this element on this edge
                   for (unsigned other_node_this_edge = 0; other_node_this_edge != n_nodes; other_node_this_edge++)
-                    if ( (elem->is_node_on_edge(other_node_this_edge, current_edge)) && // On the current edge
-                         (elem->node_id(other_node_this_edge) != global_id))               // But not the original node
-                      {
-                        // We've found a nodal neighbor!  Save a pointer to it..
-                        node_to_save = elem->node_ptr(other_node_this_edge);
-                        break;
-                      }
-
-                  // Make sure we found something
-                  libmesh_assert(node_to_save != nullptr);
-
-                  neighbor_set.insert(node_to_save);
+                    {
+                      const bool both_vertices = elem->is_vertex(local_node_number) &&
+                                                 elem->is_vertex(other_node_this_edge);
+                      if ( elem->is_node_on_edge(other_node_this_edge, current_edge) && // On the current edge
+                           elem->node_id(other_node_this_edge) != global_id          && // But not the original node
+                            // vertex nodes on the same edge of higher order elements are not nodal neighbors
+                          (elem_order == 1 || !both_vertices))
+                        {
+                          // We've found a nodal neighbor!  Save a pointer to it..
+                          node_to_save = elem->node_ptr(other_node_this_edge);
+
+                          // Make sure we found something
+                          libmesh_assert(node_to_save != nullptr);
+
+                          neighbor_set.insert(node_to_save);
+                        }
+                    }
                 }
 
               // Keep looking for edges, node may be on more than one edge
diff --git a/tests/mesh/nodal_neighbors.C b/tests/mesh/nodal_neighbors.C
index 9ac570a5127..8ece7047bae 100644
--- a/tests/mesh/nodal_neighbors.C
+++ b/tests/mesh/nodal_neighbors.C
@@ -4,6 +4,9 @@
 #include <libmesh/mesh_tools.h>
 #include <libmesh/replicated_mesh.h>
 #include <libmesh/elem.h>
+#include <libmesh/utility.h>
+#include <libmesh/reference_elem.h>
+#include <unordered_map>
 
 #include "test_comm.h"
 #include "libmesh_cppunit.h"
@@ -15,7 +18,7 @@ class NodalNeighborsTest : public CppUnit::TestCase
 {
   /**
    * The goal of this test is to ensure that MeshTools::find_nodal_neighbors()
-   * works in 1D.  If the numbering of MeshGeneration::build_line()
+   * works as expected.  If the numbering of MeshGeneration::build_line()
    * ever changes, this test will break, as it compares hand-checked
    * hard-coded "validation" data with the results of
    * MeshTools::find_nodal_neighbors().  We also use a ReplicatedMesh
@@ -38,53 +41,291 @@ public:
   CPPUNIT_TEST( testEdge4 );
   CPPUNIT_TEST( testOrientation );
 
+  CPPUNIT_TEST( testTri3 );
+  CPPUNIT_TEST( testTri6 );
+  CPPUNIT_TEST( testTri7 );
+
+  CPPUNIT_TEST( testQuad4 );
+  CPPUNIT_TEST( testQuad8 );
+  CPPUNIT_TEST( testQuad9 );
+
+  CPPUNIT_TEST( testTet4 );
+  CPPUNIT_TEST( testTet10 );
+  CPPUNIT_TEST( testTet14 );
+
+  CPPUNIT_TEST( testPyramid5 );
+  CPPUNIT_TEST( testPyramid13 );
+  CPPUNIT_TEST( testPyramid14 );
+  CPPUNIT_TEST( testPyramid18 );
+
+  CPPUNIT_TEST( testPrism6 );
+  CPPUNIT_TEST( testPrism18 );
+  CPPUNIT_TEST( testPrism20 );
+  CPPUNIT_TEST( testPrism21 );
+
+  CPPUNIT_TEST( testHex8 );
+  CPPUNIT_TEST( testHex20 );
+  CPPUNIT_TEST( testHex27 );
+
   CPPUNIT_TEST_SUITE_END();
 
+private:
+
+  // Hard-coded neighbor relationships for each element type
+  using Map = std::map<ElemType, std::map<dof_id_type, std::set<dof_id_type>>>;
+
+  static Map build_elem_type_to_neighbor_map()
+  {
+    Map m;
+    //m.reserve(10);
+
+    m[EDGE2][0].insert(1);
+    m[EDGE2][1].insert(0);
+
+    m[EDGE3][0].insert(2);
+    m[EDGE3][1].insert(2);
+    m[EDGE3][2].insert({0, 1});
+
+    m[EDGE4][0].insert(2);
+    m[EDGE4][1].insert(3);
+    m[EDGE4][2].insert({0, 3});
+    m[EDGE4][3].insert({2, 1});
+
+    m[TRI3][0].insert({1, 2});
+    m[TRI3][1].insert({0, 2});
+    m[TRI3][2].insert({1, 0});
+
+    m[TRI6][0].insert({3, 5});
+    m[TRI6][1].insert({3, 4});
+    m[TRI6][2].insert({5, 4});
+    m[TRI6][3].insert({0, 1});
+    m[TRI6][4].insert({1, 2});
+    m[TRI6][5].insert({2, 0});
+
+    m[TRI7] = m[TRI6];
+    // Node 6 of TRI7 is a face node, has no neighbors because it is not on an edge
+
+    m[QUAD4][0].insert({1, 3});
+    m[QUAD4][1].insert({0, 2});
+    m[QUAD4][2].insert({1, 3});
+    m[QUAD4][3].insert({0, 2});
+
+    m[QUAD8][0].insert({4, 7});
+    m[QUAD8][1].insert({4, 5});
+    m[QUAD8][2].insert({5, 6});
+    m[QUAD8][3].insert({6, 7});
+    m[QUAD8][4].insert({0, 1});
+    m[QUAD8][5].insert({1, 2});
+    m[QUAD8][6].insert({2, 3});
+    m[QUAD8][7].insert({0, 3});
+
+    m[QUAD9] = m[QUAD8];
+    // Node 8 of QUAD9 is a face node, has no neighbors because it is not on an edge
+
+    m[TET4][0].insert({1, 2, 3});
+    m[TET4][1].insert({0, 2, 3});
+    m[TET4][2].insert({0, 1, 3});
+    m[TET4][3].insert({0, 1, 2});
+
+    m[TET10][0].insert({4, 6, 7});
+    m[TET10][1].insert({4, 5, 8});
+    m[TET10][2].insert({5, 6, 9});
+    m[TET10][3].insert({7, 8, 9});
+    m[TET10][4].insert({0, 1});
+    m[TET10][5].insert({1, 2});
+    m[TET10][6].insert({0, 2});
+    m[TET10][7].insert({0, 3});
+    m[TET10][8].insert({1, 3});
+    m[TET10][9].insert({2, 3});
+
+    m[TET14] = m[TET10];
+    // Nodes 10-13 are face nodes and have no neighbors
+
+    m[PYRAMID5][0].insert({1, 3, 4});
+    m[PYRAMID5][1].insert({0, 2, 4});
+    m[PYRAMID5][2].insert({1, 3, 4});
+    m[PYRAMID5][3].insert({0, 2, 4});
+    m[PYRAMID5][4].insert({0, 1, 2, 3});
+
+    m[PYRAMID13][0].insert({5, 8, 9});
+    m[PYRAMID13][1].insert({5, 6, 10});
+    m[PYRAMID13][2].insert({6, 7, 11});
+    m[PYRAMID13][3].insert({7, 8, 12});
+    m[PYRAMID13][4].insert({9, 10, 11, 12});
+    m[PYRAMID13][5].insert({0, 1});
+    m[PYRAMID13][6].insert({1, 2});
+    m[PYRAMID13][7].insert({2, 3});
+    m[PYRAMID13][8].insert({0, 3});
+    m[PYRAMID13][9].insert({0, 4});
+    m[PYRAMID13][10].insert({1, 4});
+    m[PYRAMID13][11].insert({2, 4});
+    m[PYRAMID13][12].insert({3, 4});
+
+    m[PYRAMID14] = m[PYRAMID13];
+    // Node 13 is a face node and has no neighbors
+
+    m[PYRAMID18] = m[PYRAMID14];
+    // Nodes 14-17 are face nodes and have no neighbors
+
+    m[PRISM6][0].insert({1, 2, 3});
+    m[PRISM6][1].insert({0, 2, 4});
+    m[PRISM6][2].insert({0, 1, 5});
+    m[PRISM6][3].insert({0, 4, 5});
+    m[PRISM6][4].insert({1, 3, 5});
+    m[PRISM6][5].insert({2, 3, 4});
+
+    m[PRISM18][0].insert({6, 8, 9});
+    m[PRISM18][1].insert({6, 7, 10});
+    m[PRISM18][2].insert({7, 8, 11});
+    m[PRISM18][3].insert({9, 12, 14});
+    m[PRISM18][4].insert({10, 12, 13});
+    m[PRISM18][5].insert({11, 13, 14});
+    m[PRISM18][6].insert({0, 1});
+    m[PRISM18][7].insert({1, 2});
+    m[PRISM18][8].insert({0, 2});
+    m[PRISM18][9].insert({0, 3});
+    m[PRISM18][10].insert({1, 4});
+    m[PRISM18][11].insert({2, 5});
+    m[PRISM18][12].insert({3, 4});
+    m[PRISM18][13].insert({4, 5});
+    m[PRISM18][14].insert({3, 5});
+    // Nodes 15-17 are face nodes and have no neighbors
+
+    m[PRISM20] = m[PRISM18];
+    // Nodes 18-19 are face nodes and have no neighbors
+
+    m[PRISM21] = m[PRISM20];
+    // Node 20 is an interior node and has no neighbors
+
+    m[HEX8][0].insert({1, 3, 4});
+    m[HEX8][1].insert({0, 2, 5});
+    m[HEX8][2].insert({1, 3, 6});
+    m[HEX8][3].insert({0, 2, 7});
+    m[HEX8][4].insert({0, 5, 7});
+    m[HEX8][5].insert({1, 4, 6});
+    m[HEX8][6].insert({2, 5, 7});
+    m[HEX8][7].insert({3, 4, 6});
+
+    m[HEX20][0].insert({8, 11, 12});
+    m[HEX20][1].insert({8, 9, 13});
+    m[HEX20][2].insert({9, 10, 14});
+    m[HEX20][3].insert({10, 11, 15});
+    m[HEX20][4].insert({12, 16, 19});
+    m[HEX20][5].insert({13, 16, 17});
+    m[HEX20][6].insert({14, 17, 18});
+    m[HEX20][7].insert({15, 18, 19});
+    m[HEX20][8].insert({0, 1});
+    m[HEX20][9].insert({1, 2});
+    m[HEX20][10].insert({2, 3});
+    m[HEX20][11].insert({0, 3});
+    m[HEX20][12].insert({0, 4});
+    m[HEX20][13].insert({1, 5});
+    m[HEX20][14].insert({2, 6});
+    m[HEX20][15].insert({3, 7});
+    m[HEX20][16].insert({4, 5});
+    m[HEX20][17].insert({5, 6});
+    m[HEX20][18].insert({6, 7});
+    m[HEX20][19].insert({4, 7});
+
+    m[HEX27] = m[HEX20];
+    // Nodes 20-26 are face nodes and have no neighbors
+
+    return m;
+  };
+
+  inline static const Map elem_type_to_neighbor_map = build_elem_type_to_neighbor_map();
+
 protected:
 
   // Builds a 1D mesh with the specified ElemType and number of elements
-  void do_test(unsigned n_elem,
-               ElemType elem_type,
-               dof_id_type * validation_data)
+  void do_test(ElemType elem_type)
   {
-    ReplicatedMesh mesh(*TestCommWorld, /*dim=*/1);
+    const auto * ref_elem = &(ReferenceElem::get(elem_type));
+    const auto dim = ref_elem->dim();
 
-    MeshTools::Generation::build_line(mesh,
-                                      n_elem,
-                                      /*xmin=*/0.,
-                                      /*xmax=*/1.,
-                                      elem_type);
+    ReplicatedMesh mesh(*TestCommWorld, dim);
 
-    // find_nodal_neighbors() needs a data structure which is prepared by another function
-    std::vector<std::vector<const Elem *>> nodes_to_elem_map;
-    MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
+    unsigned int n_elems_per_side = 4 / Elem::type_to_default_order_map[elem_type];
 
-    // Loop over the nodes and call find_nodal_neighbors()
-    {
-      std::vector<const Node*> neighbor_nodes;
+    switch (dim)
+      {
+        case 1:
+          MeshTools::Generation::build_line(mesh, n_elems_per_side, 0., 1., elem_type);
+          break;
+        case 2:
+          MeshTools::Generation::build_square(
+              mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., elem_type);
+          break;
+
+        case 3:
+          MeshTools::Generation::build_cube(mesh,
+                                            n_elems_per_side,
+                                            n_elems_per_side,
+                                            n_elems_per_side,
+                                            0.,
+                                            1.,
+                                            0.,
+                                            1.,
+                                            0.,
+                                            1.,
+                                            elem_type);
+          break;
+
+        default:
+          libmesh_error_msg("Unsupported dimension " << dim);
+      }
 
-      unsigned ctr = 0;
-      for (const auto & node : mesh.node_ptr_range())
-        {
-          MeshTools::find_nodal_neighbors(mesh, *node, nodes_to_elem_map, neighbor_nodes);
+    const auto & neighbor_map = libmesh_map_find(elem_type_to_neighbor_map, elem_type);
 
-          // The entries in neighbor_nodes are just sorted according
-          // to memory address, which is somewhat arbitrary, so create
-          // a vector sorted by IDs for test purposes.
-          std::vector<dof_id_type> neighbor_node_ids(neighbor_nodes.size());
-          for (std::size_t i=0; i<neighbor_nodes.size(); ++i)
-            neighbor_node_ids[i] = neighbor_nodes[i]->id();
-          std::sort(neighbor_node_ids.begin(), neighbor_node_ids.end());
+    // find_nodal_neighbors() needs a data structure which is prepared by another function
+    std::unordered_map<dof_id_type, std::vector<const Elem *>> nodes_to_elem_map;
+    MeshTools::build_nodes_to_elem_map(mesh, nodes_to_elem_map);
 
-          // Compare to validation_data
-          for (std::size_t j=0; j<neighbor_node_ids.size(); ++j)
-            {
-              CPPUNIT_ASSERT_EQUAL( validation_data[2*ctr + j], neighbor_node_ids[j] );
-            }
+    // Loop over the nodes and call find_nodal_neighbors()
+    std::vector<const Node*> neighbor_nodes;
 
-          ++ctr;
-        }
-    }
+    std::set<dof_id_type> node_ids_checked;
+    for (const auto * elem : mesh.element_ptr_range())
+      {
+        for (const auto & node : elem->node_ref_range())
+          {
+            // Have we already checked this node?
+            if (node_ids_checked.find(node.id()) != node_ids_checked.end())
+              continue;
+
+            MeshTools::find_nodal_neighbors(mesh, node, nodes_to_elem_map, neighbor_nodes);
+
+            for (const auto * neigh : neighbor_nodes)
+              {
+                // Find the common elements between node and neigh
+                const auto & elems_containing_node = libmesh_map_find(nodes_to_elem_map, node.id());
+                const auto & elems_containing_neigh = libmesh_map_find(nodes_to_elem_map, neigh->id());
+                const Elem * common_elem = nullptr;
+                for (const auto * neigh_elem : elems_containing_neigh)
+                  {
+                    if ((std::find(elems_containing_node.begin(), elems_containing_node.end(), neigh_elem) !=
+                         elems_containing_node.end()))
+                      common_elem = neigh_elem;
+                    else
+                      continue;
+
+                    // Which local elem node numbers are node and neighbor?
+                    unsigned node_id_local = common_elem->local_node(node.id());
+                    unsigned neigh_id_local = common_elem->local_node(neigh->id());
+
+                    // For a given element type, we know the hard-coded neighbor relationship
+                    // between local node ids. Make sure the neigbor from find_nodal_neighbors
+                    // was intended to be a neighbor. Note that this test does not currently check
+                    // for neighbors that find_nodal_neighbors may not have found.
+                    const auto & correct_neighbor_ids = libmesh_map_find(neighbor_map, node_id_local);
+                    CPPUNIT_ASSERT(correct_neighbor_ids.find(neigh_id_local) != correct_neighbor_ids.end());
+
+                  }
+              }
+
+            node_ids_checked.insert(node.id());
+          }
+      }
   }
 
 public:
@@ -95,71 +336,141 @@ public:
   void testEdge2()
   {
     LOG_UNIT_TEST;
+    do_test(EDGE2);
+  }
 
-    // 11 nodes, 2 neighbor entries per node
-    dof_id_type validation_data[22] =
-      {
-        1, DofObject::invalid_id,
-        0, 2,
-        1, 3,
-        2, 4,
-        3, 5,
-        4, 6,
-        5, 7,
-        6, 8,
-        7, 9,
-        8, 10,
-        9, DofObject::invalid_id
-      };
 
-    do_test(/*n_elem=*/10, EDGE2, validation_data);
+  void testEdge3()
+  {
+    LOG_UNIT_TEST;
+    do_test(EDGE3);
   }
 
 
-  void testEdge3()
+  void testEdge4()
   {
     LOG_UNIT_TEST;
+    do_test(EDGE4);
+  }
 
-    // 11 nodes, 2 neighbor entries per node
-    dof_id_type validation_data[22] =
-      {
-        2, DofObject::invalid_id,
-        2, 4,
-        0, 1,
-        4, 6,
-        1, 3,
-        6, 8,
-        3, 5,
-        8, 10,
-        5, 7,
-        10, DofObject::invalid_id,
-        7, 9
-      };
+  void testTri3()
+  {
+    LOG_UNIT_TEST;
+    do_test(TRI3);
+  }
 
-    do_test(/*n_elem=*/5, EDGE3, validation_data);
+  void testTri6()
+  {
+    LOG_UNIT_TEST;
+    do_test(TRI6);
   }
 
+  void testTri7()
+  {
+    LOG_UNIT_TEST;
+    do_test(TRI7);
+  }
 
-  void testEdge4()
+  void testQuad4()
   {
     LOG_UNIT_TEST;
+    do_test(QUAD4);
+  }
 
-    // 10 nodes, 2 neighbor entries per node
-    dof_id_type validation_data[20] =
-      {
-        2, DofObject::invalid_id,
-        3, 5,
-        0, 3,
-        1, 2,
-        6, 8,
-        1, 6,
-        4, 5,
-        9, DofObject::invalid_id,
-        4, 9,
-        7, 8
-      };
+  void testQuad8()
+  {
+    LOG_UNIT_TEST;
+    do_test(QUAD8);
+  }
+
+  void testQuad9()
+  {
+    LOG_UNIT_TEST;
+    do_test(QUAD9);
+  }
 
-    do_test(/*n_elem=*/3, EDGE4, validation_data);
+  void testTet4()
+  {
+    LOG_UNIT_TEST;
+    do_test(TET4);
+  }
+
+  void testTet10()
+  {
+    LOG_UNIT_TEST;
+    do_test(TET10);
+  }
+
+  void testTet14()
+  {
+    LOG_UNIT_TEST;
+    do_test(TET14);
+  }
+
+  void testPyramid5()
+  {
+    LOG_UNIT_TEST;
+    do_test(PYRAMID5);
+  }
+
+  void testPyramid13()
+  {
+    LOG_UNIT_TEST;
+    do_test(PYRAMID13);
+  }
+
+  void testPyramid14()
+  {
+    LOG_UNIT_TEST;
+    do_test(PYRAMID14);
+  }
+
+  void testPyramid18()
+  {
+    LOG_UNIT_TEST;
+    do_test(PYRAMID18);
+  }
+
+  void testPrism6()
+  {
+    LOG_UNIT_TEST;
+    do_test(PRISM6);
+  }
+
+  void testPrism18()
+  {
+    LOG_UNIT_TEST;
+    do_test(PRISM18);
+  }
+
+  void testPrism20()
+  {
+    LOG_UNIT_TEST;
+    do_test(PRISM20);
+  }
+
+  void testPrism21()
+  {
+    LOG_UNIT_TEST;
+    do_test(PRISM21);
+  }
+
+  void testHex8()
+  {
+    LOG_UNIT_TEST;
+    do_test(HEX8);
+  }
+
+  void testHex20()
+  {
+    LOG_UNIT_TEST;
+    do_test(HEX20);
+  }
+
+  void testHex27()
+  {
+    LOG_UNIT_TEST;
+    do_test(HEX27);
   }
 
   void testOrientation()