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

Fixed: Using direction-dependent Gauss quadrature on face integrals

Browse Source
This commit is contained in:
Knut Morten Okstad
2021-03-11 17:54:18 +01:00
parent 9421df2838
commit ff60831796
4 changed files with 163 additions and 68 deletions
Download Patch File Download Diff File Expand all files Collapse all files

75
src/ASM/ASMs3D.C
View File

@@ -2478,28 +2478,43 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
const int t2 = 1 + t1%3; // second tangent direction
// Get Gaussian quadrature points and weights
// For now, use the largest polynomial order of the two tangent directions
int nG1 = this->getNoGaussPt(std::max(svol->order(t1-1),svol->order(t2-1)),true);
int nGP = integrand.getBouIntegrationPoints(nG1);
const double* xg = GaussQuadrature::getCoord(nGP);
const double* wg = GaussQuadrature::getWeight(nGP);
if (!xg || !wg) return false;
// Compute parameter values of the Gauss points over the whole patch face
// and compute parameter values of the Gauss points over the whole patch face
std::array<int,3> ng;
std::array<const double*,3> xg, wg;
std::array<Matrix,3> gpar;
for (int d = 0; d < 3; d++)
if (-1-d == faceDir)
{
ng[d] = 1;
xg[d] = nullptr;
wg[d] = nullptr;
gpar[d].resize(1,1);
gpar[d].fill(svol->startparam(d));
}
else if (1+d == faceDir)
{
ng[d] = 1;
xg[d] = nullptr;
wg[d] = nullptr;
gpar[d].resize(1,1);
gpar[d].fill(svol->endparam(d));
}
else
this->getGaussPointParameters(gpar[d],d,nGP,xg);
{
int n = this->getNoGaussPt(svol->order(d),true);
ng[d] = integrand.getBouIntegrationPoints(n);
xg[d] = GaussQuadrature::getCoord(ng[d]);
wg[d] = GaussQuadrature::getWeight(ng[d]);
if (xg[d] && wg[d])
this->getGaussPointParameters(gpar[d],d,ng[d],xg[d]);
else
return false;
}
const int tt1 = t1 > t2 ? t2-1 : t1-1;
const int tt2 = t1 > t2 ? t1-1 : t2-1;
const int nG1 = ng[tt1];
const int nG2 = ng[tt2];
// Extract the Neumann order flag (1 or higher) for the integrand
integrand.setNeumannOrder(1 + lIndex/10);
@@ -2607,48 +2622,54 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
}
// Define some loop control variables depending on which face we are on
int nf1, j1, j2;
int nf1 = 0, j1 = 0, j2 = 0;
switch (abs(faceDir))
{
case 1: nf1 = nel2; j2 = i3-p3; j1 = i2-p2; break;
case 2: nf1 = nel1; j2 = i3-p3; j1 = i1-p1; break;
case 3: nf1 = nel1; j2 = i2-p2; j1 = i1-p1; break;
default: nf1 = j1 = j2 = 0;
}
// --- Integration loop over all Gauss points in each direction --------
int k1, k2, k3;
int ip = (j2*nGP*nf1 + j1)*nGP;
int jp = (j2*nf1 + j1)*nGP*nGP;
int ip = (j2*nf1*nG2 + j1)*nG1;
int jp = (j2*nf1 + j1)*nG1*nG2;
fe.iGP = firstp + jp; // Global integration point counter
for (int j = 0; j < nGP; j++, ip += nGP*(nf1-1))
for (int i = 0; i < nGP; i++, ip++, fe.iGP++)
for (int j = 0; j < nG2; j++, ip += nG1*(nf1-1))
for (int i = 0; i < nG1; i++, ip++, fe.iGP++)
{
#if SP_DEBUG > 4
std::cout <<"Elem "<< iel <<": "<< i1-p1 <<" "<< i2-p2 <<" "<< i3-p3
<<", Face "<< faceDir <<": "<< j1 <<" "<< j2
<<", Point "<< ip <<": "<< i <<" "<< j <<" "
<< spline.size() << std::endl;
#endif
// Local element coordinates and parameter values
// of current integration point
int k1 = 0, k2 = 0, k3 = 0;
switch (abs(faceDir))
{
case 1: k2 = i; k3 = j; k1 = 0; break;
case 2: k1 = i; k3 = j; k2 = 0; break;
case 3: k1 = i; k2 = j; k3 = 0; break;
default: k1 = k2 = k3 = 0;
case 1: k2 = i; k3 = j; break;
case 2: k1 = i; k3 = j; break;
case 3: k1 = i; k2 = j; break;
}
if (gpar[0].size() > 1)
if (xg[0])
{
fe.xi = xg[k1];
fe.xi = xg[0][k1];
fe.u = param[0] = gpar[0](k1+1,i1-p1+1);
}
if (gpar[1].size() > 1)
if (xg[1])
{
fe.eta = xg[k2];
fe.eta = xg[1][k2];
fe.v = param[1] = gpar[1](k2+1,i2-p2+1);
}
if (gpar[2].size() > 1)
if (xg[2])
{
fe.zeta = xg[k3];
fe.zeta = xg[2][k3];
fe.w = param[2] = gpar[2](k3+1,i3-p3+1);
}
@@ -2675,7 +2696,7 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
X.t = time.t;
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dA*wg[i]*wg[j];
fe.detJxW *= dA*wg[tt1][i]*wg[tt2][j];
if (!integrand.evalBou(*A,fe,time,X,normal))
ok = false;
}

