From a6debc90af2f1c684b5ef0e10d8c470b2f3e671b 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 ASMs2DmxLag --- src/ASM/ASMs2DmxLag.C | 96 +++++++++++++++++++++++++++++++------------ src/ASM/ASMs2DmxLag.h | 29 +++++++++++++ 2 files changed, 98 insertions(+), 27 deletions(-) diff --git a/src/ASM/ASMs2DmxLag.C b/src/ASM/ASMs2DmxLag.C index 5bc280c8..f539c140 100644 --- a/src/ASM/ASMs2DmxLag.C +++ b/src/ASM/ASMs2DmxLag.C @@ -238,17 +238,26 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, { if (this->empty()) return true; // silently ignore empty patches + if (myCache.empty()) { + myCache.emplace_back(std::make_unique(*this, cachePolicy, 1)); + const BasisFunctionCache& front = static_cast(*myCache.front()); + for (size_t b = 2; b <= this->getNoBasis(); ++b) + myCache.emplace_back(std::make_unique(front, b)); + } + + for (std::unique_ptr& cache : myCache) { + cache->setIntegrand(&integrand); + cache->init(1); + } + + BasisFunctionCache& cache = static_cast(*myCache.front()); + // Get Gaussian quadrature points and weights - const double* xg = GaussQuadrature::getCoord(nGauss); - const double* wg = GaussQuadrature::getWeight(nGauss); - if (!xg || !wg) return false; + const std::array& ng = cache.nGauss(); + const std::array& xg = cache.coord(); + const std::array& wg = cache.weight(); - // Get parametric coordinates of the elements - RealArray upar, vpar; - this->getGridParameters(upar,0,1); - this->getGridParameters(vpar,1,1); - - const int nelx = upar.size() - 1; + const int nelx = cache.noElms()[0]; // === Assembly loop over all elements in the patch ========================== @@ -260,7 +269,6 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, for (size_t t = 0; t < threadGroups[g].size(); t++) { MxFiniteElement fe(elem_size); - Matrices dNxdu(nxx.size()); Matrix Xnod, Jac; Vec4 X(nullptr,time.t); for (size_t e = 0; e < threadGroups[g][t].size() && ok; ++e) @@ -288,37 +296,37 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, // --- Integration loop over all Gauss points in each direction -------- - int jp = (i2*nelx + i1)*nGauss*nGauss; + size_t ip = 0; + int jp = (i2*nelx + i1)*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++) + for (int j = 0; j < ng[1]; j++) + for (int i = 0; i < ng[0]; i++, fe.iGP++, ++ip) { // Parameter value of current integration point - fe.u = 0.5*(upar[i1]*(1.0-xg[i]) + upar[i1+1]*(1.0+xg[i])); - fe.v = 0.5*(vpar[i2]*(1.0-xg[j]) + vpar[i2+1]*(1.0+xg[j])); + fe.u = cache.getParam(0,i1,i); + fe.v = cache.getParam(1,i2,j); // Local coordinates of current integration point - fe.xi = xg[i]; - fe.eta = xg[j]; + fe.xi = xg[0][i]; + fe.eta = xg[1][j]; - // Compute basis function derivatives at current integration point - // using tensor product of one-dimensional Lagrange polynomials - for (size_t b = 0; b < nxx.size(); ++b) - if (!Lagrange::computeBasis(fe.basis(b+1),dNxdu[b], - elem_sizes[b][0],xg[i], - elem_sizes[b][1],xg[j])) - ok = false; + // Fetch basis function derivatives at current integration point + std::vector bfs(this->getNoBasis()); + for (size_t b = 0; b < this->getNoBasis(); ++b) { + bfs[b] = &myCache[b]->getVals(iel-1,ip); + fe.basis(b+1) = bfs[b]->N; + } // Compute Jacobian inverse of coordinate mapping and derivatives - if (!fe.Jacobian(Jac,Xnod,dNxdu,geoBasis)) + if (!fe.Jacobian(Jac,Xnod,geoBasis,&bfs)) continue; // skip singular points // Cartesian coordinates of current integration point X.assign(Xnod * fe.basis(geoBasis)); // Evaluate the integrand and accumulate element contributions - fe.detJxW *= wg[i]*wg[j]; + fe.detJxW *= wg[0][i]*wg[1][j]; if (!integrand.evalIntMx(*A,fe,time,X)) ok = false; } @@ -525,7 +533,7 @@ bool ASMs2DmxLag::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 // Now evaluate the solution field @@ -547,3 +555,37 @@ bool ASMs2DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand, return true; } + + +ASMs2DmxLag::BasisFunctionCache::BasisFunctionCache (const ASMs2DLag& pch, + ASM::CachePolicy plcy, + int b) : + ASMs2DLag::BasisFunctionCache(pch,plcy,b) +{ +} + + +ASMs2DmxLag::BasisFunctionCache::BasisFunctionCache (const BasisFunctionCache& cache, + int b) : + ASMs2DLag::BasisFunctionCache(cache,b) +{ +} + + +BasisFunctionVals ASMs2DmxLag::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 ASMs2DmxLag& pch = static_cast(patch); + + BasisFunctionVals result; + if (nderiv == 1) + Lagrange::computeBasis(result.N,result.dNdu, + pch.elem_sizes[basis-1][0],q.xg[0][gpIdx[0]], + pch.elem_sizes[basis-1][1],q.xg[1][gpIdx[1]]); + + return result; +} diff --git a/src/ASM/ASMs2DmxLag.h b/src/ASM/ASMs2DmxLag.h index 771e6679..f860f9fa 100644 --- a/src/ASM/ASMs2DmxLag.h +++ b/src/ASM/ASMs2DmxLag.h @@ -28,6 +28,33 @@ class ASMs2DmxLag : public ASMs2DLag, private ASMmxBase { +protected: + //! \brief Implementation of basis function cache. + class BasisFunctionCache : public ASMs2DLag::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; + + protected: + //! \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; + }; + public: //! \brief The constructor initializes the dimension of each basis. ASMs2DmxLag(unsigned char n_s, const CharVec& n_f); @@ -50,6 +77,8 @@ public: //! This is used to reinitialize the patch after it has been refined. virtual void clear(bool retainGeometry); + //! \brief Returns the number of bases. + virtual size_t getNoBasis() const { return 2; } //! \brief Returns the total number of nodes in this patch. virtual size_t getNoNodes(int basis) const; //! \brief Returns the number of solution fields.