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

Added: Default Gauss point selection for LRSplines

Browse Source
This commit is contained in:
Knut Morten Okstad 2019-03-12 14:32:39 +01:00
parent d5d3987b3c
commit b54d19d19c
4 changed files with 150 additions and 282 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -1058,15 +1058,19 @@ bool ASMu2D::integrate (Integrand& integrand,
bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES; bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES;
bool use3rdDer = integrand.getIntegrandType() & Integrand::THIRD_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 // Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss); const int nGP = this->getNoGaussPt(p1 > p2 ? p1 : p2);
const double* wg = GaussQuadrature::getWeight(nGauss); const double* xg = GaussQuadrature::getCoord(nGP);
const double* wg = GaussQuadrature::getWeight(nGP);
if (!xg || !wg) return false; if (!xg || !wg) return false;
// Get the reduced integration quadrature points, if needed // Get the reduced integration quadrature points, if needed
const double* xr = nullptr; const double* xr = nullptr;
const double* wr = nullptr; const double* wr = nullptr;
int nRed = integrand.getReducedIntegration(nGauss); int nRed = integrand.getReducedIntegration(nGP);
if (nRed > 0) if (nRed > 0)
{ {
xr = GaussQuadrature::getCoord(nRed); xr = GaussQuadrature::getCoord(nRed);
@ -1074,7 +1078,7 @@ bool ASMu2D::integrate (Integrand& integrand,
if (!xr || !wr) return false; if (!xr || !wr) return false;
} }
else if (nRed < 0) 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. // Evaluate basis function values and derivatives at all integration points.
// We do this before the integration point loop to exploit multi-threading // We do this before the integration point loop to exploit multi-threading
@ -1085,11 +1089,11 @@ bool ASMu2D::integrate (Integrand& integrand,
std::vector<Go::BasisDerivsSf3> spline3; std::vector<Go::BasisDerivsSf3> spline3;
if (use3rdDer) if (use3rdDer)
spline3.resize(nel*nGauss*nGauss); spline3.resize(nel*nGP*nGP);
else if (use2ndDer) else if (use2ndDer)
spline2.resize(nel*nGauss*nGauss); spline2.resize(nel*nGP*nGP);
else else
spline1.resize(nel*nGauss*nGauss); spline1.resize(nel*nGP*nGP);
if (xr) if (xr)
splineRed.resize(nel*nRed*nRed); splineRed.resize(nel*nRed*nRed);
@ -1097,10 +1101,10 @@ bool ASMu2D::integrate (Integrand& integrand,
for (iel = jp = rp = 0; iel < nel; iel++) for (iel = jp = rp = 0; iel < nel; iel++)
{ {
RealArray u, v; RealArray u, v;
this->getGaussPointParameters(u,0,nGauss,1+iel,xg); this->getGaussPointParameters(u,0,nGP,1+iel,xg);
this->getGaussPointParameters(v,1,nGauss,1+iel,xg); this->getGaussPointParameters(v,1,nGP,1+iel,xg);
for (int j = 0; j < nGauss; j++) for (int j = 0; j < nGP; j++)
for (int i = 0; i < nGauss; i++, jp++) for (int i = 0; i < nGP; i++, jp++)
if (use3rdDer) if (use3rdDer)
this->computeBasis(u[i],v[j],spline3[jp],iel); this->computeBasis(u[i],v[j],spline3[jp],iel);
else if (use2ndDer) else if (use2ndDer)
@ -1141,8 +1145,8 @@ bool ASMu2D::integrate (Integrand& integrand,
FiniteElement fe; FiniteElement fe;
fe.iel = MLGE[iel-1]; fe.iel = MLGE[iel-1];
fe.p = lrspline->order(0) - 1; fe.p = p1 - 1;
fe.q = lrspline->order(1) - 1; fe.q = p2 - 1;
Matrix dNdu, Xnod, Jac; Matrix dNdu, Xnod, Jac;
Matrix3D d2Ndu2, Hess; Matrix3D d2Ndu2, Hess;
Matrix4D d3Ndu3; Matrix4D d3Ndu3;
@ -1169,7 +1173,7 @@ bool ASMu2D::integrate (Integrand& integrand,
std::array<RealArray,2> gpar, redpar; std::array<RealArray,2> gpar, redpar;
for (int d = 0; d < 2; d++) 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) if (xr)
this->getGaussPointParameters(redpar[d],d,nRed,iel,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 ---------- // --- 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 fe.iGP = firstIp + jp; // Global integration point counter
for (int j = 0; j < nGauss; j++) for (int j = 0; j < nGP; j++)
for (int i = 0; i < nGauss; i++, fe.iGP++) for (int i = 0; i < nGP; i++, fe.iGP++)
{ {
// Local element coordinates of current integration point // Local element coordinates of current integration point
fe.xi = xg[i]; fe.xi = xg[i];
@ -1357,7 +1361,7 @@ bool ASMu2D::integrate (Integrand& integrand,
{ {
if (!lrspline) return true; // silently ignore empty patches if (!lrspline) return true; // silently ignore empty patches
if (integrand.getReducedIntegration(nGauss) != 0) if (integrand.getReducedIntegration(2) != 0)
{ {
std::cerr <<" *** ASMu2D::integrate(Integrand&,GlobalIntegral&," std::cerr <<" *** ASMu2D::integrate(Integrand&,GlobalIntegral&,"
<<"const TimeDomain&,const Real3DMat&): Available for standard" <<"const TimeDomain&,const Real3DMat&): Available for standard"
@ -1546,8 +1550,12 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
PROFILE2("ASMu2D::integrate(B)"); PROFILE2("ASMu2D::integrate(B)");
const int p1 = lrspline->order(0);
const int p2 = lrspline->order(1);
// Get Gaussian quadrature points and weights // 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* xg = GaussQuadrature::getCoord(nGP);
const double* wg = GaussQuadrature::getWeight(nGP); const double* wg = GaussQuadrature::getWeight(nGP);
if (!xg || !wg) return false; 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; size_t firstp = iit == firstBp.end() ? 0 : iit->second;
FiniteElement fe; FiniteElement fe;
fe.p = lrspline->order(0) - 1; fe.p = p1 - 1;
fe.q = lrspline->order(1) - 1; fe.q = p2 - 1;
fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0; fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0;
double param[3] = { 0.0, 0.0, 0.0 }; double param[3] = { 0.0, 0.0, 0.0 };

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

@ -106,13 +106,14 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
const int p1 = projBasis->order(0); const int p1 = projBasis->order(0);
const int p2 = projBasis->order(1); const int p2 = projBasis->order(1);
const int pm = p1 > p2 ? p1 : p2;
// Get Gaussian quadrature points // Get Gaussian quadrature points
const int ng1 = continuous ? nGauss : p1 - 1; const int ng1 = continuous ? this->getNoGaussPt(pm,true) : p1 - 1;
const int ng2 = continuous ? nGauss : p2 - 1; const int ng2 = continuous ? ng1 : p2 - 1;
const double* xg = GaussQuadrature::getCoord(ng1); const double* xg = GaussQuadrature::getCoord(ng1);
const double* yg = GaussQuadrature::getCoord(ng2); 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 (!xg || !yg) return false;
if (continuous && !wg) return false; if (continuous && !wg) return false;
@ -463,11 +464,6 @@ bool ASMu2D::edgeL2projection (const DirichletEdge& edge,
StdVector B(n*m); StdVector B(n*m);
A.resize(n,n); 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 // find the normal and tangent direction for the edge
int edgeDir, t1, t2; int edgeDir, t1, t2;
switch (edge.edg) switch (edge.edg)
@ -479,16 +475,22 @@ bool ASMu2D::edgeL2projection (const DirichletEdge& edge,
default: return false; 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<Vector,2> gpar; std::array<Vector,2> gpar;
for (int d = 0; d < 2; d++) for (int d = 0; d < 2; d++)
if (-1-d == edgeDir) if (-1-d == edgeDir)
{ {
gpar[d].resize(nGauss); gpar[d].resize(nGP);
gpar[d].fill(d == 0 ? lrspline->startparam(0) : lrspline->startparam(1)); gpar[d].fill(d == 0 ? lrspline->startparam(0) : lrspline->startparam(1));
} }
else if (1+d == edgeDir) else if (1+d == edgeDir)
{ {
gpar[d].resize(nGauss); gpar[d].resize(nGP);
gpar[d].fill(d == 0 ? lrspline->endparam(0) : lrspline->endparam(1)); 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; if (!this->getElementCoordinates(Xnod,iel)) return false;
// Get integration gauss points over this element // 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 --------------- // --- 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 // Parameter values of current integration point
double u = param[0] = gpar[0][j]; double u = param[0] = gpar[0][j];

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

@ -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, void ASMu3D::evaluateBasis (int iel, int basis, double u, double v, double w,
Vector& N, Matrix& dNdu) const Vector& N, Matrix& dNdu) const
{ {
PROFILE2("Spline evaluation"); PROFILE3("ASMu3D::evalBasis(1)");
std::vector<RealArray> result; std::vector<RealArray> result;
this->getBasis(basis)->computeBasis(u, v, w, result, 1, iel); 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, void ASMu3D::evaluateBasis (FiniteElement& fe, Matrix& dNdu,
const Matrix& C, const Matrix& B, int basis) const const Matrix& C, const Matrix& B, int basis) const
{ {
PROFILE2("BeSpline evaluation"); PROFILE3("ASMu3D::evalBasis(BE)");
Matrix N = C*B; Matrix N = C*B;
dNdu.resize(N.rows(),3); dNdu.resize(N.rows(),3);
@ -763,7 +763,7 @@ void ASMu3D::evaluateBasis (int iel, FiniteElement& fe,
Matrix& dNdu, Matrix3D& d2Ndu2, Matrix& dNdu, Matrix3D& d2Ndu2,
int basis) const int basis) const
{ {
PROFILE2("Spline evaluation"); PROFILE3("ASMu3D::evalBasis(2)");
std::vector<RealArray> result; std::vector<RealArray> result;
this->getBasis(basis)->computeBasis(fe.u, fe.v, fe.w, result, 2, iel); 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, void ASMu3D::evaluateBasis (int iel, FiniteElement& fe, int derivs,
int basis) const int basis) const
{ {
PROFILE2("Spline evaluation"); PROFILE3("ASMu2D::evalBasis");
std::vector<RealArray> result; std::vector<RealArray> result;
this->getBasis(basis)->computeBasis(fe.u, fe.v, fe.w, result, derivs, iel); 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)"); 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 p1 = lrspline->order(0);
int p2 = lrspline->order(1); int p2 = lrspline->order(1);
int p3 = lrspline->order(2); 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]; double u[2*p1], v[2*p2], w[2*p3];
Go::BsplineBasis basis1 = getBezierBasis(p1); Go::BsplineBasis basis1 = getBezierBasis(p1);
Go::BsplineBasis basis2 = getBezierBasis(p2); Go::BsplineBasis basis2 = getBezierBasis(p2);
Go::BsplineBasis basis3 = getBezierBasis(p3); Go::BsplineBasis basis3 = getBezierBasis(p3);
Matrix BN( p1*p2*p3, nGauss*nGauss*nGauss), rBN; Matrix BN( p1*p2*p3, nGP*nGP*nGP), rBN;
Matrix BdNdu(p1*p2*p3, nGauss*nGauss*nGauss), rBdNdu; Matrix BdNdu(p1*p2*p3, nGP*nGP*nGP), rBdNdu;
Matrix BdNdv(p1*p2*p3, nGauss*nGauss*nGauss), rBdNdv; Matrix BdNdv(p1*p2*p3, nGP*nGP*nGP), rBdNdv;
Matrix BdNdw(p1*p2*p3, nGauss*nGauss*nGauss), rBdNdw; Matrix BdNdw(p1*p2*p3, nGP*nGP*nGP), rBdNdw;
int ig = 1; // gauss point iterator int ig = 1; // gauss point iterator
for (int zeta = 0; zeta < nGauss; zeta++) for (int zeta = 0; zeta < nGP; zeta++)
for (int eta = 0; eta < nGauss; eta++) for (int eta = 0; eta < nGP; eta++)
for (int xi = 0; xi < nGauss; xi++, ig++) { for (int xi = 0; xi < nGP; xi++, ig++) {
basis1.computeBasisValues(xg[xi], u, 1); basis1.computeBasisValues(xg[xi], u, 1);
basis2.computeBasisValues(xg[eta], v, 1); basis2.computeBasisValues(xg[eta], v, 1);
basis3.computeBasisValues(xg[zeta], w, 1); basis3.computeBasisValues(xg[zeta], w, 1);
@ -869,7 +872,7 @@ bool ASMu3D::integrate (Integrand& integrand,
// Get the reduced integration quadrature points, if needed // Get the reduced integration quadrature points, if needed
const double* xr = nullptr; const double* xr = nullptr;
const double* wr = nullptr; const double* wr = nullptr;
int nRed = integrand.getReducedIntegration(nGauss); int nRed = integrand.getReducedIntegration(nGP);
if (nRed > 0) if (nRed > 0)
{ {
xr = GaussQuadrature::getCoord(nRed); xr = GaussQuadrature::getCoord(nRed);
@ -899,7 +902,7 @@ bool ASMu3D::integrate (Integrand& integrand,
} }
} }
else if (nRed < 0) else if (nRed < 0)
nRed = nGauss; // The integrand needs to know nGauss nRed = nGP; // The integrand needs to know nGauss
ThreadGroups oneGroup; ThreadGroups oneGroup;
if (glInt.threadSafe()) oneGroup.oneGroup(nel); if (glInt.threadSafe()) oneGroup.oneGroup(nel);
@ -958,7 +961,7 @@ bool ASMu3D::integrate (Integrand& integrand,
std::array<RealArray,3> gpar, redpar; std::array<RealArray,3> gpar, redpar;
for (int d = 0; d < 3; d++) 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) if (xr)
this->getGaussPointParameters(redpar[d],d,nRed,iel,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 ---------- // --- Integration loop over all Gauss points in each direction ----------
int ig = 1; 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 fe.iGP = firstIp + jp; // Global integration point counter
for (int k = 0; k < nGauss; k++) for (int k = 0; k < nGP; k++)
for (int j = 0; j < nGauss; j++) for (int j = 0; j < nGP; j++)
for (int i = 0; i < nGauss; i++, fe.iGP++, ig++) for (int i = 0; i < nGP; i++, fe.iGP++, ig++)
{ {
// Local element coordinates of current integration point // Local element coordinates of current integration point
fe.xi = xg[i]; fe.xi = xg[i];
@ -1150,18 +1153,21 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
PROFILE2("ASMu3D::integrate(B)"); 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} // 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 faceDir = (lIndex%10+1)/((lIndex%2) ? -2 : 2);
const int t1 = 1 + abs(faceDir)%3; // first tangent direction const int t1 = 1 + abs(faceDir)%3; // first tangent direction
const int t2 = 1 + t1%3; // second 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 // Extract the Neumann order flag (1 or higher) for the integrand
integrand.setNeumannOrder(1 + lIndex/10); integrand.setNeumannOrder(1 + lIndex/10);
@ -1230,7 +1236,7 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
double dXidu[3]; double dXidu[3];
// Get element face area in the parameter space // 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) if (dA < 0.0) // topology error (probably logic error)
{ {
ok = false; ok = false;
@ -1311,12 +1317,17 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
if (integrand.getIntegrandType() & Integrand::G_MATRIX) if (integrand.getIntegrandType() & Integrand::G_MATRIX)
utl::getGmat(Jac,dXidu,fe.G); 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 // Cartesian coordinates of current integration point
X.assign(Xnod * fe.N); X.assign(Xnod * fe.N);
X.t = time.t; X.t = time.t;
// Evaluate the integrand and accumulate element contributions // 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)) if (!integrand.evalBou(*A,fe,time,X,normal))
ok = false; ok = false;
} }
@ -1347,177 +1358,6 @@ bool ASMu3D::integrateEdge (Integrand& integrand, int lEdge,
{ {
std::cerr << "ASMu3D::integrateEdge(...) is not properly implemented yet :(" << std::endl; std::cerr << "ASMu3D::integrateEdge(...) is not properly implemented yet :(" << std::endl;
exit(776654); 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<Matrix,3> 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<Go::BasisDerivs> 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<char,size_t>::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; return false;
} }
@ -1678,6 +1518,8 @@ bool ASMu3D::evalSolution (Matrix& sField, const Vector& locSol,
bool ASMu3D::evalSolution (Matrix& sField, const Vector& locSol, bool ASMu3D::evalSolution (Matrix& sField, const Vector& locSol,
const RealArray* gpar, bool, int deriv, int) const const RealArray* gpar, bool, int deriv, int) const
{ {
PROFILE2("ASMu3D::evalSol(P)");
size_t nComp = locSol.size() / this->getNoNodes(); size_t nComp = locSol.size() / this->getNoNodes();
if (nComp*this->getNoNodes() != locSol.size()) if (nComp*this->getNoNodes() != locSol.size())
return false; return false;
@ -1693,6 +1535,7 @@ bool ASMu3D::evalSolution (Matrix& sField, const Vector& locSol,
Go::BasisPts spline0; Go::BasisPts spline0;
Go::BasisDerivs spline1; Go::BasisDerivs spline1;
Go::BasisDerivs2 spline2; Go::BasisDerivs2 spline2;
int lel = -1;
// Evaluate the primary solution field at each point // Evaluate the primary solution field at each point
sField.resize(nComp,nPoints); sField.resize(nComp,nPoints);
@ -1708,9 +1551,12 @@ bool ASMu3D::evalSolution (Matrix& sField, const Vector& locSol,
} }
utl::gather(MNPC[iel],nComp,locSol,eSol); utl::gather(MNPC[iel],nComp,locSol,eSol);
// Set up control point (nodal) coordinates for current element if (iel != lel && deriv > 0)
if (deriv > 0 && !this->getElementCoordinates(Xnod,iel+1)) {
return false; 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 // Evaluate basis function values/derivatives at current parametric point
// and multiply with control point values to get the point-wise solution // 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, bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
const RealArray* gpar, bool) const const RealArray* gpar, bool) const
{ {
#ifdef SP_DEBUG
std::cout <<"ASMu3D::evalSolution(Matrix&,const IntegrandBase&,const RealArray*,bool)\n";
#endif
sField.resize(0,0); 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 // TODO: investigate the possibility of doing "regular" refinement by
// uniform tesselation grid and ignoring LR mesh lines // uniform tesselation grid and ignoring LR mesh lines
size_t nPoints = gpar[0].size();
bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES; bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES;
if (nPoints != gpar[1].size() || nPoints != gpar[2].size())
return false;
Vector solPt; Vector solPt;
Matrix dNdu, Jac, Xnod; 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 Matrix B(p1*p2*p3, 4); // Bezier evaluation points and derivatives
// Evaluate the secondary solution field at each point // Evaluate the secondary solution field at each point
int lel = -1;
for (size_t i = 0; i < nPoints; i++) for (size_t i = 0; i < nPoints; i++)
{ {
// Fetch element containing evaluation point // Fetch element containing evaluation point
@ -1884,8 +1729,12 @@ bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
else else
this->evaluateBasis(fe, dNdu, bezierExtract[iel], B); this->evaluateBasis(fe, dNdu, bezierExtract[iel], B);
// Set up control point (nodal) coordinates for current element if (iel != lel)
if (!this->getElementCoordinates(Xnod,iel+1)) return false; {
lel = iel; // Set up control point (nodal) coordinates for current element
if (!this->getElementCoordinates(Xnod,iel+1))
return false;
}
// Compute the Jacobian inverse // Compute the Jacobian inverse
fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); 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)) if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
continue; continue;
#if SP_DEBUG > 4
if (1+iel == dbgElm || dbgElm == 0)
std::cout <<"\n"<< fe;
#endif
// Now evaluate the solution field // Now evaluate the solution field
if (!integrand.evalSol(solPt,fe,Xnod*fe.N,MNPC[iel])) if (!integrand.evalSol(solPt,fe,Xnod*fe.N,MNPC[iel]))
return false; return false;
@ -2289,7 +2143,7 @@ bool ASMu3D::refine (const RealFunc& refC, double refTol)
{ {
Go::Point X0; Go::Point X0;
int iel = 0; int iel = 0;
std::vector<int> elements; IntVec elements;
for (const LR::Element* elm : lrspline->getAllElements()) for (const LR::Element* elm : lrspline->getAllElements())
{ {
double u0 = 0.5*(elm->umin() + elm->umax()); double u0 = 0.5*(elm->umin() + elm->umax());

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

@ -30,9 +30,11 @@ bool ASMu3D::getGrevilleParameters (RealArray& prm, int dir, int basisNum) const
if (!this->getBasis(basisNum) || dir < 0 || dir > 2) return false; if (!this->getBasis(basisNum) || dir < 0 || dir > 2) return false;
const LR::LRSpline* lrspline = this->getBasis(basisNum); const LR::LRSpline* lrspline = this->getBasis(basisNum);
prm.clear(); prm.clear();
prm.reserve(lrspline->nBasisFunctions()); prm.reserve(lrspline->nBasisFunctions());
for (const LR::Basisfunction *b : lrspline->getAllBasisfunctions())
for (const LR::Basisfunction* b : lrspline->getAllBasisfunctions())
prm.push_back(b->getGrevilleParameter()[dir]); prm.push_back(b->getGrevilleParameter()[dir]);
return true; return true;
@ -109,15 +111,16 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
const int p1 = lrspline->order(0); const int p1 = lrspline->order(0);
const int p2 = lrspline->order(1); const int p2 = lrspline->order(1);
const int p3 = lrspline->order(2); const int p3 = lrspline->order(2);
const int pm = std::max(std::max(p1,p2),p3);
// Get Gaussian quadrature points // Get Gaussian quadrature points
const int ng1 = continuous ? nGauss : p1 - 1; const int ng1 = continuous ? this->getNoGaussPt(pm,true) : p1 - 1;
const int ng2 = continuous ? nGauss : p2 - 1; const int ng2 = continuous ? ng1 : p2 - 1;
const int ng3 = continuous ? nGauss : p3 - 1; const int ng3 = continuous ? ng1 : p3 - 1;
const double* xg = GaussQuadrature::getCoord(ng1); const double* xg = GaussQuadrature::getCoord(ng1);
const double* yg = GaussQuadrature::getCoord(ng2); const double* yg = GaussQuadrature::getCoord(ng2);
const double* zg = GaussQuadrature::getCoord(ng3); 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 (!xg || !yg || !zg) return false;
if (continuous && !wg) 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[0],0,ng1,iel,xg);
this->getGaussPointParameters(gpar[1],1,ng2,iel,yg); this->getGaussPointParameters(gpar[1],1,ng2,iel,yg);
this->getGaussPointParameters(gpar[2],2,ng3,iel,zg); this->getGaussPointParameters(gpar[2],2,ng3,iel,zg);
// convert to unstructured mesh representation
expandTensorGrid(gpar.data(),unstrGpar.data()); expandTensorGrid(gpar.data(),unstrGpar.data());
// Evaluate the secondary solution at all integration points // 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); eB[r-1].add(phi,sField(r,ip+1)*dJw);
} }
for (size_t i = 0; i < MNPC[iel-1].size(); ++i) { for (size_t i = 0; i < eA.rows(); ++i) {
for (size_t j = 0; j < MNPC[iel-1].size(); ++j) 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); A(MNPC[iel-1][i]+1, MNPC[iel-1][j]+1) += eA(i+1, j+1);
int jp = MNPC[iel-1][i]+1; int jp = MNPC[iel-1][i]+1;
@ -269,7 +270,7 @@ LR::LRSplineVolume* ASMu3D::scRecovery (const IntegrandBase& integrand) const
size_t ip = 0; size_t ip = 0;
std::vector<LR::Element*>::const_iterator elStart, elEnd, el; std::vector<LR::Element*>::const_iterator elStart, elEnd, el;
std::vector<LR::Element*> supportElements; std::vector<LR::Element*> supportElements;
for (LR::Basisfunction *b : lrspline->getAllBasisfunctions()) for (LR::Basisfunction* b : lrspline->getAllBasisfunctions())
{ {
#if SP_DEBUG > 2 #if SP_DEBUG > 2
std::cout <<"Basis: "<< *b <<"\n ng1 ="<< ng1 <<"\n ng2 ="<< ng2 <<"\n ng3 ="<< ng3 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) for (el = elStart; el != elEnd; ++el)
{ {
int iel = (**el).getId()+1; int iel = (**el).getId()+1;
#if SP_DEBUG > 2
std::cout <<"Element "<< **el << std::endl;
#endif
// evaluate all gauss points for this element // evaluate all gauss points for this element
std::array<RealArray,3> gaussPt, unstrGauss; std::array<RealArray,3> gaussPt, unstrGauss;
this->getGaussPointParameters(gaussPt[0],0,ng1,iel,xg); this->getGaussPointParameters(gaussPt[0],0,ng1,iel,xg);
this->getGaussPointParameters(gaussPt[1],1,ng2,iel,yg); this->getGaussPointParameters(gaussPt[1],1,ng2,iel,yg);
this->getGaussPointParameters(gaussPt[2],2,ng3,iel,zg); 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()); expandTensorGrid(gaussPt.data(),unstrGauss.data());
// Evaluate the secondary solution at all Gauss points // Evaluate the secondary solution at all Gauss points
@ -471,11 +469,6 @@ bool ASMu3D::faceL2projection (const DirichletFace& face,
StdVector B(n*m); StdVector B(n*m);
A.resize(n,n); 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 // find the normal direction for the face
int faceDir; int faceDir;
switch (face.edg) 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 t1 = 1 + abs(faceDir)%3; // first tangent direction
const int t2 = 1 + t1%3; // second 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; Vector N;
Matrix dNdu, dNdX, Xnod, Jac; 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 ===================== // === Assembly loop over all elements on the patch face =====================
for (size_t ie = 0; ie < face.MLGE.size(); ie++) // for all face elements for (size_t ie = 0; ie < face.MLGE.size(); ie++) // for all face elements
{ {
int iel = 1 + face.MLGE[ie]; int iel = 1 + face.MLGE[ie];
@ -513,7 +516,7 @@ bool ASMu3D::faceL2projection (const DirichletFace& face,
gpar[d].fill(lrspline->endparam(d)); gpar[d].fill(lrspline->endparam(d));
} }
else 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 // Get element face area in the parameter space
double dA = 0.25*this->getParametricArea(iel,abs(faceDir)); 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 // Set up control point coordinates for current element
if (!this->getElementCoordinates(Xnod,iel)) return false; if (!this->getElementCoordinates(Xnod,iel)) return false;
double u = gpar[0].front(); double u = param[0] = gpar[0].front();
double v = gpar[1].front(); double v = param[1] = gpar[1].front();
double w = gpar[2].front(); double w = param[2] = gpar[2].front();
// --- Integration loop over all Gauss points over the face ------------- // --- Integration loop over all Gauss points over the face ----------------
for (int j = 0; j < nGauss; j++)
for (int i = 0; i < nGauss; i++) for (int j = 0; j < nGP; j++)
for (int i = 0; i < nGP; i++)
{ {
// Parameter values of current integration point // Parameter values of current integration point
int k1, k2, k3; int k1, k2, k3;
@ -539,9 +543,9 @@ bool ASMu3D::faceL2projection (const DirichletFace& face,
case 3: k1 = i; k2 = j; k3 = 0; break; case 3: k1 = i; k2 = j; k3 = 0; break;
default: k1 = k2 = k3 = 0; default: k1 = k2 = k3 = 0;
} }
if (gpar[0].size() > 1) u = gpar[0](k1+1); if (gpar[0].size() > 1) u = param[0] = gpar[0](k1+1);
if (gpar[1].size() > 1) v = gpar[1](k2+1); if (gpar[1].size() > 1) v = param[1] = gpar[1](k2+1);
if (gpar[2].size() > 1) w = gpar[2](k3+1); if (gpar[2].size() > 1) w = param[2] = gpar[2](k3+1);
// Evaluate basis function derivatives at integration points // Evaluate basis function derivatives at integration points
this->evaluateBasis(iel-1, myGeoBasis, u, v, w, N, dNdu); 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 if (dJxW == 0.0) continue; // skip singular points
// Cartesian coordinates of current integration point // Cartesian coordinates of current integration point
X = Xnod * N; X.assign(Xnod * N);
X.t = time; X.t = time;
// For mixed basis, we need to compute functions separate from geometry // For mixed basis, we need to compute functions separate from geometry