diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C
index 7a708642..0ad6ad0c 100644
--- a/src/ASM/ASMbase.C
+++ b/src/ASM/ASMbase.C
@@ -1630,3 +1630,14 @@ bool ASMbase::writeLagBasis (std::ostream& os, const char* type) const
 
   return true;
 }
+
+
+void ASMbase::getBoundaryElms (int lIndex, IntVec& elms,
+                               int orient, bool local) const
+{
+  this->findBoundaryElms(elms,lIndex,orient);
+  if (local) return;
+
+  for (int& elm : elms)
+    elm = MLGE[elm]-1;
+}
diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h
index db5c0dba..f365885e 100644
--- a/src/ASM/ASMbase.h
+++ b/src/ASM/ASMbase.h
@@ -229,14 +229,15 @@ public:
   virtual Vec3 getCoord(size_t inod) const = 0;
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X nsd\f$\times\f$n-matrix, where \a n is the number of nodes
-  //! \param[in] geo If true returns coordinates for geometry basis
-  //! in the patch
+  //! \param[in] geo If \e true, coordinates for the geometry basis are returned
+  //! otherwise the integration basis coordinates are returned
   virtual void getNodalCoordinates(Matrix& X, bool geo = false) const = 0;
   //! \brief Returns a matrix with nodal coordinates for an element.
-  //! \param[in] iel 1-based element index local to current patch
-  //! \param[in] forceItg If true force returning integration basis coordinates
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in one element
+  //! \param[in] iel 1-based element index local to current patch
+  //! \param[in] forceItg If \e true, return the integration basis coordinates
+  //! otherwise the geometry basis coordinates are returned
   virtual bool getElementCoordinates(Matrix& X, int iel,
                                      bool forceItg = false) const = 0;
 
@@ -245,7 +246,7 @@ public:
   //! \param nodes Array of node numbers
   //! \param[in] basis Which basis to grab nodes for (for mixed methods)
   //! \param[in] thick Thickness of connection
-  //! \param[in] orient Local orientation of the boundary face/edge
+  //! \param[in] orient Local orientation flag (for LR splines only)
   //! \param[in] local If \e true, return patch-local numbers
   virtual void getBoundaryNodes(int lIndex, IntVec& nodes,
                                 int basis = 0, int thick = 1,
@@ -262,11 +263,13 @@ public:
                                  int basis = 0, int orient = -1,
                                  bool local = false, bool open = false) const {}
 
-  //! \brief Finds the global (or patch-local) node numbers on a patch boundary.
-  //! \param[in] lIndex Local index of the boundary face/edge
-  //! \param[in] orient Local orientation of the boundary face/edge
+  //! \brief Finds the global (or patch-local) element numbers on a boundary.
+  //! \param[in] lIndex Local index of the boundary face/edge/vertex
   //! \param[out] elms Array of element numbers
-  virtual void getBoundaryElms(int lIndex, int orient, IntVec& elms) const = 0;
+  //! \param[in] orient Local orientation flag (for LR splines only)
+  //! \param[in] local If \e true, return patch-local element numbers
+  void getBoundaryElms(int lIndex, IntVec& elms,
+                       int orient = -1, bool local = false) const;
 
   //! \brief Returns (1-based) index of a predefined node set in the patch.
   virtual int getNodeSetIdx(const std::string&) const { return 0; }
@@ -821,6 +824,12 @@ protected:
   int searchCtrlPt(RealArray::const_iterator cit, RealArray::const_iterator end,
                    const Vec3& X, int dimension, double tol = 0.001) const;
 
+  //! \brief Finds the patch-local element numbers on a patch boundary.
+  //! \param[out] elms Array of element numbers
+  //! \param[in] lIndex Local index of the boundary face/edge/vertex
+  //! \param[in] orient Local orientation flag (for LR splines only)
+  virtual void findBoundaryElms(IntVec& elms, int lIndex, int orient) const = 0;
+
   //! \brief Assembles L2-projection matrices for the secondary solution.
   //! \param[out] A Left-hand-side matrix
   //! \param[out] B Right-hand-side vectors
