diff --git a/src/ASM/ASMs2DLag.C b/src/ASM/ASMs2DLag.C
index eb7adcc6..899790bc 100644
--- a/src/ASM/ASMs2DLag.C
+++ b/src/ASM/ASMs2DLag.C
@@ -789,6 +789,28 @@ void ASMs2DLag::generateThreadGroups (const Integrand&, bool, bool)
 }
 
 
+void ASMs2DLag::findBoundaryElms (IntVec& elms, int lIndex, int) const
+{
+  const int N1m = (nx-1)/(p1-1);
+  const int N2m = (ny-1)/(p2-1);
+
+  elms.clear();
+  switch (lIndex) {
+  case 1:
+  case 2:
+    elms.reserve(N2m);
+    for (int i = 0; i < N2m; ++i)
+      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(i + (lIndex-3)*N1m*(N2m-1));
+  }
+}
+
+
 bool ASMs2DLag::write(std::ostream& os, int) const
 {
   return this->writeLagBasis(os, "quad");
diff --git a/src/ASM/ASMs2DLag.h b/src/ASM/ASMs2DLag.h
index 29b1c625..a72daf29 100644
--- a/src/ASM/ASMs2DLag.h
+++ b/src/ASM/ASMs2DLag.h
@@ -128,6 +128,11 @@ protected:
   //! \param[in] Xnod Coordinates of the node
   void setCoord(size_t inod, const Vec3& Xnod);
 
+  //! \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 Find element for parameter, and optionally calculate local coordinates.
   int findElement(double u, double v,
                   double* xi = nullptr, double* eta = nullptr) const;
diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C
index 2f3ab533..1b2f3c45 100644
--- a/src/ASM/ASMs3D.C
+++ b/src/ASM/ASMs3D.C
@@ -3507,59 +3507,50 @@ void ASMs3D::generateThreadGroups (char lIndex, bool silence, bool)
 {
   if (threadGroupsFace.find(lIndex) != threadGroupsFace.end()) return;
 
-  const int p1 = svol->order(0);
-  const int p2 = svol->order(1);
-  const int p3 = svol->order(2);
+  const int p1 = svol->order(0) - 1;
+  const int p2 = svol->order(1) - 1;
+  const int p3 = svol->order(2) - 1;
   const int n1 = svol->numCoefs(0);
   const int n2 = svol->numCoefs(1);
   const int n3 = svol->numCoefs(2);
 
-  // Find elements that are on the boundary face 'lIndex'
-  IntVec map; map.reserve(this->getNoBoundaryElms(lIndex,2));
-  int iel = 0;
-  for (int i3 = p3; i3 <= n3; i3++)
-    for (int i2 = p2; i2 <= n2; i2++)
-      for (int i1 = p1; i1 <= n1; i1++, iel++)
-	switch (lIndex)
-          {
-	  case 1: if (i1 == p1) map.push_back(iel); break;
-	  case 2: if (i1 == n1) map.push_back(iel); break;
-	  case 3: if (i2 == p2) map.push_back(iel); break;
-	  case 4: if (i2 == n2) map.push_back(iel); break;
-	  case 5: if (i3 == p3) map.push_back(iel); break;
-	  case 6: if (i3 == n3) map.push_back(iel); break;
-          }
-
+  // Flag the non-zero knot-spans
   std::vector<bool> el1, el2, el3;
-  el1.reserve(n1 - p1 + 1);
-  el2.reserve(n2 - p2 + 1);
-  el3.reserve(n3 - p3 + 1);
-
-  if (lIndex > 2)
-    for (int i = p1-1; i < n1; i++)
+  if (lIndex > 2) {
+    el1.reserve(n1-p1);
+    for (int i = p1; i < n1; i++)
       el1.push_back(svol->knotSpan(0,i) > 0.0);
-  if (lIndex < 3 || lIndex > 4)
-    for (int i = p2-1; i < n2; i++)
+  }
+  if (lIndex < 3 || lIndex > 4) {
+    el2.reserve(n2-p2);
+    for (int i = p2; i < n2; i++)
       el2.push_back(svol->knotSpan(1,i) > 0.0);
-  if (lIndex < 6)
-    for (int i = p3-1; i < n3; i++)
+  }
+  if (lIndex < 6) {
+    el3.reserve(n3-p3);
+    for (int i = p3; i < n3; i++)
       el3.push_back(svol->knotSpan(2,i) > 0.0);
+  }
 
   ThreadGroups& fGrp = threadGroupsFace[lIndex];
   switch (lIndex)
     {
     case 1:
     case 2:
-      fGrp.calcGroups(el2,el3,p2-1,p3-1);
+      fGrp.calcGroups(el2,el3,p2,p3);
       break;
     case 3:
     case 4:
-      fGrp.calcGroups(el1,el3,p1-1,p3-1);
+      fGrp.calcGroups(el1,el3,p1,p3);
       break;
     default:
-      fGrp.calcGroups(el1,el2,p1-1,p2-1);
+      fGrp.calcGroups(el1,el2,p1,p2);
     }
 
+  // Find elements that are on the boundary face 'lIndex'
+  IntVec map;
+  this->findBoundaryElms(map,lIndex);
+
   fGrp.applyMap(map);
 
   if (silence || fGrp.size() < 2) return;
diff --git a/src/ASM/ASMs3DLag.C b/src/ASM/ASMs3DLag.C
index cabb311c..9fead94e 100644
--- a/src/ASM/ASMs3DLag.C
+++ b/src/ASM/ASMs3DLag.C
@@ -1079,45 +1079,58 @@ void ASMs3DLag::generateThreadGroups (char lIndex, bool, bool)
 {
   if (threadGroupsFace.find(lIndex) != threadGroupsFace.end()) return;
 
-  const int n1 = (nx-1)/(p1-1);
-  const int n2 = (ny-1)/(p2-1);
-  const int n3 = (nz-1)/(p3-1);
-
-  // Find elements that are on the boundary face 'lIndex'
-  IntVec map; map.reserve(this->getNoBoundaryElms(lIndex,2));
-  int d1, d2, iel = 0;
-  for (int i3 = 1; i3 <= n3; i3++)
-    for (int i2 = 1; i2 <= n2; i2++)
-      for (int i1 = 1; i1 <= n1; i1++, iel++)
-        switch (lIndex)
-        {
-          case 1: if (i1 ==  1) map.push_back(iel); break;
-          case 2: if (i1 == n1) map.push_back(iel); break;
-          case 3: if (i2 ==  1) map.push_back(iel); break;
-          case 4: if (i2 == n2) map.push_back(iel); break;
-          case 5: if (i3 ==  1) map.push_back(iel); break;
-          case 6: if (i3 == n3) map.push_back(iel); break;
-        }
-
+  ThreadGroups& fGrp = threadGroupsFace[lIndex];
   switch (lIndex)
   {
     case 1:
     case 2:
-      d1 = n2;
-      d2 = n3;
+      fGrp.calcGroups((ny-1)/(p2-1), (nz-1)/(p3-1), 1);
       break;
     case 3:
     case 4:
-      d1 = n1;
-      d2 = n3;
+      fGrp.calcGroups((nx-1)/(p1-1), (nz-1)/(p3-1), 1);
       break;
     default:
-      d1 = n1;
-      d2 = n2;
+      fGrp.calcGroups((nx-1)/(p1-1), (ny-1)/(p2-1), 1);
   }
 
-  threadGroupsFace[lIndex].calcGroups(d1,d2,1);
-  threadGroupsFace[lIndex].applyMap(map);
+  // Find elements that are on the boundary face 'lIndex'
+  IntVec map;
+  this->findBoundaryElms(map,lIndex);
+
+  fGrp.applyMap(map);
+}
+
+
+void ASMs3DLag::findBoundaryElms (IntVec& elms, int lIndex, int) const
+{
+  const int N1m = (nx-1)/(p1-1);
+  const int N2m = (ny-1)/(p2-1);
+  const int N3m = (nz-1)/(p3-1);
+
+  elms.clear();
+  switch (lIndex) {
+  case 1:
+  case 2:
+    elms.reserve(N2m*N3m);
+    for (int k = 0; k < N3m; ++k)
+      for (int j = 0; j < N2m; ++j)
+        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(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(i + j*N1m + (lIndex-5)*(N1m*N2m*(N3m-1)));
+  }
 }
 
 
diff --git a/src/ASM/ASMs3DLag.h b/src/ASM/ASMs3DLag.h
index 27bf9737..8428080d 100644
--- a/src/ASM/ASMs3DLag.h
+++ b/src/ASM/ASMs3DLag.h
@@ -128,17 +128,23 @@ public:
   //! \param[in] nSegSpan Number of visualization segments over each knot-span
   virtual bool getGridParameters(RealArray& prm, int dir, int nSegSpan) const;
 
-  //! \brief Find element for parameter, and optionally calculate local coordinates.
-  virtual int findElement(double u, double v, double w, double* xi = nullptr,
-                          double* eta = nullptr, double* zeta = nullptr) const;
-
 protected:
   //! \brief Assigned global coordinates for the given node.
   //! \param[in] inod 1-based node index local to current patch
   //! \param[in] Xnod Coordinates of the node
   void setCoord(size_t inod, const Vec3& Xnod);
 
+  //! \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;
+
 public:
+  //! \brief Find element for parameter, and optionally calculate local coordinates.
+  virtual int findElement(double u, double v, double w,
+                          double* xi = nullptr, double* eta = nullptr,
+                          double* zeta = nullptr) const;
+
   //! \brief Updates the nodal coordinates for this patch.
   //! \param[in] displ Incremental displacements to update the coordinates with
   virtual bool updateCoords(const Vector& displ);