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

added: support for Integrand::AVERAGE flag in LR integrands

Browse Source
This commit is contained in:
Arne Morten Kvarving 2021-12-21 19:10:50 +01:00
parent 4a347ed412
commit a4b38fff9d
2 changed files with 58 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
src/ASM/LR/ASMu2D.C
View File

@ -1187,6 +1187,32 @@ bool ASMu2D::integrate (Integrand& integrand,
dXidu[1] = el->vmax() - el->vmin();
}
if (integrand.getIntegrandType() & Integrand::AVERAGE)
{
// --- Compute average value of basis functions over the element -----
fe.Navg.resize(el->support().size(), true);
double area = 0.0;
int ip = (iel-1)*nGP*nGP;
for (int j = 0; j < nGP; j++)
for (int i = 0; i < nGP; i++, ip++)
{
SplineUtils::extractBasis(spline1[ip],fe.N,dNdu);
// Compute Jacobian determinant of coordinate mapping
// and multiply by weight of current integration point
double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false);
double weight = dA*wg[i]*wg[j];
// Numerical quadrature
fe.Navg.add(fe.N,detJac*weight);
area += detJac*weight;
}
// Divide by element area
fe.Navg /= area;
}
// Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(el->support().size(),fe.iel);
if (!integrand.initElement(MNPC[iel-1],fe,X,nRed*nRed,*A))

32
src/ASM/LR/ASMu3D.C
View File

@ -1089,6 +1089,38 @@ bool ASMu3D::integrate (Integrand& integrand,
dXidu[2] = el->getParmin(2);
}
if (integrand.getIntegrandType() & Integrand::AVERAGE)
{
// --- Compute average value of basis functions over the element -----
int ip = (iel-1)*nGP*nGP*nGP + firstIp;
fe.Navg.resize(el->support().size(),true);
double vol = 0.0;
int ig = 1;
for (int k = 0; k < nGP; k++)
for (int j = 0; j < nGP; j++)
for (int i = 0; i < nGP; i++, ++ip, ++ig)
{
B.fillColumn(1, BN.getColumn(ig));
B.fillColumn(2, BdNdu.getColumn(ig)*2.0/du);
B.fillColumn(3, BdNdv.getColumn(ig)*2.0/dv);
B.fillColumn(4, BdNdw.getColumn(ig)*2.0/dw);
this->evaluateBasis(fe, dNdu, C, B);
// Compute Jacobian determinant of coordinate mapping
// and multiply by weight of current integration point
double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false);
double weight = dV*wg[i]*wg[j]*wg[k];
// Numerical quadrature
fe.Navg.add(fe.N,detJac*weight);
vol += detJac*weight;
}
// Divide by element volume
fe.Navg /= vol;
}
// Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(el->support().size(),fe.iel);
if (!integrand.initElement(MNPC[iel-1],fe,X,nRed*nRed*nRed,*A))