diff --git a/src/ASM/ASMs1D.C b/src/ASM/ASMs1D.C
index d80781f7..224a2644 100644
--- a/src/ASM/ASMs1D.C
+++ b/src/ASM/ASMs1D.C
@@ -1893,10 +1893,7 @@ void ASMs1D::getElmConnectivities (IntMat& neigh) const
 }
 
 
-void ASMs1D::getBoundaryElms (int lIndex, int, IntVec& elms) const
+void ASMs1D::findBoundaryElms (IntVec& elms, int lIndex, int) const
 {
-  if (lIndex == 1)
-    elms = {0};
-  else
-    elms = {MLGE[nel-1]-1};
+  elms = { lIndex == 1 ? 0 : static_cast<int>(nel)-1 };
 }
diff --git a/src/ASM/ASMs1D.h b/src/ASM/ASMs1D.h
index efbd4a50..6ebf8c3f 100644
--- a/src/ASM/ASMs1D.h
+++ b/src/ASM/ASMs1D.h
@@ -125,11 +125,6 @@ public:
   virtual void getBoundaryNodes(int lIndex, IntVec& nodes,
                                 int, int thick, int, bool local) const;
 
-  //! \brief Finds the global (or patch-local) node number on a patch end.
-  //! \param[in] lIndex Local index of the end point
-  //! \param[out] elms Array of global element numbers
-  virtual void getBoundaryElms(int lIndex, int, IntVec& elms) const;
-
   //! \brief Finds the node that is closest to the given point.
   //! \param[in] X Global coordinates of point to search for
   //! \return 1-based nodal index and distance to to the found node
@@ -344,6 +339,11 @@ protected:
   // Internal utility methods
   // ========================
 
+  //! \brief Finds the path-local element numbers on a patch boundary.
+  //! \param[out] elms Array of element numbers
+  //! \param[in] lIndex Local index of the end point
+  virtual void findBoundaryElms(IntVec& elms, int lIndex, int = 0) const;
+
   //! \brief Assembles L2-projection matrices for the secondary solution.
   //! \param[out] A Left-hand-side matrix
   //! \param[out] B Right-hand-side vectors
diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C
index fcd9243d..c238d12b 100644
--- a/src/ASM/ASMs2D.C
+++ b/src/ASM/ASMs2D.C
@@ -3239,7 +3239,7 @@ void ASMs2D::getElmConnectivities (IntMat& neigh) const
 }
 
 
-void ASMs2D::getBoundaryElms (int lIndex, int, IntVec& elms) const
+void ASMs2D::findBoundaryElms (IntVec& elms, int lIndex, int) const
 {
   int N1m = surf->numCoefs_u() - surf->order_u() + 1;
   int N2m = surf->numCoefs_v() - surf->order_v() + 1;
@@ -3250,13 +3250,13 @@ void ASMs2D::getBoundaryElms (int lIndex, int, IntVec& elms) const
   case 2:
     elms.reserve(N2m);
     for (int i = 0; i < N2m; ++i)
-      elms.push_back(MLGE[i*N1m + (lIndex-1)*(N1m-1)] - 1);
+      elms.push_back(i*N1m + (lIndex-1)*(N1m-1));
     break;
   case 3:
   case 4:
     elms.reserve(N1m);
     for (int i = 0; i < N1m; ++i)
-      elms.push_back(MLGE[i + (lIndex-3)*N1m*(N2m-1)] - 1);
+      elms.push_back(i + (lIndex-3)*N1m*(N2m-1));
   }
 }
 
diff --git a/src/ASM/ASMs2D.h b/src/ASM/ASMs2D.h
index f773e95c..fc9e947f 100644
--- a/src/ASM/ASMs2D.h
+++ b/src/ASM/ASMs2D.h
@@ -80,8 +80,7 @@ protected:
     virtual ~BasisFunctionCache() = default;
 
     //! \brief Returns number of elements in each direction.
-    const std::array<size_t,2>& noElms() const
-    { return nel; }
+    const std::array<size_t,2>& noElms() const { return nel; }
 
   protected:
     //! \brief Implementation specific initialization.
@@ -105,22 +104,21 @@ protected:
     //! \param reduced True to return index for reduced quadrature
     size_t index(size_t el, size_t gp, bool reduced) const override;
 
