diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index c0c18425..7cc16492 100644 --- a/src/ASM/LR/ASMu2D.C +++ b/src/ASM/LR/ASMu2D.C @@ -1058,15 +1058,19 @@ bool ASMu2D::integrate (Integrand& integrand, bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES; bool use3rdDer = integrand.getIntegrandType() & Integrand::THIRD_DERIVATIVES; + const int p1 = lrspline->order(0); + const int p2 = lrspline->order(1); + // Get Gaussian quadrature points and weights - const double* xg = GaussQuadrature::getCoord(nGauss); - const double* wg = GaussQuadrature::getWeight(nGauss); + const int nGP = this->getNoGaussPt(p1 > p2 ? p1 : p2); + const double* xg = GaussQuadrature::getCoord(nGP); + const double* wg = GaussQuadrature::getWeight(nGP); if (!xg || !wg) return false; // Get the reduced integration quadrature points, if needed const double* xr = nullptr; const double* wr = nullptr; - int nRed = integrand.getReducedIntegration(nGauss); + int nRed = integrand.getReducedIntegration(nGP); if (nRed > 0) { xr = GaussQuadrature::getCoord(nRed); @@ -1074,7 +1078,7 @@ bool ASMu2D::integrate (Integrand& integrand, if (!xr || !wr) return false; } else if (nRed < 0) - nRed = nGauss; // The integrand needs to know nGauss + nRed = nGP; // The integrand needs to know nGauss // Evaluate basis function values and derivatives at all integration points. // We do this before the integration point loop to exploit multi-threading @@ -1085,11 +1089,11 @@ bool ASMu2D::integrate (Integrand& integrand, std::vector spline3; if (use3rdDer) - spline3.resize(nel*nGauss*nGauss); + spline3.resize(nel*nGP*nGP); else if (use2ndDer) - spline2.resize(nel*nGauss*nGauss); + spline2.resize(nel*nGP*nGP); else - spline1.resize(nel*nGauss*nGauss); + spline1.resize(nel*nGP*nGP); if (xr) splineRed.resize(nel*nRed*nRed); @@ -1097,10 +1101,10 @@ bool ASMu2D::integrate (Integrand& integrand, for (iel = jp = rp = 0; iel < nel; iel++) { RealArray u, v; - this->getGaussPointParameters(u,0,nGauss,1+iel,xg); - this->getGaussPointParameters(v,1,nGauss,1+iel,xg); - for (int j = 0; j < nGauss; j++) - for (int i = 0; i < nGauss; i++, jp++) + this->getGaussPointParameters(u,0,nGP,1+iel,xg); + this->getGaussPointParameters(v,1,nGP,1+iel,xg); + for (int j = 0; j < nGP; j++) + for (int i = 0; i < nGP; i++, jp++) if (use3rdDer) this->computeBasis(u[i],v[j],spline3[jp],iel); else if (use2ndDer) @@ -1141,8 +1145,8 @@ bool ASMu2D::integrate (Integrand& integrand, FiniteElement fe; fe.iel = MLGE[iel-1]; - fe.p = lrspline->order(0) - 1; - fe.q = lrspline->order(1) - 1; + fe.p = p1 - 1; + fe.q = p2 - 1; Matrix dNdu, Xnod, Jac; Matrix3D d2Ndu2, Hess; Matrix4D d3Ndu3; @@ -1169,7 +1173,7 @@ bool ASMu2D::integrate (Integrand& integrand, std::array gpar, redpar; for (int d = 0; d < 2; d++) { - this->getGaussPointParameters(gpar[d],d,nGauss,iel,xg); + this->getGaussPointParameters(gpar[d],d,nGP,iel,xg); if (xr) this->getGaussPointParameters(redpar[d],d,nRed,iel,xr); } @@ -1251,11 +1255,11 @@ bool ASMu2D::integrate (Integrand& integrand, // --- Integration loop over all Gauss points in each direction ---------- - int jp = (iel-1)*nGauss*nGauss; + int jp = (iel-1)*nGP*nGP; fe.iGP = firstIp + jp; // Global integration point counter - for (int j = 0; j < nGauss; j++) - for (int i = 0; i < nGauss; i++, fe.iGP++) + for (int j = 0; j < nGP; j++) + for (int i = 0; i < nGP; i++, fe.iGP++) { // Local element coordinates of current integration point fe.xi = xg[i]; @@ -1357,7 +1361,7 @@ bool ASMu2D::integrate (Integrand& integrand, { if (!lrspline) return true; // silently ignore empty patches - if (integrand.getReducedIntegration(nGauss) != 0) + if (integrand.getReducedIntegration(2) != 0) { std::cerr <<" *** ASMu2D::integrate(Integrand&,GlobalIntegral&," <<"const TimeDomain&,const Real3DMat&): Available for standard" @@ -1546,8 +1550,12 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex, PROFILE2("ASMu2D::integrate(B)"); + const int p1 = lrspline->order(0); + const int p2 = lrspline->order(1); + // Get Gaussian quadrature points and weights - int nGP = integrand.getBouIntegrationPoints(nGauss); + int nG1 = this->getNoGaussPt(lIndex%10 < 3 ? p1 : p2, true); + int nGP = integrand.getBouIntegrationPoints(nG1); const double* xg = GaussQuadrature::getCoord(nGP); const double* wg = GaussQuadrature::getWeight(nGP); if (!xg || !wg) return false; @@ -1578,8 +1586,8 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex, size_t firstp = iit == firstBp.end() ? 0 : iit->second; FiniteElement fe; - fe.p = lrspline->order(0) - 1; - fe.q = lrspline->order(1) - 1; + fe.p = p1 - 1; + fe.q = p2 - 1; fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0; double param[3] = { 0.0, 0.0, 0.0 }; diff --git a/src/ASM/LR/ASMu2Drecovery.C b/src/ASM/LR/ASMu2Drecovery.C index 727da8f4..6107889f 100644 --- a/src/ASM/LR/ASMu2Drecovery.C +++ b/src/ASM/LR/ASMu2Drecovery.C @@ -106,13 +106,14 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B, const int p1 = projBasis->order(0); const int p2 = projBasis->order(1); + const int pm = p1 > p2 ? p1 : p2; // Get Gaussian quadrature points - const int ng1 = continuous ? nGauss : p1 - 1; - const int ng2 = continuous ? nGauss : p2 - 1; + const int ng1 = continuous ? this->getNoGaussPt(pm,true) : p1 - 1; + const int ng2 = continuous ? ng1 : p2 - 1; const double* xg = GaussQuadrature::getCoord(ng1); const double* yg = GaussQuadrature::getCoord(ng2); - const double* wg = continuous ? GaussQuadrature::getWeight(nGauss) : nullptr; + const double* wg = continuous ? GaussQuadrature::getWeight(ng1) : nullptr; if (!xg || !yg) return false; if (continuous && !wg) return false; @@ -463,11 +464,6 @@ bool ASMu2D::edgeL2projection (const DirichletEdge& edge, StdVector B(n*m); A.resize(n,n); - // Get Gaussian quadrature points and weights - const double* xg = GaussQuadrature::getCoord(nGauss); - const double* wg = GaussQuadrature::getWeight(nGauss); - if (!xg || !wg) return false; - // find the normal and tangent direction for the edge int edgeDir, t1, t2; switch (edge.edg) @@ -479,16 +475,22 @@ bool ASMu2D::edgeL2projection (const DirichletEdge& edge, default: return false; } + // Get Gaussian quadrature points and weights + const int nGP = this->getNoGaussPt(lrspline->order(t2),true); + const double* xg = GaussQuadrature::getCoord(nGP); + const double* wg = GaussQuadrature::getWeight(nGP); + if (!xg || !wg) return false; + std::array gpar; for (int d = 0; d < 2; d++) if (-1-d == edgeDir) { - gpar[d].resize(nGauss); + gpar[d].resize(nGP); gpar[d].fill(d == 0 ? lrspline->startparam(0) : lrspline->startparam(1)); } else if (1+d == edgeDir) { - gpar[d].resize(nGauss); + gpar[d].resize(nGP); gpar[d].fill(d == 0 ? lrspline->endparam(0) : lrspline->endparam(1)); } @@ -511,11 +513,11 @@ bool ASMu2D::edgeL2projection (const DirichletEdge& edge, if (!this->getElementCoordinates(Xnod,iel)) return false; // Get integration gauss points over this element - this->getGaussPointParameters(gpar[t2-1],t2-1,nGauss,iel,xg); + this->getGaussPointParameters(gpar[t2-1],t2-1,nGP,iel,xg); // --- Integration loop over all Gauss points along the edge --------------- - for (int j = 0; j < nGauss; j++) + for (int j = 0; j < nGP; j++) { // Parameter values of current integration point double u = param[0] = gpar[0][j]; diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C index 7f82f228..5fdf2731 100644 --- a/src/ASM/LR/ASMu3D.C +++ b/src/ASM/LR/ASMu3D.C @@ -724,7 +724,7 @@ double ASMu3D::getElementCorners (int iEl, Vec3Vec& XC) const void ASMu3D::evaluateBasis (int iel, int basis, double u, double v, double w, Vector& N, Matrix& dNdu) const { - PROFILE2("Spline evaluation"); + PROFILE3("ASMu3D::evalBasis(1)"); std::vector result; this->getBasis(basis)->computeBasis(u, v, w, result, 1, iel); @@ -749,7 +749,7 @@ void ASMu3D::evaluateBasis (int iel, FiniteElement& fe, Matrix& dNdu, void ASMu3D::evaluateBasis (FiniteElement& fe, Matrix& dNdu, const Matrix& C, const Matrix& B, int basis) const { - PROFILE2("BeSpline evaluation"); + PROFILE3("ASMu3D::evalBasis(BE)"); Matrix N = C*B; dNdu.resize(N.rows(),3); @@ -763,7 +763,7 @@ void ASMu3D::evaluateBasis (int iel, FiniteElement& fe, Matrix& dNdu, Matrix3D& d2Ndu2, int basis) const { - PROFILE2("Spline evaluation"); + PROFILE3("ASMu3D::evalBasis(2)"); std::vector result; this->getBasis(basis)->computeBasis(fe.u, fe.v, fe.w, result, 2, iel); @@ -790,7 +790,7 @@ void ASMu3D::evaluateBasis (int iel, FiniteElement& fe, void ASMu3D::evaluateBasis (int iel, FiniteElement& fe, int derivs, int basis) const { - PROFILE2("Spline evaluation"); + PROFILE3("ASMu2D::evalBasis"); std::vector result; this->getBasis(basis)->computeBasis(fe.u, fe.v, fe.w, result, derivs, iel); @@ -831,27 +831,30 @@ bool ASMu3D::integrate (Integrand& integrand, PROFILE2("ASMu3D::integrate(I)"); - // Get Gaussian quadrature points and weights - const double* xg = GaussQuadrature::getCoord(nGauss); - const double* wg = GaussQuadrature::getWeight(nGauss); - if (!xg || !wg) return false; - - // Evaluate all gauss points on the bezier patch (-1, 1) int p1 = lrspline->order(0); int p2 = lrspline->order(1); int p3 = lrspline->order(2); + int pm = std::max(std::max(p1,p2),p3); + + // Get Gaussian quadrature points and weights + int nGP = this->getNoGaussPt(pm); + const double* xg = GaussQuadrature::getCoord(nGP); + const double* wg = GaussQuadrature::getWeight(nGP); + if (!xg || !wg) return false; + + // Evaluate all gauss points on the bezier patch (-1, 1) double u[2*p1], v[2*p2], w[2*p3]; Go::BsplineBasis basis1 = getBezierBasis(p1); Go::BsplineBasis basis2 = getBezierBasis(p2); Go::BsplineBasis basis3 = getBezierBasis(p3); - Matrix BN( p1*p2*p3, nGauss*nGauss*nGauss), rBN; - Matrix BdNdu(p1*p2*p3, nGauss*nGauss*nGauss), rBdNdu; - Matrix BdNdv(p1*p2*p3, nGauss*nGauss*nGauss), rBdNdv; - Matrix BdNdw(p1*p2*p3, nGauss*nGauss*nGauss), rBdNdw; + Matrix BN( p1*p2*p3, nGP*nGP*nGP), rBN; + Matrix BdNdu(p1*p2*p3, nGP*nGP*nGP), rBdNdu; + Matrix BdNdv(p1*p2*p3, nGP*nGP*nGP), rBdNdv; + Matrix BdNdw(p1*p2*p3, nGP*nGP*nGP), rBdNdw; int ig = 1; // gauss point iterator - for (int zeta = 0; zeta < nGauss; zeta++) - for (int eta = 0; eta < nGauss; eta++) - for (int xi = 0; xi < nGauss; xi++, ig++) { + for (int zeta = 0; zeta < nGP; zeta++) + for (int eta = 0; eta < nGP; eta++) + for (int xi = 0; xi < nGP; xi++, ig++) { basis1.computeBasisValues(xg[xi], u, 1); basis2.computeBasisValues(xg[eta], v, 1); basis3.computeBasisValues(xg[zeta], w, 1); @@ -869,7 +872,7 @@ bool ASMu3D::integrate (Integrand& integrand, // Get the reduced integration quadrature points, if needed const double* xr = nullptr; const double* wr = nullptr; - int nRed = integrand.getReducedIntegration(nGauss); + int nRed = integrand.getReducedIntegration(nGP); if (nRed > 0) { xr = GaussQuadrature::getCoord(nRed); @@ -899,7 +902,7 @@ bool ASMu3D::integrate (Integrand& integrand, } } else if (nRed < 0) - nRed = nGauss; // The integrand needs to know nGauss + nRed = nGP; // The integrand needs to know nGauss ThreadGroups oneGroup; if (glInt.threadSafe()) oneGroup.oneGroup(nel); @@ -958,7 +961,7 @@ bool ASMu3D::integrate (Integrand& integrand, std::array gpar, redpar; for (int d = 0; d < 3; d++) { - this->getGaussPointParameters(gpar[d],d,nGauss,iel,xg); + this->getGaussPointParameters(gpar[d],d,nGP,iel,xg); if (xr) this->getGaussPointParameters(redpar[d],d,nRed,iel,xr); } @@ -1040,12 +1043,12 @@ bool ASMu3D::integrate (Integrand& integrand, // --- Integration loop over all Gauss points in each direction ---------- int ig = 1; - int jp = (iel-1)*nGauss*nGauss*nGauss; + int jp = (iel-1)*nGP*nGP*nGP; fe.iGP = firstIp + jp; // Global integration point counter - for (int k = 0; k < nGauss; k++) - for (int j = 0; j < nGauss; j++) - for (int i = 0; i < nGauss; i++, fe.iGP++, ig++) + for (int k = 0; k < nGP; k++) + for (int j = 0; j < nGP; j++) + for (int i = 0; i < nGP; i++, fe.iGP++, ig++) { // Local element coordinates of current integration point fe.xi = xg[i]; @@ -1150,18 +1153,21 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex, PROFILE2("ASMu3D::integrate(B)"); - // Get Gaussian quadrature points and weights - int nGP = integrand.getBouIntegrationPoints(nGauss); - const double* xg = GaussQuadrature::getCoord(nGP); - const double* wg = GaussQuadrature::getWeight(nGP); - if (!xg || !wg) return false; - // Find the parametric direction of the face normal {-3,-2,-1, 1, 2, 3} const int faceDir = (lIndex%10+1)/((lIndex%2) ? -2 : 2); const int t1 = 1 + abs(faceDir)%3; // first tangent direction const int t2 = 1 + t1%3; // second tangent direction + // Get Gaussian quadrature points and weights + // Use the largest polynomial order of the two tangent directions + const int pmax = std::max(lrspline->order(t1),lrspline->order(t2)); + const int nG1 = this->getNoGaussPt(pmax,true); + const int nGP = integrand.getBouIntegrationPoints(nG1); + const double* xg = GaussQuadrature::getCoord(nGP); + const double* wg = GaussQuadrature::getWeight(nGP); + if (!xg || !wg) return false; + // Extract the Neumann order flag (1 or higher) for the integrand integrand.setNeumannOrder(1 + lIndex/10); @@ -1230,7 +1236,7 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex, double dXidu[3]; // Get element face area in the parameter space - double dA = this->getParametricArea(iEl+1,abs(faceDir)); + double dA = 0.25*this->getParametricArea(iEl+1,abs(faceDir)); if (dA < 0.0) // topology error (probably logic error) { ok = false; @@ -1311,12 +1317,17 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex, if (integrand.getIntegrandType() & Integrand::G_MATRIX) utl::getGmat(Jac,dXidu,fe.G); +#if SP_DEBUG > 4 + if (iEl+1 == dbgElm || iEl+1 == -dbgElm || dbgElm == 0) + std::cout <<"\n"<< fe; +#endif + // Cartesian coordinates of current integration point X.assign(Xnod * fe.N); X.t = time.t; // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.25*dA*wg[i]*wg[j]; + fe.detJxW *= dA*wg[i]*wg[j]; if (!integrand.evalBou(*A,fe,time,X,normal)) ok = false; } @@ -1347,177 +1358,6 @@ bool ASMu3D::integrateEdge (Integrand& integrand, int lEdge, { std::cerr << "ASMu3D::integrateEdge(...) is not properly implemented yet :(" << std::endl; exit(776654); -#if 0 - if (!lrspline) return true; // silently ignore empty patches - - PROFILE2("ASMu3D::integrate(E)"); - - // Get Gaussian quadrature points and weights - const double* xg = GaussQuadrature::getCoord(nGauss); - const double* wg = GaussQuadrature::getWeight(nGauss); - if (!xg || !wg) return false; - - // Compute parameter values of the Gauss points along the whole patch edge - std::array gpar; - for (int d = 0; d < 3; d++) - if (lEdge < d*4+1 || lEdge >= d*4+5) - { - gpar[d].resize(1,1); - if (lEdge%4 == 1) - gpar[d].fill(lrspline->startparam(d)); - else if (lEdge%4 == 0) - gpar[d].fill(lrspline->endparam(d)); - else if (lEdge == 6 || lEdge == 10) - gpar[d].fill(d == 0 ? lrspline->endparam(d) : lrspline->startparam(d)); - else if (lEdge == 2 || lEdge == 11) - gpar[d].fill(d == 1 ? lrspline->endparam(d) : lrspline->startparam(d)); - else if (lEdge == 3 || lEdge == 7) - gpar[d].fill(d == 2 ? lrspline->endparam(d) : lrspline->startparam(d)); - } - else - { - int pm1 = lrspline->order(d) - 1; - RealArray::const_iterator uit = lrspline->basis(d).begin() + pm1; - double ucurr, uprev = *(uit++); - int nCol = lrspline->numCoefs(d) - pm1; - gpar[d].resize(nGauss,nCol); - for (int j = 1; j <= nCol; uit++, j++) - { - ucurr = *uit; - for (int i = 1; i <= nGauss; i++) - gpar[d](i,j) = 0.5*((ucurr-uprev)*xg[i-1] + ucurr+uprev); - uprev = ucurr; - } - } - - // Evaluate basis function derivatives at all integration points - std::vector spline; - { - PROFILE2("Spline evaluation"); - lrspline->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline); - } - - const int n1 = lrspline->numCoefs(0); - const int n2 = lrspline->numCoefs(1); - const int n3 = lrspline->numCoefs(2); - - const int p1 = lrspline->order(0); - const int p2 = lrspline->order(1); - const int p3 = lrspline->order(2); - - std::map::const_iterator iit = firstBp.find(lEdge); - size_t firstp = iit == firstBp.end() ? 0 : iit->second; - - FiniteElement fe(p1*p2*p3); - fe.u = gpar[0](1,1); - fe.v = gpar[1](1,1); - fe.w = gpar[2](1,1); - if (gpar[0].size() == 1) fe.xi = fe.u == lrspline->startparam(0) ? -1.0 : 1.0; - if (gpar[1].size() == 1) fe.eta = fe.v == lrspline->startparam(1) ? -1.0 : 1.0; - if (gpar[2].size() == 1) fe.zeta = fe.w == lrspline->startparam(2) ? -1.0 : 1.0; - - Matrix dNdu, Xnod, Jac; - Vec4 X; - Vec3 tang; - - - // === Assembly loop over all elements on the patch edge ===================== - - int iel = 1; - for (int i3 = p3; i3 <= n3; i3++) - for (int i2 = p2; i2 <= n2; i2++) - for (int i1 = p1; i1 <= n1; i1++, iel++) - { - fe.iel = MLGE[iel-1]; - if (fe.iel < 1) continue; // zero-volume element - - // Skip elements that are not on current boundary edge - bool skipMe = false; - switch (lEdge) - { - case 1: if (i2 > p2 || i3 > p3) skipMe = true; break; - case 2: if (i2 < n2 || i3 > p3) skipMe = true; break; - case 3: if (i2 > p2 || i3 < n3) skipMe = true; break; - case 4: if (i2 < n2 || i3 < n3) skipMe = true; break; - case 5: if (i1 > p1 || i3 > p3) skipMe = true; break; - case 6: if (i1 < n1 || i3 > p3) skipMe = true; break; - case 7: if (i1 > p1 || i3 < n3) skipMe = true; break; - case 8: if (i1 < n1 || i3 < n3) skipMe = true; break; - case 9: if (i1 > p1 || i2 > p2) skipMe = true; break; - case 10: if (i1 < n1 || i2 > p2) skipMe = true; break; - case 11: if (i1 > p1 || i2 < n2) skipMe = true; break; - case 12: if (i1 < n1 || i2 < n2) skipMe = true; break; - } - if (skipMe) continue; - - // Get element edge length in the parameter space - double dS = 0.0; - int ip = MNPC[iel-1].back(); -#ifdef INDEX_CHECK - if (ip < 0 || (size_t)ip >= nodeInd.size()) return false; -#endif - if (lEdge < 5) - { - dS = lrspline->knotSpan(0,nodeInd[ip].I); - ip = (i1-p1)*nGauss; - } - else if (lEdge < 9) - { - dS = lrspline->knotSpan(1,nodeInd[ip].J); - ip = (i2-p2)*nGauss; - } - else if (lEdge < 13) - { - dS = lrspline->knotSpan(2,nodeInd[ip].K); - ip = (i3-p3)*nGauss; - } - - // Set up control point coordinates for current element - if (!this->getElementCoordinates(Xnod,iel)) return false; - - // Initialize element quantities - LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel,true); - bool ok = integrand.initElementBou(MNPC[iel-1],*A); - - - // --- Integration loop over all Gauss points along the edge ----------------- - - fe.iGP = firstp + ip; // Global integration point counter - - for (int i = 0; i < nGauss && ok; i++, ip++, fe.iGP++) - { - // Parameter values of current integration point - if (gpar[0].size() > 1) fe.u = gpar[0](i+1,i1-p1+1); - if (gpar[1].size() > 1) fe.v = gpar[1](i+1,i2-p2+1); - if (gpar[2].size() > 1) fe.w = gpar[2](i+1,i3-p3+1); - - // Fetch basis function derivatives at current integration point - this->evaluateBasis(iel, fe, dNdu); - - // Compute basis function derivatives and the edge tang - fe.detJxW = utl::Jacobian(Jac,tang,fe.dNdX,Xnod,dNdu,1+(lEdge-1)/4); - if (fe.detJxW == 0.0) continue; // skip singular points - - // Cartesian coordinates of current integration point - X = Xnod * fe.N; - X.t = time.t; - - // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.5*dS*wg[i]; - ok = integrand.evalBou(*A,fe,time,X,tang); - } - - // Assembly of global system integral - if (ok && !glInt.assemble(A->ref(),fe.iel)) - ok = false; - - A->destruct(); - - if (!ok) return false; - } - - return true; -#endif return false; } @@ -1678,6 +1518,8 @@ bool ASMu3D::evalSolution (Matrix& sField, const Vector& locSol, bool ASMu3D::evalSolution (Matrix& sField, const Vector& locSol, const RealArray* gpar, bool, int deriv, int) const { + PROFILE2("ASMu3D::evalSol(P)"); + size_t nComp = locSol.size() / this->getNoNodes(); if (nComp*this->getNoNodes() != locSol.size()) return false; @@ -1693,6 +1535,7 @@ bool ASMu3D::evalSolution (Matrix& sField, const Vector& locSol, Go::BasisPts spline0; Go::BasisDerivs spline1; Go::BasisDerivs2 spline2; + int lel = -1; // Evaluate the primary solution field at each point sField.resize(nComp,nPoints); @@ -1708,9 +1551,12 @@ bool ASMu3D::evalSolution (Matrix& sField, const Vector& locSol, } utl::gather(MNPC[iel],nComp,locSol,eSol); - // Set up control point (nodal) coordinates for current element - if (deriv > 0 && !this->getElementCoordinates(Xnod,iel+1)) - return false; + if (iel != lel && deriv > 0) + { + lel = iel; // Set up control point (nodal) coordinates for current element + if (!this->getElementCoordinates(Xnod,iel+1)) + return false; + } // Evaluate basis function values/derivatives at current parametric point // and multiply with control point values to get the point-wise solution @@ -1815,19 +1661,17 @@ bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand, bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand, const RealArray* gpar, bool) const { -#ifdef SP_DEBUG - std::cout <<"ASMu3D::evalSolution(Matrix&,const IntegrandBase&,const RealArray*,bool)\n"; -#endif - sField.resize(0,0); + size_t nPoints = gpar[0].size(); + if (nPoints != gpar[1].size() || nPoints != gpar[2].size()) + return false; + + PROFILE2("ASMu3D::evalSol(S)"); // TODO: investigate the possibility of doing "regular" refinement by // uniform tesselation grid and ignoring LR mesh lines - size_t nPoints = gpar[0].size(); bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES; - if (nPoints != gpar[1].size() || nPoints != gpar[2].size()) - return false; Vector solPt; Matrix dNdu, Jac, Xnod; @@ -1839,6 +1683,7 @@ bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand, Matrix B(p1*p2*p3, 4); // Bezier evaluation points and derivatives // Evaluate the secondary solution field at each point + int lel = -1; for (size_t i = 0; i < nPoints; i++) { // Fetch element containing evaluation point @@ -1884,8 +1729,12 @@ bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand, else this->evaluateBasis(fe, dNdu, bezierExtract[iel], B); - // Set up control point (nodal) coordinates for current element - if (!this->getElementCoordinates(Xnod,iel+1)) return false; + if (iel != lel) + { + lel = iel; // Set up control point (nodal) coordinates for current element + if (!this->getElementCoordinates(Xnod,iel+1)) + return false; + } // Compute the Jacobian inverse fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); @@ -1896,6 +1745,11 @@ bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand, if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) continue; +#if SP_DEBUG > 4 + if (1+iel == dbgElm || dbgElm == 0) + std::cout <<"\n"<< fe; +#endif + // Now evaluate the solution field if (!integrand.evalSol(solPt,fe,Xnod*fe.N,MNPC[iel])) return false; @@ -2289,7 +2143,7 @@ bool ASMu3D::refine (const RealFunc& refC, double refTol) { Go::Point X0; int iel = 0; - std::vector elements; + IntVec elements; for (const LR::Element* elm : lrspline->getAllElements()) { double u0 = 0.5*(elm->umin() + elm->umax()); diff --git a/src/ASM/LR/ASMu3Drecovery.C b/src/ASM/LR/ASMu3Drecovery.C index 9e93f55d..627aa4de 100644 --- a/src/ASM/LR/ASMu3Drecovery.C +++ b/src/ASM/LR/ASMu3Drecovery.C @@ -30,9 +30,11 @@ bool ASMu3D::getGrevilleParameters (RealArray& prm, int dir, int basisNum) const if (!this->getBasis(basisNum) || dir < 0 || dir > 2) return false; const LR::LRSpline* lrspline = this->getBasis(basisNum); + prm.clear(); prm.reserve(lrspline->nBasisFunctions()); - for (const LR::Basisfunction *b : lrspline->getAllBasisfunctions()) + + for (const LR::Basisfunction* b : lrspline->getAllBasisfunctions()) prm.push_back(b->getGrevilleParameter()[dir]); return true; @@ -109,15 +111,16 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B, const int p1 = lrspline->order(0); const int p2 = lrspline->order(1); const int p3 = lrspline->order(2); + const int pm = std::max(std::max(p1,p2),p3); // Get Gaussian quadrature points - const int ng1 = continuous ? nGauss : p1 - 1; - const int ng2 = continuous ? nGauss : p2 - 1; - const int ng3 = continuous ? nGauss : p3 - 1; + const int ng1 = continuous ? this->getNoGaussPt(pm,true) : p1 - 1; + const int ng2 = continuous ? ng1 : p2 - 1; + const int ng3 = continuous ? ng1 : p3 - 1; const double* xg = GaussQuadrature::getCoord(ng1); const double* yg = GaussQuadrature::getCoord(ng2); const double* zg = GaussQuadrature::getCoord(ng3); - const double* wg = continuous ? GaussQuadrature::getWeight(nGauss) : nullptr; + const double* wg = continuous ? GaussQuadrature::getWeight(ng1) : nullptr; if (!xg || !yg || !zg) return false; if (continuous && !wg) return false; @@ -147,8 +150,6 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B, this->getGaussPointParameters(gpar[0],0,ng1,iel,xg); this->getGaussPointParameters(gpar[1],1,ng2,iel,yg); this->getGaussPointParameters(gpar[2],2,ng3,iel,zg); - - // convert to unstructured mesh representation expandTensorGrid(gpar.data(),unstrGpar.data()); // Evaluate the secondary solution at all integration points @@ -194,8 +195,8 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B, eB[r-1].add(phi,sField(r,ip+1)*dJw); } - for (size_t i = 0; i < MNPC[iel-1].size(); ++i) { - for (size_t j = 0; j < MNPC[iel-1].size(); ++j) + for (size_t i = 0; i < eA.rows(); ++i) { + for (size_t j = 0; j < eA.cols(); ++j) A(MNPC[iel-1][i]+1, MNPC[iel-1][j]+1) += eA(i+1, j+1); int jp = MNPC[iel-1][i]+1; @@ -269,7 +270,7 @@ LR::LRSplineVolume* ASMu3D::scRecovery (const IntegrandBase& integrand) const size_t ip = 0; std::vector::const_iterator elStart, elEnd, el; std::vector supportElements; - for (LR::Basisfunction *b : lrspline->getAllBasisfunctions()) + for (LR::Basisfunction* b : lrspline->getAllBasisfunctions()) { #if SP_DEBUG > 2 std::cout <<"Basis: "<< *b <<"\n ng1 ="<< ng1 <<"\n ng2 ="<< ng2 <<"\n ng3 ="<< ng3 @@ -316,18 +317,15 @@ LR::LRSplineVolume* ASMu3D::scRecovery (const IntegrandBase& integrand) const for (el = elStart; el != elEnd; ++el) { int iel = (**el).getId()+1; +#if SP_DEBUG > 2 + std::cout <<"Element "<< **el << std::endl; +#endif // evaluate all gauss points for this element std::array gaussPt, unstrGauss; this->getGaussPointParameters(gaussPt[0],0,ng1,iel,xg); this->getGaussPointParameters(gaussPt[1],1,ng2,iel,yg); this->getGaussPointParameters(gaussPt[2],2,ng3,iel,zg); - -#if SP_DEBUG > 2 - std::cout << "Element " << **el << std::endl; -#endif - - // convert to unstructured mesh representation expandTensorGrid(gaussPt.data(),unstrGauss.data()); // Evaluate the secondary solution at all Gauss points @@ -471,11 +469,6 @@ bool ASMu3D::faceL2projection (const DirichletFace& face, StdVector B(n*m); A.resize(n,n); - // Get Gaussian quadrature points and weights - const double* xg = GaussQuadrature::getCoord(nGauss); - const double* wg = GaussQuadrature::getWeight(nGauss); - if (!xg || !wg) return false; - // find the normal direction for the face int faceDir; switch (face.edg) @@ -491,11 +484,21 @@ bool ASMu3D::faceL2projection (const DirichletFace& face, const int t1 = 1 + abs(faceDir)%3; // first tangent direction const int t2 = 1 + t1%3; // second tangent direction + // Get Gaussian quadrature points and weights + // Use the largest polynomial order of the two tangent directions + const int pmax = std::max(lrspline->order(t1),lrspline->order(t2)); + const int nGP = this->getNoGaussPt(pmax,true); + const double* xg = GaussQuadrature::getCoord(nGP); + const double* wg = GaussQuadrature::getWeight(nGP); + if (!xg || !wg) return false; + Vector N; Matrix dNdu, dNdX, Xnod, Jac; - Vec4 X; + double param[3] = { 0.0, 0.0, 0.0 }; + Vec4 X(param); // === Assembly loop over all elements on the patch face ===================== + for (size_t ie = 0; ie < face.MLGE.size(); ie++) // for all face elements { int iel = 1 + face.MLGE[ie]; @@ -513,7 +516,7 @@ bool ASMu3D::faceL2projection (const DirichletFace& face, gpar[d].fill(lrspline->endparam(d)); } else - this->getGaussPointParameters(gpar[d],d,nGauss,iel,xg); + this->getGaussPointParameters(gpar[d],d,nGP,iel,xg); // Get element face area in the parameter space double dA = 0.25*this->getParametricArea(iel,abs(faceDir)); @@ -522,13 +525,14 @@ bool ASMu3D::faceL2projection (const DirichletFace& face, // Set up control point coordinates for current element if (!this->getElementCoordinates(Xnod,iel)) return false; - double u = gpar[0].front(); - double v = gpar[1].front(); - double w = gpar[2].front(); + double u = param[0] = gpar[0].front(); + double v = param[1] = gpar[1].front(); + double w = param[2] = gpar[2].front(); - // --- Integration loop over all Gauss points over the face ------------- - for (int j = 0; j < nGauss; j++) - for (int i = 0; i < nGauss; i++) + // --- Integration loop over all Gauss points over the face ---------------- + + for (int j = 0; j < nGP; j++) + for (int i = 0; i < nGP; i++) { // Parameter values of current integration point int k1, k2, k3; @@ -539,9 +543,9 @@ bool ASMu3D::faceL2projection (const DirichletFace& face, case 3: k1 = i; k2 = j; k3 = 0; break; default: k1 = k2 = k3 = 0; } - if (gpar[0].size() > 1) u = gpar[0](k1+1); - if (gpar[1].size() > 1) v = gpar[1](k2+1); - if (gpar[2].size() > 1) w = gpar[2](k3+1); + if (gpar[0].size() > 1) u = param[0] = gpar[0](k1+1); + if (gpar[1].size() > 1) v = param[1] = gpar[1](k2+1); + if (gpar[2].size() > 1) w = param[2] = gpar[2](k3+1); // Evaluate basis function derivatives at integration points this->evaluateBasis(iel-1, myGeoBasis, u, v, w, N, dNdu); @@ -551,7 +555,7 @@ bool ASMu3D::faceL2projection (const DirichletFace& face, if (dJxW == 0.0) continue; // skip singular points // Cartesian coordinates of current integration point - X = Xnod * N; + X.assign(Xnod * N); X.t = time; // For mixed basis, we need to compute functions separate from geometry