From ab100c290d97e41a1405806ec279e56e6511b4cf Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 28 Aug 2023 13:37:54 +0200 Subject: [PATCH] changed: use getBasis in ASMxxD::assembleL2matrices make it clear that we want the integration and projection bases. --- src/ASM/ASMs2Drecovery.C | 45 ++++++------- src/ASM/ASMs3Drecovery.C | 126 +++++++++++++++++------------------- src/ASM/LR/ASMu2Drecovery.C | 42 ++++++------ src/ASM/LR/ASMu3Drecovery.C | 42 ++++++------ 4 files changed, 123 insertions(+), 132 deletions(-) diff --git a/src/ASM/ASMs2Drecovery.C b/src/ASM/ASMs2Drecovery.C index 4d648bb1..23f97478 100644 --- a/src/ASM/ASMs2Drecovery.C +++ b/src/ASM/ASMs2Drecovery.C @@ -166,15 +166,17 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B, { const size_t nnod = this->getNoProjectionNodes(); - const Go::SplineSurface* proj = static_cast(projB); + const Go::SplineSurface* itg = this->getBasis(ASM::INTEGRATION_BASIS); + const Go::SplineSurface* proj = this->getBasis(ASM::PROJECTION_BASIS); + const bool separateProjBasis = proj != itg; - const int g1 = surf->order_u(); - const int g2 = surf->order_v(); + const int g1 = itg->order_u(); + const int g2 = itg->order_v(); const int p1 = proj->order_u(); const int p2 = proj->order_v(); const int n1 = proj->numCoefs_u(); - const int nel1 = surf->numCoefs_u() - g1 + 1; - const int nel2 = surf->numCoefs_v() - g2 + 1; + const int nel1 = itg->numCoefs_u() - g1 + 1; + const int nel2 = itg->numCoefs_v() - g2 + 1; const int pmax = p1 > p2 ? p1 : p2; // Get Gaussian quadrature point coordinates (and weights if continuous) @@ -193,15 +195,13 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B, gpar[1] = this->getGaussPointParameters(gp,1,ng2,yg); // Evaluate basis functions at all integration points - std::vector spl0; - std::vector spl1, spl2; + std::vector spl1; + std::vector spl2; if (continuous) - { + itg->computeBasisGrid(gpar[0],gpar[1],spl2); + + if (!continuous || separateProjBasis) proj->computeBasisGrid(gpar[0],gpar[1],spl1); - surf->computeBasisGrid(gpar[0],gpar[1],spl2); - } - else - proj->computeBasisGrid(gpar[0],gpar[1],spl0); // Evaluate the secondary solution at all integration points Matrix sField; @@ -213,7 +213,7 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B, } double dA = 1.0; - Vector phi(p1*p2), phi2(g1*g2); + Vector phi(p1*p2); Matrix dNdu, Xnod, J; @@ -236,20 +236,17 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B, int ip = (i2*ng1*nel1 + i1)*ng2; IntVec lmnpc; - if (proj != surf) + if (separateProjBasis) { // Establish nodal point correspondance for the projection element int i, j, vidx; lmnpc.reserve(phi.size()); - if (continuous) - vidx = (spl1[ip].left_idx[1]-p1+1)*n1 + (spl1[ip].left_idx[0]-p1+1); - else - vidx = (spl0[ip].left_idx[1]-p1+1)*n1 + (spl0[ip].left_idx[0]-p1+1); + vidx = (spl1[ip].left_idx[1]-p1+1)*n1 + (spl1[ip].left_idx[0]-p1+1); for (j = 0; j < p2; j++, vidx += n1) for (i = 0; i < p1; i++) lmnpc.push_back(vidx+i); } - const IntVec& mnpc = proj == surf ? MNPC[iel] : lmnpc; + const IntVec& mnpc = separateProjBasis ? lmnpc : MNPC[iel]; // --- Integration loop over all Gauss points in each direction ---------- @@ -259,12 +256,10 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B, for (int i = 0; i < ng1; i++, ip++) { if (continuous) - { - SplineUtils::extractBasis(spl1[ip],phi,dNdu); - SplineUtils::extractBasis(spl2[ip],phi2,dNdu); - } - else - phi = spl0[ip].basisValues; + SplineUtils::extractBasis(spl2[ip],phi,dNdu); + + if (!continuous || separateProjBasis) + phi = spl1[ip].basisValues; // Compute the Jacobian inverse and derivatives double dJw = 1.0; diff --git a/src/ASM/ASMs3Drecovery.C b/src/ASM/ASMs3Drecovery.C index 04d2f9f1..eb0f4d18 100644 --- a/src/ASM/ASMs3Drecovery.C +++ b/src/ASM/ASMs3Drecovery.C @@ -63,7 +63,7 @@ bool ASMs3D::getQuasiInterplParameters (RealArray& prm, int dir) const bool ASMs3D::evaluate (const ASMbase* basis, const Vector& locVec, - RealArray& vec, int basisNum) const + RealArray& vec, int basisNum) const { const ASMs3D* pch = dynamic_cast(basis); if (!pch) return false; @@ -97,13 +97,13 @@ bool ASMs3D::evaluate (const ASMbase* basis, const Vector& locVec, const Vector& vec2 = sValues; Go::SplineVolume* vol_new = Go::VolumeInterpolator::regularInterpolation(svol->basis(0), - svol->basis(1), - svol->basis(2), - gpar[0], gpar[1], gpar[2], - const_cast(vec2), - sValues.rows(), - svol->rational(), - weights); + svol->basis(1), + svol->basis(2), + gpar[0], gpar[1], gpar[2], + const_cast(vec2), + sValues.rows(), + svol->rational(), + weights); vec.assign(vol_new->coefs_begin(),vol_new->coefs_end()); delete vol_new; @@ -166,19 +166,21 @@ bool ASMs3D::assembleL2matrices (SparseMatrix& A, StdVector& B, { const size_t nnod = this->getNoProjectionNodes(); - const Go::SplineVolume* proj = static_cast(projB); + const Go::SplineVolume* itg = this->getBasis(ASM::INTEGRATION_BASIS); + const Go::SplineVolume* proj = this->getBasis(ASM::PROJECTION_BASIS); + const bool separateProjBasis = proj != itg; - const int g1 = svol->order(0); - const int g2 = svol->order(1); - const int g3 = svol->order(2); + const int g1 = itg->order(0); + const int g2 = itg->order(1); + const int g3 = itg->order(2); const int p1 = proj->order(0); const int p2 = proj->order(1); const int p3 = proj->order(2); const int n1 = proj->numCoefs(0); const int n2 = proj->numCoefs(1); - const int nel1 = svol->numCoefs(0) - g1 + 1; - const int nel2 = svol->numCoefs(1) - g2 + 1; - const int nel3 = svol->numCoefs(2) - g3 + 1; + const int nel1 = itg->numCoefs(0) - g1 + 1; + const int nel2 = itg->numCoefs(1) - g2 + 1; + const int nel3 = itg->numCoefs(2) - g3 + 1; int pmax = p1 > p2 ? p1 : p2; if (pmax < p3) pmax = p3; @@ -202,15 +204,13 @@ bool ASMs3D::assembleL2matrices (SparseMatrix& A, StdVector& B, gpar[2] = this->getGaussPointParameters(gp,2,ng3,zg); // Evaluate basis functions at all integration points - std::vector spl0; - std::vector spl1, spl2; + std::vector spl1; + std::vector spl2; if (continuous) - { + itg->computeBasisGrid(gpar[0],gpar[1],gpar[2],spl2); + + if (!continuous || separateProjBasis) proj->computeBasisGrid(gpar[0],gpar[1],gpar[2],spl1); - svol->computeBasisGrid(gpar[0],gpar[1],gpar[2],spl2); - } - else - proj->computeBasisGrid(gpar[0],gpar[1],gpar[2],spl0); // Evaluate the secondary solution at all integration points Matrix sField; @@ -234,62 +234,56 @@ bool ASMs3D::assembleL2matrices (SparseMatrix& A, StdVector& B, for (int i2 = 0; i2 < nel2; i2++) for (int i1 = 0; i1 < nel1; i1++, iel++) { - if (MLGE[iel] < 1) continue; // zero-volume element + if (MLGE[iel] < 1) continue; // zero-volume element - if (continuous) - { - // Set up control point (nodal) coordinates for current element - if (!this->getElementCoordinates(Xnod,1+iel)) - return false; - else if ((dV = 0.125*this->getParametricVolume(1+iel)) < 0.0) - return false; // topology error (probably logic error) - } + if (continuous) + { + // Set up control point (nodal) coordinates for current element + if (!this->getElementCoordinates(Xnod,1+iel)) + return false; + else if ((dV = 0.125*this->getParametricVolume(1+iel)) < 0.0) + return false; // topology error (probably logic error) + } int ip = ((i3*ng2*nel2 + i2)*ng1*nel1 + i1)*ng3; IntVec lmnpc; - if (proj != svol) + if (separateProjBasis) { // Establish nodal point correspondance for the projection element int i, j, k, vidx, widx; lmnpc.reserve(phi.size()); - if (continuous) - widx = ((spl1[ip].left_idx[2]-p3+1)*n1*n2 + - (spl1[ip].left_idx[1]-p2+1)*n1 + - (spl1[ip].left_idx[0]-p1+1)); - else - widx = ((spl0[ip].left_idx[2]-p3+1)*n1*n2 + - (spl0[ip].left_idx[1]-p2+1)*n1 + - (spl0[ip].left_idx[0]-p1+1)); + widx = ((spl1[ip].left_idx[2]-p3+1)*n1*n2 + + (spl1[ip].left_idx[1]-p2+1)*n1 + + (spl1[ip].left_idx[0]-p1+1)); + for (k = 0; k < p3; k++, widx += n1*n2) for (j = vidx = 0; j < p2; j++, vidx += n1) for (i = 0; i < p1; i++) lmnpc.push_back(widx+vidx+i); } - const IntVec& mnpc = proj == svol ? MNPC[iel] : lmnpc; + const IntVec& mnpc = separateProjBasis ? lmnpc : MNPC[iel]; - // --- Integration loop over all Gauss points in each direction -------- + // --- Integration loop over all Gauss points in each direction -------- Matrix eA(p1*p2*p3, p1*p2*p3); Vectors eB(sField.rows(), Vector(p1*p2*p3)); for (int k = 0; k < ng3; k++, ip += ng2*(nel2-1)*ng1*nel1) for (int j = 0; j < ng2; j++, ip += ng1*(nel1-1)) for (int i = 0; i < ng1; i++, ip++) - { - if (continuous) - { - SplineUtils::extractBasis(spl1[ip],phi,dNdu); - SplineUtils::extractBasis(spl2[ip],phi2,dNdu); - } - else - phi = spl0[ip].basisValues; + { + if (continuous) + SplineUtils::extractBasis(spl2[ip],phi,dNdu); - // Compute the Jacobian inverse and derivatives - double dJw = dV; - if (continuous) - { - dJw *= wg[i]*wg[j]*wg[k]*utl::Jacobian(J,dNdu,Xnod,dNdu,false); - if (dJw == 0.0) continue; // skip singular points - } + if (!continuous || separateProjBasis) + phi = spl1[ip].basisValues; + + // Compute the Jacobian inverse and derivatives + double dJw = dV; + if (continuous) + { + dJw *= wg[i]*wg[j]*wg[k]*utl::Jacobian(J,dNdu,Xnod,dNdu,false); + if (dJw == 0.0) continue; // skip singular points + } // Integrate the mass matrix eA.outer_product(phi, phi, true, dJw); @@ -297,7 +291,7 @@ bool ASMs3D::assembleL2matrices (SparseMatrix& A, StdVector& B, // Integrate the rhs vector B for (size_t r = 1; r <= sField.rows(); r++) eB[r-1].add(phi,sField(r,ip+1)*dJw); - } + } for (int i = 0; i < p1*p2*p3; ++i) { for (int j = 0; j < p1*p2*p3; ++j) @@ -307,7 +301,7 @@ bool ASMs3D::assembleL2matrices (SparseMatrix& A, StdVector& B, for (size_t r = 0; r < sField.rows(); r++, jp += nnod) B(jp) += eB[r](1+i); } - } + } return true; } @@ -432,13 +426,13 @@ Go::SplineVolume* ASMs3D::projectSolutionLocal (const IntegrandBase& integrand) svol->getWeights(weights); return quasiInterpolation(svol->basis(0), - svol->basis(1), - svol->basis(2), - gpar[0], gpar[1], gpar[2], - sValues, - sValues.rows(), - svol->rational(), - weights); + svol->basis(1), + svol->basis(2), + gpar[0], gpar[1], gpar[2], + sValues, + sValues.rows(), + svol->rational(), + weights); } diff --git a/src/ASM/LR/ASMu2Drecovery.C b/src/ASM/LR/ASMu2Drecovery.C index b4a93ae2..c2ed92b4 100644 --- a/src/ASM/LR/ASMu2Drecovery.C +++ b/src/ASM/LR/ASMu2Drecovery.C @@ -108,8 +108,12 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B, { size_t nnod = this->getNoProjectionNodes(); - const int p1 = projB->order(0); - const int p2 = projB->order(1); + const LR::LRSplineSurface* itg = this->getBasis(ASM::INTEGRATION_BASIS); + const LR::LRSplineSurface* proj = this->getBasis(ASM::PROJECTION_BASIS); + const bool separateProjBasis = proj != itg; + + const int p1 = proj->order(0); + const int p2 = proj->order(1); const int pm = p1 > p2 ? p1 : p2; // Get Gaussian quadrature points @@ -121,17 +125,16 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B, if (!xg || !yg) return false; if (continuous && !wg) return false; - bool singleBasis = (this->getNoBasis() == 1 && projB == lrspline); IntMat lmnpc; - const IntMat& gmnpc = singleBasis ? MNPC : lmnpc; - if (!singleBasis) { - lmnpc.resize(projB->nElements()); - for (const LR::Element* elm : projB->getAllElements()) { + if (separateProjBasis) { + lmnpc.resize(proj->nElements()); + for (const LR::Element* elm : proj->getAllElements()) { lmnpc[elm->getId()].reserve(elm->nBasisFunctions()); for (const LR::Basisfunction* f : elm->support()) lmnpc[elm->getId()].push_back(f->getId()); } } + const IntMat& gmnpc = separateProjBasis ? lmnpc : MNPC; A.preAssemble(gmnpc, gmnpc.size()); // === Assembly loop over all elements in the patch ========================== @@ -142,12 +145,12 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B, for (size_t e = 0; e < group[t].size(); e++) { double dA = 0.0; - Vector phi, phi2; + Vector phi; Matrix dNdu, Xnod, Jac; - Go::BasisPtsSf spl0; - Go::BasisDerivsSf spl1, spl2; + Go::BasisPtsSf spl1; + Go::BasisDerivsSf spl2; int ielp = group[t][e]; - const LR::Element* elm = projB->getElement(ielp); + const LR::Element* elm = proj->getElement(ielp); int iel = lrspline->getElementContaining(elm->midpoint())+1; if (continuous) @@ -178,7 +181,7 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B, // Set up basis function size (for extractBasis subroutine) size_t nbf = elm->nBasisFunctions(); - const IntVec& mnpc = singleBasis ? gmnpc[iel-1] : gmnpc[ielp]; + const IntVec& mnpc = separateProjBasis ? gmnpc[ielp] : gmnpc[iel-1]; // --- Integration loop over all Gauss points in each direction ------------ Matrix eA(nbf, nbf); @@ -189,17 +192,14 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B, { if (continuous) { - this->computeBasis(gpar[0][i],gpar[1][j],spl1,ielp, - static_cast(projB.get())); - SplineUtils::extractBasis(spl1,phi,dNdu); - this->computeBasis(gpar[0][i],gpar[1][j],spl2,iel-1); - SplineUtils::extractBasis(spl2,phi2,dNdu); + this->computeBasis(gpar[0][i],gpar[1][j],spl2,iel-1,itg); + SplineUtils::extractBasis(spl2,phi,dNdu); } - else + + if (!continuous || separateProjBasis) { - this->computeBasis(gpar[0][i],gpar[1][j],spl0,ielp, - static_cast(projB.get())); - phi = spl0.basisValues; + this->computeBasis(gpar[0][i],gpar[1][j],spl1,ielp,proj); + phi = spl1.basisValues; } // Compute the Jacobian inverse and derivatives diff --git a/src/ASM/LR/ASMu3Drecovery.C b/src/ASM/LR/ASMu3Drecovery.C index 48315817..f9b0e666 100644 --- a/src/ASM/LR/ASMu3Drecovery.C +++ b/src/ASM/LR/ASMu3Drecovery.C @@ -110,9 +110,13 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B, { size_t nnod = this->getNoProjectionNodes(); - const int p1 = projB->order(0); - const int p2 = projB->order(1); - const int p3 = projB->order(2); + const LR::LRSplineVolume* itg = this->getBasis(ASM::INTEGRATION_BASIS); + const LR::LRSplineVolume* proj = this->getBasis(ASM::PROJECTION_BASIS); + const bool separateProjBasis = proj != itg; + + const int p1 = proj->order(0); + const int p2 = proj->order(1); + const int p3 = proj->order(2); const int pm = std::max(std::max(p1,p2),p3); // Get Gaussian quadrature points @@ -126,17 +130,16 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B, if (!xg || !yg || !zg) return false; if (continuous && !wg) return false; - bool singleBasis = (this->getNoBasis() == 1 && projB == lrspline); IntMat lmnpc; - const IntMat& gmnpc = singleBasis ? MNPC : lmnpc; - if (!singleBasis) { - lmnpc.resize(projB->nElements()); - for (const LR::Element* elm : projB->getAllElements()) { + if (separateProjBasis) { + lmnpc.resize(proj->nElements()); + for (const LR::Element* elm : proj->getAllElements()) { lmnpc[elm->getId()].reserve(elm->nBasisFunctions()); for (const LR::Basisfunction* f : elm->support()) lmnpc[elm->getId()].push_back(f->getId()); } } + const IntMat& gmnpc = separateProjBasis ? lmnpc : MNPC; A.preAssemble(gmnpc, gmnpc.size()); // === Assembly loop over all elements in the patch ========================== @@ -147,12 +150,12 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B, for (size_t e = 0; e < group[t].size(); e++) { double dV = 0.0; - Vector phi, phi2; + Vector phi; Matrix dNdu, Xnod, Jac; - Go::BasisPts spl0; - Go::BasisDerivs spl1, spl2; + Go::BasisPts spl1; + Go::BasisDerivs spl2; int ielp = group[t][e]; - const LR::Element* elm = projB->getElement(ielp); + const LR::Element* elm = proj->getElement(ielp); int iel = lrspline->getElementContaining(elm->midpoint())+1; if (continuous) @@ -183,7 +186,7 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B, // Set up basis function size (for extractBasis subroutine) size_t nbf = elm->nBasisFunctions(); - const IntVec& mnpc = singleBasis ? gmnpc[iel-1] : gmnpc[ielp]; + const IntVec& mnpc = separateProjBasis ? gmnpc[ielp] : gmnpc[iel-1]; // --- Integration loop over all Gauss points in each direction ------------ @@ -196,15 +199,14 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B, { if (continuous) { - static_cast(projB.get())->computeBasis(gpar[0][i],gpar[1][j],gpar[2][k],spl1,ielp); - SplineUtils::extractBasis(spl1,phi,dNdu); - lrspline->computeBasis(gpar[0][i],gpar[1][j],gpar[2][k],spl2,iel-1); - SplineUtils::extractBasis(spl2,phi2,dNdu); + itg->computeBasis(gpar[0][i],gpar[1][j],gpar[2][k],spl2,iel-1); + SplineUtils::extractBasis(spl2,phi,dNdu); } - else + + if (!continuous || separateProjBasis) { - static_cast(projB.get())->computeBasis(gpar[0][i],gpar[1][j],gpar[2][k],spl0,ielp); - phi = spl0.basisValues; + proj->computeBasis(gpar[0][i],gpar[1][j],gpar[2][k],spl1,ielp); + phi = spl1.basisValues; } // Compute the Jacobian inverse and derivatives