16
src/ASM/ASMs3D.h
View File

@@ -671,13 +671,6 @@ protected:
void getCornerPoints(int i1, int i2, int i3,
std::vector<utl::Point>& XC) const;
//! \brief Generates element groups for multi-threading of interior integrals.
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] silence If \e true, suppress threading group outprint
//! \param[in] ignoreGlobalLM Sanity check option
virtual void generateThreadGroups(const Integrand& integrand, bool silence,
bool ignoreGlobalLM);
//! \brief Generates element groups for multi-threading of interior integrals.
//! \param[in] strip1 Strip width in first direction
//! \param[in] strip2 Strip width in second direction
@@ -687,6 +680,14 @@ protected:
void generateThreadGroups(size_t strip1, size_t strip2, size_t strip3,
bool silence, bool ignoreGlobalLM);
public:
//! \brief Generates element groups for multi-threading of interior integrals.
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] silence If \e true, suppress threading group outprint
//! \param[in] ignoreGlobalLM Sanity check option
virtual void generateThreadGroups(const Integrand& integrand, bool silence,
bool ignoreGlobalLM);
//! \brief Generates element groups for multi-threading of boundary integrals.
//! \param[in] lIndex Local index [1,6] of the boundary face
//! \param[in] silence If \e true, suppress threading group outprint
@@ -695,7 +696,6 @@ protected:
//! \brief Generates element groups from a partition.
virtual void generateThreadGroupsFromElms(const IntVec& elms);
public:
//! \brief Auxilliary function for computation of basis function indices.
static void scatterInd(int n1, int n2, int n3, int p1, int p2, int p3,
const int* start, IntVec& index);

98
src/ASM/ASMs3DLag.C
View File

