diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C index 2b226af9..07521629 100644 --- a/src/ASM/ASMbase.C +++ b/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& madof, int basis) const { - if (madof.empty()) - { + if (madof.empty()) + { std::cerr <<" *** ASMbase::injectNodalVec: Empty MADOF array."<< std::endl; return false; } diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index 07736746..abe5ae5a 100644 --- a/src/ASM/ASMbase.h +++ b/src/ASM/ASMbase.h @@ -654,6 +654,8 @@ protected: 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: diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C index 1ac5dfdb..3bc79520 100644 --- a/src/ASM/ASMs2D.C +++ b/src/ASM/ASMs2D.C @@ -688,7 +688,11 @@ 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) @@ -696,8 +700,6 @@ bool ASMs2D::connectBasis (int edge, ASMs2D& neighbor, int nedge, bool revers, << 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++) diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C index 8534e89c..c24f89ad 100644 --- a/src/ASM/ASMs3D.C +++ b/src/ASM/ASMs3D.C @@ -756,7 +756,11 @@ 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) @@ -764,8 +768,6 @@ bool ASMs3D::connectBasis (int face, ASMs3D& neighbor, int nface, int norient, << 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 ASMs3D::getEdge(int lEdge, bool open, int basis) const +std::vector ASMs3D::getEdge (int lEdge, bool open, int basis) const { std::vector 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); diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index 465c87c0..97365fec 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/src/ASM/ASMs3Dmx.C @@ -448,15 +448,14 @@ bool ASMs3Dmx::integrate (Integrand& integrand, // Evaluate basis function derivatives at all integration points std::vector> splinex(m_basis.size()); std::vector> splinex2(m_basis.size()); - if (use2ndDer) { + if (use2ndDer) #pragma omp parallel for schedule(static) for (size_t i=0;icomputeBasisGrid(gpar[0],gpar[1],gpar[2],splinex2[i]); - } else { + else #pragma omp parallel for schedule(static) for (size_t i=0;icomputeBasisGrid(gpar[0],gpar[1],gpar[2],splinex[i]); - } std::vector 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 diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index 510ad9ee..d5aebc0b 100644 --- a/src/ASM/LR/ASMu2D.C +++ b/src/ASM/LR/ASMu2D.C @@ -637,13 +637,13 @@ void ASMu2D::constrainEdge (int dir, bool open, int dof, int code, char basis) // get all elements and functions on this edge std::vector edgeFunctions; std::vector edgeElements; - lr->getEdgeFunctions(edgeFunctions, edge); - lr->getEdgeElements( edgeElements, edge); + lr->getEdgeFunctions(edgeFunctions,edge); + lr->getEdgeElements (edgeElements ,edge); // find the corners since these are not to be included in the L2-fitting // of the inhomogenuous dirichlet boundaries; corners are interpolatory. - // Optimization note: loop over the "edge"-container to manually pick up the - // end nodes. LRspine::getEdgeFunctions() does a global search. + // Optimization note: loop over the "edge"-container to manually pick up + // the end nodes. LRspine::getEdgeFunctions() does a global search. std::vector c1, c2; switch (edge) { @@ -676,7 +676,7 @@ void ASMu2D::constrainEdge (int dir, bool open, int dof, int code, char basis) int bcode = abs(code); int j = 0; - for (auto b : edgeFunctions ) + for (auto b : edgeFunctions) { de.MLGN[j++] = b->getId(); // skip corners for open boundaries @@ -993,6 +993,11 @@ bool ASMu2D::integrate (Integrand& integrand, std::vector::iterator el = lrspline->elementBegin(); for (int iel = 1; el != lrspline->elementEnd(); el++, iel++) { +#ifdef SP_DEBUG + if (dbgElm < 0 && iel != -dbgElm) + continue; // Skipping all elements, except for -dbgElm +#endif + FiniteElement fe(MNPC[iel-1].size()); fe.iel = MLGE[iel-1]; @@ -1091,13 +1096,18 @@ bool ASMu2D::integrate (Integrand& integrand, lrspline->computeBasis(fe.u,fe.v,spline, iel-1); SplineUtils::extractBasis(spline,fe.N,dNdu); #if SP_DEBUG > 4 - std::cout <<"\nBasis functions at a integration point " - <<" : (u,v) = "<< spline.param[0] <<" "<< spline.param[1] - <<" left_idx = "<< spline.left_idx[0] <<" "<< spline.left_idx[1]; - for (size_t ii = 0; ii < spline.basisValues.size(); ii++) - std::cout <<'\n'<< 1+ii <<'\t' << spline.basisValues[ii] <<'\t' - << spline.basisDerivs_u[ii] <<'\t'<< spline.basisDerivs_v[ii]; - std::cout << std::endl; + if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) + { + std::cout <<"\nBasis functions at a integration point " + <<" : (u,v) = "<< spline.param[0] <<" "<< spline.param[1] + <<" left_idx = "<< spline.left_idx[0] + <<" "<< spline.left_idx[1]; + for (size_t ii = 0; ii < spline.basisValues.size(); ii++) + std::cout <<'\n'<< 1+ii <<'\t' << spline.basisValues[ii] <<'\t' + << spline.basisDerivs_u[ii] <<'\t' + << spline.basisDerivs_v[ii]; + std::cout << std::endl; + } #endif } @@ -1111,7 +1121,8 @@ bool ASMu2D::integrate (Integrand& integrand, return false; #if SP_DEBUG > 4 - std::cout <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; + if (iel == dbgElm || iel == -dbgElm || dbgElm == 0) + std::cout <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX; #endif // Cartesian coordinates of current integration point @@ -1134,6 +1145,10 @@ bool ASMu2D::integrate (Integrand& integrand, return false; A->destruct(); + +#ifdef SP_DEBUG + if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm +#endif } return true; @@ -1178,6 +1193,11 @@ bool ASMu2D::integrate (Integrand& integrand, double dA = this->getParametricArea(iel); if (dA < 0.0) return false; // topology error (probably logic error) +#ifdef SP_DEBUG + if (dbgElm < 0 && iel != -dbgElm) + continue; // Skipping all elements, except for -dbgElm +#endif + // Set up control point (nodal) coordinates for current element if (!this->getElementCoordinates(Xnod,iel)) return false; @@ -1230,8 +1250,9 @@ bool ASMu2D::integrate (Integrand& integrand, return false; #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; #endif // Cartesian coordinates of current integration point @@ -1254,6 +1275,10 @@ bool ASMu2D::integrate (Integrand& integrand, return false; A->destruct(); + +#ifdef SP_DEBUG + if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm +#endif } return true; @@ -1306,6 +1331,11 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex, std::vector::iterator el = lrspline->elementBegin(); for (int iel = 1; el != lrspline->elementEnd(); el++, iel++) { +#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) @@ -1383,6 +1413,11 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex, return false; A->destruct(); + +#ifdef SP_DEBUG + if (dbgElm < 0 && iel == -dbgElm) + break; // Skipping all elements, except for -dbgElm +#endif } return true; @@ -1870,9 +1905,9 @@ bool ASMu2D::updateDirichlet (const std::map& func, (*mit)->setSlaveCoeff(edgeControlmatrix[dof-1][nit->first]); else //scalar condition (*mit)->setSlaveCoeff(edgeControlpoints[nit->first]); - #if SP_DEBUG > 1 +#if SP_DEBUG > 1 std::cout <<"Updated constraint: "<< **mit; - #endif +#endif } } }