OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

changed: use getBasis in ASMxxD::assembleL2matrices

Browse Source 
make it clear that we want the integration and projection bases.
This commit is contained in:
Arne Morten Kvarving 2023-08-28 13:37:54 +02:00
parent ad7ff2b19e
commit ab100c290d
4 changed files with 123 additions and 132 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
src/ASM/ASMs2Drecovery.C
View File

@ -166,15 +166,17 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
{
const size_t nnod = this->getNoProjectionNodes();
const Go::SplineSurface* proj = static_cast<const Go::SplineSurface*>(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<Go::BasisPtsSf> spl0;
std::vector<Go::BasisDerivsSf> spl1, spl2;
std::vector<Go::BasisPtsSf> spl1;
std::vector<Go::BasisDerivsSf> 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;

126
src/ASM/ASMs3Drecovery.C
View File

@ -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<const ASMs3D*>(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<Vector&>(vec2),
sValues.rows(),
svol->rational(),
weights);
svol->basis(1),
svol->basis(2),
gpar[0], gpar[1], gpar[2],
const_cast<Vector&>(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<const Go::SplineVolume*>(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<Go::BasisPts> spl0;
std::vector<Go::BasisDerivs> spl1, spl2;
std::vector<Go::BasisPts> spl1;
std::vector<Go::BasisDerivs> 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);
}

42
src/ASM/LR/ASMu2Drecovery.C
View File

@ -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<const LR::LRSplineSurface*>(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<const LR::LRSplineSurface*>(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

42
src/ASM/LR/ASMu3Drecovery.C
View File

@ -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<const LR::LRSplineVolume*>(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<const LR::LRSplineVolume*>(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