-  protected:
-    //! \brief Configure quadratures.
-    bool setupQuadrature();
-
     //! \brief Setup integration point parameters.
     virtual void setupParameters();
 
     const ASMs2D& patch; //!< Reference to patch cache is for
 
     int basis; //!< Basis to use
-    std::array<size_t,2> nel{0,0}; //!< Number of elements in each direction
+    std::array<size_t,2> nel; //!< Number of elements in each direction
 
   private:
     //! \brief Obtain structured element indices.
     //! \param el Global element index
     std::array<size_t,2> elmIndex(size_t el) const;
+
+    //! \brief Configure quadratures.
+    bool setupQuadrature();
   };
 
 public:
@@ -238,16 +236,17 @@ public:
   virtual int getNodeID(size_t inod, bool noAddedNodes = false) const;
 
   //! \brief Returns a matrix with nodal coordinates for an element.
-  //! \param[in] iel Element index
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in one element
-  //! \param[in] forceItg If true return integration basis element coordinates
+  //! \param[in] iel 1-based element index local to current patch
+  //! \param[in] forceItg If \e true, return the integration basis coordinates
+  //! otherwise the geometry basis coordinates are returned
   virtual bool getElementCoordinates(Matrix& X, int iel, bool forceItg = false) const;
 
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
-  //! \param[in] geo If true return coordinates of geometry basis
-  //! in the patch
+  //! \param[in] geo If \e true, coordinates for the geometry basis are returned
+  //! otherwise the integration basis coordinates are returned
   virtual void getNodalCoordinates(Matrix& X, bool geo = false) const;
 
   //! \brief Returns the global coordinates for the given node.
@@ -268,11 +267,6 @@ public:
                                 int basis, int thick = 1,
                                 int = 0, bool local = false) const;
 
-  //! \brief Finds the global (or patch-local) node numbers on a patch boundary.
-  //! \param[in] lIndex Local index of the boundary face/edge
-  //! \param[out] elms Array of element numbers
-  virtual void getBoundaryElms(int lIndex, int, IntVec& elms) const;
-
   //! \brief Returns the node index for a given corner.
   //! \param[in] I -1 or +1 for either umin or umax corner
   //! \param[in] J -1 or +1 for either vmin or vmax corner
@@ -317,6 +311,7 @@ public:
   //! \brief Checks if a separate projection basis is used for this patch.
   virtual bool separateProjectionBasis() const;
 
+
   // Various methods for preprocessing of boundary conditions and patch topology
   // ===========================================================================
 
@@ -589,6 +584,11 @@ protected:
   // Internal utility methods
   // ========================
 
+  //! \brief Finds the patch-local element numbers on a patch boundary.
+  //! \param[out] elms Array of element numbers
+  //! \param[in] lIndex Local index of the boundary edge
+  virtual void findBoundaryElms(IntVec& elms, int lIndex, int = 0) const;
+
   //! \brief Assembles L2-projection matrices for the secondary solution.
   //! \param[out] A Left-hand-side matrix
   //! \param[out] B Right-hand-side vectors
diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C
index baa81f06..2f3ab533 100644
--- a/src/ASM/ASMs3D.C
+++ b/src/ASM/ASMs3D.C
@@ -3776,7 +3776,7 @@ void ASMs3D::getElmConnectivities (IntMat& neigh) const
 }
 
 
