From a2472712416996d4e162c9ab4e9c010961a61264 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 ASMs3D --- src/ASM/ASMs3D.C | 282 ++++++++++++++++++++++++++++++++++------------- src/ASM/ASMs3D.h | 70 +++++++++++- 2 files changed, 276 insertions(+), 76 deletions(-) diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C index 5e081f30..547634d0 100644 --- a/src/ASM/ASMs3D.C +++ b/src/ASM/ASMs3D.C @@ -160,6 +160,8 @@ void ASMs3D::clear (bool retainGeometry) threadGroupsVol[0].clear(); threadGroupsVol[1].clear(); threadGroupsFace.clear(); + + myCache.clear(); } @@ -1980,52 +1982,22 @@ bool ASMs3D::integrate (Integrand& integrand, bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES; bool useElmVtx = integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS; + if (myCache.empty()) + myCache.emplace_back(std::make_unique(*this, cachePolicy, 1)); + + BasisFunctionCache& cache = *myCache.front(); + cache.setIntegrand(&integrand); + if (!cache.init(use2ndDer ? 2 : 1)) + return false; + // Get Gaussian quadrature points and weights - std::array ng; - std::array xg, wg; - for (int d = 0; d < 3; d++) - { - ng[d] = this->getNoGaussPt(svol->order(d)); - 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 - - // Compute parameter values of the Gauss points over the whole patch - std::array gpar, redpar; - for (int d = 0; d < 3; d++) - { - this->getGaussPointParameters(gpar[d],d,ng[d],xg[d]); - if (xr) - this->getGaussPointParameters(redpar[d],d,nRed,xr); - } - - // Evaluate basis function derivatives at all integration points - std::vector spline; - std::vector spline2; - std::vector splineRed; - { - PROFILE2("Spline evaluation"); - if (use2ndDer) - svol->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline2); - else - svol->computeBasisGrid(gpar[0],gpar[1],gpar[2],spline); - if (xr) - svol->computeBasisGrid(redpar[0],redpar[1],redpar[2],splineRed); - } + const double* xr = cache.coord(true)[0]; + const double* wr = cache.weight(true)[0]; const int n1 = svol->numCoefs(0); const int n2 = svol->numCoefs(1); @@ -2102,21 +2074,21 @@ bool ASMs3D::integrate (Integrand& integrand, fe.Navg.resize(p1*p2*p3,true); double vol = 0.0; - int ip = (((i3-p3)*ng[2]*nel2 + i2-p2)*ng[1]*nel1 + i1-p1)*ng[0]; - for (int k = 0; k < ng[2]; k++, ip += ng[1]*(nel2-1)*ng[0]*nel1) - for (int j = 0; j < ng[1]; j++, ip += ng[0]*(nel1-1)) + size_t ip = 0; + for (int k = 0; k < ng[2]; k++) + for (int j = 0; j < ng[1]; j++) for (int i = 0; i < ng[0]; i++, ip++) { // Fetch basis function derivatives at current integration point - SplineUtils::extractBasis(spline[ip],fe.N,dNdu); + const BasisFunctionVals& bfs = cache.getVals(iel-1,ip); // Compute Jacobian determinant of coordinate mapping // and multiply by weight of current integration point - double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false); + double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,bfs.dNdu,false); double weight = dV*wg[0][i]*wg[1][j]*wg[2][k]; // Numerical quadrature - fe.Navg.add(fe.N,detJac*weight); + fe.Navg.add(bfs.N,detJac*weight); vol += detJac*weight; } @@ -2127,9 +2099,9 @@ bool ASMs3D::integrate (Integrand& integrand, if (integrand.getIntegrandType() & Integrand::ELEMENT_CENTER) { // Compute the element center - param[0] = 0.5*(gpar[0](1,i1-p1+1) + gpar[0](ng[0],i1-p1+1)); - param[1] = 0.5*(gpar[1](1,i2-p2+1) + gpar[1](ng[1],i2-p2+1)); - param[2] = 0.5*(gpar[2](1,i3-p3+1) + gpar[2](ng[2],i3-p3+1)); + param[0] = 0.5*(cache.getParam(0,i1-p1,0) + cache.getParam(0,i1-p1,ng[0]-1)); + param[1] = 0.5*(cache.getParam(1,i2-p2,0) + cache.getParam(1,i2-p2,ng[1]-1)); + param[2] = 0.5*(cache.getParam(2,i3-p3,0) + cache.getParam(2,i3-p3,ng[2]-1)); SplineUtils::point(X,param[0],param[1],param[2],svol); if (!useElmVtx) { @@ -2143,6 +2115,7 @@ bool ASMs3D::integrate (Integrand& integrand, // Initialize element quantities LocalIntegral* A = integrand.getLocalIntegral(fe.N.size(),fe.iel); + int nRed = cache.nGauss(true)[0]; if (!integrand.initElement(MNPC[iel-1],fe,X,nRed*nRed*nRed,*A)) { A->destruct(); @@ -2154,9 +2127,9 @@ bool ASMs3D::integrate (Integrand& integrand, { // --- Selective reduced integration loop ---------------------------- - int ip = (((i3-p3)*nRed*nel2 + i2-p2)*nRed*nel1 + i1-p1)*nRed; - for (int k = 0; k < nRed; k++, ip += nRed*(nel2-1)*nRed*nel1) - for (int j = 0; j < nRed; j++, ip += nRed*(nel1-1)) + size_t ip = 0; + for (int k = 0; k < nRed; k++) + for (int j = 0; j < nRed; j++) for (int i = 0; i < nRed; i++, ip++) { // Local element coordinates of current integration point @@ -2165,15 +2138,15 @@ bool ASMs3D::integrate (Integrand& integrand, fe.zeta = xr[k]; // Parameter values of current integration point - fe.u = param[0] = redpar[0](i+1,i1-p1+1); - fe.v = param[1] = redpar[1](j+1,i2-p2+1); - fe.w = param[2] = redpar[2](k+1,i3-p3+1); + fe.u = param[0] = cache.getParam(0,i1-p1,i,true); + fe.v = param[1] = cache.getParam(1,i2-p2,j,true); + fe.w = param[2] = cache.getParam(2,i3-p3,k,true); - // Fetch basis function derivatives at current point - SplineUtils::extractBasis(splineRed[ip],fe.N,dNdu); + const BasisFunctionVals& bfs = cache.getVals(iel-1,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); @@ -2188,13 +2161,12 @@ bool ASMs3D::integrate (Integrand& integrand, // --- Integration loop over all Gauss points in each direction -------- - int ip = (((i3-p3)*ng[2]*nel2 + i2-p2)*ng[1]*nel1 + i1-p1)*ng[0]; + size_t ip = 0; int jp = (((i3-p3)*nel2 + i2-p2)*nel1 + i1-p1)*ng[0]*ng[1]*ng[2]; fe.iGP = firstIp + jp; // Global integration point counter - - for (int k = 0; k < ng[2]; k++, ip += ng[1]*(nel2-1)*ng[0]*nel1) - for (int j = 0; j < ng[1]; j++, ip += ng[0]*(nel1-1)) - for (int i = 0; i < ng[0]; i++, ip++, fe.iGP++) + for (int k = 0; k < ng[2]; k++) + for (int j = 0; j < ng[1]; j++) + for (int i = 0; i < ng[0]; i++, ++ip, ++fe.iGP) { // Local element coordinates of current integration point fe.xi = xg[0][i]; @@ -2202,23 +2174,21 @@ bool ASMs3D::integrate (Integrand& integrand, fe.zeta = xg[2][k]; // Parameter values of current integration point - fe.u = param[0] = gpar[0](i+1,i1-p1+1); - fe.v = param[1] = gpar[1](j+1,i2-p2+1); - fe.w = param[2] = gpar[2](k+1,i3-p3+1); + fe.u = param[0] = cache.getParam(0,i1-p1,i); + fe.v = param[1] = cache.getParam(1,i2-p2,j); + fe.w = param[2] = cache.getParam(2,i3-p3,k); // Fetch basis function derivatives at current integration point - if (use2ndDer) - SplineUtils::extractBasis(spline2[ip],fe.N,dNdu,d2Ndu2); - else - SplineUtils::extractBasis(spline[ip],fe.N,dNdu); + const BasisFunctionVals& bfs = cache.getVals(iel-1,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 // Compute Hessian of coordinate mapping and 2nd order derivatives if (use2ndDer) - if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,fe.dNdX)) + if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,bfs.d2Ndu2,fe.dNdX)) ok = false; // Compute G-matrix @@ -2259,6 +2229,7 @@ bool ASMs3D::integrate (Integrand& integrand, } } + cache.finalizeAssembly(); return ok; } @@ -3793,3 +3764,164 @@ bool ASMs3D::addRigidCpl (int lindx, int ldim, int basis, return this->ASMstruct::addRigidCpl(lindx,ldim,basis,gMaster,Xmaster,extraPt); } + + +ASMs3D::BasisFunctionCache::BasisFunctionCache (const ASMs3D& pch, + ASM::CachePolicy plcy, + int b) : + ::BasisFunctionCache<3>(plcy), + patch(pch), + basis(b) +{ + for (size_t d = 0; d < 3 && patch.svol; ++d) + nel[d] = patch.svol->numCoefs(d) - patch.svol->order(d) + 1; +} + + +ASMs3D::BasisFunctionCache::BasisFunctionCache (const BasisFunctionCache& cache, + int b) : + ::BasisFunctionCache<3>(cache), + patch(cache.patch), + basis(b), + nel(cache.nel) +{ +} + + +bool ASMs3D::BasisFunctionCache::internalInit () +{ + if (!mainQ->xg[0]) + this->setupQuadrature(); + + nTotal = nel[0]*nel[1]*nel[2]*mainQ->ng[0]*mainQ->ng[1]*mainQ->ng[2]; + if (reducedQ->xg[0]) + nTotalRed = nel[0]*nel[1]*nel[2]*reducedQ->ng[0]*reducedQ->ng[1]*reducedQ->ng[2]; + + return true; +} + + +void ASMs3D::BasisFunctionCache::internalCleanup () +{ + if (basis == 1) { + mainQ->reset(); + reducedQ->reset(); + } +} + + +bool ASMs3D::BasisFunctionCache::setupQuadrature () +{ + // Get Gaussian quadrature points and weights + for (int d = 0; d < 3; d++) + { + mainQ->ng[d] = patch.getNoGaussPt(patch.svol ? patch.svol->order(d) : 2); + mainQ->xg[d] = GaussQuadrature::getCoord(mainQ->ng[d]); + mainQ->wg[d] = GaussQuadrature::getWeight(mainQ->ng[d]); + if (!mainQ->xg[d] || !mainQ->wg[d]) return false; + } + + // Get the reduced integration quadrature points, if needed + int nRed = integrand ? integrand->getReducedIntegration(mainQ->ng[0]) : 0; + if (nRed > 0) + { + reducedQ->xg[0] = reducedQ->xg[1] = reducedQ->xg[2] = GaussQuadrature::getCoord(nRed); + reducedQ->wg[0] = reducedQ->wg[1] = reducedQ->wg[2] = GaussQuadrature::getWeight(nRed); + if (!reducedQ->xg[0] || !reducedQ->wg[0]) return false; + } else if (nRed < 0) + nRed = mainQ->ng[0]; + + reducedQ->ng[0] = reducedQ->ng[1] = reducedQ->ng[2] = nRed; + + this->setupParameters(); + return true; +} + + +void ASMs3D::BasisFunctionCache::setupParameters () +{ + // Compute parameter values of the Gauss points over the whole patch + for (int d = 0; d < 3; d++) { + patch.getGaussPointParameters(mainQ->gpar[d],d,mainQ->ng[d],mainQ->xg[d]); + if (reducedQ->xg[0]) + patch.getGaussPointParameters(reducedQ->gpar[d],d,reducedQ->ng[d],reducedQ->xg[d]); + } +} + + +BasisFunctionVals ASMs3D::BasisFunctionCache::calculatePt (size_t el, + size_t gp, + bool reduced) const +{ + PROFILE2("Spline evaluation"); + std::array gpIdx = this->gpIndex(gp,reduced); + std::array elIdx = this->elmIndex(el); + + BasisFunctionVals result; + if (nderiv == 1) { + Go::BasisDerivs spline; +#pragma omp critical + patch.getBasis(basis)->computeBasis(this->getParam(0,elIdx[0],gpIdx[0],reduced), + this->getParam(1,elIdx[1],gpIdx[1],reduced), + this->getParam(2,elIdx[2],gpIdx[2],reduced), spline); + SplineUtils::extractBasis(spline,result.N,result.dNdu); + } else if (nderiv == 2) { + Go::BasisDerivs2 spline; +#pragma omp critical + patch.getBasis(basis)->computeBasis(this->getParam(0,elIdx[0],gpIdx[0],reduced), + this->getParam(1,elIdx[1],gpIdx[1],reduced), + this->getParam(2,elIdx[2],gpIdx[2],reduced), spline); + SplineUtils::extractBasis(spline,result.N,result.dNdu,result.d2Ndu2); + } + + return result; +} + + +size_t ASMs3D::BasisFunctionCache::index (size_t el, size_t gp, bool reduced) const +{ + const Quadrature& q = reduced ? *reducedQ : *mainQ; + std::array elIdx = this->elmIndex(el); + std::array gpIdx = this->gpIndex(gp,reduced); + + return ((elIdx[2]*q.ng[2]*nel[1] + elIdx[1])*q.ng[1]*nel[0] + elIdx[0])*q.ng[0] + + (gpIdx[2]*q.ng[1]*nel[1] + gpIdx[1])*q.ng[0]*nel[0] + gpIdx[0]; +} + + +std::array ASMs3D::BasisFunctionCache::elmIndex (size_t el) const +{ + return { el % nel[0], (el/nel[0]) % nel[1], el / (nel[0]*nel[1]) }; +} + + +void ASMs3D::BasisFunctionCache::calculateAll () +{ + PROFILE2("Spline evaluation"); + auto&& extract1 = [this](const Matrix& par1, + const Matrix& par2, + const Matrix& par3, + std::vector& values) + { + std::vector spline; + patch.getBasis(basis)->computeBasisGrid(par1,par2,par3,spline); + size_t idx = 0; + for (const Go::BasisDerivs& spl : spline) { + SplineUtils::extractBasis(spl,values[idx].N,values[idx].dNdu); + ++idx; + } + }; + if (nderiv == 1) + extract1(mainQ->gpar[0],mainQ->gpar[1],mainQ->gpar[2],values); + else if (nderiv == 2) { + std::vector spline; + patch.getBasis(basis)->computeBasisGrid(mainQ->gpar[0],mainQ->gpar[1],mainQ->gpar[2],spline); + size_t idx = 0; + for (const Go::BasisDerivs2& spl : spline) { + SplineUtils::extractBasis(spl,values[idx].N,values[idx].dNdu,values[idx].d2Ndu2); + ++idx; + } + } + if (!reducedQ->gpar[0].empty()) + extract1(reducedQ->gpar[0],reducedQ->gpar[1],reducedQ->gpar[2],valuesRed); +} diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h index aedee099..0a0385bd 100644 --- a/src/ASM/ASMs3D.h +++ b/src/ASM/ASMs3D.h @@ -16,9 +16,12 @@ #include "ASMstruct.h" #include "ASM3D.h" +#include "BasisFunctionCache.h" #include "Interface.h" #include "ThreadGroups.h" +#include + namespace utl { class Point; } @@ -73,6 +76,68 @@ class ASMs3D : public ASMstruct, public ASM3D int next(); }; +protected: + //! \brief Implementation of basis function cache. + class BasisFunctionCache : public ::BasisFunctionCache<3> + { + 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 ASMs3D& 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 Returns number of elements in each direction. + const std::array& noElms() const + { return nel; } + + protected: + //! \brief Implementation specific initialization. + bool internalInit() override; + + //! \brief Implementation specific cleanup. + void internalCleanup() override; + + //! \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 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; + + //! \brief Setup integration point parameters. + virtual void setupParameters(); + + const ASMs3D& patch; //!< Reference to patch cache is for + + int basis; //!< Basis to use + std::array nel; //!< Number of elements in each direction + +private: + //! \brief Obtain structured element indices. + //! \param el Global element index + std::array elmIndex(size_t el) const; + + //! \brief Configure quadratures. + bool setupQuadrature(); + }; + public: //! \brief Struct with data for definition of global node numbers of a patch. struct BlockNodes @@ -620,7 +685,7 @@ protected: //! \param[in] xi Dimensionless Gauss point coordinates [-1,1] //! \return The parameter value matrix casted into a one-dimensional vector const Vector& getGaussPointParameters(Matrix& uGP, int dir, int nGauss, - const double* xi) const; + const double* xi) const; //! \brief Calculates parameter values for the Greville points. //! \param[out] prm Parameter values in given direction for all points @@ -821,6 +886,9 @@ protected: ThreadGroups threadGroupsVol; //! Element groups for multi-threaded face assembly std::map threadGroupsFace; + + //! Basis function cache + std::vector> myCache; }; #endif