From eff2b08e037ccc75d865722d8e89f248823ad88b Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 8 Jun 2022 13:29:06 +0200 Subject: [PATCH] added: BasisFunctionCache in ASMs2DLag --- src/ASM/ASMs2DLag.C | 156 +++++++++++++++++++++++++++++++------------- src/ASM/ASMs2DLag.h | 47 +++++++++++++ 2 files changed, 157 insertions(+), 46 deletions(-) diff --git a/src/ASM/ASMs2DLag.C b/src/ASM/ASMs2DLag.C index 9c645d15..6c233cbd 100644 --- a/src/ASM/ASMs2DLag.C +++ b/src/ASM/ASMs2DLag.C @@ -305,37 +305,25 @@ bool ASMs2DLag::integrate (Integrand& integrand, { if (this->empty()) return true; // silently ignore empty patches + if (myCache.empty()) + myCache.emplace_back(std::make_unique(*this, cachePolicy, 1)); + + BasisFunctionCache& cache = static_cast(*myCache.front()); + cache.setIntegrand(&integrand); + if (!cache.init(1)) + return false; + // Get Gaussian quadrature points and weights - std::array ng; - std::array xg, wg; - for (int d = 0; d < 2; d++) - { - ng[d] = this->getNoGaussPt(d == 0 ? p1 : p2); - xg[d] = GaussQuadrature::getCoord(ng[d]); - wg[d] = GaussQuadrature::getWeight(ng[d]); - if (!xg[d] || !wg[d]) return false; - } + const std::array& ng = cache.nGauss(); + const std::array& xg = cache.coord(); + const std::array& wg = cache.weight(); // Get the reduced integration quadrature points, if needed - const double* xr = nullptr; - const double* wr = nullptr; - int nRed = integrand.getReducedIntegration(ng[0]); - if (nRed > 0) - { - xr = GaussQuadrature::getCoord(nRed); - wr = GaussQuadrature::getWeight(nRed); - if (!xr || !wr) return false; - } - else if (nRed < 0) - nRed = ng[0]; // The integrand needs to know nGauss - - // Get parametric coordinates of the elements - RealArray upar, vpar; - this->getGridParameters(upar,0,1); - this->getGridParameters(vpar,1,1); + const double* xr = cache.coord(true)[0]; + const double* wr = cache.weight(true)[0]; // Number of elements in each direction - const int nelx = upar.empty() ? 0 : upar.size() - 1; + const int nelx = cache.noElms()[0]; // === Assembly loop over all elements in the patch ========================== @@ -376,6 +364,7 @@ bool ASMs2DLag::integrate (Integrand& integrand, // Initialize element quantities fe.iel = MLGE[iel]; LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel); + int nRed = cache.nGauss(true)[0]; if (!integrand.initElement(MNPC[iel],fe,X,nRed*nRed,*A)) { A->destruct(); @@ -387,26 +376,23 @@ bool ASMs2DLag::integrate (Integrand& integrand, { // --- Selective reduced integration loop ---------------------------- + size_t ip = 0; for (int j = 0; j < nRed; j++) - for (int i = 0; i < nRed; i++) + for (int i = 0; i < nRed; i++, ++ip) { // Local element coordinates of current integration point fe.xi = xr[i]; fe.eta = xr[j]; // Parameter value of current integration point - if (!upar.empty()) - fe.u = 0.5*(upar[i1]*(1.0-xr[i]) + upar[i1+1]*(1.0+xr[i])); - if (!vpar.empty()) - fe.v = 0.5*(vpar[i2]*(1.0-xr[j]) + vpar[i2+1]*(1.0+xr[j])); + fe.u = cache.getParam(0,i1,i,true); + fe.v = cache.getParam(1,i2,j,true); - // Compute basis function derivatives at current point - // using tensor product of one-dimensional Lagrange polynomials - if (!Lagrange::computeBasis(fe.N,dNdu,p1,xr[i],p2,xr[j])) - ok = false; + const BasisFunctionVals& bfs = cache.getVals(iel,ip,true); + fe.N = bfs.N; // Compute Jacobian inverse and derivatives - fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,bfs.dNdu); // Cartesian coordinates of current integration point X.assign(Xnod * fe.N); @@ -421,29 +407,26 @@ bool ASMs2DLag::integrate (Integrand& integrand, // --- Integration loop over all Gauss points in each direction -------- + size_t ip = 0; int jp = iel*ng[0]*ng[1]; fe.iGP = firstIp + jp; // Global integration point counter for (int j = 0; j < ng[1]; j++) - for (int i = 0; i < ng[0]; i++, fe.iGP++) + for (int i = 0; i < ng[0]; i++, fe.iGP++, ++ip) { // Local element coordinates of current integration point fe.xi = xg[0][i]; fe.eta = xg[1][j]; // Parameter value of current integration point - if (!upar.empty()) - fe.u = 0.5*(upar[i1]*(1.0-xg[0][i]) + upar[i1+1]*(1.0+xg[0][i])); - if (!vpar.empty()) - fe.v = 0.5*(vpar[i2]*(1.0-xg[1][j]) + vpar[i2+1]*(1.0+xg[1][j])); + fe.u = cache.getParam(0,i1,i); + fe.v = cache.getParam(1,i2,j); - // Compute basis function derivatives at current integration point - // using tensor product of one-dimensional Lagrange polynomials - if (!Lagrange::computeBasis(fe.N,dNdu,p1,xg[0][i],p2,xg[1][j])) - ok = false; + const BasisFunctionVals& bfs = cache.getVals(iel,ip); + fe.N = bfs.N; // Compute Jacobian inverse of coordinate mapping and derivatives - fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); + fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,bfs.dNdu); if (fe.detJxW == 0.0) continue; // skip singular points // Cartesian coordinates of current integration point @@ -468,6 +451,7 @@ bool ASMs2DLag::integrate (Integrand& integrand, } } + cache.finalizeAssembly(); return ok; } @@ -848,3 +832,83 @@ bool ASMs2DLag::evaluate (const ASMbase* basis, const Vector& locVec, vec = sValues; return true; } + + +ASMs2DLag::BasisFunctionCache::BasisFunctionCache (const ASMs2DLag& pch, + ASM::CachePolicy plcy, + int b) : + ASMs2D::BasisFunctionCache(pch,plcy,b) +{ + nel[0] = (pch.nx-1) / (pch.p1-1); + nel[1] = (pch.ny-1) / (pch.p2-1); +} + + +ASMs2DLag::BasisFunctionCache::BasisFunctionCache (const BasisFunctionCache& cache, + int b) : + ASMs2D::BasisFunctionCache(cache,b) +{ +} + + +void ASMs2DLag::BasisFunctionCache::setupParameters () +{ + // Compute parameter values of the Gauss points over the whole patch + for (int d = 0; d < 2; d++) { + RealArray par; + patch.getGridParameters(par,d,1); + mainQ->gpar[d].resize(par.size(), 1); + mainQ->gpar[d].fillColumn(1,par.data()); + if (reducedQ->xg[0]) + reducedQ->gpar[d] = mainQ->gpar[d]; + } +} + + +double ASMs2DLag::BasisFunctionCache::getParam (int dir, size_t el, + size_t gp, bool reduced) const +{ + const Quadrature& q = reduced ? *reducedQ : *mainQ; + return 0.5*(q.gpar[dir](el+1,1)*(1.0-q.xg[dir][gp]) + + q.gpar[dir](el+2,1)*(1.0+q.xg[dir][gp])); +} + + +BasisFunctionVals ASMs2DLag::BasisFunctionCache::calculatePt (size_t el, + size_t gp, + bool reduced) const +{ + std::array gpIdx = this->gpIndex(gp,reduced); + const Quadrature& q = reduced ? *reducedQ : *mainQ; + + const ASMs2DLag& pch = static_cast(patch); + + BasisFunctionVals result; + if (nderiv == 1) + Lagrange::computeBasis(result.N,result.dNdu, + pch.p1,q.xg[0][gpIdx[0]], + pch.p2,q.xg[1][gpIdx[1]]); + + return result; +} + + +void ASMs2DLag::BasisFunctionCache::calculateAll () +{ + // Evaluate basis function values and derivatives at all integration points. + // We do this before the integration point loop to exploit multi-threading + // in the integrand evaluations, which may be the computational bottleneck. + size_t iel, jp, rp; + const ASMs2DLag& pch = static_cast(patch); + for (iel = jp = rp = 0; iel < pch.nel; iel++) + { + for (int j = 0; j < mainQ->ng[1]; j++) + for (int i = 0; i < mainQ->ng[0]; i++, jp++) + values[jp] = this->calculatePt(iel,j*mainQ->ng[0]+i,false); + + if (reducedQ->xg[0]) + for (int j = 0; j < reducedQ->ng[1]; j++) + for (int i = 0; i < reducedQ->ng[0]; i++, rp++) + valuesRed[rp] = this->calculatePt(iel,j*reducedQ->ng[0]+i,true); + } +} diff --git a/src/ASM/ASMs2DLag.h b/src/ASM/ASMs2DLag.h index 50a700eb..14cdf748 100644 --- a/src/ASM/ASMs2DLag.h +++ b/src/ASM/ASMs2DLag.h @@ -25,6 +25,53 @@ class ASMs2DLag : public ASMs2D { +protected: + //! \brief Implementation of basis function cache. + class BasisFunctionCache : public ASMs2D::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 ASMs2DLag& pch, ASM::CachePolicy plcy, int 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); + + //! \brief Empty destructor. + virtual ~BasisFunctionCache() = default; + + //! \brief Obtain a single integration point parameter. + //! \param dir Direction of for integration point + //! \param el Element number in given direction + //! \param gp Integration point in given direction + //! \param reduced True to return parameter for reduced quadrature + double getParam(int dir, size_t el, size_t gp, bool reduced = false) const override; + + protected: + //! \brief Obtain global integration point index. + //! \param el Element of integration point (0-indexed) + //! \param gp Integration point on element (0-indexed) + //! \param reduced True to return index for reduced quadrature + size_t index(size_t el, size_t gp, bool reduced) const override + { return ::BasisFunctionCache<2>::index(el,gp,reduced); } + + //! \brief Calculates basis function info in a single integration point. + //! \param el Element of integration point (0-indexed) + //! \param gp Integratin 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; + + //! \brief Calculates basis function info in all integration points. + void calculateAll() override; + + //! \brief Configure quadrature points. + void setupParameters() override; + }; + public: //! \brief Default constructor. ASMs2DLag(unsigned char n_s = 2, unsigned char n_f = 2);