-void ASMs3D::getBoundaryElms (int lIndex, int, IntVec& elms) const
+void ASMs3D::findBoundaryElms (IntVec& elms, int lIndex, int) const
 {
   int N1m = svol->numCoefs(0) - svol->order(0) + 1;
   int N2m = svol->numCoefs(1) - svol->order(1) + 1;
@@ -3789,21 +3789,21 @@ void ASMs3D::getBoundaryElms (int lIndex, int, IntVec& elms) const
     elms.reserve(N2m*N3m);
     for (int k = 0; k < N3m; ++k)
       for (int j = 0; j < N2m; ++j)
-        elms.push_back(MLGE[j*N1m + k*N1m*N2m + (lIndex-1)*(N1m-1)] - 1);
+        elms.push_back(j*N1m + k*N1m*N2m + (lIndex-1)*(N1m-1));
     break;
   case 3:
   case 4:
     elms.reserve(N1m*N3m);
     for (int k = 0; k < N3m; ++k)
       for (int i = 0; i < N1m; ++i)
-        elms.push_back(MLGE[i + k*N1m*N2m + (lIndex-3)*(N1m*(N2m-1))] - 1);
+        elms.push_back(i + k*N1m*N2m + (lIndex-3)*(N1m*(N2m-1)));
     break;
   case 5:
   case 6:
     elms.reserve(N1m*N2m);
     for (int j = 0; j < N2m; ++j)
       for (int i = 0; i < N1m; ++i)
-        elms.push_back(MLGE[i + j*N1m + (lIndex-5)*(N1m*N2m*(N3m-1))] - 1);
+        elms.push_back(i + j*N1m + (lIndex-5)*(N1m*N2m*(N3m-1)));
   }
 }
 
diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h
index a1792c2c..9171b667 100644
--- a/src/ASM/ASMs3D.h
+++ b/src/ASM/ASMs3D.h
@@ -96,10 +96,9 @@ protected:
     virtual ~BasisFunctionCache() = default;
 
     //! \brief Returns number of elements in each direction.
-    const std::array<size_t,3>& noElms() const
-    { return nel; }
+    const std::array<size_t,3>& noElms() const { return nel; }
 
- protected:
+  protected:
     //! \brief Implementation specific initialization.
     bool internalInit() override;
 
@@ -108,7 +107,7 @@ protected:
 
     //! \brief Calculates basis function info in a single integration point.
     //! \param el Element of integration point (0-indexed)
-    //! \param gp Integratin point on element (0-indexed)
+    //! \param gp Integration point on element (0-indexed)
     //! \param reduced If true, returns values for reduced integration scheme
     BasisFunctionVals calculatePt(size_t el, size_t gp, bool reduced) const override;
 
@@ -129,7 +128,7 @@ protected:
     int basis; //!< Basis to use
     std::array<size_t,3> nel; //!< Number of elements in each direction
 
-private:
+  private:
     //! \brief Obtain structured element indices.
     //! \param el Global element index
     std::array<size_t,3> elmIndex(size_t el) const;
@@ -210,9 +209,9 @@ public:
   //! \param[in] dir Parameter direction defining which boundary to return
   virtual Go::SplineSurface* getBoundary(int dir, int = 1);
   //! \brief Returns the spline volume representing a basis of this patch.
-  virtual const Go::SplineVolume* getBasis(int basis = 1) const;
-  //! \brief Returns the spline volume representing a basis of this patch.
   virtual Go::SplineVolume* getBasis(int basis = 1);
+  //! \brief Returns the spline volume representing a basis of this patch.
+  virtual const Go::SplineVolume* getBasis(int basis = 1) const;
   //! \brief Copies the parameter domain from the \a other patch.
   virtual void copyParameterDomain(const ASMbase* other);
 
@@ -253,16 +252,17 @@ public:
   virtual int getNodeID(size_t inod, bool noAddedNodes = false) const;
 
   //! \brief Returns a matrix with nodal coordinates for an element.
-  //! \param[in] iel Element index
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
   //! in one element
-  //! \param[in] forceItg If true return integration basis element coordinates
+  //! \param[in] iel 1-based element index local to current patch
+  //! \param[in] forceItg If \e true, return the integration basis coordinates
+  //! otherwise the geometry basis coordinates are returned
   virtual bool getElementCoordinates(Matrix& X, int iel, bool forceItg = false) const;
 
   //! \brief Returns a matrix with all nodal coordinates within the patch.
   //! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
-  //! \param[in] geo If true returns coordinates for geometry basis
-  //! in the patch
+  //! \param[in] geo If \e true, coordinates for the geometry basis are returned
+  //! otherwise the integration basis coordinates are returned
   virtual void getNodalCoordinates(Matrix& X, bool geo = false) const;
 
   //! \brief Returns the global coordinates for the given node.
