diff --git a/src/ASM/DomainDecomposition.C b/src/ASM/DomainDecomposition.C
index 9c8ebfce..798a1d28 100644
--- a/src/ASM/DomainDecomposition.C
+++ b/src/ASM/DomainDecomposition.C
@@ -14,12 +14,15 @@
 #include "DomainDecomposition.h"
 #include "ASMbase.h"
 #include "ASMstruct.h"
+#include "ASM2D.h"
+#include "ASM3D.h"
 #include "LinSolParams.h"
 #include "ProcessAdm.h"
 #include "SAMpatch.h"
 #include "SIMbase.h"
 #include "IFEM.h"
 #include "Utilities.h"
+#include "Vec3.h"
 #include <cassert>
 #include <numeric>
 #include <iostream>
@@ -33,12 +36,34 @@ public:
   //! \param orient Orientation of boundary on slave
   //! \param lIdx Index of boundary on slave
   //! \param basis Basis to use for boundary on slave
-  NodeIterator(const ASMstruct* pch, int orient, int lIdx, int basis)
+  //! \param dim Dimension to iterate over
+  NodeIterator(const ASMstruct* pch, int orient, int lIdx, int basis, int dim = 2)
   {
+    if (dim == 0) {
+      nodes.resize(1);
+      return;
+    }
+
     int n1, n2, n3;
     pch->getSize(n1,n2,n3,basis);
     int nsd = pch->getNoSpaceDim();
     if (nsd == 3) {
+      if (dim == 1) {
+        if (lIdx <= 4)
+          nodes.resize(n1);
+        else if (lIdx >= 5 && lIdx <= 8)
+          nodes.resize(n2);
+        else
+          nodes.resize(n3);
+
+        if (orient == 0)
+          std::iota(nodes.begin(), nodes.end(), 0);
+        else
+          std::iota(nodes.rbegin(), nodes.rend(), 0);
+
+        return;
+      }
+
       int dim1, dim2;
       if (lIdx == 1 || lIdx == 2)
         dim1 = n2, dim2 = n3;
@@ -270,6 +295,58 @@ std::vector<IntVec> DomainDecomposition::calcSubdomains3D(size_t nel1, size_t ne
 }
 
 
+void DomainDecomposition::setupNodeNumbers(int basis, IntVec& lNodes,
+                                           std::set<int>& cbasis,
+                                           const ASMbase* pch,
+                                           int dim, int lidx)
+{
+  if (basis != 0) // specified base
+    cbasis = utl::getDigits(basis);
+  else if (dim == 0 || (dim == 1 && pch->getNoSpaceDim() == 3)) {
+    // need to expand to all bases for corners and edges
+    for (size_t b = 1; b <= pch->getNoBasis(); ++b)
+      cbasis.insert(b);
+  } else // directly add nodes, cbasis remains empty
+    pch->getBoundaryNodes(lidx, lNodes);
+
+  const ASM2D* pch2D = dynamic_cast<const ASM2D*>(pch);
+  const ASM3D* pch3D = dynamic_cast<const ASM3D*>(pch);
+  for (auto& it2 : cbasis) {
+    if (dim == 0) {
+      int node = 0;
+      if (pch2D) {
+        static std::map<int, std::array<int,2>> vmap2D
+         {{1, {-1, -1}},
+          {2, { 1, -1}},
+          {3, {-1,  1}},
+          {4, { 1,  1}}};
+        auto itc = vmap2D.find(lidx);
+        node = pch2D->getCorner(itc->second[0], itc->second[1], it2);
+      } else if (pch3D) {
+        static std::map<int, std::array<int,3>> vmap3D
+         {{1, {-1, -1, -1}},
+          {2, { 1, -1, -1}},
+          {3, {-1,  1, -1}},
+          {4, { 1,  1, -1}},
+          {5, {-1, -1,  1}},
+          {6, { 1, -1,  1}},
+          {7, {-1,  1,  1}},
+          {8, { 1,  1,  1}}};
+        auto itc = vmap3D.find(lidx);
+        node = pch3D->getCorner(itc->second[0], itc->second[1],
+                                          itc->second[2], it2);
+      }
+      lNodes.push_back(pch->getNodeID(node));
+    } else if (dim == 1 && pch3D) {
+      std::vector<int> eNodes = pch3D->getEdge(lidx, false, it2);
+      for (const int& it : eNodes)
+        lNodes.push_back(pch->getNodeID(it));
+    } else
+      pch->getBoundaryNodes(lidx, lNodes, it2);
+  }
+}
+
+
 bool DomainDecomposition::calcGlobalNodeNumbers(const ProcessAdm& adm,
                                                 const SIMbase& sim)
 {
@@ -323,12 +400,8 @@ bool DomainDecomposition::calcGlobalNodeNumbers(const ProcessAdm& adm,
 
     std::set<int> cbasis;
     IntVec lNodes;
-    if (it.basis != 0) {
-      cbasis = utl::getDigits(it.basis);
-      for (auto& it2 : cbasis)
-        sim.getPatch(sidx)->getBoundaryNodes(it.sidx, lNodes, it2);
-    } else
-      sim.getPatch(sidx)->getBoundaryNodes(it.sidx, lNodes);
+
+    setupNodeNumbers(it.basis, lNodes, cbasis, sim.getPatch(sidx), it.dim, it.sidx);
 
     int nRecv;
     adm.receive(nRecv, getPatchOwner(it.master));
@@ -344,7 +417,7 @@ bool DomainDecomposition::calcGlobalNodeNumbers(const ProcessAdm& adm,
     for (size_t b = 1; b <= sim.getPatch(sidx)->getNoBasis(); ++b) {
       if (cbasis.empty() || cbasis.find(b) != cbasis.end()) {
         NodeIterator iter(dynamic_cast<const ASMstruct*>(sim.getPatch(sidx)),
-                          it.orient, it.sidx, b);
+                          it.orient, it.sidx, b, it.dim);
         auto it_n = iter.begin();
         for (size_t i = 0; i < iter.size(); ++i, ++it_n) {
           int node = MLGN[lNodes[i+ofs]-1];
@@ -389,13 +462,8 @@ bool DomainDecomposition::calcGlobalNodeNumbers(const ProcessAdm& adm,
       continue;
 
     std::set<int> cbasis;
-    if (it.basis != 0)
-      cbasis = utl::getDigits(it.basis);
-
     std::vector<int> glbNodes;
-    for (size_t b = 1; b <= sim.getPatch(midx)->getNoBasis(); ++b)
-      if (cbasis.empty() || cbasis.find(b) != cbasis.end())
-        sim.getPatch(midx)->getBoundaryNodes(it.midx, glbNodes, b);
+    setupNodeNumbers(it.basis, glbNodes, cbasis, sim.getPatch(midx), it.dim, it.midx);
 
     for (size_t i = 0; i < glbNodes.size(); ++i)
       glbNodes[i] = MLGN[glbNodes[i]-1];
@@ -410,7 +478,9 @@ bool DomainDecomposition::calcGlobalNodeNumbers(const ProcessAdm& adm,
 
 
 std::vector<int> DomainDecomposition::setupEquationNumbers(const SIMbase& sim,
-                                                           int pidx, int lidx, const std::set<int>& cbasis)
+                                                           int pidx, int lidx,
+                                                           const std::set<int>& cbasis,
+                                                           int dim)
 {
   std::vector<IntVec> lNodes(sim.getPatch(pidx)->getNoBasis());
   std::vector<int> result;
@@ -427,8 +497,10 @@ std::vector<int> DomainDecomposition::setupEquationNumbers(const SIMbase& sim,
       if (!cbasis.empty() && cbasis.find(basis) == cbasis.end())
         continue;
 
-      if (lNodes[basis-1].empty())
-        sim.getPatch(pidx)->getBoundaryNodes(lidx, lNodes[basis-1], basis);
+      if (lNodes[basis-1].empty()) {
+        IntSet dummy;
+        setupNodeNumbers(basis, lNodes[basis-1], dummy, sim.getPatch(pidx), dim, lidx);
+      }
 
       std::set<int> components;
       if (block == 0 || sim.getSolParams()->getBlock(block-1).comps == 0) {
@@ -552,7 +624,7 @@ bool DomainDecomposition::calcGlobalEqNumbers(const ProcessAdm& adm,
     if (it.basis != 0)
       cbasis = utl::getDigits(it.basis);
 
-    IntVec locEqs = setupEquationNumbers(sim, sidx, it.sidx, cbasis);
+    IntVec locEqs = setupEquationNumbers(sim, sidx, it.sidx, cbasis, it.dim);
 
     int nRecv;
     adm.receive(nRecv, getPatchOwner(it.master));
@@ -587,7 +659,7 @@ bool DomainDecomposition::calcGlobalEqNumbers(const ProcessAdm& adm,
           components = utl::getDigits(sim.getSolParams()->getBlock(block-1).comps);
 
         NodeIterator iter(dynamic_cast<const ASMstruct*>(sim.getPatch(sidx)),
-                          it.orient, it.sidx, basis);
+                          it.orient, it.sidx, basis, it.dim);
         auto it_n = iter.begin();
         int nodeDofs = components.size();
         for (size_t i = 0; i < iter.size(); ++i, ++it_n) {
@@ -667,7 +739,7 @@ bool DomainDecomposition::calcGlobalEqNumbers(const ProcessAdm& adm,
     if (it.basis != 0)
       cbasis = utl::getDigits(it.basis);
 
-    IntVec glbEqs = setupEquationNumbers(sim, midx, it.midx, cbasis);
+    IntVec glbEqs = setupEquationNumbers(sim, midx, it.midx, cbasis, it.dim);
     adm.send(int(glbEqs.size()), getPatchOwner(it.slave));
     adm.send(glbEqs, getPatchOwner(it.slave));
   }
@@ -696,10 +768,96 @@ int DomainDecomposition::getGlobalEq(int lEq, size_t idx) const
 }
 
 
+bool DomainDecomposition::sanityCheckCorners(const SIMbase& sim)
+{
+#ifdef HAVE_MPI
+  const ProcessAdm& adm = sim.getProcessAdm();
+  if (!adm.isParallel())
+    return true;
+
+  std::vector<int> sizes(adm.getNoProcs());
+  for (int i = 1; i <= sim.getNoPatches(); ++i)
+    sizes[adm.dd.getPatchOwner(i)] += 5*sim.getPatch(1)->getNoBasis() *
+                                      pow(2,sim.getNoSpaceDim());
+
+  std::vector<int> displ(adm.getNoProcs());
+  for (int i = 1; i < adm.getNoProcs(); ++i)
+    displ[i] = displ[i-1] + sizes[i-1];
+
+  std::vector<double> loc_data;
+  for (int i = 1; i <= sim.getNoPatches(); ++i) {
+    int pIdx;
+    if ((pIdx=sim.getLocalPatchIndex(i)) > 0) {
+      const ASM2D* pch2D = nullptr;
+      if (sim.getNoSpaceDim() == 2)
+        pch2D = dynamic_cast<const ASM2D*>(sim.getPatch(pIdx));
+      const ASM3D* pch3D = nullptr;
+      if (sim.getNoSpaceDim() == 3)
+        pch3D = dynamic_cast<const ASM3D*>(sim.getPatch(pIdx));
+      if (pch2D || pch3D) {
+        for (size_t c = 0; c < pow(2,sim.getPatch(pIdx)->getNoSpaceDim()); ++c) {
+          for (size_t b = 1; b <= sim.getPatch(pIdx)->getNoBasis(); ++b) {
+            int node;
+            if (pch2D)
+              node = pch2D->getCorner((c == 1 || c == 3) ? 1 : -1,
+                                                 c >= 2  ? 1 : -1, b);
+            else if (pch3D)
+              node = pch3D->getCorner((c == 1 || c == 3 || c == 5 || c == 7) ? 1 : -1,
+                                      (c == 2 || c == 3 || c == 6 || c == 7) ? 1 : -1,
+                                                                     c >= 4  ? 1 : -1, b);
+
+            Vec3 pos = sim.getPatch(pIdx)->getCoord(node);
+            loc_data.push_back(pos[0]);
+            loc_data.push_back(pos[1]);
+            loc_data.push_back(pos[2]);
+            loc_data.push_back(b);
+            loc_data.push_back(adm.dd.getMLGN()[sim.getPatch(pIdx)->getNodeID(node)-1]);
+          }
+        }
+      }
+    }
+  }
+
+  std::vector<double> glob_data(displ.back()+sizes.back());
+  MPI_Allgatherv(loc_data.data(), loc_data.size(), MPI_DOUBLE,
+                 glob_data.data(), sizes.data(), displ.data(), MPI_DOUBLE,
+                 *adm.getCommunicator());
+
+  // unpack data
+  std::vector<std::array<double,5>> corners;
+  for (size_t i = 0; i < glob_data.size(); i += 5)
+    corners.push_back({glob_data[i], glob_data[i+1],
+                        glob_data[i+2], glob_data[i+3], glob_data[i+4]});
+
+  for (const auto& c : corners) {
+    auto fail = std::find_if(corners.begin(), corners.end(),
+                            [c](const std::array<double,5>& C)
+                            {
+                              return std::fabs(C[0]-c[0]) < 1e-6 &&
+                                     std::fabs(C[1]-c[1]) < 1e-6 &&
+                                     std::fabs(C[2]-c[2]) < 1e-6 &&
+                                    (int)c[3] == (int)C[3] &&
+                                    (int)c[4] != (int)C[4];
+                            });
+    if (fail != corners.end()) {
+      std::cerr << "** DomainDecomposition::setup ** Corner node " << (int)c[4]
+                << " with coordinates (" << c[0] << "," << c[1] << "," << c[2] << ")"
+                << " found with different global ID " << (int)(*fail)[4]
+                << " for basis " << (int)c[3] << "."
+                << " You have to add vertex / edge connections." << std::endl;
+      return false;
+    }
+  }
+#endif
+
+  return true;
+}
+
+
 bool DomainDecomposition::setup(const ProcessAdm& adm, const SIMbase& sim)
 {
 #ifdef HAVE_MPI
-  if (adm.isParallel()) {
+  if (!adm.isParallel())
     IFEM::cout << "Establishing domain decomposition" << std::endl;
 
 #if SP_DEBUG > 1
@@ -709,8 +867,8 @@ bool DomainDecomposition::setup(const ProcessAdm& adm, const SIMbase& sim)
                   ", slave/idx=" << it.slave << "/" << it.sidx <<
                   ", orient=" << it.orient << ", dim=" << it.dim <<
                   ", basis=" << it.basis << std::endl;
-#endif
   }
+#endif
 #endif
 
   sam = dynamic_cast<const SAMpatch*>(sim.getSAM());
@@ -736,7 +894,12 @@ bool DomainDecomposition::setup(const ProcessAdm& adm, const SIMbase& sim)
   if (ok < adm.getNoProcs())
     return false;
 
+  // sanity check all corners of the patches
+  if (!sanityCheckCorners(sim))
+    return false;
+
   ok = 1;
+
   // Establish local equation mappings for each block.
   if (sim.getSolParams() && sim.getSolParams()->getNoBlocks() > 1) {
     IFEM::cout << "\tEstablishing local block equation numbers" << std::endl;
diff --git a/src/ASM/DomainDecomposition.h b/src/ASM/DomainDecomposition.h
index 3e19c9b0..537a79f6 100644
--- a/src/ASM/DomainDecomposition.h
+++ b/src/ASM/DomainDecomposition.h
@@ -19,6 +19,7 @@
 #include <vector>
 #include <cstddef>
 
+class ASMbase;
 class ProcessAdm;
 class SAMpatch;
 class SIMbase;
@@ -64,6 +65,10 @@ public:
         if (A.slave != B.slave)
           return A.slave < B.slave;
 
+        // then by dim
+        if (A.dim != B.dim)
+          return A.dim < B.dim;
+
         // finally by index on master
         return A.midx < B.midx;
       }
@@ -181,7 +186,20 @@ private:
   //! \param cbasis If non-empty, bases to connect
   std::vector<int> setupEquationNumbers(const SIMbase& sim,
                                         int pidx, int lidx,
-                                        const std::set<int>& cbasis);
+                                        const std::set<int>& cbasis,
+                                        int dim);
+
+  //! \brief Setup node numbers for all bases on a boundary.
+  //! \param basis Bases to grab nodes for
+  //! \param lNodes Resulting nodes
+  //! \param cbasis Bases to connect
+  //! \param pch Patch to obtain nodes for
+  //! \param dim Dimension of boundary to obtain nodes for
+  //! \param lidx Local index of boundary to obtain nodes for
+  void setupNodeNumbers(int basis, std::vector<int>& lNodes,
+                        std::set<int>& cbasis,
+                        const ASMbase* pch,
+                        int dim, int lidx);
 
   //! \brief Calculate the global node numbers for given finite element model.
   bool calcGlobalNodeNumbers(const ProcessAdm& adm, const SIMbase& sim);
@@ -189,6 +207,11 @@ private:
   //! \brief Calculate the global equation numbers for given finite element model.
   bool calcGlobalEqNumbers(const ProcessAdm& adm, const SIMbase& sim);
 
+  //! \brief Sanity check model.
+  //! \details Collects the corners of all patches in the model and make sure nodes
+  //!          with matching coordinates have the same global node ID.
+  bool sanityCheckCorners(const SIMbase& sim);
+
   std::map<int,int> patchOwner; //!< Process that owns a particular patch
 
   //! \brief Struct with information per matrix block.
diff --git a/src/ASM/Test/MPI/TestDomainDecomposition.C b/src/ASM/Test/MPI/TestDomainDecomposition.C
index be0158ab..9cfb1da1 100644
--- a/src/ASM/Test/MPI/TestDomainDecomposition.C
+++ b/src/ASM/Test/MPI/TestDomainDecomposition.C
@@ -13,6 +13,7 @@
 #include "DomainDecomposition.h"
 #include "SAM.h"
 #include "ASMbase.h"
+#include "IntegrandBase.h"
 #include "SIM2D.h"
 #include "SIM3D.h"
 #include "IFEM.h"
@@ -44,6 +45,14 @@ protected:
 };
 
 
+template<class Dim>
+class DummySIM : public Dim {
+public:
+  class DummyIntegrand : public IntegrandBase {};
+  DummySIM() : Dim(1) { Dim::myProblem = new DummyIntegrand; }
+};
+
+
 typedef std::vector<int> IntVec;
 
 static IntVec readIntVector(const std::string& file)
@@ -84,6 +93,29 @@ class TestDomainDecomposition3D : public testing::Test,
 };
 
 
+TEST_P(TestDomainDecomposition2D, Corner)
+{
+  DummySIM<SIM2D> sim;
+  std::stringstream str;
+  str << "src/ASM/Test/refdata/DomainDecomposition_MPI_2D_corner";
+  if (GetParam() == 1)
+    str << "_fail";
+  str << ".xinp";
+  sim.read(str.str().c_str());
+
+  if (GetParam() == 0) {
+    ASSERT_TRUE(sim.preprocess());
+    const ProcessAdm& adm = sim.getProcessAdm();
+    str.str("");
+    str << "src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes";
+    str << adm.getProcId() << ".ref";
+    IntVec B = readIntVector(str.str());
+    check_intvectors_equal(adm.dd.getMLGN(), B);
+  } else
+    ASSERT_FALSE(sim.preprocess());
+}
+
+
 TEST_P(TestDomainDecomposition2D, SetupSingleBasis)
 {
   SIM2D sim(2);
@@ -468,6 +500,32 @@ TEST_P(TestDomainDecomposition3D, SetupMixedBasisPeriodicLM)
 }
 
 
+TEST_P(TestDomainDecomposition3D, Corner)
+{
+  if (GetParam() > 1)
+    return;
+
+  DummySIM<SIM3D> sim;
+  std::stringstream str;
+  str << "src/ASM/Test/refdata/DomainDecomposition_MPI_3D_corner";
+  if (GetParam() == 1)
+    str << "_fail";
+  str << ".xinp";
+  sim.read(str.str().c_str());
+
+  if (GetParam() == 0) {
+    ASSERT_TRUE(sim.preprocess());
+    const ProcessAdm& adm = sim.getProcessAdm();
+    str.str("");
+    str << "src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes";
+    str << adm.getProcId() << ".ref";
+    IntVec B = readIntVector(str.str());
+    check_intvectors_equal(adm.dd.getMLGN(), B);
+  } else
+    ASSERT_FALSE(sim.preprocess());
+}
+
+
 const std::vector<int> orientations2D = {0,1};
 INSTANTIATE_TEST_CASE_P(TestDomainDecomposition2D, TestDomainDecomposition2D, testing::ValuesIn(orientations2D));
 const std::vector<int> orientations3D = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes0.ref b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes0.ref
new file mode 100644
index 00000000..b1bbd6ae
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes0.ref
@@ -0,0 +1,2 @@
+9
+1 2 3 4 5 6 7 8 9
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes1.ref b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes1.ref
new file mode 100644
index 00000000..57842ee7
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes1.ref
@@ -0,0 +1,2 @@
+9
+3 10 11 6 12 13 9 14 15
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes2.ref b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes2.ref
new file mode 100644
index 00000000..80d6d113
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes2.ref
@@ -0,0 +1,2 @@
+9
+16 17 18 19 20 21 11 22 23
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes3.ref b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes3.ref
new file mode 100644
index 00000000..5e85f98a
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_4_corner_nodes3.ref
@@ -0,0 +1,2 @@
+9
+11 22 23 13 24 25 15 26 27
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_corner.xinp b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_corner.xinp
new file mode 100644
index 00000000..516f399f
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_corner.xinp
@@ -0,0 +1,22 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<simulation>
+
+  <console>
+    <logging output_prefix="foo"/>
+  </console>
+
+  <!-- General - geometry definitions !-->
+  <geometry>
+    <partitioning procs="4" nperproc="1"/>
+    <patchfile>src/ASM/Test/refdata/backstep2D-4.g2</patchfile>
+    <refine lowerpatch="1" upperpatch="4" u="1" v="1"/>
+    <topology>
+      <connection master="1" medge="2" slave="2" sedge="1"/>
+      <connection master="2" medge="2" slave="4" sedge="1"/>
+      <connection master="3" medge="4" slave="4" sedge="3"/>
+      <connection master="2" midx="2"  slave="3" sidx="3" dim="0"/>
+    </topology>
+  </geometry>
+
+</simulation>
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_corner_fail.xinp b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_corner_fail.xinp
new file mode 100644
index 00000000..66ef035b
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_2D_corner_fail.xinp
@@ -0,0 +1,21 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<simulation>
+
+  <console>
+    <logging output_prefix="foo"/>
+  </console>
+
+  <!-- General - geometry definitions !-->
+  <geometry>
+    <partitioning procs="4" nperproc="1"/>
+    <patchfile>src/ASM/Test/refdata/backstep2D-4.g2</patchfile>
+    <refine lowerpatch="1" upperpatch="4" u="1" v="1"/>
+    <topology>
+      <connection master="1" medge="2" slave="2" sedge="1"/>
+      <connection master="2" medge="2" slave="4" sedge="1"/>
+      <connection master="3" medge="4" slave="4" sedge="3"/>
+    </topology>
+  </geometry>
+
+</simulation>
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes0.ref b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes0.ref
new file mode 100644
index 00000000..ec4447b9
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes0.ref
@@ -0,0 +1,2 @@
+27
+1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes1.ref b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes1.ref
new file mode 100644
index 00000000..ff6a3e10
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes1.ref
@@ -0,0 +1,2 @@
+27
+3 28 29 6 30 31 9 32 33 12 34 35 15 36 37 18 38 39 21 40 41 24 42 43 27 44 45
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes2.ref b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes2.ref
new file mode 100644
index 00000000..394d124d
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes2.ref
@@ -0,0 +1,2 @@
+27
+46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 29 64 65 31 66 67 33 68 69
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes3.ref b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes3.ref
new file mode 100644
index 00000000..0302c620
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_4_corner_nodes3.ref
@@ -0,0 +1,2 @@
+27
+29 64 65 31 66 67 33 68 69 35 70 71 37 72 73 39 74 75 41 76 77 43 78 79 45 80 81
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_corner.xinp b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_corner.xinp
new file mode 100644
index 00000000..f773293b
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_corner.xinp
@@ -0,0 +1,22 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<simulation>
+
+  <console>
+    <logging output_prefix="foo"/>
+  </console>
+
+  <!-- General - geometry definitions !-->
+  <geometry>
+    <partitioning procs="4" nperproc="1"/>
+    <patchfile>src/ASM/Test/refdata/backstep-4.g2</patchfile>
+    <refine lowerpatch="1" upperpatch="4" u="1" v="1" w="1"/>
+    <topology>
+      <connection master="3" mface="6" slave="4" sface="5"/>
+      <connection master="2" mface="2" slave="4" sface="1"/>
+      <connection master="1" mface="2" slave="2" sface="1"/>
+      <connection master="2" midx="6" slave="3" sidx="7" dim="1"/>
+    </topology>
+  </geometry>
+
+</simulation>
diff --git a/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_corner_fail.xinp b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_corner_fail.xinp
new file mode 100644
index 00000000..8af61c64
--- /dev/null
+++ b/src/ASM/Test/refdata/DomainDecomposition_MPI_3D_corner_fail.xinp
@@ -0,0 +1,21 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+
+<simulation>
+
+  <console>
+    <logging output_prefix="foo"/>
+  </console>
+
+  <!-- General - geometry definitions !-->
+  <geometry>
+    <partitioning procs="4" nperproc="1"/>
+    <patchfile>src/ASM/Test/refdata/backstep-4.g2</patchfile>
+    <refine lowerpatch="1" upperpatch="4" u="1" v="1" w="1"/>
+    <topology>
+      <connection master="3" mface="6" slave="4" sface="5"/>
+      <connection master="2" mface="2" slave="4" sface="1"/>
+      <connection master="1" mface="2" slave="2" sface="1"/>
+    </topology>
+  </geometry>
+
+</simulation>
diff --git a/src/ASM/Test/refdata/backstep-4.g2 b/src/ASM/Test/refdata/backstep-4.g2
new file mode 100644
index 00000000..6536a4cd
--- /dev/null
+++ b/src/ASM/Test/refdata/backstep-4.g2
@@ -0,0 +1,67 @@
+700 1 0 0
+3 0
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+0 0 0
+1 0 0
+0 1 0
+1 1 0
+0 0 1
+1 0 1
+0 1 1
+1 1 1 
+
+700 1 0 0
+3 0
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+1 0 0
+2 0 0
+1 1 0
+2 1 0
+1 0 1
+2 0 1
+1 1 1
+2 1 1 
+
+700 1 0 0
+3 0
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+2 0 -1
+3 0 -1
+2 1 -1
+3 1 -1
+2 0 0
+3 0 0
+2 1 0
+3 1 0 
+
+700 1 0 0
+3 0
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+2 0 0
+3 0 0
+2 1 0
+3 1 0
+2 0 1
+3 0 1
+2 1 1
+3 1 1 
diff --git a/src/ASM/Test/refdata/backstep2D-4.g2 b/src/ASM/Test/refdata/backstep2D-4.g2
new file mode 100644
index 00000000..87268d1f
--- /dev/null
+++ b/src/ASM/Test/refdata/backstep2D-4.g2
@@ -0,0 +1,43 @@
+200 1 0 0
+3 0
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+0 0 0
+1 0 0
+0 1 0
+1 1 0
+
+200 1 0 0
+3 0
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+1 0 0
+2 0 0
+1 1 0
+2 1 0
+
+200 1 0 0
+3 0
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+2 -1 0
+3 -1 0
+2 0 0
+3 0 0
+
+200 1 0 0
+3 0
+2 2
+0 0 1 1
+2 2
+0 0 1 1
+2 0 0
+3 0 0
+2 1 0
+3 1 0