@@ -553,12 +553,44 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
const int t2 = 1 + t1%3; // second tangent direction of the face
// Get Gaussian quadrature points and weights
// For now, use the largest polynomial order of the two tangent directions
int nG1 = this->getNoGaussPt(std::max(svol->order(t1-1),svol->order(t2-1)),true);
int nGP = integrand.getBouIntegrationPoints(nG1);
const double* xg = GaussQuadrature::getCoord(nGP);
const double* wg = GaussQuadrature::getWeight(nGP);
if (!xg || !wg) return false;
// and compute parameter values of the Gauss points over the whole patch face
std::array<int,3> ng;
std::array<const double*,3> xg, wg;
std::array<Matrix,3> gpar;
for (int d = 0; d < 3; d++)
if (-1-d == faceDir)
{
ng[d] = 1;
xg[d] = nullptr;
wg[d] = nullptr;
gpar[d].resize(1,1);
gpar[d].fill(svol->startparam(d));
}
else if (1+d == faceDir)
{
ng[d] = 1;
xg[d] = nullptr;
wg[d] = nullptr;
gpar[d].resize(1,1);
gpar[d].fill(svol->endparam(d));
}
else
{
int n = this->getNoGaussPt(svol->order(d),true);
ng[d] = integrand.getBouIntegrationPoints(n);
xg[d] = GaussQuadrature::getCoord(ng[d]);
wg[d] = GaussQuadrature::getWeight(ng[d]);
if (xg[d] && wg[d])
this->getGaussPointParameters(gpar[d],d,ng[d],xg[d]);
else
return false;
}
const int tt0 = t0-1;
const int tt1 = t1 > t2 ? t2-1 : t1-1;
const int tt2 = t1 > t2 ? t1-1 : t2-1;
const int nG1 = ng[tt1];
const int nG2 = ng[tt2];
// Number of elements in each direction
const int nel1 = (nx-1)/(p1-1);
@@ -594,6 +626,12 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
std::map<char,size_t>::const_iterator iit = firstBp.find(lIndex%10);
size_t firstp = iit == firstBp.end() ? 0 : iit->second;
// Lambda function for linear interpolation between two values
auto&& linearInt = [](double xi, double v1, double v2)
{
return 0.5*(v1*(1.0-xi) + v2*(1.0+xi));
};
// === Assembly loop over all elements on the patch face =====================
@@ -638,48 +676,43 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
}
// Define some loop control variables depending on which face we are on
int nf1, j1, j2;
int nf1 = 0, j1 = 0, j2 = 0;
switch (abs(faceDir))
{
case 1: nf1 = nel2; j2 = i3; j1 = i2; break;
case 2: nf1 = nel1; j2 = i3; j1 = i1; break;
case 3: nf1 = nel1; j2 = i2; j1 = i1; break;
default: nf1 = j1 = j2 = 0;
}
// --- Integration loop over all Gauss points in each direction --------
int k1, k2, k3;
int jp = (j2*nf1 + j1)*nGP*nGP;
int jp = (j2*nf1 + j1)*nG1*nG2;
fe.iGP = firstp + jp; // Global integration point counter
for (int j = 0; j < nGP; j++)
for (int i = 0; i < nGP; i++, fe.iGP++)
for (int j = 0; j < nG2; j++)
for (int i = 0; i < nG1; i++, fe.iGP++)
{
// Local element coordinates of current integration point
xi[t0-1] = faceDir < 0 ? -1.0 : 1.0;
xi[t1-1] = xg[i];
xi[t2-1] = xg[j];
xi[tt0] = faceDir < 0 ? -1.0 : 1.0;
xi[tt1] = xg[tt1][i];
xi[tt2] = xg[tt2][j];
fe.xi = xi[0];
fe.eta = xi[1];
fe.zeta = xi[2];
// Local element coordinates and parameter values
// of current integration point
switch (abs(faceDir))
{
case 1: k2 = i; k3 = j; k1 = -1; break;
case 2: k1 = i; k3 = j; k2 = -1; break;
case 3: k1 = i; k2 = j; k3 = -1; break;
default: k1 = k2 = k3 = -1;
}
if (upar.size() > 1)
fe.u = 0.5*(upar[i1]*(1.0-xg[k1]) + upar[i1+1]*(1.0+xg[k1]));
if (vpar.size() > 1)
fe.v = 0.5*(vpar[i2]*(1.0-xg[k2]) + vpar[i2+1]*(1.0+xg[k2]));
if (wpar.size() > 1)
fe.w = 0.5*(wpar[i3]*(1.0-xg[k3]) + wpar[i3+1]*(1.0+xg[k3]));
int k1 = -1, k2 = -1, k3 = -1;
switch (abs(faceDir))
{
case 1: k2 = i; k3 = j; break;
case 2: k1 = i; k3 = j; break;
case 3: k1 = i; k2 = j; break;
}
if (xg[0]) fe.u = linearInt(xg[0][k1],upar[i1],upar[i1+1]);
if (xg[1]) fe.v = linearInt(xg[1][k2],vpar[i2],vpar[i2+1]);
if (xg[2]) fe.w = linearInt(xg[2][k3],wpar[i3],wpar[i3+1]);
// Compute the basis functions and their derivatives, using
// tensor product of one-dimensional Lagrange polynomials
@@ -697,7 +730,7 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
X.t = time.t;
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i]*wg[j];
fe.detJxW *= wg[tt1][i]*wg[tt2][j];
if (!integrand.evalBou(*A,fe,time,X,normal))
ok = false;
}
@@ -967,6 +1000,7 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
const int nel = this->getNoElms(true);
for (int iel = 1; iel <= nel; iel++)
{
fe.iel = MLGE[iel-1];
const IntVec& mnpc = MNPC[iel-1];
this->getElementCoordinates(Xnod,iel);
@@ -981,8 +1015,6 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
if (!Lagrange::computeBasis(fe.N,dNdu,p1,fe.xi,p2,fe.eta,p3,fe.zeta))
return false;
fe.iel = MLGE[iel-1];
// Compute the Jacobian inverse
fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
@@ -1094,7 +1126,7 @@ bool ASMs3DLag::getGridParameters (RealArray& prm, int dir, int nSegPerSpan) con
}
bool ASMs3DLag::write(std::ostream& os, int) const
bool ASMs3DLag::write (std::ostream& os, int) const
{
return this->writeLagBasis(os, "hexahedron");
return this->writeLagBasis(os,"hexahedron");
}