@@ -291,11 +291,6 @@ public:
   virtual void getBoundary1Nodes(int lEdge, IntVec& nodes, int basis, int = 0,
                                  bool local = false, bool open = false) const;
 
-  //! \brief Finds the global (or patch-local) node numbers on a patch boundary.
-  //! \param[in] lIndex Local index of the boundary face/edge
-  //! \param[out] elms Array of element numbers
-  virtual void getBoundaryElms(int lIndex, int, IntVec& elms) const;
-
   //! \brief Returns the node index for a given corner.
   //! \param[in] I -1 or +1 for either umin or umax corner
   //! \param[in] J -1 or +1 for either vmin or vmax corner
@@ -660,6 +655,11 @@ protected:
   // Internal utility methods
   // ========================
 
+  //! \brief Finds the patch-local element numbers on a patch boundary.
+  //! \param[out] elms Array of element numbers
+  //! \param[in] lIndex Local index of the boundary face
+  virtual void findBoundaryElms(IntVec& elms, int lIndex, int = 0) const;
+
   //! \brief Assembles L2-projection matrices for the secondary solution.
   //! \param[out] A Left-hand-side matrix
   //! \param[out] B Right-hand-side vectors
diff --git a/src/ASM/ASMsupel.h b/src/ASM/ASMsupel.h
index a1c56906..2b36d06a 100644
--- a/src/ASM/ASMsupel.h
+++ b/src/ASM/ASMsupel.h
@@ -76,7 +76,7 @@ public:
                                 int, int, int, bool local) const;
 
   //! \brief Dummy method doing nothing.
-  virtual void getBoundaryElms(int, int, IntVec&) const {}
+  virtual void findBoundaryElms(IntVec&, int, int) const {}
   //! \brief Dummy method doing nothing.
   virtual bool getParameterDomain(Real2DMat&, IntVec*) const { return false; }
   //! \brief Dummy method doing nothing.
diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C
index 131bd6ae..a652738b 100644
--- a/src/ASM/LR/ASMu2D.C
+++ b/src/ASM/LR/ASMu2D.C
@@ -2928,7 +2928,7 @@ void ASMu2D::getElmConnectivities (IntMat& neigh) const
 }
 
 
-void ASMu2D::getBoundaryElms (int lIndex, int orient, IntVec& elms) const
+void ASMu2D::findBoundaryElms (IntVec& elms, int lIndex, int orient) const
 {
   std::vector<LR::Element*> elements;
   switch (lIndex) {
@@ -2939,18 +2939,18 @@ void ASMu2D::getBoundaryElms (int lIndex, int orient, IntVec& elms) const
   default: return;
   }
 
-  // Lambda function for sorting wrt. element centre coordinate
-  auto&& onMidPoint = [orient,lIndex](LR::Element* a, LR::Element* b)
-  {
-    int index = lIndex < 3 ? 1 : 0;
-    double am = a->midpoint()[index];
-    double bm = b->midpoint()[index];
-    return orient == 1 ? bm < am : am < bm;
-  };
+  if (orient >= 0)
+    std::sort(elements.begin(), elements.end(),
+              [orient,lIndex](LR::Element* a, LR::Element* b)
+              {
+                int index = lIndex < 3 ? 1 : 0;
+                double am = a->midpoint()[index];
+                double bm = b->midpoint()[index];
+                return orient == 1 ? bm < am : am < bm;
+              });
 
-  std::sort(elements.begin(),elements.end(),onMidPoint);
   for (const LR::Element* elem : elements)
-    elms.push_back(MLGE[elem->getId()]-1);
+    elms.push_back(elem->getId());
 }
 
 
diff --git a/src/ASM/LR/ASMu2D.h b/src/ASM/LR/ASMu2D.h
index 1e063253..5b884d95 100644
--- a/src/ASM/LR/ASMu2D.h
+++ b/src/ASM/LR/ASMu2D.h
@@ -190,12 +190,6 @@ public:
   virtual void getBoundaryNodes(int lIndex, IntVec& nodes, int basis, int = 1,
                                 int orient = 0, bool local = false) const;
 
