diff --git a/src/ASM/ASMs2Dmx.C b/src/ASM/ASMs2Dmx.C index 3e12400d5..01903a0ac 100644 --- a/src/ASM/ASMs2Dmx.C +++ b/src/ASM/ASMs2Dmx.C @@ -215,8 +215,11 @@ bool ASMs2Dmx::generateFEMTopology () else if (ASMmxBase::Type == ASMmxBase::SUBGRID) { projB = proj = m_basis.front()->clone(); altProjBasis = ASMmxBase::raiseBasis(surf); - } else - projB = proj = m_basis[2-ASMmxBase::geoBasis]->clone(); + } + else if (geoBasis < 3) + projB = proj = m_basis[2-geoBasis]->clone(); + else + return false; // Logic error } delete surf; geomB = surf = m_basis[geoBasis-1]->clone(); @@ -657,25 +660,12 @@ bool ASMs2Dmx::integrate (Integrand& integrand, // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod, - dNxdu[geoBasis-1]); - 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) - fe.grad(b+1).multiply(dNxdu[b],Jac); + if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives - if (use2ndDer) { - if (!utl::Hessian(Hess,fe.hess(geoBasis),Jac,Xnod, - d2Nxdu2[geoBasis-1],fe.grad(geoBasis),true)) - ok = false; - for (size_t b = 0; b < m_basis.size() && ok; ++b) - if ((int)b != geoBasis) - if (!utl::Hessian(Hess,fe.hess(b+1),Jac,Xnod, - d2Nxdu2[b],fe.grad(b+1),false)) - ok = false; - } + if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis)) + ok = false; // Compute G-matrix if (integrand.getIntegrandType() & Integrand::G_MATRIX) @@ -1197,12 +1187,8 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xtmp,dNxdu[geoBasis-1]); - if (fe.detJxW == 0.0) continue; // skip singular points - - for (size_t b = 1; b <= m_basis.size(); b++) - if (b != (size_t)geoBasis) - fe.grad(b).multiply(dNxdu[b-1],Jac); + if (!fe.Jacobian(Jac,Xtmp,dNxdu,geoBasis)) + continue; // skip singular points // Cartesian coordinates of current integration point utl::Point X4(Xtmp*fe.basis(geoBasis),{fe.u,fe.v}); diff --git a/src/ASM/ASMs2DmxLag.C b/src/ASM/ASMs2DmxLag.C index b06581152..5bc280c8e 100644 --- a/src/ASM/ASMs2DmxLag.C +++ b/src/ASM/ASMs2DmxLag.C @@ -311,12 +311,8 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, ok = false; // Compute Jacobian inverse of coordinate mapping and derivatives - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod, - dNxdu[geoBasis-1]); - if (fe.detJxW == 0.0) continue; // skip singular points - for (size_t b = 0; b < nxx.size(); ++b) - if (b != (size_t)geoBasis-1) - fe.grad(b+1).multiply(dNxdu[b],Jac); + if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + continue; // skip singular points // Cartesian coordinates of current integration point X.assign(Xnod * fe.basis(geoBasis)); @@ -423,10 +419,11 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex, elem_sizes[b][1],xi[1])) ok = false; - // Compute basis function derivatives and the edge normal - fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod, + // Compute basis function derivatives and the edge normal + fe.detJxW = utl::Jacobian(Jac,normal,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 < nxx.size(); ++b) if (b != (size_t)geoBasis-1) fe.grad(b+1).multiply(dNxdu[b],Jac); @@ -526,13 +523,10 @@ bool ASMs2DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand, elem_sizes[b][1],eta)) return false; - // Compute the Jacobian inverse - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1]); - if (fe.detJxW == 0.0) continue; // skip singular points - - for (size_t b = 1; b <= nxx.size(); b++) - if (b != (size_t)geoBasis) - fe.grad(b).multiply(dNxdu[b-1],Jac); + // Compute Jacobian inverse of the coordinate mapping and + // basis function derivatives w.r.t. Cartesian coordinates + if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + continue; // skip singular points // Now evaluate the solution field if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis), diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index 1bd6337d7..ade7392f9 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/src/ASM/ASMs3Dmx.C @@ -216,8 +216,11 @@ bool ASMs3Dmx::generateFEMTopology () else if (ASMmxBase::Type == ASMmxBase::SUBGRID) { projB = proj = m_basis.front()->clone(); altProjBasis = ASMmxBase::raiseBasis(svol); - } else - projB = proj = m_basis[2-ASMmxBase::geoBasis]->clone(); + } + else if (geoBasis < 3) + projB = proj = m_basis[2-geoBasis]->clone(); + else + return false; // Logic error } delete svol; geomB = svol = m_basis[geoBasis-1]->clone(); @@ -620,25 +623,12 @@ bool ASMs3Dmx::integrate (Integrand& integrand, // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod, - dNxdu[geoBasis-1]); - 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) - fe.grad(b+1).multiply(dNxdu[b],Jac); + if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives - if (use2ndDer) { - if (!utl::Hessian(Hess,fe.hess(geoBasis),Jac,Xnod, - d2Nxdu2[geoBasis-1],fe.grad(geoBasis),true)) - ok = false; - for (size_t b = 0; b < m_basis.size() && ok; ++b) - if ((int)b != geoBasis) - if (!utl::Hessian(Hess,fe.hess(b+1),Jac,Xnod, - d2Nxdu2[b],fe.grad(b+1),false)) - ok = false; - } + if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis)) + ok = false; // Compute G-matrix if (integrand.getIntegrandType() & Integrand::G_MATRIX) @@ -1246,13 +1236,9 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, SplineUtils::extractBasis(splinex[b][i],fe.basis(b+1),dNxdu[b]); // Compute Jacobian inverse of the coordinate mapping and - // basis function derivatives w.r.t. Cartesian coordinates - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xtmp,dNxdu[geoBasis-1]); - if (fe.detJxW == 0.0) continue; // skip singular points - - for (size_t b = 1; b <= m_basis.size(); b++) - if (b != (size_t)geoBasis) - fe.grad(b).multiply(dNxdu[b-1],Jac); + // basis function derivatives w.r.t. Cartesian coordinate + if (!fe.Jacobian(Jac,Xtmp,dNxdu,geoBasis)) + continue; // skip singular points // Cartesian coordinates of current integration point utl::Point X4(Xtmp * fe.basis(geoBasis),{fe.u,fe.v,fe.w}); diff --git a/src/ASM/ASMs3DmxLag.C b/src/ASM/ASMs3DmxLag.C index bc784602d..7ca9da763 100644 --- a/src/ASM/ASMs3DmxLag.C +++ b/src/ASM/ASMs3DmxLag.C @@ -348,12 +348,8 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, ok = false; // Compute Jacobian inverse of coordinate mapping and derivatives - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod, - dNxdu[geoBasis-1]); - if (fe.detJxW == 0.0) continue; // skip singular points - for (size_t b = 0; b < nxx.size(); ++b) - if (b != (size_t)geoBasis-1) - fe.grad(b+1).multiply(dNxdu[b],Jac); + if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + continue; // skip singular points // Cartesian coordinates of current integration point X.assign(Xnod * fe.basis(geoBasis)); @@ -494,6 +490,7 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex, fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(geoBasis),Xnod, dNxdu[geoBasis-1],t1,t2); if (fe.detJxW == 0.0) continue; // skip singular points + for (size_t b = 0; b < nxx.size(); ++b) if (b != (size_t)geoBasis-1) fe.grad(b+1).multiply(dNxdu[b],Jac); @@ -597,14 +594,10 @@ bool ASMs3DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand, elem_sizes[b][2],fe.zeta)) return false; - // Compute the Jacobian inverse - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod, - dNxdu[geoBasis-1]); - if (fe.detJxW == 0.0) continue; // skip singular points - - for (size_t b = 1; b <= nxx.size(); b++) - if (b != (size_t)geoBasis) - fe.grad(b).multiply(dNxdu[b-1],Jac); + // Compute Jacobian inverse of the coordinate mapping and + // basis function derivatives w.r.t. Cartesian coordinates + if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + continue; // skip singular points // Now evaluate the solution field if (!integrand.evalSol(solPt,fe,Xnod*fe.basis(geoBasis), diff --git a/src/ASM/FiniteElement.C b/src/ASM/FiniteElement.C index 848869405..e97c83d92 100644 --- a/src/ASM/FiniteElement.C +++ b/src/ASM/FiniteElement.C @@ -12,6 +12,7 @@ //============================================================================== #include "FiniteElement.h" +#include "CoordinateMapping.h" #include "Vec3Oper.h" @@ -65,3 +66,58 @@ std::ostream& MxFiniteElement::write (std::ostream& os) const } return os; } + + +/*! + This method also calculates the first-derivatives of the basis functions + with respect to the Cartesian coordinates, using the same geometry mapping + for all bases. +*/ + +bool MxFiniteElement::Jacobian (Matrix& Jac, const Matrix& Xnod, + const std::vector& dNxdu, + unsigned short int gBasis) +{ + detJxW = utl::Jacobian(Jac, + gBasis > 1 ? dMdX[gBasis-2] : dNdX, + Xnod, dNxdu[gBasis-1]); + if (detJxW == 0.0) return false; // singular point + + if (gBasis > 1) + dNdX.multiply(dNxdu.front(),Jac); + + for (size_t basis = 2; basis <= dNxdu.size(); basis++) + if (basis != gBasis) + dMdX[basis-2].multiply(dNxdu[basis-1],Jac); + + return true; +} + + +/*! + This method also calculates the second-derivatives of the basis functions + with respect to the Cartesian coordinates, using the same geometry mapping + for all bases. +*/ + +bool MxFiniteElement::Hessian (Matrix3D& Hess, const Matrix& Jac, + const Matrix& Xnod, + const std::vector& d2Nxdu2, + unsigned short int gBasis) +{ + bool ok = utl::Hessian(Hess, + gBasis > 1 ? d2MdX2[gBasis-2] : d2NdX2, + Jac, Xnod, d2Nxdu2[gBasis-1], + gBasis > 1 ? dMdX[gBasis-2] : dNdX); + + if (ok && gBasis > 1) + ok = utl::Hessian(Hess, d2NdX2, Jac, Xnod, + d2Nxdu2.front(), dNdX, false); + + for (size_t basis = 2; basis <= d2Nxdu2.size() && ok; basis++) + if (basis != gBasis) + ok = utl::Hessian(Hess, d2MdX2[basis-2], Jac, Xnod, + d2Nxdu2[basis-1], dMdX[basis-2], false); + + return ok; +} diff --git a/src/ASM/FiniteElement.h b/src/ASM/FiniteElement.h index 469b2505c..847d00d33 100644 --- a/src/ASM/FiniteElement.h +++ b/src/ASM/FiniteElement.h @@ -49,13 +49,8 @@ class FiniteElement : public ItgPoint //! \brief Returns a const reference to the basis function 2nd-derivatives. virtual const Matrix3D& hess(char) const { return d2NdX2; } - //! \brief Returns a reference to the basis function 2nd-derivatives. - virtual Matrix3D& hess(char) { return d2NdX2; } - //! \brief Returns a const reference to the basis function 3nd-derivatives. virtual const Matrix4D& hess2(char) const { return d3NdX3; } - //! \brief Returns a reference to the basis function 3nd-derivatives. - virtual Matrix4D& hess2(char) { return d3NdX3; } protected: //! \brief Writes the finite element object to the given output stream. @@ -118,13 +113,25 @@ class MxFiniteElement : public FiniteElement //! \brief Returns a const reference to the basis function 2nd-derivatives. virtual const Matrix3D& hess(char b) const { return b < 2 ? d2NdX2 : d2MdX2[b-2]; } - //! \brief Returns a reference to the basis function 2nd-derivatives. - virtual Matrix3D& hess(char b) { return b < 2 ? d2NdX2 : d2MdX2[b-2]; } - //! \brief Returns a const reference to the basis function 3rd-derivatives. virtual const Matrix4D& hess2(char b) const { return b < 2 ? d3NdX3 : d3MdX3[b-2]; } - //! \brief Returns a reference to the basis function 3rd-derivatives. - virtual Matrix4D& hess2(char b) { return b < 2 ? d3NdX3 : d3MdX3[b-2]; } + + //! \brief Sets up the Jacobian matrix of the coordinate mapping. + //! \param[out] Jac The inverse of the Jacobian matrix + //! \param[in] Xnod Matrix of element nodal coordinates + //! \param[in] dNxdu First order derivatives of basis functions + //! \param[in] gBasis 1-based index of basis representing the geometry + bool Jacobian(Matrix& Jac, const Matrix& Xnod, + const std::vector& dNxdu, unsigned short int gBasis); + + //! \brief Sets up the Hessian matrix of the coordinate mapping. + //! \param[out] Hess The Hessian matrix + //! \param[in] Jac The inverse of the Jacobian matrix + //! \param[in] Xnod Matrix of element nodal coordinates + //! \param[in] d2Nxdu2 Second order derivatives of basis functions + //! \param[in] gBasis 1-based index of basis representing the geometry + bool Hessian(Matrix3D& Hess, const Matrix& Jac, const Matrix& Xnod, + const std::vector& d2Nxdu2, unsigned short int gBasis); protected: //! \brief Writes the finite element object to the given output stream. diff --git a/src/ASM/LR/ASMu2Dmx.C b/src/ASM/LR/ASMu2Dmx.C index 1e89f6001..4119aa032 100644 --- a/src/ASM/LR/ASMu2Dmx.C +++ b/src/ASM/LR/ASMu2Dmx.C @@ -109,7 +109,7 @@ void ASMu2Dmx::clear (bool retainGeometry) { if (!retainGeometry) { // Erase the spline data - for (auto& patch : m_basis) + for (SplinePtr& patch : m_basis) patch.reset(); m_basis.clear(); @@ -258,7 +258,7 @@ bool ASMu2Dmx::generateFEMTopology () nb.clear(); nb.reserve(m_basis.size()); - for (const auto& it : m_basis) { + for (const SplinePtr& it : m_basis) { nb.push_back(it->nBasisFunctions()); #ifdef SP_DEBUG std::cout <<"Basis "<< nb.size() @@ -275,23 +275,21 @@ bool ASMu2Dmx::generateFEMTopology () myMLGE.resize(nel,0); myMLGN.resize(nnod); myMNPC.resize(nel); - for (auto&& it : m_basis) + for (SplinePtr& it : m_basis) it->generateIDs(); size_t iel = 0; for (const LR::Element* el1 : m_basis[geoBasis-1]->getAllElements()) { size_t nfunc = 0; - for (const auto& it : m_basis) { - auto el_it2 = it->elementBegin() + it->getElementContaining(el1->midpoint()); - nfunc += (*el_it2)->nBasisFunctions(); - } + for (const SplinePtr& it : m_basis) + nfunc += it->getElement(it->getElementContaining(el1->midpoint()))->nBasisFunctions(); myMNPC[iel].resize(nfunc); size_t lnod = 0, ofs = 0; - for (const auto& it : m_basis) { - auto el_it2 = it->elementBegin() + it->getElementContaining(el1->midpoint()); - for (LR::Basisfunction* b : (*el_it2)->support()) + for (const SplinePtr& it : m_basis) { + const LR::Element* el2 = it->getElement(it->getElementContaining(el1->midpoint())); + for (LR::Basisfunction* b : el2->support()) myMNPC[iel][lnod++] = b->getId() + ofs; ofs += it->nBasisFunctions(); } @@ -327,6 +325,14 @@ bool ASMu2Dmx::integrate (Integrand& integrand, const double* xg = GaussQuadrature::getCoord(nGauss); const double* wg = GaussQuadrature::getWeight(nGauss); 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; ThreadGroups oneGroup; @@ -335,40 +341,36 @@ bool ASMu2Dmx::integrate (Integrand& integrand, // === Assembly loop over all elements in the patch ========================== + bool ok = true; for (size_t t = 0; t < groups.size() && ok; ++t) - { #pragma omp parallel for schedule(static) for (size_t e = 0; e < groups[t].size(); ++e) { if (!ok) continue; - int iel = groups[t][e] + 1; - auto el1 = threadBasis->elementBegin()+iel-1; - std::vector els; - std::vector elem_sizes; - for (size_t i=0; i < m_basis.size(); ++i) { - els.push_back(m_basis[i]->getElementContaining((*el1)->midpoint())+1); - elem_sizes.push_back(m_basis[i]->getElement(els.back()-1)->nBasisFunctions()); - } - int geoEl = els[geoBasis-1]; + std::vector els; + std::vector elem_sizes; + this->getElementsAt(threadBasis->getElement(groups[t][e])->midpoint(),els,elem_sizes); - MxFiniteElement fe(elem_sizes); - fe.iel = MLGE[geoEl-1]; - std::vector dNxdu(m_basis.size()); + MxFiniteElement fe(elem_sizes); + std::vector dNxdu(m_basis.size()); + std::vector d2Nxdu2(use2ndDer ? m_basis.size() : 0); Matrix Xnod, Jac; - double param[3] = { 0.0, 0.0, 0.0 }; - Vec4 X(param,time.t); - std::vector d2Nxdu2(m_basis.size()); Matrix3D Hess; 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 - double dA = this->getParametricArea(geoEl); - if (dA < 0.0) // topology error (probably logic error) + double dA = 0.25*this->getParametricArea(geoEl); + if (dA < 0.0) { - ok = false; + ok = false; // topology error (probably logic error) continue; } @@ -395,17 +397,17 @@ bool ASMu2Dmx::integrate (Integrand& integrand, this->getGaussPointParameters(gpar[d],d,nGauss,geoEl,xg); // Initialize element quantities - LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,false); - if (!integrand.initElement(MNPC[geoEl-1], fe, elem_sizes, nb, *A)) + LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel); + if (!integrand.initElement(MNPC[geoEl-1],fe,elem_sizes,nb,*A)) { A->destruct(); ok = false; continue; } - // --- 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 = (geoEl-1)*nGauss*nGauss; fe.iGP = firstIp + jp; // Global integration point counter for (int j = 0; j < nGauss; j++) @@ -433,25 +435,14 @@ bool ASMu2Dmx::integrate (Integrand& integrand, SplineUtils::extractBasis(spline,fe.basis(b+1),dNxdu[b]); } - // Compute Jacobian inverse of coordinate mapping and derivatives + // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1]); - 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) - fe.grad(b+1).multiply(dNxdu[b],Jac); + if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives - if (use2ndDer) { - if (!utl::Hessian(Hess,fe.hess(geoBasis),Jac,Xnod, - d2Nxdu2[geoBasis-1],fe.grad(geoBasis),true)) - ok = false; - - for (size_t b = 0; b < m_basis.size() && ok; ++b) - if ((int)b != geoBasis-1) - utl::Hessian(Hess,fe.hess(b+1),Jac,Xnod, - d2Nxdu2[b],fe.grad(b+1),false); - } + if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis)) + ok = false; // Compute G-matrix if (integrand.getIntegrandType() & Integrand::G_MATRIX) @@ -461,7 +452,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand, X.assign(Xnod * fe.basis(geoBasis)); // 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)) ok = false; } @@ -476,7 +467,6 @@ bool ASMu2Dmx::integrate (Integrand& integrand, A->destruct(); } - } return ok; } @@ -486,7 +476,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex, GlobalIntegral& glInt, const TimeDomain& time) { - if (!m_basis[0] || !m_basis[1]) + if (m_basis.empty()) return true; // silently ignore empty patches PROFILE2("ASMu2Dmx::integrate(B)"); @@ -531,17 +521,19 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex, // === Assembly loop over all elements on the patch edge ===================== - std::vector::iterator el1 = m_basis[geoBasis-1]->elementBegin(); - for (int iel = 1; el1 != m_basis[geoBasis-1]->elementEnd(); ++el1, ++iel) + int iel = 0; + for (const LR::Element* el1 : m_basis[geoBasis-1]->getAllElements()) { + ++iel; + // Skip elements that are not on current boundary edge bool skipMe = false; switch (edgeDir) { - case -1: if ((*el1)->umin() != m_basis[geoBasis-1]->startparam(0)) skipMe = true; break; - case 1: if ((*el1)->umax() != m_basis[geoBasis-1]->endparam(0) ) skipMe = true; break; - case -2: if ((*el1)->vmin() != m_basis[geoBasis-1]->startparam(1)) skipMe = true; break; - case 2: if ((*el1)->vmax() != m_basis[geoBasis-1]->endparam(1) ) skipMe = true; break; + case -1: if (el1->umin() != m_basis[geoBasis-1]->startparam(0)) skipMe = true; break; + case 1: if (el1->umax() != m_basis[geoBasis-1]->endparam(0) ) skipMe = true; break; + case -2: if (el1->vmin() != m_basis[geoBasis-1]->startparam(1)) skipMe = true; break; + case 2: if (el1->vmax() != m_basis[geoBasis-1]->endparam(1) ) skipMe = true; break; } if (skipMe) continue; @@ -549,12 +541,9 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex, std::find(myElms.begin(), myElms.end(), iel-1) == myElms.end()) continue; - std::vector els; + std::vector els; std::vector elem_sizes; - for (size_t i=0; i < m_basis.size(); ++i) { - els.push_back(m_basis[i]->getElementContaining((*el1)->midpoint())+1); - elem_sizes.push_back((*(m_basis[i]->elementBegin()+(els.back()-1)))->nBasisFunctions()); - } + this->getElementsAt(el1->midpoint(),els,elem_sizes); int geoEl = els[geoBasis-1]; MxFiniteElement fe(elem_sizes,firstp); @@ -607,9 +596,9 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex, dNxdu[geoBasis-1],t1,t2); 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) - fe.grad(b+1).multiply(dNxdu[b],Jac); + for (size_t b = 1; b <= m_basis.size(); ++b) + if ((int)b != geoBasis) + fe.grad(b).multiply(dNxdu[b-1],Jac); if (edgeDir < 0) normal *= -1.0; @@ -644,8 +633,11 @@ bool ASMu2Dmx::integrate (Integrand& integrand, const TimeDomain& time, const ASM::InterfaceChecker& iChkgen) { - if (!geo) return true; // silently ignore empty patches - if (!(integrand.getIntegrandType() & Integrand::INTERFACE_TERMS)) return true; + if (m_basis.empty()) + return true; // silently ignore empty patches + + if (!(integrand.getIntegrandType() & Integrand::INTERFACE_TERMS)) + return true; // No interface terms PROFILE2("ASMu2Dmx::integrate(J)"); @@ -662,24 +654,22 @@ bool ASMu2Dmx::integrate (Integrand& integrand, Vec4 X(nullptr,time.t); Vec3 normal; - std::vector::iterator el1 = m_basis[0]->elementBegin(); - for (int iel = 1; el1 != m_basis[0]->elementEnd(); ++el1, ++iel) { - short int status = iChk.hasContribution(iel); + int iel = 0; + for (const LR::Element* el1 : m_basis.front()->getAllElements()) + { + short int status = iChk.hasContribution(++iel); if (!status) continue; // no interface contributions for this element if (!myElms.empty() && !glInt.threadSafe() && - std::find(myElms.begin(), myElms.end(), iel-1) == myElms.end()) { - ++iel, ++el1; + std::find(myElms.begin(), myElms.end(), iel-1) == myElms.end()) continue; - } - // first push the hosting elements - std::vector els(1,iel); - std::vector elem_sizes(1,(*el1)->nBasisFunctions()); - for (size_t i=1; i < m_basis.size(); ++i) { - els.push_back(m_basis[i]->getElementContaining((*el1)->midpoint())+1); - elem_sizes.push_back(m_basis[i]->getElement(els.back()-1)->nBasisFunctions()); - } + std::vector els; + std::vector elem_sizes; + this->getElementsAt(el1->midpoint(),els,elem_sizes); + // Replace first entry by hosting element + els.front() = iel; + elem_sizes.front() = el1->nBasisFunctions(); // Set up control point coordinates for current element if (!this->getElementCoordinates(Xnod,els[geoBasis-1])) @@ -697,8 +687,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand, const int t2 = 3-abs(edgeDir); // Tangent direction along the edge // Set up parameters - double u1 = iedge != 2 ? (*el1)->umin() : (*el1)->umax(); - double v1 = iedge < 4 ? (*el1)->vmin() : (*el1)->vmax(); + double u1 = iedge != 2 ? el1->umin() : el1->umax(); + double v1 = iedge < 4 ? el1->vmin() : el1->vmax(); double u2(u1); double v2(v1); const double epsilon = 1e-8; @@ -714,24 +704,18 @@ bool ASMu2Dmx::integrate (Integrand& integrand, std::vector intersections = iChk.getIntersections(iel, iedge); for (size_t i = 0; i < intersections.size() && ok; ++i) { - std::vector parval(2); - parval[0] = u1-epsu; - parval[1] = v1-epsv; - if (iedge == 1 || iedge == 2) v2 = intersections[i]; else u2 = intersections[i]; - int el_neigh = this->getBasis(1)->getElementContaining(parval)+1; - const LR::Element* el2 = m_basis[0]->getElement(el_neigh-1); - - std::vector els2(1,el_neigh); - std::vector elem_sizes2(1,el2->nBasisFunctions()); - for (size_t i=1; i < m_basis.size(); ++i) { - els2.push_back(m_basis[i]->getElementContaining(el2->midpoint())+1); - elem_sizes2.push_back(m_basis[i]->getElement(els2.back()-1)->nBasisFunctions()); - } + int el_neigh = this->getBasis(1)->getElementContaining({u1-epsu,v1-epsv})+1; + const LR::Element* el2 = m_basis.front()->getElement(el_neigh-1); + std::vector els2; + std::vector elem_sizes2; + this->getElementsAt(el2->midpoint(),els2,elem_sizes2); + els2.front() = el_neigh; + elem_sizes2.front() = el2->nBasisFunctions(); LocalIntegral* A_neigh = integrand.getLocalIntegral(elem_sizes2, el_neigh); ok = integrand.initElement(MNPC[els2[geoBasis-1]-1], @@ -793,16 +777,18 @@ bool ASMu2Dmx::integrate (Integrand& integrand, // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates - fe.detJxW = utl::Jacobian(Jac2, normal, - fe.grad(geoBasis+m_basis.size()), - Xnod2,dNxdu[geoBasis-1+m_basis.size()],t1,t2); - fe.detJxW = utl::Jacobian(Jac, normal, - fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1],t1,t2); + fe.detJxW = utl::Jacobian(Jac2,normal, + fe.grad(geoBasis+m_basis.size()),Xnod2, + dNxdu[geoBasis-1+m_basis.size()],t1,t2); + fe.detJxW = utl::Jacobian(Jac,normal, + fe.grad(geoBasis),Xnod, + dNxdu[geoBasis-1],t1,t2); 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) { - fe.grad(b+1).multiply(dNxdu[b],Jac); - fe.grad(b+1+m_basis.size()).multiply(dNxdu[b+m_basis.size()],Jac); + + for (size_t b = 1; b <= m_basis.size(); ++b) + if ((int)b != geoBasis) { + 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) @@ -865,31 +851,29 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const Vector& locSol, sField.resize(std::accumulate(nc.begin(), nc.end(), 0), nPoints); for (size_t i = 0; i < nPoints; i++) { - size_t ofs=0; - Vector Ztmp; + RealArray Ztmp; + const double* locPtr = locSol.data(); for (size_t j=0; j < m_basis.size(); ++j) { if (nc[j] == 0) continue; + // Fetch element containing evaluation point. // 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]); + const LR::Element* el = m_basis[j]->getElement(iel); // Evaluate basis function values/derivatives at current parametric point // and multiply with control point values to get the point-wise solution this->computeBasis(gpar[0][i],gpar[1][i],splinex[j],iel,m_basis[j].get()); - std::vector::iterator el_it = m_basis[j]->elementBegin()+iel; Matrix val1(nc[j], splinex[j].basisValues.size()); - size_t col=1; - for (auto* b : (*el_it)->support()) { - for (size_t n = 1; n <= nc[j]; ++n) - val1(n, col) = locSol(b->getId()*nc[j]+n+ofs); - ++col; - } + size_t col = 0; + for (LR::Basisfunction* b : el->support()) + val1.fillColumn(++col, locPtr+b->getId()*nc[j]); Vector Ytmp; val1.multiply(splinex[j].basisValues,Ytmp); Ztmp.insert(Ztmp.end(),Ytmp.begin(),Ytmp.end()); - ofs += nb[j]*nc[j]; + locPtr += nb[j]*nc[j]; } sField.fillColumn(i+1, Ztmp); @@ -921,18 +905,15 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, { // Fetch element containing evaluation point // sadly, points are not always ordered in the same way as the elements - std::vector els; + std::vector els; std::vector elem_sizes; - for (size_t b = 0; b < m_basis.size(); ++b) { - els.push_back(m_basis[b]->getElementContaining(gpar[0][i],gpar[1][i])+1); - elem_sizes.push_back((*(m_basis[b]->elementBegin()+els.back()-1))->nBasisFunctions()); - } + this->getElementsAt({gpar[0][i],gpar[1][i]},els,elem_sizes); // Evaluate the basis functions at current parametric point - MxFiniteElement fe(elem_sizes); - std::vector dNxdu(m_basis.size()); + MxFiniteElement fe(elem_sizes,firstIp+i); + std::vector dNxdu(m_basis.size()); + std::vector d2Nxdu2(use2ndDer ? m_basis.size() : 0); Matrix Jac, Xnod; - std::vector d2Nxdu2(m_basis.size()); Matrix3D Hess; if (use2ndDer) for (size_t b = 0; b < m_basis.size(); ++b) { @@ -950,33 +931,22 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, // Set up control point (nodal) coordinates for current element if (!this->getElementCoordinates(Xnod,els[geoBasis-1])) return false; - // Compute the Jacobian inverse - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1]); - 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) - fe.grad(b+1).multiply(dNxdu[b],Jac); + // Compute Jacobian inverse of the coordinate mapping and + // basis function derivatives w.r.t. Cartesian coordinates + if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives - if (use2ndDer) { - if (!utl::Hessian(Hess,fe.hess(geoBasis),Jac,Xnod, - d2Nxdu2[geoBasis-1],fe.grad(geoBasis),true)) - return false; - - for (size_t b = 0; b < m_basis.size(); ++b) - if (b != (size_t)geoBasis) - utl::Hessian(Hess,fe.hess(b+1),Jac,Xnod, - d2Nxdu2[b],fe.grad(b+1),false); - } - - // Now evaluate the solution field - Vector solPt; + if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis)) + return false; // Cartesian coordinates of current integration point fe.u = gpar[0][i]; fe.v = gpar[1][i]; utl::Point X4(Xnod*fe.basis(geoBasis),{fe.u,fe.v}); + // Now evaluate the solution field + Vector solPt; if (!integrand.evalSol(solPt,fe,X4, MNPC[els[geoBasis-1]-1],elem_sizes,nb)) return false; @@ -998,13 +968,12 @@ bool ASMu2Dmx::refine (const LR::RefineData& prm, Vectors& sol) if (prm.errors.empty() && prm.elements.empty()) return true; - for (Vector& solvec : sol) { + for (Vector& solvec : sol) for (size_t j = 0; j < m_basis.size(); j++) { Vector bVec; this->extractNodeVec(solvec, bVec, 0, j+1); LR::extendControlPoints(m_basis[j].get(), bVec, nfx[j]); } - } if (doRefine(prm, refBasis.get())) { for (size_t j = 0; j < m_basis.size(); ++j) @@ -1046,10 +1015,10 @@ bool ASMu2Dmx::refine (const LR::RefineData& prm, Vectors& sol) #ifdef SP_DEBUG std::cout <<"Refined mesh: "; - for (const auto& it : m_basis) + for (const SplinePtr& it : m_basis) std::cout << it->nElements() <<" "; std::cout <<"elements "; - for (const auto& it : m_basis) + for (const SplinePtr& it : m_basis) std::cout << it->nBasisFunctions() <<" "; std::cout <<"nodes."<< std::endl; std::cout << "Projection basis: " @@ -1114,7 +1083,7 @@ void ASMu2Dmx::generateThreadGroups (const Integrand& integrand, bool silence, LR::generateThreadGroups(altProjThreadGroups,altProjBasis.get()); std::vector bases; - for (const std::shared_ptr& basis : m_basis) + for (const SplinePtr& basis : m_basis) bases.push_back(basis.get()); if (silence || threadGroups[0].size() < 2) return; @@ -1160,8 +1129,8 @@ void ASMu2Dmx::getBoundaryNodes (int lIndex, IntVec& nodes, int basis, } -void ASMu2Dmx::remapErrors(RealArray& errors, - const RealArray& origErr, bool elemErrors) const +void ASMu2Dmx::remapErrors (RealArray& errors, + const RealArray& origErr, bool elemErrors) const { const LR::LRSplineSurface* geo = this->getBasis(ASMmxBase::geoBasis); for (const LR::Element* elm : geo->getAllElements()) { @@ -1219,7 +1188,7 @@ void ASMu2Dmx::storeMesh (const std::string& fName, int fType) const }; std::string btag("basis1"); - for (const auto& patch : m_basis) + for (const SplinePtr& patch : m_basis) { writeBasis(patch, btag); ++btag.back(); @@ -1256,3 +1225,20 @@ void ASMu2Dmx::swapProjectionBasis () this->generateBezierExtraction(); } } + + +void ASMu2Dmx::getElementsAt (const RealArray& param, + std::vector& elms, + std::vector& sizes) const +{ + elms.clear(); + sizes.clear(); + elms.reserve(m_basis.size()); + sizes.reserve(m_basis.size()); + for (const SplinePtr& basis : m_basis) + { + int iel = basis->getElementContaining(param); + elms.push_back(1+iel); + sizes.push_back(basis->getElement(iel)->nBasisFunctions()); + } +} diff --git a/src/ASM/LR/ASMu2Dmx.h b/src/ASM/LR/ASMu2Dmx.h index dd2766426..ad15939aa 100644 --- a/src/ASM/LR/ASMu2Dmx.h +++ b/src/ASM/LR/ASMu2Dmx.h @@ -219,6 +219,15 @@ class ASMu2Dmx : public ASMu2D, private ASMmxBase //! \brief Swaps between the main and alternative projection basis. virtual void swapProjectionBasis(); +private: + //! \brief Finds the elements and associted sizes at given parametric point. + //! \param[in] param Parametric point to find elements at + //! \param[out] elms List of elements on each basis containign given point + //! \param[out] sizes List of element sizes (numer of element nodes) + void getElementsAt(const RealArray& param, + std::vector& elms, + std::vector& sizes) const; + protected: using ASMu2D::generateThreadGroups; //! \brief Generates element groups for multi-threading of interior integrals. @@ -229,10 +238,13 @@ class ASMu2Dmx : public ASMu2D, private ASMmxBase bool ignoreGlobalLM); private: - std::vector> m_basis; //!< All bases - LR::LRSplineSurface* threadBasis; //!< Basis for thread groups - std::shared_ptr refBasis; //!< Basis to refine based on - std::shared_ptr altProjBasis; //!< Alternative projection basis + typedef std::shared_ptr SplinePtr; //!< Pointer to spline + + std::vector m_basis; //!< All bases + LR::LRSplineSurface* threadBasis; //!< Basis for thread groups + SplinePtr refBasis; //!< Basis to refine based on + SplinePtr altProjBasis; //!< Alternative projection basis + ThreadGroups altProjThreadGroups; //!< Element groups for multi-threaded assembly - alternative projection basis }; diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C index 8affd4b6f..176583b7d 100644 --- a/src/ASM/LR/ASMu3D.C +++ b/src/ASM/LR/ASMu3D.C @@ -887,41 +887,6 @@ void ASMu3D::evaluateBasis (int iel, FiniteElement& fe, } } -void ASMu3D::evaluateBasis (int iel, FiniteElement& fe, int derivs, - int basis) const -{ - PROFILE3("ASMu3D::evalBasis"); - - std::vector result; - this->getBasis(basis)->computeBasis(fe.u, fe.v, fe.w, result, derivs, iel); - size_t n = 0, nBasis = result.size(); - Vector& N = fe.basis(basis); - Matrix& dNdu = fe.grad(basis); - Matrix3D& d2Ndu2 = fe.hess(basis); - - N.resize(nBasis); - if (derivs > 0) - dNdu.resize(nBasis,3); - if (derivs > 1) - d2Ndu2.resize(nBasis,3,3); - for (const RealArray& phi : result) { - N(++n) = phi[0]; - if (derivs > 0) { - dNdu(n,1) = phi[1]; - dNdu(n,2) = phi[2]; - dNdu(n,3) = phi[3]; - } - if (derivs > 1) { - d2Ndu2(n,1,1) = phi[4]; - d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = phi[5]; - d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = phi[6]; - d2Ndu2(n,2,2) = phi[7]; - d2Ndu2(n,2,3) = d2Ndu2(n,3,2) = phi[8]; - d2Ndu2(n,3,3) = phi[9]; - } - } -} - bool ASMu3D::integrate (Integrand& integrand, GlobalIntegral& glInt, @@ -1040,15 +1005,10 @@ bool ASMu3D::integrate (Integrand& integrand, // Get element volume in the parameter space const LR::Element* el = lrspline->getElement(iel-1); - double du = el->umax() - el->umin(); - double dv = el->vmax() - el->vmin(); - double dw = el->wmax() - el->wmin(); - double dV = el->volume(); - if (dV < 0.0) - { - ok = false; - continue; - } + double du = 0.5*(el->umax() - el->umin()); + double dv = 0.5*(el->vmax() - el->vmin()); + double dw = 0.5*(el->wmax() - el->wmin()); + double dV = du*dv*dw; // Set up control point (nodal) coordinates for current element if (!this->getElementCoordinates(Xnod,iel)) @@ -1082,12 +1042,9 @@ bool ASMu3D::integrate (Integrand& integrand, } if (integrand.getIntegrandType() & Integrand::G_MATRIX) - { // Element size in parametric space - dXidu[0] = el->getParmin(0); - dXidu[1] = el->getParmin(1); - dXidu[2] = el->getParmin(2); - } + for (int i = 0; i < 3; i++) + dXidu[i] = el->getParmax(i) - el->getParmin(i); size_t nen = el->support().size(); if (integrand.getIntegrandType() & Integrand::AVERAGE) @@ -1103,9 +1060,9 @@ bool ASMu3D::integrate (Integrand& integrand, 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); + B.fillColumn(2, BdNdu.getColumn(ig)/du); + B.fillColumn(3, BdNdv.getColumn(ig)/dv); + B.fillColumn(4, BdNdw.getColumn(ig)/dw); this->evaluateBasis(fe, dNdu, C, B); // Compute Jacobian determinant of coordinate mapping @@ -1152,9 +1109,9 @@ bool ASMu3D::integrate (Integrand& integrand, // Extract bezier basis functions B.fillColumn(1, rBN.getColumn(ig)); - B.fillColumn(2, rBdNdu.getColumn(ig)*2.0/du); - B.fillColumn(3, rBdNdv.getColumn(ig)*2.0/dv); - B.fillColumn(4, rBdNdw.getColumn(ig)*2.0/dw); + B.fillColumn(2, rBdNdu.getColumn(ig)/du); + B.fillColumn(3, rBdNdv.getColumn(ig)/dv); + B.fillColumn(4, rBdNdw.getColumn(ig)/dw); this->evaluateBasis(fe, dNdu, C, B); // Compute Jacobian inverse and derivatives @@ -1165,7 +1122,7 @@ bool ASMu3D::integrate (Integrand& integrand, X.assign(Xnod * fe.N); // Compute the reduced integration terms of the integrand - fe.detJxW *= 0.125*dV*wr[i]*wr[j]*wr[k]; + fe.detJxW *= dV*wr[i]*wr[j]*wr[k]; if (!integrand.reducedInt(*A,fe,X)) ok = false; } @@ -1200,9 +1157,9 @@ bool ASMu3D::integrate (Integrand& integrand, { // Extract bezier basis functions 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); + B.fillColumn(2, BdNdu.getColumn(ig)/du); + B.fillColumn(3, BdNdv.getColumn(ig)/dv); + B.fillColumn(4, BdNdw.getColumn(ig)/dw); this->evaluateBasis(fe, dNdu, C, B); #ifdef SP_DEBUG @@ -1250,7 +1207,7 @@ bool ASMu3D::integrate (Integrand& integrand, X.assign(Xnod * fe.N); // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.125*dV*wg[i]*wg[j]*wg[k]; + fe.detJxW *= dV*wg[i]*wg[j]*wg[k]; PROFILE3("Integrand::evalInt"); if (!integrand.evalInt(*A,fe,time,X)) ok = false; @@ -1376,12 +1333,9 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex, fe.h = this->getElementCorners(iEl+1,fe.XC); if (integrand.getIntegrandType() & Integrand::G_MATRIX) - { // Element size in parametric space - dXidu[0] = el->getParmax(0) - el->getParmin(0); - dXidu[1] = el->getParmax(1) - el->getParmin(1); - dXidu[2] = el->getParmax(2) - el->getParmin(2); - } + for (int i = 0; i < 3; i++) + dXidu[i] = el->getParmax(i) - el->getParmin(i); // Initialize element quantities size_t nen = el->support().size(); @@ -1495,7 +1449,9 @@ bool ASMu3D::diracPoint (Integrand& integrand, GlobalIntegral& glInt, fe.u = param[0]; fe.v = param[1]; fe.w = param[2]; - this->evaluateBasis(iel,fe); + Go::BasisPts pt; + this->getBasis()->computeBasis(fe.u, fe.v, fe.w, pt, iel); + fe.N = pt.basisValues; LocalIntegral* A = integrand.getLocalIntegral(MNPC[iel].size(),fe.iel,true); bool ok = integrand.evalPoint(*A,fe,pval) && glInt.assemble(A,fe.iel); diff --git a/src/ASM/LR/ASMu3D.h b/src/ASM/LR/ASMu3D.h index a5084562c..fe8daa9a8 100644 --- a/src/ASM/LR/ASMu3D.h +++ b/src/ASM/LR/ASMu3D.h @@ -574,10 +574,6 @@ class ASMu3D : public ASMLRSpline, public ASM3D void evaluateBasis(int iel, int basis, double u, double v, double w, Vector& N, Matrix& dNdu) const; - //! \brief Evaluate all basis functions and \a derivs number of derivatives on one element - void evaluateBasis(int iel, FiniteElement& fe, - int derivs = 0, int basis = 1) const; - //! \brief Evaluate all basis functions and first derivatives on one element void evaluateBasis(FiniteElement& fe, Matrix& dNdu, const Matrix& C, const Matrix& B, int basis = 1) const; diff --git a/src/ASM/LR/ASMu3Dmx.C b/src/ASM/LR/ASMu3Dmx.C index 46f4f8b12..841efc0bf 100644 --- a/src/ASM/LR/ASMu3Dmx.C +++ b/src/ASM/LR/ASMu3Dmx.C @@ -40,6 +40,7 @@ ASMu3Dmx::ASMu3Dmx (const CharVec& n_f) : ASMu3D(std::accumulate(n_f.begin(), n_f.end(), 0)), ASMmxBase(n_f), bezierExtractmx(myBezierExtractmx) { + threadBasis = nullptr; myGeoBasis = ASMmxBase::geoBasis; } @@ -49,6 +50,7 @@ ASMu3Dmx::ASMu3Dmx (const ASMu3Dmx& patch, const CharVec& n_f) m_basis(patch.m_basis), bezierExtractmx(patch.myBezierExtractmx) { + threadBasis = patch.threadBasis; nfx = patch.nfx; nb = patch.nb; myGeoBasis = ASMmxBase::geoBasis; @@ -110,7 +112,7 @@ void ASMu3Dmx::clear (bool retainGeometry) { if (!retainGeometry) { // Erase the spline data - for (auto& patch : m_basis) + for (SplinePtr& patch : m_basis) patch.reset(); m_basis.clear(); @@ -199,11 +201,14 @@ bool ASMu3Dmx::getSolution (Matrix& sField, const Vector& locSol, bool ASMu3Dmx::generateFEMTopology () { + if (!myMLGN.empty()) + return true; + if (m_basis.empty()) { - auto vec = ASMmxBase::establishBases(tensorspline, ASMmxBase::Type); - m_basis.resize(vec.size()); - for (size_t i=0;inBasisFunctions(); - - if (shareFE == 'F') return true; - + nb.clear(); + nb.reserve(m_basis.size()); + for (const SplinePtr& it : m_basis) { + nb.push_back(it->nBasisFunctions()); #ifdef SP_DEBUG - size_t nbasis=0; - for (auto& it : m_basis) { - std::cout << "Basis " << ++nbasis << ":\n"; - std::cout <<"numCoefs: "<< it->nBasisFunctions(); - std::cout <<"\norder: "<< it->order(0) <<" "<< + std::cout <<"Basis "<< nb.size() + <<":\nnumCoefs: "<< nb.back() + <<"\norder: "<< it->order(0) <<" "<< it->order(1) <<" "<< it->order(2) << std::endl; - } #endif + } - nel = m_basis[geoBasis-1]->nElements(); + if (shareFE == 'F') return true; + nel = m_basis[geoBasis-1]->nElements(); nnod = std::accumulate(nb.begin(), nb.end(), 0); myMLGE.resize(nel,0); myMLGN.resize(nnod); myMNPC.resize(nel); - for (auto& it : m_basis) + for (SplinePtr& it : m_basis) it->generateIDs(); - std::vector::iterator el_it1 = m_basis[geoBasis-1]->elementBegin(); - for (size_t iel=0; ielgetAllElements()) { size_t nfunc = 0; - for (size_t i=0; ielementBegin() + - m_basis[i]->getElementContaining((*el_it1)->midpoint()); - nfunc += (*el_it2)->nBasisFunctions(); - } - myMLGE[iel] = ++gEl; // global element number over all patches + for (const SplinePtr& it : m_basis) + nfunc += it->getElement(it->getElementContaining(el1->midpoint()))->nBasisFunctions(); myMNPC[iel].resize(nfunc); - int lnod = 0; - size_t ofs=0; - for (size_t i=0; ielementBegin() + - m_basis[i]->getElementContaining((*el_it1)->midpoint()); - for (LR::Basisfunction *b : (*el_it2)->support()) - myMNPC[iel][lnod++] = b->getId()+ofs; - ofs += nb[i]; + size_t lnod = 0, ofs = 0; + for (const SplinePtr& it : m_basis) { + const LR::Element* el2 = it->getElement(it->getElementContaining(el1->midpoint())); + for (LR::Basisfunction* b : el2->support()) + myMNPC[iel][lnod++] = b->getId() + ofs; + ofs += it->nBasisFunctions(); } + + myMLGE[iel++] = ++gEl; // global element number over all patches } + size_t b = 0; myBezierExtractmx.resize(m_basis.size()); - for (size_t b = 1; b <= m_basis.size(); ++b) { - myBezierExtractmx[b-1].resize(this->getBasis(b)->nElements()); - std::vector::const_iterator eit = this->getBasis(b)->elementBegin(); - for (int iel = 0; iel < this->getBasis(b)->nElements(); iel++, ++eit) + for (const SplinePtr& it : m_basis) { + myBezierExtractmx[b].resize(it->nElements()); + for (int iel = 0; iel < it->nElements(); iel++) { PROFILE("Bezier extraction"); - int p1 = this->getBasis(b)->order(0); - int p2 = this->getBasis(b)->order(1); - int p3 = this->getBasis(b)->order(2); - // Get bezier extraction matrix RealArray extrMat; - this->getBasis(b)->getBezierExtraction(iel,extrMat); - myBezierExtractmx[b-1][iel].resize((*eit)->nBasisFunctions(),p1*p2*p3); - myBezierExtractmx[b-1][iel].fill(extrMat.data(),extrMat.size()); + it->getBezierExtraction(iel,extrMat); + myBezierExtractmx[b][iel].resize(it->getElement(iel)->nBasisFunctions(), + it->order(0)*it->order(1)*it->order(2)); + myBezierExtractmx[b][iel].fill(extrMat.data(),extrMat.size()); } + ++b; } for (size_t inod = 0; inod < nnod; ++inod) @@ -323,6 +319,9 @@ bool ASMu3Dmx::integrate (Integrand& integrand, GlobalIntegral& glInt, const TimeDomain& time) { + if (m_basis.empty()) + return true; // silently ignore empty patches + PROFILE2("ASMu3Dmx::integrate(I)"); // Get Gaussian quadrature points and weights @@ -330,6 +329,15 @@ bool ASMu3Dmx::integrate (Integrand& integrand, const double* wg = GaussQuadrature::getWeight(nGauss); 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) std::vector BN(m_basis.size()); std::vector BdNdu(m_basis.size()); @@ -340,98 +348,70 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int p1 = lrspline->order(0); int p2 = lrspline->order(1); int p3 = lrspline->order(2); + double u[2*p1], v[2*p2], w[2*p3]; Go::BsplineBasis basis1 = getBezierBasis(p1); Go::BsplineBasis basis2 = getBezierBasis(p2); 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); BdNdv[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 - for(int zeta=0; zeta 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; if (glInt.threadSafe()) oneGroup.oneGroup(nel); - const IntMat& group = glInt.threadSafe() ? oneGroup[0] : threadGroups[0]; + const IntMat& groups = glInt.threadSafe() ? oneGroup[0] : threadGroups[0]; // === Assembly loop over all elements in the patch ========================== bool ok = true; - for (size_t t = 0; t < group.size() && ok; t++) + for (size_t t = 0; t < groups.size() && ok; t++) #pragma omp parallel for schedule(static) - for (size_t e = 0; e < group[t].size(); e++) + for (size_t e = 0; e < groups[t].size(); e++) { if (!ok) continue; - int iel = group[t][e] + 1; - const LR::Element* el = threadBasis->getElement(iel-1); - std::vector els; - std::vector elem_sizes; - for (size_t i = 0; i < m_basis.size(); ++i) { - els.push_back(m_basis[i]->getElementContaining(el->midpoint())+1); - elem_sizes.push_back((*(m_basis[i]->elementBegin()+els.back()-1))->nBasisFunctions()); - } - int iEl = el->getId(); - MxFiniteElement fe(elem_sizes); - fe.iel = MLGE[iEl]; - std::vector dNxdu(m_basis.size()); + const LR::Element* el = threadBasis->getElement(groups[t][e]); + + std::vector els; + std::vector elem_sizes; + this->getElementsAt(el->midpoint(),els,elem_sizes); + MxFiniteElement fe(elem_sizes); + std::vector dNxdu(m_basis.size()); + std::vector d2Nxdu2(use2ndDer ? m_basis.size() : 0); Matrix Xnod, Jac; - std::vector d2Nxdu2(m_basis.size()); Matrix3D Hess; double dXidu[3]; double param[3] = { 0.0, 0.0, 0.0 }; Vec4 X(param,time.t); + + int iEl = el->getId(); + fe.iel = MLGE[iEl]; + // Get element volume in the parameter space - double du = el->umax() - el->umin(); - double dv = el->vmax() - el->vmin(); - double dw = el->wmax() - el->wmin(); - double vol = el->volume(); - if (vol < 0.0) - { - ok = false; // topology error (probably logic error) - continue; - } + double du = 0.5*(el->umax() - el->umin()); + double dv = 0.5*(el->vmax() - el->vmin()); + double dw = 0.5*(el->wmax() - el->wmin()); + double dV = du*dv*dw; // Set up control point (nodal) coordinates for current element if (!this->getElementCoordinates(Xnod,iEl+1)) @@ -441,25 +421,18 @@ bool ASMu3Dmx::integrate (Integrand& integrand, } // Compute parameter values of the Gauss points over the whole element - std::array gpar, redpar; + std::array gpar; for (int d = 0; d < 3; d++) - { 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) fe.h = this->getElementCorners(iEl+1, fe.XC); if (integrand.getIntegrandType() & Integrand::G_MATRIX) - { // Element size in parametric space - dXidu[0] = el->getParmin(0); - dXidu[1] = el->getParmin(1); - dXidu[2] = el->getParmin(2); - } + for (int i = 0; i < 3; i++) + dXidu[i] = el->getParmax(i) - el->getParmin(i); + else if (integrand.getIntegrandType() & Integrand::AVERAGE) { // --- Compute average value of basis functions over the element ----- @@ -478,7 +451,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, // and multiply by weight of current integration point double detJac = utl::Jacobian(Jac,fe.grad(geoBasis), Xnod,dNxdu[geoBasis-1],false); - double weight = 0.125*vol*wg[i]*wg[j]*wg[k]; + double weight = dV*wg[i]*wg[j]*wg[k]; // Numerical quadrature fe.Navg.add(fe.N,detJac*weight); @@ -509,15 +482,10 @@ bool ASMu3Dmx::integrate (Integrand& integrand, 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 ---------- - fe.iGP = iEl*nGauss*nGauss*nGauss; // Global integration point counter + int jp = iEl*nGauss*nGauss*nGauss; + fe.iGP = firstIp + jp; // Global integration point counter std::vector B(m_basis.size()); size_t ig = 1; @@ -539,40 +507,25 @@ bool ASMu3Dmx::integrate (Integrand& integrand, for (size_t b = 0; b < m_basis.size(); ++b) { Matrix B(BN[b].rows(), 4); B.fillColumn(1, BN[b].getColumn(ig)); - B.fillColumn(2, BdNdu[b].getColumn(ig)*2.0/du); - B.fillColumn(3, BdNdv[b].getColumn(ig)*2.0/dv); - B.fillColumn(4, BdNdw[b].getColumn(ig)*2.0/dw); + B.fillColumn(2, BdNdu[b].getColumn(ig)/du); + B.fillColumn(3, BdNdv[b].getColumn(ig)/dv); + B.fillColumn(4, BdNdw[b].getColumn(ig)/dw); // 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); else this->evaluateBasis(fe, dNxdu[b], bezierExtractmx[b][els[b]-1], B, b+1); } - // Compute Jacobian inverse of coordinate mapping and derivatives - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1]); - 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) - fe.grad(b+1).multiply(dNxdu[b],Jac); + // Compute Jacobian inverse of the coordinate mapping and + // basis function derivatives w.r.t. Cartesian coordinates + if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives - if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES) { - if (!utl::Hessian(Hess,fe.hess(geoBasis),Jac,Xnod, - d2Nxdu2[geoBasis-1],dNxdu[geoBasis-1])) { - ok = false; - continue; - } - - for (size_t b = 0; b < m_basis.size() && ok; ++b) - if ((int)b != geoBasis) - if (!utl::Hessian(Hess,fe.hess(b+1),Jac,Xnod, - d2Nxdu2[b],fe.grad(b+1),false)) { - ok = false; - break; - } - } + if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis)) + ok = false; // Compute G-matrix if (integrand.getIntegrandType() & Integrand::G_MATRIX) @@ -582,14 +535,13 @@ bool ASMu3Dmx::integrate (Integrand& integrand, X.assign(Xnod * fe.basis(geoBasis)); // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.125*vol*wg[i]*wg[j]*wg[k]; + fe.detJxW *= dV*wg[i]*wg[j]*wg[k]; if (!integrand.evalIntMx(*A,fe,time,X)) ok = false; - - } // end gauss integrand + } // Finalize the element quantities - if (ok && !integrand.finalizeElement(*A,fe,time)) + if (ok && !integrand.finalizeElement(*A,fe,time,firstIp+jp)) ok = false; // Assembly of global system integral @@ -607,6 +559,9 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex, GlobalIntegral& glInt, const TimeDomain& time) { + if (m_basis.empty()) + return true; // silently ignore empty patches + PROFILE2("ASMu3Dmx::integrate(B)"); // Get Gaussian quadrature points and weights @@ -624,36 +579,37 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex, std::map::const_iterator iit = firstBp.find(lIndex); size_t firstp = iit == firstBp.end() ? 0 : iit->second; - LR::parameterEdge edge; - switch(lIndex) + // Lambda function mapping face index to LR enum value + auto&& getFaceEnum = [](int faceIndex) -> LR::parameterEdge { - case 1: edge = LR::WEST; break; - case 2: edge = LR::EAST; break; - case 3: edge = LR::SOUTH; break; - case 4: edge = LR::NORTH; break; - case 5: edge = LR::BOTTOM; break; - case 6: edge = LR::TOP; break; - default:edge = LR::NONE; - } + switch (faceIndex) { + case 1: return LR::WEST; + case 2: return LR::EAST; + case 3: return LR::SOUTH; + case 4: return LR::NORTH; + case 5: return LR::BOTTOM; + case 6: return LR::TOP; + } + return LR::NONE; + }; - // fetch all elements along the chosen edge + // Fetch all elements on the chosen face std::vector edgeElms; - this->getBasis(geoBasis)->getEdgeElements(edgeElms, (LR::parameterEdge) edge); + this->getBasis(geoBasis)->getEdgeElements(edgeElms,getFaceEnum(lIndex)); - // iterate over all edge elements - bool ok = true; - for(LR::Element *el : edgeElms) { + + // === Assembly loop over all elements on the patch face ===================== + + for (LR::Element* el : edgeElms) + { int iEl = el->getId(); if (!myElms.empty() && !glInt.threadSafe() && std::find(myElms.begin(), myElms.end(), iEl) == myElms.end()) continue; - std::vector els; + std::vector els; std::vector elem_sizes; - for (size_t i=0; i < m_basis.size(); ++i) { - els.push_back(m_basis[i]->getElementContaining(el->midpoint())+1); - elem_sizes.push_back((*(m_basis[i]->elementBegin()+els.back()-1))->nBasisFunctions()); - } + this->getElementsAt(el->midpoint(),els,elem_sizes); MxFiniteElement fe(elem_sizes); fe.iel = MLGE[iEl]; @@ -687,49 +643,37 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex, double dXidu[3]; // Get element face area in the parameter space - double dA = this->getParametricArea(iEl+1,abs(faceDir)); - if (dA < 0.0) // topology error (probably logic error) - { - ok = false; - break; - } + double dA = 0.25*this->getParametricArea(iEl+1,abs(faceDir)); + if (dA < 0.0) return false; // topology error (probably logic error) // Set up control point coordinates for current element if (!this->getElementCoordinates(Xnod,iEl+1)) - { - ok = false; - break; - } + return false; if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS) fe.h = this->getElementCorners(iEl+1,fe.XC); if (integrand.getIntegrandType() & Integrand::G_MATRIX) - { // Element size in parametric space - dXidu[0] = el->getParmax(0) - el->getParmin(0); - dXidu[1] = el->getParmax(1) - el->getParmin(1); - dXidu[2] = el->getParmax(2) - el->getParmin(2); - } + for (int i = 0; i < 3; i++) + dXidu[i] = el->getParmax(i) - el->getParmin(i); // Initialize element quantities LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,true); - if (!integrand.initElementBou(MNPC[iEl],elem_sizes,nb,*A)) - { - A->destruct(); - ok = false; - break; - } + bool ok = integrand.initElementBou(MNPC[iEl],elem_sizes,nb,*A); + - // --- Integration loop over all Gauss points in each direction -------- + // --- Integration loop over all Gauss points in each direction ------------ fe.iGP = firstp; // Global integration point counter - int k1,k2,k3; + firstp += nGP*nGP; + for (int j = 0; j < nGP; j++) - for (int i = 0; i < nGP; i++, fe.iGP++) + for (int i = 0; i < nGP && ok; i++, fe.iGP++) { // Local element coordinates and parameter values // of current integration point + int k1, k2, k3; switch (abs(faceDir)) { case 1: k2 = i; k3 = j; k1 = 0; break; @@ -761,9 +705,10 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex, fe.detJxW = utl::Jacobian(Jac, normal, fe.grad(geoBasis), Xnod, dNxdu[geoBasis-1], t1, t2); 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) - fe.grad(b+1).multiply(dNxdu[b],Jac); + + for (size_t b = 1; b <= m_basis.size(); ++b) + if ((int)b != geoBasis) + fe.grad(b).multiply(dNxdu[b-1],Jac); if (faceDir < 0) normal *= -1.0; @@ -775,10 +720,9 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex, X.assign(Xnod * fe.basis(geoBasis)); // Evaluate the integrand and accumulate element contributions - fe.detJxW *= 0.25*dA*wg[i]*wg[j]; - if (!integrand.evalBouMx(*A,fe,time,X,normal)) - ok = false; - } + fe.detJxW *= dA*wg[i]*wg[j]; + ok = integrand.evalBouMx(*A,fe,time,X,normal); + } // Finalize the element quantities if (ok && !integrand.finalizeElementBou(*A,fe,time)) @@ -790,10 +734,10 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex, A->destruct(); - firstp += nGP*nGP*nGP; + if (!ok) return false; } - return ok; + return true; } @@ -821,31 +765,29 @@ bool ASMu3Dmx::evalSolution (Matrix& sField, const Vector& locSol, sField.resize(std::accumulate(nc.begin(), nc.end(), 0), nPoints); for (size_t i = 0; i < nPoints; i++) { - size_t ofs=0; - Vector Ztmp; + RealArray Ztmp; + const double* locPtr = locSol.data(); for (size_t j=0; j < m_basis.size(); ++j) { if (nc[j] == 0) continue; + // Fetch element containing evaluation point. // 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]); + const LR::Element* el = m_basis[j]->getElement(iel); // Evaluate basis function values/derivatives at current parametric point // and multiply with control point values to get the point-wise solution m_basis[j]->computeBasis(gpar[0][i],gpar[1][i],gpar[2][i],splinex[j],iel); - std::vector::iterator el_it = m_basis[j]->elementBegin()+iel; Matrix val1(nc[j], splinex[j].basisValues.size()); - size_t col=1; - for (auto* b : (*el_it)->support()) { - for (size_t n = 1; n <= nc[j]; ++n) - val1(n, col) = locSol(b->getId()*nc[j]+n+ofs); - ++col; - } + size_t col = 0; + for (LR::Basisfunction* b : el->support()) + val1.fillColumn(++col, locPtr+b->getId()*nc[j]); Vector Ytmp; val1.multiply(splinex[j].basisValues,Ytmp); Ztmp.insert(Ztmp.end(),Ytmp.begin(),Ytmp.end()); - ofs += nb[j]*nc[j]; + locPtr += nb[j]*nc[j]; } sField.fillColumn(i+1, Ztmp); @@ -877,18 +819,15 @@ bool ASMu3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, { // Fetch element containing evaluation point // sadly, points are not always ordered in the same way as the elements - std::vector els; + std::vector els; std::vector elem_sizes; - for (size_t b = 0; b < m_basis.size(); ++b) { - els.push_back(m_basis[b]->getElementContaining(gpar[0][i],gpar[1][i],gpar[2][i])+1); - elem_sizes.push_back(m_basis[b]->getElement(els.back()-1)->nBasisFunctions()); - } + this->getElementsAt({gpar[0][i],gpar[1][i],gpar[2][i]},els,elem_sizes); // Evaluate the basis functions at current parametric point - MxFiniteElement fe(elem_sizes); - std::vector dNxdu(m_basis.size()); + MxFiniteElement fe(elem_sizes,firstIp+i); + std::vector dNxdu(m_basis.size()); + std::vector d2Nxdu2(use2ndDer ? m_basis.size() : 0); Matrix Jac, Xnod; - std::vector d2Nxdu2(m_basis.size()); Matrix3D Hess; if (use2ndDer) for (size_t b = 0; b < m_basis.size(); ++b) { @@ -906,24 +845,14 @@ bool ASMu3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, // Set up control point (nodal) coordinates for current element if (!this->getElementCoordinates(Xnod,els[geoBasis-1])) return false; - // Compute the Jacobian inverse - fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1]); - 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) - fe.grad(b+1).multiply(dNxdu[b],Jac); + // Compute Jacobian inverse of the coordinate mapping and + // basis function derivatives w.r.t. Cartesian coordinates + if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives - if (use2ndDer) { - if (!utl::Hessian(Hess,fe.hess(geoBasis),Jac,Xnod, - d2Nxdu2[geoBasis-1],fe.grad(geoBasis),true)) - return false; - - for (size_t b = 0; b < m_basis.size(); ++b) - if (b != (size_t)geoBasis) - utl::Hessian(Hess,fe.hess(b+1),Jac,Xnod, - d2Nxdu2[b],fe.grad(b+1),false); - } + if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis)) + return false; // Cartesian coordinates of current integration point fe.u = gpar[0][i]; @@ -954,13 +883,12 @@ bool ASMu3Dmx::refine (const LR::RefineData& prm, Vectors& sol) if (prm.errors.empty() && prm.elements.empty()) return true; - for (Vector& solvec : sol) { - size_t ofs = 0; + for (Vector& solvec : sol) for (size_t j = 0; j < m_basis.size(); j++) { - LR::extendControlPoints(m_basis[j].get(), solvec, nfx[j], ofs); - ofs += nfx[j]*nb[j]; + Vector bVec; + this->extractNodeVec(solvec, bVec, 0, j+1); + LR::extendControlPoints(m_basis[j].get(), bVec, nfx[j]); } - } if (doRefine(prm, refBasis.get())) { for (size_t j = 0; j < m_basis.size(); ++j) @@ -993,7 +921,7 @@ bool ASMu3Dmx::refine (const LR::RefineData& prm, Vectors& sol) } 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) { sol[i].resize(len); LR::contractControlPoints(m_basis[j].get(), sol[i], nfx[j], ofs); @@ -1002,10 +930,10 @@ bool ASMu3Dmx::refine (const LR::RefineData& prm, Vectors& sol) #ifdef SP_DEBUG std::cout <<"Refined mesh: "; - for (const auto& it : m_basis) + for (const SplinePtr& it : m_basis) std::cout << it->nElements() <<" "; std::cout <<"elements "; - for (const auto& it : m_basis) + for (const SplinePtr& it : m_basis) std::cout << it->nBasisFunctions() <<" "; std::cout <<"nodes."<< std::endl; std::cout << "Projection basis: " @@ -1042,48 +970,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, bool ignoreGlobalLM) { @@ -1112,7 +998,7 @@ void ASMu3Dmx::generateThreadGroups (const Integrand& integrand, bool silence, LR::generateThreadGroups(altProjThreadGroups,altProjBasis.get()); std::vector bases; - for (const std::shared_ptr& basis : m_basis) + for (const SplinePtr& basis : m_basis) bases.push_back(basis.get()); if (silence || threadGroups[0].size() < 2) return; @@ -1131,6 +1017,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 () { if (altProjBasis) { @@ -1139,3 +1068,20 @@ void ASMu3Dmx::swapProjectionBasis () std::swap(projThreadGroups, altProjThreadGroups); } } + + +void ASMu3Dmx::getElementsAt (const RealArray& param, + std::vector& elms, + std::vector& sizes) const +{ + elms.clear(); + sizes.clear(); + elms.reserve(m_basis.size()); + sizes.reserve(m_basis.size()); + for (const SplinePtr& basis : m_basis) + { + int iel = basis->getElementContaining(param); + elms.push_back(1+iel); + sizes.push_back(basis->getElement(iel)->nBasisFunctions()); + } +} diff --git a/src/ASM/LR/ASMu3Dmx.h b/src/ASM/LR/ASMu3Dmx.h index c372df742..a5b14c12e 100644 --- a/src/ASM/LR/ASMu3Dmx.h +++ b/src/ASM/LR/ASMu3Dmx.h @@ -184,6 +184,15 @@ class ASMu3Dmx : public ASMu3D, private ASMmxBase //! \brief Swaps between the main and alternative projection basis. virtual void swapProjectionBasis(); +private: + //! \brief Finds the elements and associted sizes at given parametric point. + //! \param[in] param Parametric point to find elements at + //! \param[out] elms List of elements on each basis containign given point + //! \param[out] sizes List of element sizes (numer of element nodes) + void getElementsAt(const RealArray& param, + std::vector& elms, + std::vector& sizes) const; + protected: //! \brief Generates element groups for multi-threading of interior integrals. //! \param[in] integrand Object with problem-specific data and methods @@ -193,12 +202,16 @@ class ASMu3Dmx : public ASMu3D, private ASMmxBase bool ignoreGlobalLM); private: - std::vector> m_basis; //!< All bases - LR::LRSplineVolume* threadBasis; //!< Basis for thread groups - std::shared_ptr refBasis; //!< Basis to refine based on - std::shared_ptr altProjBasis; //!< Alternative projection basis + typedef std::shared_ptr SplinePtr; //!< Pointer to spline + + std::vector m_basis; //!< All bases + LR::LRSplineVolume* threadBasis; //!< Basis for thread groups + SplinePtr refBasis; //!< Basis to refine based on + SplinePtr altProjBasis; //!< Alternative projection basis + const std::vector& bezierExtractmx; //!< Bezier extraction matrices std::vector myBezierExtractmx; //!< Bezier extraction matrices + ThreadGroups altProjThreadGroups; //!< Element groups for multi-threaded assembly - alternative projection basis };