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

Added: Error exit if integrand with reduced integration for mixed LR.

Browse Source 
Fixed: Solution transfer for mixed in 3D.
This commit is contained in:
Knut Morten Okstad 2018-09-20 17:43:26 +02:00
parent d8ac636fee
commit 58a1bf78d1
2 changed files with 150 additions and 150 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

93
src/ASM/LR/ASMu2Dmx.C
View File

@ -325,6 +325,14 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
const double* xg = GaussQuadrature::getCoord(nGauss); const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss); const double* wg = GaussQuadrature::getWeight(nGauss);
if (!xg || !wg) return false; if (!xg || !wg) return false;
if (integrand.getReducedIntegration(nGauss))
{
std::cerr <<" *** Reduced integration not available for mixed LR splines"
<< std::endl;
return false;
}
bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES; bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES;
ThreadGroups oneGroup; ThreadGroups oneGroup;
@ -333,9 +341,9 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
// === Assembly loop over all elements in the patch ========================== // === Assembly loop over all elements in the patch ==========================
bool ok = true; bool ok = true;
for (size_t t = 0; t < groups.size() && ok; ++t) for (size_t t = 0; t < groups.size() && ok; ++t)
{
#pragma omp parallel for schedule(static) #pragma omp parallel for schedule(static)
for (size_t e = 0; e < groups[t].size(); ++e) for (size_t e = 0; e < groups[t].size(); ++e)
{ {
@ -346,23 +354,23 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
std::vector<size_t> elem_sizes; std::vector<size_t> elem_sizes;
this->getElementsAt(threadBasis->getElement(groups[t][e])->midpoint(),els,elem_sizes); this->getElementsAt(threadBasis->getElement(groups[t][e])->midpoint(),els,elem_sizes);
int geoEl = els[geoBasis-1]; MxFiniteElement fe(elem_sizes);
std::vector<Matrix> dNxdu(m_basis.size());
MxFiniteElement fe(elem_sizes); std::vector<Matrix3D> d2Nxdu2(use2ndDer ? m_basis.size() : 0);
fe.iel = MLGE[geoEl-1];
std::vector<Matrix> dNxdu(m_basis.size());
Matrix Xnod, Jac; Matrix Xnod, Jac;
double param[3] = { 0.0, 0.0, 0.0 };
Vec4 X(param,time.t);
std::vector<Matrix3D> d2Nxdu2(m_basis.size());
Matrix3D Hess; Matrix3D Hess;
double dXidu[2]; double dXidu[2];
double param[3] = { 0.0, 0.0, 0.0 };
Vec4 X(param,time.t);
int geoEl = els[geoBasis-1];
fe.iel = MLGE[geoEl-1];
// Get element area in the parameter space // Get element area in the parameter space
double dA = this->getParametricArea(geoEl); double dA = 0.25*this->getParametricArea(geoEl);
if (dA < 0.0) // topology error (probably logic error) if (dA < 0.0)
{ {
ok = false; ok = false; // topology error (probably logic error)
continue; continue;
} }
@ -389,8 +397,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
this->getGaussPointParameters(gpar[d],d,nGauss,geoEl,xg); this->getGaussPointParameters(gpar[d],d,nGauss,geoEl,xg);
// Initialize element quantities // Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,false); LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel);
if (!integrand.initElement(MNPC[geoEl-1], fe, elem_sizes, nb, *A)) if (!integrand.initElement(MNPC[geoEl-1],fe,elem_sizes,nb,*A))
{ {
A->destruct(); A->destruct();
ok = false; ok = false;
@ -444,7 +452,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
X.assign(Xnod * fe.basis(geoBasis)); X.assign(Xnod * fe.basis(geoBasis));
// 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.evalIntMx(*A,fe,time,X)) if (!integrand.evalIntMx(*A,fe,time,X))
ok = false; ok = false;
} }
@ -459,7 +467,6 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
A->destruct(); A->destruct();
} }
}
return ok; return ok;
} }
@ -469,7 +476,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt, GlobalIntegral& glInt,
const TimeDomain& time) const TimeDomain& time)
{ {
if (!m_basis[0] || !m_basis[1]) if (m_basis.empty())
return true; // silently ignore empty patches return true; // silently ignore empty patches
PROFILE2("ASMu2Dmx::integrate(B)"); PROFILE2("ASMu2Dmx::integrate(B)");
@ -589,9 +596,9 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
dNxdu[geoBasis-1],t1,t2); dNxdu[geoBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < m_basis.size(); ++b) for (size_t b = 1; b <= m_basis.size(); ++b)
if (b != (size_t)geoBasis-1) if ((int)b != geoBasis)
fe.grad(b+1).multiply(dNxdu[b],Jac); fe.grad(b).multiply(dNxdu[b-1],Jac);
if (edgeDir < 0) if (edgeDir < 0)
normal *= -1.0; normal *= -1.0;
@ -626,8 +633,11 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
const TimeDomain& time, const TimeDomain& time,
const ASM::InterfaceChecker& iChkgen) const ASM::InterfaceChecker& iChkgen)
{ {
if (!geo) return true; // silently ignore empty patches if (m_basis.empty())
if (!(integrand.getIntegrandType() & Integrand::INTERFACE_TERMS)) return true; return true; // silently ignore empty patches
if (!(integrand.getIntegrandType() & Integrand::INTERFACE_TERMS))
return true; // No interface terms
PROFILE2("ASMu2Dmx::integrate(J)"); PROFILE2("ASMu2Dmx::integrate(J)");
@ -767,16 +777,18 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
// Compute Jacobian inverse of the coordinate mapping and // Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates // basis function derivatives w.r.t. Cartesian coordinates
fe.detJxW = utl::Jacobian(Jac2, normal, fe.detJxW = utl::Jacobian(Jac2,normal,
fe.grad(geoBasis+m_basis.size()), fe.grad(geoBasis+m_basis.size()),Xnod2,
Xnod2,dNxdu[geoBasis-1+m_basis.size()],t1,t2); dNxdu[geoBasis-1+m_basis.size()],t1,t2);
fe.detJxW = utl::Jacobian(Jac, normal, fe.detJxW = utl::Jacobian(Jac,normal,
fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1],t1,t2); fe.grad(geoBasis),Xnod,
dNxdu[geoBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < m_basis.size(); ++b)
if (b != (size_t)geoBasis-1) { for (size_t b = 1; b <= m_basis.size(); ++b)
fe.grad(b+1).multiply(dNxdu[b],Jac); if ((int)b != geoBasis) {
fe.grad(b+1+m_basis.size()).multiply(dNxdu[b+m_basis.size()],Jac); fe.grad(b).multiply(dNxdu[b-1],Jac);
fe.grad(b+m_basis.size()).multiply(dNxdu[b-1+m_basis.size()],Jac);
} }
if (edgeDir < 0) if (edgeDir < 0)
@ -844,6 +856,7 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const Vector& locSol,
for (size_t j=0; j < m_basis.size(); ++j) { for (size_t j=0; j < m_basis.size(); ++j) {
if (nc[j] == 0) if (nc[j] == 0)
continue; continue;
// Fetch element containing evaluation point. // Fetch element containing evaluation point.
// Sadly, points are not always ordered in the same way as the elements. // Sadly, points are not always ordered in the same way as the elements.
int iel = m_basis[j]->getElementContaining(gpar[0][i],gpar[1][i]); int iel = m_basis[j]->getElementContaining(gpar[0][i],gpar[1][i]);
@ -897,10 +910,10 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
this->getElementsAt({gpar[0][i],gpar[1][i]},els,elem_sizes); this->getElementsAt({gpar[0][i],gpar[1][i]},els,elem_sizes);
// Evaluate the basis functions at current parametric point // Evaluate the basis functions at current parametric point
MxFiniteElement fe(elem_sizes,firstIp+i); MxFiniteElement fe(elem_sizes,firstIp+i);
std::vector<Matrix> dNxdu(m_basis.size()); std::vector<Matrix> dNxdu(m_basis.size());
std::vector<Matrix3D> d2Nxdu2(use2ndDer ? m_basis.size() : 0);
Matrix Jac, Xnod; Matrix Jac, Xnod;
std::vector<Matrix3D> d2Nxdu2(m_basis.size());
Matrix3D Hess; Matrix3D Hess;
if (use2ndDer) if (use2ndDer)
for (size_t b = 0; b < m_basis.size(); ++b) { for (size_t b = 0; b < m_basis.size(); ++b) {
@ -927,14 +940,13 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis)) if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis))
return false; return false;
// Now evaluate the solution field
Vector solPt;
// Cartesian coordinates of current integration point // Cartesian coordinates of current integration point
fe.u = gpar[0][i]; fe.u = gpar[0][i];
fe.v = gpar[1][i]; fe.v = gpar[1][i];
utl::Point X4(Xnod*fe.basis(geoBasis),{fe.u,fe.v}); utl::Point X4(Xnod*fe.basis(geoBasis),{fe.u,fe.v});
// Now evaluate the solution field
Vector solPt;
if (!integrand.evalSol(solPt,fe,X4, if (!integrand.evalSol(solPt,fe,X4,
MNPC[els[geoBasis-1]-1],elem_sizes,nb)) MNPC[els[geoBasis-1]-1],elem_sizes,nb))
return false; return false;
@ -956,13 +968,12 @@ bool ASMu2Dmx::refine (const LR::RefineData& prm, Vectors& sol)
if (prm.errors.empty() && prm.elements.empty()) if (prm.errors.empty() && prm.elements.empty())
return true; return true;
for (Vector& solvec : sol) { for (Vector& solvec : sol)
for (size_t j = 0; j < m_basis.size(); j++) { for (size_t j = 0; j < m_basis.size(); j++) {
Vector bVec; Vector bVec;
this->extractNodeVec(solvec, bVec, 0, j+1); this->extractNodeVec(solvec, bVec, 0, j+1);
LR::extendControlPoints(m_basis[j].get(), bVec, nfx[j]); LR::extendControlPoints(m_basis[j].get(), bVec, nfx[j]);
} }
}
if (doRefine(prm, refBasis.get())) { if (doRefine(prm, refBasis.get())) {
for (size_t j = 0; j < m_basis.size(); ++j) for (size_t j = 0; j < m_basis.size(); ++j)
@ -1118,8 +1129,8 @@ void ASMu2Dmx::getBoundaryNodes (int lIndex, IntVec& nodes, int basis,
} }
void ASMu2Dmx::remapErrors(RealArray& errors, void ASMu2Dmx::remapErrors (RealArray& errors,
const RealArray& origErr, bool elemErrors) const const RealArray& origErr, bool elemErrors) const
{ {
const LR::LRSplineSurface* geo = this->getBasis(ASMmxBase::geoBasis); const LR::LRSplineSurface* geo = this->getBasis(ASMmxBase::geoBasis);
for (const LR::Element* elm : geo->getAllElements()) { for (const LR::Element* elm : geo->getAllElements()) {

207
src/ASM/LR/ASMu3Dmx.C
View File

@ -40,6 +40,7 @@ ASMu3Dmx::ASMu3Dmx (const CharVec& n_f)
: ASMu3D(std::accumulate(n_f.begin(), n_f.end(), 0)), ASMmxBase(n_f), : ASMu3D(std::accumulate(n_f.begin(), n_f.end(), 0)), ASMmxBase(n_f),
bezierExtractmx(myBezierExtractmx) bezierExtractmx(myBezierExtractmx)
{ {
threadBasis = nullptr;
myGeoBasis = ASMmxBase::geoBasis; myGeoBasis = ASMmxBase::geoBasis;
} }
@ -49,6 +50,7 @@ ASMu3Dmx::ASMu3Dmx (const ASMu3Dmx& patch, const CharVec& n_f)
m_basis(patch.m_basis), m_basis(patch.m_basis),
bezierExtractmx(patch.myBezierExtractmx) bezierExtractmx(patch.myBezierExtractmx)
{ {
threadBasis = patch.threadBasis;
nfx = patch.nfx; nfx = patch.nfx;
nb = patch.nb; nb = patch.nb;
myGeoBasis = ASMmxBase::geoBasis; myGeoBasis = ASMmxBase::geoBasis;
@ -199,6 +201,9 @@ bool ASMu3Dmx::getSolution (Matrix& sField, const Vector& locSol,
bool ASMu3Dmx::generateFEMTopology () bool ASMu3Dmx::generateFEMTopology ()
{ {
if (!myMLGN.empty())
return true;
if (m_basis.empty()) { if (m_basis.empty()) {
VolumeVec vvec = ASMmxBase::establishBases(tensorspline, ASMmxBase::Type); VolumeVec vvec = ASMmxBase::establishBases(tensorspline, ASMmxBase::Type);
m_basis.resize(vvec.size()); m_basis.resize(vvec.size());
@ -314,6 +319,9 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
GlobalIntegral& glInt, GlobalIntegral& glInt,
const TimeDomain& time) const TimeDomain& time)
{ {
if (m_basis.empty())
return true; // silently ignore empty patches
PROFILE2("ASMu3Dmx::integrate(I)"); PROFILE2("ASMu3Dmx::integrate(I)");
// Get Gaussian quadrature points and weights // Get Gaussian quadrature points and weights
@ -321,6 +329,15 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
const double* wg = GaussQuadrature::getWeight(nGauss); const double* wg = GaussQuadrature::getWeight(nGauss);
if (!xg || !wg) return false; if (!xg || !wg) return false;
if (integrand.getReducedIntegration(nGauss))
{
std::cerr <<" *** Reduced integration not available for mixed LR splines"
<< std::endl;
return false;
}
bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES;
// evaluate all gauss points on the bezier patch (-1, 1) // evaluate all gauss points on the bezier patch (-1, 1)
std::vector<Matrix> BN(m_basis.size()); std::vector<Matrix> BN(m_basis.size());
std::vector<Matrix> BdNdu(m_basis.size()); std::vector<Matrix> BdNdu(m_basis.size());
@ -331,55 +348,33 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
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);
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);
BN [b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss);
BN[b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss);
BdNdu[b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss); BdNdu[b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss);
BdNdv[b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss); BdNdv[b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss);
BdNdw[b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss); BdNdw[b-1].resize(p1*p2*p3, nGauss*nGauss*nGauss);
int ig=1; // gauss point iterator int ig = 1; // gauss point iterator
for(int zeta=0; zeta<nGauss; zeta++) { for (int zeta = 0; zeta < nGauss; zeta++)
for(int eta=0; eta<nGauss; eta++) { for (int eta = 0; eta < nGauss; eta++)
for(int xi=0; xi<nGauss; xi++, ig++) { for (int xi = 0; xi < nGauss; xi++, ig++) {
double u[2*p1];
double v[2*p2];
double w[2*p3];
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);
int ib=1; // basis function iterator int ib = 1; // basis function iterator
double sum = 0; for (int k = 0; k < p3; k++)
for(int k=0; k<p3; k++) { for (int j = 0; j < p2; j++)
for(int j=0; j<p2; j++) { for (int i = 0; i < p1; i++, ib++) {
for(int i=0; i<p1; i++, ib++) {
BN[b-1](ib,ig) = u[2*i ]*v[2*j ]*w[2*k ]; BN[b-1](ib,ig) = u[2*i ]*v[2*j ]*w[2*k ];
BdNdu[b-1](ib,ig) = u[2*i+1]*v[2*j ]*w[2*k ]; BdNdu[b-1](ib,ig) = u[2*i+1]*v[2*j ]*w[2*k ];
BdNdv[b-1](ib,ig) = u[2*i ]*v[2*j+1]*w[2*k ]; BdNdv[b-1](ib,ig) = u[2*i ]*v[2*j+1]*w[2*k ];
BdNdw[b-1](ib,ig) = u[2*i ]*v[2*j ]*w[2*k+1]; BdNdw[b-1](ib,ig) = u[2*i ]*v[2*j ]*w[2*k+1];
sum += BN[b-1](ib,ig);
} }
}
}
} }
}
}
} }
// Get the reduced integration quadrature points, if needed
const double* xr = nullptr;
const double* wr = nullptr;
int nRed = integrand.getReducedIntegration(nGauss);
if (nRed > 0)
{
xr = GaussQuadrature::getCoord(nRed);
wr = GaussQuadrature::getWeight(nRed);
if (!xr || !wr) return false;
}
else if (nRed < 0)
nRed = nGauss; // The integrand needs to know nGauss
ThreadGroups oneGroup; ThreadGroups oneGroup;
if (glInt.threadSafe()) oneGroup.oneGroup(nel); if (glInt.threadSafe()) oneGroup.oneGroup(nel);
const IntMat& groups = glInt.threadSafe() ? oneGroup[0] : threadGroups[0]; const IntMat& groups = glInt.threadSafe() ? oneGroup[0] : threadGroups[0];
@ -400,17 +395,18 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
std::vector<int> els; std::vector<int> els;
std::vector<size_t> elem_sizes; std::vector<size_t> elem_sizes;
this->getElementsAt(el->midpoint(),els,elem_sizes); this->getElementsAt(el->midpoint(),els,elem_sizes);
int iEl = el->getId(); MxFiniteElement fe(elem_sizes);
MxFiniteElement fe(elem_sizes); std::vector<Matrix> dNxdu(m_basis.size());
fe.iel = MLGE[iEl]; std::vector<Matrix3D> d2Nxdu2(use2ndDer ? m_basis.size() : 0);
std::vector<Matrix> dNxdu(m_basis.size());
Matrix Xnod, Jac; Matrix Xnod, Jac;
std::vector<Matrix3D> d2Nxdu2(m_basis.size());
Matrix3D Hess; Matrix3D Hess;
double dXidu[3]; double dXidu[3];
double param[3] = { 0.0, 0.0, 0.0 }; double param[3] = { 0.0, 0.0, 0.0 };
Vec4 X(param,time.t); Vec4 X(param,time.t);
int iEl = el->getId();
fe.iel = MLGE[iEl];
// Get element volume in the parameter space // Get element volume in the parameter space
double du = el->umax() - el->umin(); double du = el->umax() - el->umin();
double dv = el->vmax() - el->vmin(); double dv = el->vmax() - el->vmin();
@ -430,14 +426,9 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
} }
// Compute parameter values of the Gauss points over the whole element // Compute parameter values of the Gauss points over the whole element
std::array<RealArray,3> gpar, redpar; std::array<RealArray,3> gpar;
for (int d = 0; d < 3; d++) for (int d = 0; d < 3; d++)
{
this->getGaussPointParameters(gpar[d],d,nGauss,iEl+1,xg); this->getGaussPointParameters(gpar[d],d,nGauss,iEl+1,xg);
if (xr)
this->getGaussPointParameters(redpar[d],d,nRed,iEl+1,xr);
}
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS) if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
fe.h = this->getElementCorners(iEl+1, fe.XC); fe.h = this->getElementCorners(iEl+1, fe.XC);
@ -498,12 +489,6 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
continue; continue;
} }
if (xr)
{
std::cerr << "Haven't really figured out what this part does yet\n";
exit(42142);
}
// --- Integration loop over all Gauss points in each direction ---------- // --- Integration loop over all Gauss points in each direction ----------
int jp = iEl*nGauss*nGauss*nGauss; int jp = iEl*nGauss*nGauss*nGauss;
@ -534,7 +519,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
B.fillColumn(4, BdNdw[b].getColumn(ig)*2.0/dw); B.fillColumn(4, BdNdw[b].getColumn(ig)*2.0/dw);
// Fetch basis function derivatives at current integration point // Fetch basis function derivatives at current integration point
if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES) if (use2ndDer)
this->evaluateBasis(els[b]-1, fe, dNxdu[b], d2Nxdu2[b], b+1); this->evaluateBasis(els[b]-1, fe, dNxdu[b], d2Nxdu2[b], b+1);
else else
this->evaluateBasis(fe, dNxdu[b], bezierExtractmx[b][els[b]-1], B, b+1); this->evaluateBasis(fe, dNxdu[b], bezierExtractmx[b][els[b]-1], B, b+1);
@ -560,8 +545,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
fe.detJxW *= 0.125*vol*wg[i]*wg[j]*wg[k]; fe.detJxW *= 0.125*vol*wg[i]*wg[j]*wg[k];
if (!integrand.evalIntMx(*A,fe,time,X)) if (!integrand.evalIntMx(*A,fe,time,X))
ok = false; ok = false;
}
} // end gauss integrand
// Finalize the element quantities // Finalize the element quantities
if (ok && !integrand.finalizeElement(*A,fe,time,firstIp+jp)) if (ok && !integrand.finalizeElement(*A,fe,time,firstIp+jp))
@ -582,6 +566,9 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
GlobalIntegral& glInt, GlobalIntegral& glInt,
const TimeDomain& time) const TimeDomain& time)
{ {
if (m_basis.empty())
return true; // silently ignore empty patches
PROFILE2("ASMu3Dmx::integrate(B)"); PROFILE2("ASMu3Dmx::integrate(B)");
// Get Gaussian quadrature points and weights // Get Gaussian quadrature points and weights
@ -663,7 +650,7 @@ bool ASMu3Dmx::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) return false; // topology error (probably logic error) if (dA < 0.0) return false; // topology error (probably logic error)
// Set up control point coordinates for current element // Set up control point coordinates for current element
@ -728,9 +715,10 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
fe.detJxW = utl::Jacobian(Jac, normal, fe.grad(geoBasis), Xnod, fe.detJxW = utl::Jacobian(Jac, normal, fe.grad(geoBasis), Xnod,
dNxdu[geoBasis-1], t1, t2); dNxdu[geoBasis-1], t1, t2);
if (fe.detJxW == 0.0) continue; // skip singular points if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < m_basis.size(); ++b)
if (b != (size_t)geoBasis-1) for (size_t b = 1; b <= m_basis.size(); ++b)
fe.grad(b+1).multiply(dNxdu[b],Jac); if ((int)b != geoBasis)
fe.grad(b).multiply(dNxdu[b-1],Jac);
if (faceDir < 0) normal *= -1.0; if (faceDir < 0) normal *= -1.0;
@ -742,7 +730,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
X.assign(Xnod * fe.basis(geoBasis)); X.assign(Xnod * fe.basis(geoBasis));
// 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];
ok = integrand.evalBouMx(*A,fe,time,X,normal); ok = integrand.evalBouMx(*A,fe,time,X,normal);
} }
@ -792,6 +780,7 @@ bool ASMu3Dmx::evalSolution (Matrix& sField, const Vector& locSol,
for (size_t j=0; j < m_basis.size(); ++j) { for (size_t j=0; j < m_basis.size(); ++j) {
if (nc[j] == 0) if (nc[j] == 0)
continue; continue;
// Fetch element containing evaluation point. // Fetch element containing evaluation point.
// Sadly, points are not always ordered in the same way as the elements. // Sadly, points are not always ordered in the same way as the elements.
int iel = m_basis[j]->getElementContaining(gpar[0][i],gpar[1][i],gpar[2][i]); int iel = m_basis[j]->getElementContaining(gpar[0][i],gpar[1][i],gpar[2][i]);
@ -845,10 +834,10 @@ bool ASMu3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
this->getElementsAt({gpar[0][i],gpar[1][i],gpar[2][i]},els,elem_sizes); this->getElementsAt({gpar[0][i],gpar[1][i],gpar[2][i]},els,elem_sizes);
// Evaluate the basis functions at current parametric point // Evaluate the basis functions at current parametric point
MxFiniteElement fe(elem_sizes,firstIp+i); MxFiniteElement fe(elem_sizes,firstIp+i);
std::vector<Matrix> dNxdu(m_basis.size()); std::vector<Matrix> dNxdu(m_basis.size());
std::vector<Matrix3D> d2Nxdu2(use2ndDer ? m_basis.size() : 0);
Matrix Jac, Xnod; Matrix Jac, Xnod;
std::vector<Matrix3D> d2Nxdu2(m_basis.size());
Matrix3D Hess; Matrix3D Hess;
if (use2ndDer) if (use2ndDer)
for (size_t b = 0; b < m_basis.size(); ++b) { for (size_t b = 0; b < m_basis.size(); ++b) {
@ -904,13 +893,12 @@ bool ASMu3Dmx::refine (const LR::RefineData& prm, Vectors& sol)
if (prm.errors.empty() && prm.elements.empty()) if (prm.errors.empty() && prm.elements.empty())
return true; return true;
for (Vector& solvec : sol) { for (Vector& solvec : sol)
size_t ofs = 0;
for (size_t j = 0; j < m_basis.size(); j++) { for (size_t j = 0; j < m_basis.size(); j++) {
LR::extendControlPoints(m_basis[j].get(), solvec, nfx[j], ofs); Vector bVec;
ofs += nfx[j]*nb[j]; this->extractNodeVec(solvec, bVec, 0, j+1);
LR::extendControlPoints(m_basis[j].get(), bVec, nfx[j]);
} }
}
if (doRefine(prm, refBasis.get())) { if (doRefine(prm, refBasis.get())) {
for (size_t j = 0; j < m_basis.size(); ++j) for (size_t j = 0; j < m_basis.size(); ++j)
@ -943,7 +931,7 @@ bool ASMu3Dmx::refine (const LR::RefineData& prm, Vectors& sol)
} }
size_t ofs = 0; size_t ofs = 0;
for (int i = sol.size()-1; i > 0; i--) for (int i = sol.size()-1; i >= 0; i--)
for (size_t j = 0; j < m_basis.size(); ++j) { for (size_t j = 0; j < m_basis.size(); ++j) {
sol[i].resize(len); sol[i].resize(len);
LR::contractControlPoints(m_basis[j].get(), sol[i], nfx[j], ofs); LR::contractControlPoints(m_basis[j].get(), sol[i], nfx[j], ofs);
@ -992,48 +980,6 @@ Vec3 ASMu3Dmx::getCoord (size_t inod) const
} }
void ASMu3Dmx::remapErrors (RealArray& errors,
const RealArray& origErr, bool elemErrors) const
{
const LR::LRSplineVolume* geo = this->getBasis(ASMmxBase::geoBasis);
for (const LR::Element* elm : geo->getAllElements()) {
int rEl = refBasis->getElementContaining((elm->umin()+elm->umax())/2.0,
(elm->vmin()+elm->vmax())/2.0,
(elm->wmin()+elm->wmax())/2.0);
if (elemErrors)
errors[rEl] += origErr[elm->getId()];
else
for (LR::Basisfunction* b : refBasis->getElement(rEl)->support())
errors[b->getId()] += origErr[elm->getId()];
}
}
size_t ASMu3Dmx::getNoRefineNodes() const
{
return refBasis->nBasisFunctions();
}
size_t ASMu3Dmx::getNoRefineElms() const
{
return refBasis->nElements();
}
void ASMu3Dmx::copyRefinement(LR::LRSplineVolume* basis, int multiplicity) const
{
for (const LR::MeshRectangle* rect : refBasis->getAllMeshRectangles()) {
int mult = rect->multiplicity_ > 1 ? basis->order(rect->constDirection())
: multiplicity;
LR::MeshRectangle* newRect = rect->copy();
newRect->multiplicity_ = mult;
basis->insert_line(newRect);
}
}
void ASMu3Dmx::generateThreadGroups (const Integrand& integrand, bool silence, void ASMu3Dmx::generateThreadGroups (const Integrand& integrand, bool silence,
bool ignoreGlobalLM) bool ignoreGlobalLM)
{ {
@ -1081,6 +1027,49 @@ void ASMu3Dmx::generateThreadGroups (const Integrand& integrand, bool silence,
} }
void ASMu3Dmx::remapErrors (RealArray& errors,
const RealArray& origErr, bool elemErrors) const
{
const LR::LRSplineVolume* geo = this->getBasis(ASMmxBase::geoBasis);
for (const LR::Element* elm : geo->getAllElements()) {
int rEl = refBasis->getElementContaining((elm->umin()+elm->umax())/2.0,
(elm->vmin()+elm->vmax())/2.0,
(elm->wmin()+elm->wmax())/2.0);
if (elemErrors)
errors[rEl] += origErr[elm->getId()];
else
for (LR::Basisfunction* b : refBasis->getElement(rEl)->support())
errors[b->getId()] += origErr[elm->getId()];
}
}
size_t ASMu3Dmx::getNoRefineNodes() const
{
return refBasis->nBasisFunctions();
}
size_t ASMu3Dmx::getNoRefineElms() const
{
return refBasis->nElements();
}
void ASMu3Dmx::copyRefinement (LR::LRSplineVolume* basis,
int multiplicity) const
{
for (const LR::MeshRectangle* rect : refBasis->getAllMeshRectangles()) {
int mult = rect->multiplicity_ > 1 ? basis->order(rect->constDirection())
: multiplicity;
LR::MeshRectangle* newRect = rect->copy();
newRect->multiplicity_ = mult;
basis->insert_line(newRect);
}
}
void ASMu3Dmx::swapProjectionBasis () void ASMu3Dmx::swapProjectionBasis ()
{ {
if (altProjBasis) { if (altProjBasis) {