-  //! \brief Finds the global (or patch-local) node numbers on a patch boundary.
-  //! \param[in] lIndex Local index of the boundary face/edge
-  //! \param[in] orient Orientation of boundary (used for sorting)
-  //! \param[out] elms Array of element numbers
-  virtual void getBoundaryElms(int lIndex, int orient, IntVec& elms) const;
-
   //! \brief Returns the polynomial order in each parameter direction.
   //! \param[out] p1 Order in first (u) direction
   //! \param[out] p2 Order in second (v) direction
@@ -560,6 +554,12 @@ protected:
   // Internal utility methods
   // ========================
 
+  //! \brief Finds the patch-local element numbers on a patch boundary.
+  //! \param[out] elms Array of element numbers
+  //! \param[in] lIndex Local index of the boundary edge
+  //! \param[in] orient Orientation of boundary (used for sorting)
+  virtual void findBoundaryElms(IntVec& elms, int lIndex, int orient) const;
+
   //! \brief Assembles L2-projection matrices for the secondary solution.
   //! \param[out] A Left-hand-side matrix
   //! \param[out] B Right-hand-side vectors
diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C
index 28242790..3b480360 100644
--- a/src/ASM/LR/ASMu3D.C
+++ b/src/ASM/LR/ASMu3D.C
@@ -2355,33 +2355,32 @@ void ASMu3D::getElmConnectivities (IntMat& neigh) const
 }
 
 
-void ASMu3D::getBoundaryElms (int lIndex, int orient, IntVec& elms) const
+void ASMu3D::findBoundaryElms (IntVec& elms, int lIndex, int orient) const
 {
   std::vector<LR::Element*> elements;
   this->getBasis(1)->getEdgeElements(elements,getFaceEnum(lIndex));
 
-  std::sort(elements.begin(), elements.end(),
-            [lIndex,orient](const LR::Element* a, const LR::Element* b)
-            {
-              int dir = (lIndex - 1) / 2;
-              int u = dir == 0 ? 1 : 0;
-              int v = 1 + (dir != 2 ? 1 : 0);
-              int idx = (orient & 4) ? v : u;
-              auto A = a->midpoint();
-              auto B = b->midpoint();
-              if (A[idx] != B[idx])
-                return (orient & 2) ? A[idx] > B[idx] : A[idx] < B[idx];
+  if (orient >= 0)
+    std::sort(elements.begin(), elements.end(),
+              [lIndex,orient](const LR::Element* a, const LR::Element* b)
+              {
+                int u = lIndex <= 2 ? 1 : 0;
+                int v = lIndex >= 5 ? 1 : 2;
+                int idx = (orient & 4) ? v : u;
+                std::vector<double> A = a->midpoint();
+                std::vector<double> B = b->midpoint();
+                if (A[idx] != B[idx])
+                  return (orient & 2) ? A[idx] > B[idx] : A[idx] < B[idx];
 
-              idx = (orient & 4) ? u : v;
-              if (A[idx] != B[idx])
-                return (orient & 1) ? A[idx] > B[idx] : A[idx] < B[idx];
-
-              return false;
-            });
+                idx = (orient & 4) ? u : v;
+                if (A[idx] != B[idx])
+                  return (orient & 1) ? A[idx] > B[idx] : A[idx] < B[idx];
 
+                return false;
+              });
 
   for (const LR::Element* elem : elements)
-    elms.push_back(MLGE[elem->getId()]-1);
+    elms.push_back(elem->getId());
 }
 
 
diff --git a/src/ASM/LR/ASMu3D.h b/src/ASM/LR/ASMu3D.h
index 7cf3846f..b1a3e925 100644
--- a/src/ASM/LR/ASMu3D.h
+++ b/src/ASM/LR/ASMu3D.h
@@ -189,12 +189,6 @@ public:
                                  int orient, bool local,
                                  bool open = false) const;
 
-  //! \brief Finds the global (or patch-local) node numbers on a patch boundary.
-  //! \param[in] lIndex Local index of the boundary face/edge
-  //! \param[in] orient Orientation of boundary (used for sorting)
-  //! \param[out] elms Array of element numbers
-  virtual void getBoundaryElms(int lIndex, int orient, IntVec& elms) const;
-
   //! \brief Returns the node index for a given corner.
   //! \param[in] I -1 or +1 for either umin or umax corner
   //! \param[in] J -1 or +1 for either vmin or vmax corner
