From 2f676bf6b395a1968e25a8f6bcd3a42cbe1ff891 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 9 Jun 2022 14:35:45 +0200 Subject: [PATCH] added: BasisFunctionCache in ASMu2Dmx --- src/ASM/LR/ASMu2Dmx.C | 108 +++++++++++++++++++++++++++--------------- src/ASM/LR/ASMu2Dmx.h | 27 +++++++++++ 2 files changed, 96 insertions(+), 39 deletions(-) diff --git a/src/ASM/LR/ASMu2Dmx.C b/src/ASM/LR/ASMu2Dmx.C index 4119aa03..79da3f2f 100644 --- a/src/ASM/LR/ASMu2Dmx.C +++ b/src/ASM/LR/ASMu2Dmx.C @@ -321,11 +321,6 @@ bool ASMu2Dmx::integrate (Integrand& integrand, PROFILE2("ASMu2Dmx::integrate(I)"); - // Get Gaussian quadrature points and weights - const double* xg = GaussQuadrature::getCoord(nGauss); - const double* wg = GaussQuadrature::getWeight(nGauss); - if (!xg || !wg) return false; - if (integrand.getReducedIntegration(nGauss)) { std::cerr <<" *** Reduced integration not available for mixed LR splines" @@ -335,6 +330,25 @@ bool ASMu2Dmx::integrate (Integrand& integrand, bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES; + if (myCache.empty()) { + myCache.emplace_back(std::make_unique(*this, cachePolicy, 1)); + for (size_t b = 2; b <= this->getNoBasis(); ++b) { + const BasisFunctionCache& c = static_cast(*myCache.front()); + myCache.emplace_back(std::make_unique(c,b)); + } + } + + for (std::unique_ptr& cache : myCache) { + cache->setIntegrand(&integrand); + cache->init(use2ndDer ? 2 : 1); + } + + ASMu2D::BasisFunctionCache& cache = *myCache.front(); + + const std::array& ng = cache.nGauss(); + const std::array& xg = cache.coord(); + const std::array& wg = cache.weight(); + ThreadGroups oneGroup; if (glInt.threadSafe()) oneGroup.oneGroup(nel); const IntMat& groups = glInt.threadSafe() ? oneGroup[0] : threadGroups[0]; @@ -354,9 +368,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand, std::vector elem_sizes; this->getElementsAt(threadBasis->getElement(groups[t][e])->midpoint(),els,elem_sizes); - MxFiniteElement fe(elem_sizes); - std::vector dNxdu(m_basis.size()); - std::vector d2Nxdu2(use2ndDer ? m_basis.size() : 0); + MxFiniteElement fe(elem_sizes); Matrix Xnod, Jac; Matrix3D Hess; double dXidu[2]; @@ -391,11 +403,6 @@ bool ASMu2Dmx::integrate (Integrand& integrand, dXidu[1] = geo->getElement(geoEl-1)->vmax()-geo->getElement(geoEl-1)->vmin(); } - // Compute parameter values of the Gauss points over this element - std::array gpar; - for (int d = 0; d < 2; d++) - this->getGaussPointParameters(gpar[d],d,nGauss,geoEl,xg); - // Initialize element quantities LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel); if (!integrand.initElement(MNPC[geoEl-1],fe,elem_sizes,nb,*A)) @@ -407,41 +414,34 @@ bool ASMu2Dmx::integrate (Integrand& integrand, // --- Integration loop over all Gauss points in each direction ---------- - int jp = (geoEl-1)*nGauss*nGauss; + int jp = (geoEl-1)*ng[0]*ng[1]; fe.iGP = firstIp + jp; // Global integration point counter - for (int j = 0; j < nGauss; j++) - for (int i = 0; i < nGauss; i++, fe.iGP++) + size_t ig = 0; + for (int j = 0; j < ng[1]; j++) + for (int i = 0; i < ng[0]; i++, fe.iGP++, ig++) { // Local element coordinates of current integration point - fe.xi = xg[i]; - fe.eta = xg[j]; + fe.xi = xg[0][i]; + fe.eta = xg[1][j]; // Parameter values of current integration point - fe.u = param[0] = gpar[0][i]; - fe.v = param[1] = gpar[1][j]; + fe.u = param[0] = cache.getParam(0,geoEl-1,i); + fe.v = param[1] = cache.getParam(1,geoEl-1,j); - // Compute basis function derivatives at current integration point - if (use2ndDer) - for (size_t b = 0; b < m_basis.size(); ++b) { - Go::BasisDerivsSf2 spline; - this->computeBasis(fe.u,fe.v,spline,els[b]-1,m_basis[b].get()); - SplineUtils::extractBasis(spline,fe.basis(b+1),dNxdu[b],d2Nxdu2[b]); - } - else - for (size_t b=0; b < m_basis.size(); ++b) { - Go::BasisDerivsSf spline; - this->computeBasis(fe.u, fe.v, spline, els[b]-1, m_basis[b].get()); - SplineUtils::extractBasis(spline,fe.basis(b+1),dNxdu[b]); - } + std::vector bfs(this->getNoBasis()); + for (size_t b = 0; b < m_basis.size(); ++b) { + bfs[b] = &myCache[b]->getVals(geoEl-1,ig); + fe.basis(b+1) = bfs[b]->N; + } // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates - if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + if (!fe.Jacobian(Jac,Xnod,geoBasis,&bfs)) continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives - if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis)) + if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,geoBasis,&bfs)) ok = false; // Compute G-matrix @@ -452,7 +452,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand, X.assign(Xnod * fe.basis(geoBasis)); // Evaluate the integrand and accumulate element contributions - fe.detJxW *= dA*wg[i]*wg[j]; + fe.detJxW *= dA*wg[0][i]*wg[1][j]; if (!integrand.evalIntMx(*A,fe,time,X)) ok = false; } @@ -468,6 +468,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand, A->destruct(); } + for (std::unique_ptr& cache : myCache) + cache->finalizeAssembly(); return ok; } @@ -545,8 +547,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex, std::vector elem_sizes; this->getElementsAt(el1->midpoint(),els,elem_sizes); - int geoEl = els[geoBasis-1]; MxFiniteElement fe(elem_sizes,firstp); + int geoEl = els[geoBasis-1]; fe.iel = MLGE[geoEl-1]; fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0; firstp += nGP; // Global integration point counter @@ -933,11 +935,11 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand, // Compute Jacobian inverse of the coordinate mapping and // basis function derivatives w.r.t. Cartesian coordinates - if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + if (!fe.Jacobian(Jac,Xnod,geoBasis,nullptr,&dNxdu)) continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives - if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,d2Nxdu2,geoBasis)) + if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,geoBasis,nullptr,&d2Nxdu2)) return false; // Cartesian coordinates of current integration point @@ -1242,3 +1244,31 @@ void ASMu2Dmx::getElementsAt (const RealArray& param, sizes.push_back(basis->getElement(iel)->nBasisFunctions()); } } + + +BasisFunctionVals ASMu2Dmx::BasisFunctionCache::calculatePt (size_t el, + size_t gp, + bool reduced) const +{ + const std::array gpIdx = this->gpIndex(gp,reduced); + double u = this->getParam(0,el,gpIdx[0],reduced); + double v = this->getParam(1,el,gpIdx[1],reduced); + + const ASMu2Dmx& pch = static_cast(patch); + + const LR::Element* el1 = pch.getBasis(ASMmxBase::geoBasis)->getElement(el); + size_t el_b = patch.getBasis(basis)->getElementContaining(el1->midpoint()); + + BasisFunctionVals result; + if (nderiv == 1 || reduced) { + Go::BasisDerivsSf spline; + pch.computeBasis(u,v,spline,el_b,patch.getBasis(basis)); + SplineUtils::extractBasis(spline,result.N,result.dNdu); + } else if (nderiv == 2) { + Go::BasisDerivsSf2 spline; + pch.computeBasis(u,v,spline,el_b,patch.getBasis(basis)); + SplineUtils::extractBasis(spline,result.N,result.dNdu,result.d2Ndu2); + } + + return result; +} diff --git a/src/ASM/LR/ASMu2Dmx.h b/src/ASM/LR/ASMu2Dmx.h index ad15939a..2a39c76f 100644 --- a/src/ASM/LR/ASMu2Dmx.h +++ b/src/ASM/LR/ASMu2Dmx.h @@ -30,6 +30,33 @@ class ASMu2Dmx : public ASMu2D, private ASMmxBase { +protected: + //! \brief Implementation of basis function cache. + class BasisFunctionCache : public ASMu2D::BasisFunctionCache + { + public: + //! \brief The constructor initializes the class. + //! \param pch Patch the cache is for + //! \param plcy Cache policy to use + //! \param b Basis to use + BasisFunctionCache(const ASMu2D& pch, + ASM::CachePolicy plcy, int b) + : ASMu2D::BasisFunctionCache(pch,plcy,b) {} + + //! \brief Constructor reusing quadrature info from another instance. + //! \param cache Instance holding quadrature information + //! \param b Basis to use + BasisFunctionCache(const BasisFunctionCache& cache, int b) + : ASMu2D::BasisFunctionCache(cache,b) {} + + protected: + //! \brief Calculates basis function info in a single integration point. + //! \param el Element of integration point (0-indexed) + //! \param gp Integration point on element (0-indexed) + //! \param reduced If true, returns values for reduced integration scheme + BasisFunctionVals calculatePt(size_t el, size_t gp, bool reduced) const override; + }; + public: //! \brief The constructor initializes the dimension of each basis. ASMu2Dmx(unsigned char n_s, const CharVec& n_f);