added: BasisFunctionCache in ASMs3D

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<BasisFunctionCache>(*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<int,3> ng;
+  const std::array<int,3>& ng = cache.nGauss();
-  std::array<const double*,3> xg, wg;
+  const std::array<const double*,3>& xg = cache.coord();
-  for (int d = 0; d < 3; d++)
+  const std::array<const double*,3>& wg = cache.weight();
-  {
-    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;
-  }
 
   // Get the reduced integration quadrature points, if needed
-  const double* xr = nullptr;
+  const double* xr = cache.coord(true)[0];
-  const double* wr = nullptr;
+  const double* wr = cache.weight(true)[0];
-  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<Matrix,3> 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<Go::BasisDerivs>  spline;
-  std::vector<Go::BasisDerivs2> spline2;
-  std::vector<Go::BasisDerivs>  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 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];
+          size_t ip = 0;
-          for (int k = 0; k < ng[2]; k++, ip += ng[1]*(nel2-1)*ng[0]*nel1)
+          for (int k = 0; k < ng[2]; k++)
-            for (int j = 0; j < ng[1]; j++, ip += ng[0]*(nel1-1))
+            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[0] = 0.5*(cache.getParam(0,i1-p1,0) + cache.getParam(0,i1-p1,ng[0]-1));
-          param[1] = 0.5*(gpar[1](1,i2-p2+1) + gpar[1](ng[1],i2-p2+1));
+          param[1] = 0.5*(cache.getParam(1,i2-p2,0) + cache.getParam(1,i2-p2,ng[1]-1));
-          param[2] = 0.5*(gpar[2](1,i3-p3+1) + gpar[2](ng[2],i3-p3+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;
+          size_t ip = 0;
-          for (int k = 0; k < nRed; k++, ip += nRed*(nel2-1)*nRed*nel1)
+          for (int k = 0; k < nRed; k++)
-            for (int j = 0; j < nRed; j++, ip += nRed*(nel1-1))
+            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.u = param[0] = cache.getParam(0,i1-p1,i,true);
-                fe.v = param[1] = redpar[1](j+1,i2-p2+1);
+                fe.v = param[1] = cache.getParam(1,i2-p2,j,true);
-                fe.w = param[2] = redpar[2](k+1,i3-p3+1);
+                fe.w = param[2] = cache.getParam(2,i3-p3,k,true);
 
-                // Fetch basis function derivatives at current point
+                const BasisFunctionVals& bfs = cache.getVals(iel-1,ip,true);
-                SplineUtils::extractBasis(splineRed[ip],fe.N,dNdu);
+                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++)
-        for (int k = 0; k < ng[2]; k++, ip += ng[1]*(nel2-1)*ng[0]*nel1)
+          for (int j = 0; j < ng[1]; j++)
-          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 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.u = param[0] = cache.getParam(0,i1-p1,i);
-              fe.v = param[1] = gpar[1](j+1,i2-p2+1);
+              fe.v = param[1] = cache.getParam(1,i2-p2,j);
-              fe.w = param[2] = gpar[2](k+1,i3-p3+1);
+              fe.w = param[2] = cache.getParam(2,i3-p3,k);
 
               // Fetch basis function derivatives at current integration point
-              if (use2ndDer)
+              const BasisFunctionVals& bfs = cache.getVals(iel-1,ip);
-                SplineUtils::extractBasis(spline2[ip],fe.N,dNdu,d2Ndu2);
+              fe.N = bfs.N;
-              else
-                SplineUtils::extractBasis(spline[ip],fe.N,dNdu);
 
               // 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<size_t,3> gpIdx = this->gpIndex(gp,reduced);
+  std::array<size_t,3> 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<size_t,3> elIdx = this->elmIndex(el);
+  std::array<size_t,3> 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<size_t, 3> 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<BasisFunctionVals>& values)
+  {
+    std::vector<Go::BasisDerivs>  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<Go::BasisDerivs2> 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);
+}

src/ASM/ASMs3D.h

@@ -16,9 +16,12 @@
 
 #include "ASMstruct.h"
 #include "ASM3D.h"
+#include "BasisFunctionCache.h"
 #include "Interface.h"
 #include "ThreadGroups.h"
 
+#include <memory>
+
 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<size_t,3>& 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<size_t,3> nel; //!< Number of elements in each direction
+
+private:
+    //! \brief Obtain structured element indices.
+    //! \param el Global element index
+    std::array<size_t,3> 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
@@ -821,6 +886,9 @@ protected:
   ThreadGroups                threadGroupsVol;
   //! Element groups for multi-threaded face assembly
   std::map<char,ThreadGroups> threadGroupsFace;
+
+  //! Basis function cache
+  std::vector<std::unique_ptr<BasisFunctionCache>> myCache;
 };
 
 #endif