@@ -572,6 +566,12 @@ protected:
   // Internal utility methods
   // ========================
 
+  //! \brief Finds the patch-local element numbers on a patch boundary.
+  //! \param[out] elms Array of element numbers
+  //! \param[in] lIndex Local index of the boundary face
+  //! \param[in] orient Orientation of boundary (used for sorting)
+  virtual void findBoundaryElms(IntVec& elms, int lIndex, int orient) const;
+
   //! \brief Assembles L2-projection matrices for the secondary solution.
   //! \param[out] A Left-hand-side matrix
   //! \param[out] B Right-hand-side vectors
diff --git a/src/ASM/Test/TestASMs1D.C b/src/ASM/Test/TestASMs1D.C
index 975d6566..924c6ce0 100644
--- a/src/ASM/Test/TestASMs1D.C
+++ b/src/ASM/Test/TestASMs1D.C
@@ -134,7 +134,7 @@ TEST(TestASMs1D, BoundaryElements)
   IntMat neigh(3);
   std::array<IntVec,2> n;
   for (size_t i = 1; i <= 2; ++i) {
-    pch1.getBoundaryElms(i, 0, n[i-1]);
+    pch1.getBoundaryElms(i, n[i-1]);
     ASSERT_EQ(n[i-1].size(), 1U);
     EXPECT_EQ(n[i-1][0], i == 1 ? 0 : 2);
   }
diff --git a/src/ASM/Test/TestASMs2D.C b/src/ASM/Test/TestASMs2D.C
index 4955af13..92a8d536 100644
--- a/src/ASM/Test/TestASMs2D.C
+++ b/src/ASM/Test/TestASMs2D.C
@@ -50,7 +50,7 @@ TEST(TestASMs2D, BoundaryElements)
 
   std::array<IntVec,4> n;
   for (size_t i = 1; i <= 4; ++i) {
-    pch1.getBoundaryElms(i, 0, n[i-1]);
+    pch1.getBoundaryElms(i, n[i-1]);
     ASSERT_EQ(n[i-1].size(), ref[i-1].size());
     for (size_t j = 0; j < ref[i-1].size(); ++j)
       EXPECT_EQ(n[i-1][j], ref[i-1][j]);
diff --git a/src/ASM/Test/TestASMs3D.C b/src/ASM/Test/TestASMs3D.C
index 2f4e066e..9836f88b 100644
--- a/src/ASM/Test/TestASMs3D.C
+++ b/src/ASM/Test/TestASMs3D.C
@@ -64,7 +64,7 @@ TEST(TestASMs3D, BoundaryElements)
 
   std::array<IntVec,6> n;
   for (size_t i = 1; i <= 6; ++i) {
-    pch1.getBoundaryElms(i, 0, n[i-1]);
+    pch1.getBoundaryElms(i, n[i-1]);
     ASSERT_EQ(n[i-1].size(), ref[i-1].size());
     for (size_t j = 0; j < ref[i-1].size(); ++j)
       EXPECT_EQ(n[i-1][j], ref[i-1][j]);
diff --git a/src/SIM/SIMinput.C b/src/SIM/SIMinput.C
index 55cc4554..d704d1e7 100644
--- a/src/SIM/SIMinput.C
+++ b/src/SIM/SIMinput.C
@@ -1854,8 +1854,8 @@ IntMat SIMinput::getElmConnectivities () const
     if (iface.dim == static_cast<int>(nsd)-1)
     {
       IntVec sElms, mElms;
-      myModel[iface.slave-1]->getBoundaryElms(iface.sidx, iface.orient, sElms);
-      myModel[iface.master-1]->getBoundaryElms(iface.midx, 0, mElms);
+      myModel[iface.slave-1]->getBoundaryElms(iface.sidx, sElms, iface.orient);
+      myModel[iface.master-1]->getBoundaryElms(iface.midx, mElms, 0);
       DomainDecomposition::OrientIterator iter(myModel[iface.slave-1],
                                                iface.orient, iface.sidx);