Diverse småtterier...

src/ASM/ASMbase.C

@@ -23,6 +23,7 @@
 
 
 bool ASMbase::fixHomogeneousDirichlet = true;
+int  ASMbase::dbgElm = 0;
 
 
 /*!
@@ -1008,8 +1009,8 @@ bool ASMbase::extractNodalVec (const Vector& globRes, Vector& nodeVec,
 bool ASMbase::injectNodalVec (const Vector& nodeVec, Vector& globVec,
                               const std::vector<int>& madof, int basis) const
 {
-  if (madof.empty()) 
-  { 
+  if (madof.empty())
+  {
     std::cerr <<" *** ASMbase::injectNodalVec: Empty MADOF array."<< std::endl;
     return false;
   }

src/ASM/ASMbase.h

@@ -654,6 +654,8 @@ class ASMbase
 public:
   static bool fixHomogeneousDirichlet; //!< If \e true, pre-eliminate fixed DOFs
 
+  static int dbgElm; //!< One-based element index to print debugging info for
+
   size_t idx; //!< Index of this patch in the multi-patch model
 
 protected:

src/ASM/ASMs2D.C

@@ -688,16 +688,18 @@ bool ASMs2D::connectBasis (int edge, ASMs2D& neighbor, int nedge, bool revers,
   for (i = 0; i < n1; i++, node += i1)
   {
     int slave = slaveNodes[revers ? n1-i-1 : i];
-    if (coordCheck && !neighbor.getCoord(node).equal(this->getCoord(slave),xtol))
+    if (!coordCheck)
+      ASMbase::collapseNodes(neighbor,node,*this,slave);
+    else if (neighbor.getCoord(node).equal(this->getCoord(slave),xtol))
+      ASMbase::collapseNodes(neighbor,node,*this,slave);
+    else
     {
       std::cerr <<" *** ASMs2D::connectPatch: Non-matching nodes "
 		<< node <<": "<< neighbor.getCoord(node)
 		<<"\n                                          and "
 		<< slave <<": "<< this->getCoord(slave) << std::endl;
       return false;
     }
-    else
-      ASMbase::collapseNodes(neighbor,node,*this,slave);
   }
 
   return true;
@@ -1017,7 +1019,7 @@ size_t ASMs2D::constrainEdgeLocal (int dir, bool open, int dof, int code,
 void ASMs2D::constrainCorner (int I, int J, int dof, int code, char basis)
 {
   int node = this->getCorner(I, J, basis);
-  if (node > -1)
+  if (node > 0)
     this->prescribe(node,dof,code);
 }
 
@@ -1535,6 +1537,11 @@ bool ASMs2D::integrate (Integrand& integrand,
         fe.iel = MLGE[iel];
         if (fe.iel < 1) continue; // zero-area element
 
+#ifdef SP_DEBUG
+        if (dbgElm < 0 && 1+iel != -dbgElm)
+          continue; // Skipping all elements, except for -dbgElm
+#endif
+
         int i1 = p1 + iel % nel1;
         int i2 = p2 + iel / nel1;
 
@@ -1686,10 +1693,13 @@ bool ASMs2D::integrate (Integrand& integrand,
               utl::getGmat(Jac,dXidu,fe.G);
 
 #if SP_DEBUG > 4
-            std::cout <<"\niel, ip = "<< iel <<" "<< ip
-                      <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
-            if (!fe.d2NdX2.empty())
-              std::cout <<"d2NdX2 ="<< fe.d2NdX2;
+            if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
+            {
+              std::cout <<"\niel, ip = "<< iel <<" "<< ip
+                        <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
+              if (!fe.d2NdX2.empty())
+                std::cout <<"d2NdX2 ="<< fe.d2NdX2;
+            }
 #endif
 
             // Cartesian coordinates of current integration point
@@ -1714,6 +1724,10 @@ bool ASMs2D::integrate (Integrand& integrand,
           ok = false;
 
         A->destruct();
+
+#ifdef SP_DEBUG
+	if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm
+#endif
       }
     }
   }
@@ -1788,6 +1802,11 @@ bool ASMs2D::integrate (Integrand& integrand,
         fe.iel = MLGE[iel];
         if (fe.iel < 1) continue; // zero-area element
 
+#ifdef SP_DEBUG
+        if (dbgElm < 0 && iel != -dbgElm)
+          continue; // Skipping all elements, except for -dbgElm
+#endif
+
         int i1 = p1 + iel % nel1;
         int i2 = p2 + iel / nel1;
 
@@ -1864,8 +1883,9 @@ bool ASMs2D::integrate (Integrand& integrand,
             utl::getGmat(Jac,dXidu,fe.G);
 
 #if SP_DEBUG > 4
-          std::cout <<"\niel, jp = "<< iel <<" "<< jp
-                    <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
+          if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
+            std::cout <<"\niel, jp = "<< iel <<" "<< jp
+                      <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
 #endif
 
           // Cartesian coordinates of current integration point
@@ -1890,6 +1910,10 @@ bool ASMs2D::integrate (Integrand& integrand,
           ok = false;
 
         A->destruct();
+
+#ifdef SP_DEBUG
+	if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm
+#endif
       }
     }
   }
@@ -2134,6 +2158,11 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
       fe.iel = abs(MLGE[doXelms+iel-1]);
       if (fe.iel < 1) continue; // zero-area element
 
+#ifdef SP_DEBUG
+      if (dbgElm < 0 && iel != -dbgElm)
+        continue; // Skipping all elements, except for -dbgElm
+#endif
+
       // Skip elements that are not on current boundary edge
       bool skipMe = false;
       switch (edgeDir)
@@ -2220,6 +2249,11 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
       A->destruct();
 
       if (!ok) return false;
+
+#ifdef SP_DEBUG
+      if (dbgElm < 0 && iel == -dbgElm)
+        break; // Skipping all elements, except for -dbgElm
+#endif
     }
 
   return true;
@@ -2728,7 +2762,7 @@ bool ASMs2D::evaluate (const RealFunc* func, RealArray& vec,
 }
 
 
-int ASMs2D::getCorner(int I, int J, int basis) const
+int ASMs2D::getCorner (int I, int J, int basis) const
 {
   int n1, n2, node = 1;
   for (char i = 1; i <= basis; i++)

src/ASM/ASMs3D.C

@@ -756,16 +756,18 @@ bool ASMs3D::connectBasis (int face, ASMs3D& neighbor, int nface, int norient,
 	}
 
       int slave = slaveNodes[k][l];
-      if (coordCheck && !neighbor.getCoord(node).equal(this->getCoord(slave),xtol))
+      if (!coordCheck)
+        ASMbase::collapseNodes(neighbor,node,*this,slave);
+      else if (neighbor.getCoord(node).equal(this->getCoord(slave),xtol))
+        ASMbase::collapseNodes(neighbor,node,*this,slave);
+      else
       {
 	std::cerr <<" *** ASMs3D::connectPatch: Non-matching nodes "
 		  << node <<": "<< neighbor.getCoord(node)
 		  <<"\n                                          and "
 		  << slave <<": "<< this->getCoord(slave) << std::endl;
 	return false;
       }
-      else
-	ASMbase::collapseNodes(neighbor,node,*this,slave);
     }
 
   return true;
@@ -1110,7 +1112,7 @@ size_t ASMs3D::constrainFaceLocal (int dir, bool open, int dof, int code,
 }
 
 
-std::vector<int> ASMs3D::getEdge(int lEdge, bool open, int basis) const
+std::vector<int> ASMs3D::getEdge (int lEdge, bool open, int basis) const
 {
   std::vector<int> result;
   int n1, n2, n3, n, inc = 1;
@@ -1254,10 +1256,8 @@ void ASMs3D::constrainCorner (int I, int J, int K, int dof,
                               int code, char basis)
 {
   int node = this->getCorner(I, J, K, basis);
-  if (node < 1)
-    return;
-
-  this->prescribe(node,dof,code);
+  if (node > 0)
+    this->prescribe(node,dof,code);
 }
 
 
@@ -1785,6 +1785,11 @@ bool ASMs3D::integrate (Integrand& integrand,
         fe.iel = MLGE[iel];
         if (fe.iel < 1) continue; // zero-volume element
 
+#ifdef SP_DEBUG
+        if (dbgElm < 0 && 1+iel != -dbgElm)
+          continue; // Skipping all elements, except for -dbgElm
+#endif
+
         int i1 = p1 + iel % nel1;
         int i2 = p2 + (iel / nel1) % nel2;
         int i3 = p3 + iel / (nel1*nel2);
@@ -1946,8 +1951,13 @@ bool ASMs3D::integrate (Integrand& integrand,
                 utl::getGmat(Jac,dXidu,fe.G);
 
 #if SP_DEBUG > 4
-              std::cout <<"\niel, ip = "<< iel <<" "<< ip
-                        <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
+              if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
+              {
+                std::cout <<"\niel, ip = "<< iel <<" "<< ip
+                          <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
+                if (!fe.d2NdX2.empty())
+                  std::cout <<"d2NdX2 ="<< fe.d2NdX2;
+              }
 #endif
 
               // Cartesian coordinates of current integration point
@@ -1972,6 +1982,10 @@ bool ASMs3D::integrate (Integrand& integrand,
           ok = false;
 
         A->destruct();
+
+#ifdef SP_DEBUG
+        if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm
+#endif
       }
     }
   }
@@ -2050,6 +2064,11 @@ bool ASMs3D::integrate (Integrand& integrand,
         fe.iel = MLGE[iel];
         if (fe.iel < 1) continue; // zero-volume element
 
+#ifdef SP_DEBUG
+        if (dbgElm < 0 && iel != -dbgElm)
+          continue; // Skipping all elements, except for -dbgElm
+#endif
+
         int i1 = p1 + iel % nel1;
         int i2 = p2 + (iel / nel1) % nel2;
         int i3 = p3 + iel / (nel1*nel2);
@@ -2129,8 +2148,9 @@ bool ASMs3D::integrate (Integrand& integrand,
             utl::getGmat(Jac,dXidu,fe.G);
 
 #if SP_DEBUG > 4
-          std::cout <<"\niel, jp = "<< iel <<" "<< jp
-                    <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
+          if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
+            std::cout <<"\niel, jp = "<< iel <<" "<< jp
+                      <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
 #endif
 
           // Cartesian coordinates of current integration point
@@ -2155,6 +2175,10 @@ bool ASMs3D::integrate (Integrand& integrand,
           ok = false;
 
         A->destruct();
+
+#ifdef SP_DEBUG
+        if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm
+#endif
       }
     }
   }
@@ -2278,6 +2302,11 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
         fe.iel = abs(MLGE[doXelms+iel]);
         if (fe.iel < 1) continue; // zero-volume element
 
+#ifdef SP_DEBUG
+        if (dbgElm < 0 && iel != -dbgElm)
+          continue; // Skipping all elements, except for -dbgElm
+#endif
+
         int i1 = p1 + iel % nel1;
         int i2 = p2 + (iel / nel1) % nel2;
         int i3 = p3 + iel / (nel1*nel2);
@@ -2395,6 +2424,11 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
           ok = false;
 
         A->destruct();
+
+#ifdef SP_DEBUG
+        if (dbgElm < 0 && iel == -dbgElm)
+          break; // Skipping all elements, except for -dbgElm
+#endif
       }
     }
   }
@@ -3190,7 +3224,7 @@ bool ASMs3D::evaluate (const RealFunc* func, RealArray& vec,
 }
 
 
-int ASMs3D::getCorner(int I, int J, int K, int basis) const
+int ASMs3D::getCorner (int I, int J, int K, int basis) const
 {
   int n1, n2, n3;
   int node = this->findStartNode(n1,n2,n3,basis);

src/ASM/ASMs3Dmx.C

@@ -448,15 +448,14 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
   // Evaluate basis function derivatives at all integration points
   std::vector<std::vector<Go::BasisDerivs>> splinex(m_basis.size());
   std::vector<std::vector<Go::BasisDerivs2>> splinex2(m_basis.size());
-  if (use2ndDer) {
+  if (use2ndDer)
 #pragma omp parallel for schedule(static)
     for (size_t i=0;i<m_basis.size();++i)
       m_basis[i]->computeBasisGrid(gpar[0],gpar[1],gpar[2],splinex2[i]);
-  } else {
+  else
 #pragma omp parallel for schedule(static)
     for (size_t i=0;i<m_basis.size();++i)
       m_basis[i]->computeBasisGrid(gpar[0],gpar[1],gpar[2],splinex[i]);
-  }
 
   std::vector<size_t> elem_sizes;
   for (auto& it : m_basis)
@@ -513,7 +512,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
 
         if (useElmVtx)
           this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
- 
+
         if (integrand.getIntegrandType() & Integrand::G_MATRIX)
         {
           // Element size in parametric space
@@ -553,13 +552,12 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
               fe.w = gpar[2](k+1,i3-p3+1);
 
               // Fetch basis function derivatives at current integration point
-              if (use2ndDer) {
+              if (use2ndDer)
                 for (size_t b = 0; b < m_basis.size(); ++b)
                   SplineUtils::extractBasis(splinex2[b][ip],fe.basis(b+1),dNxdu[b], d2Nxdu2[b]);
-              } else {
+              else
                 for (size_t b = 0; b < m_basis.size(); ++b)
                   SplineUtils::extractBasis(splinex[b][ip],fe.basis(b+1),dNxdu[b]);
-              }
 
               // Compute Jacobian inverse of the coordinate mapping and
               // basis function derivatives w.r.t. Cartesian coordinates