42
src/ASM/Test/TestASMs3D.C
View File

@@ -11,6 +11,8 @@
//==============================================================================
#include "ASMCube.h"
#include "IntegrandBase.h"
#include "GlobalIntegral.h"
#include "SIM3D.h"
#include <array>
@@ -326,3 +328,43 @@ TEST(TestASMs3D, Collapse)
EXPECT_FALSE(pch.collapseFace(iface,iedge));
}
}
class NoProblem : public IntegrandBase
{
public:
NoProblem() : IntegrandBase(3) {}
virtual ~NoProblem() {}
protected:
virtual bool evalBou(LocalIntegral&, const FiniteElement&,
const Vec3&, const Vec3&) const { return true; }
};
TEST(TestASMs3D, FaceIntegrate)
{
GlobalIntegral dummy;
NoProblem prb;
ASMCube patch;
ASMbase::resetNumbering();
ASSERT_TRUE(patch.raiseOrder(2,1,0));
ASSERT_TRUE(patch.generateFEMTopology());
for (int lIndex = 1; lIndex <= 6; lIndex++)
{
patch.generateThreadGroups(lIndex,false,false);
ASSERT_TRUE(patch.integrate(prb,lIndex,dummy,TimeDomain()));
}
patch.clear(true);
ASMbase::resetNumbering();
ASSERT_TRUE(patch.uniformRefine(0,3));
ASSERT_TRUE(patch.uniformRefine(1,2));
ASSERT_TRUE(patch.uniformRefine(2,1));
ASSERT_TRUE(patch.generateFEMTopology());
for (int lIndex = 1; lIndex <= 6; lIndex++)
{
patch.generateThreadGroups(lIndex,false,false);
ASSERT_TRUE(patch.integrate(prb,lIndex,dummy,TimeDomain()));
}
}