ASMu3Dmx::integrate: add support for a separate geometry basis

src/ASM/LR/ASMu3Dmx.C

@@ -303,14 +303,16 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
   }
 
   bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES;
+  const bool separateGeometry = this->getBasis(ASM::GEOMETRY_BASIS) != lrspline.get();
 
   if (myCache.empty()) {
     myCache.emplace_back(std::make_unique<BasisFunctionCache>(*this, cachePolicy, 1,
                                                               ASMmxBase::Type != SUBGRID));
-    for (size_t b = 2; b <= this->getNoBasis(); ++b) {
     const BasisFunctionCache& c = static_cast<const BasisFunctionCache&>(*myCache.front());
+    for (size_t b = 2; b <= this->getNoBasis(); ++b)
       myCache.emplace_back(std::make_unique<BasisFunctionCache>(c,b));
-    }
+    if (separateGeometry)
+      myCache.emplace_back(std::make_unique<BasisFunctionCache>(c,ASM::GEOMETRY_BASIS));
   }
 
   for (std::unique_ptr<ASMu3D::BasisFunctionCache>& cache : myCache) {
@@ -339,11 +341,9 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
       if (!ok)
         continue;
 
-      const LR::Element* el = threadBasis->getElement(groups[t][e]);
-
       std::vector<int>    els;
       std::vector<size_t> elem_sizes;
-      this->getElementsAt(el->midpoint(),els,elem_sizes);
+      this->getElementsAt(threadBasis->getElement(groups[t][e])->midpoint(),els,elem_sizes);
 
       MxFiniteElement fe(elem_sizes);
       Matrix   Xnod, Jac;
@@ -352,9 +352,11 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
       double   param[3] = { 0.0, 0.0, 0.0 };
       Vec4     X(param,time.t);
 
-      int iEl = el->getId();
+      int iEl = els[itgBasis-1]-1;
       fe.iel = MLGE[iEl];
 
+      const LR::Element* el = lrspline->getElement(iEl);
+
       // Get element volume in the parameter space
       double dV = 0.125*el->volume();
 
@@ -442,9 +444,10 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
             fe.v = param[1] = cache.getParam(1,iEl,j);
             fe.w = param[2] = cache.getParam(2,iEl,k);
 
-            std::vector<const BasisFunctionVals*> bfs(this->getNoBasis());
+            std::vector<const BasisFunctionVals*> bfs(myCache.size());
-            for (size_t b = 0; b < m_basis.size(); ++b) {
+            for (size_t b = 0; b < myCache.size(); ++b) {
               bfs[b] = &myCache[b]->getVals(iEl,ig);
+              if (b < m_basis.size())
                 fe.basis(b+1) = bfs[b]->N;
             }
 
@@ -462,7 +465,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
               utl::getGmat(Jac,dXidu,fe.G);
 
             // Cartesian coordinates of current integration point
-            X.assign(Xnod * fe.basis(itgBasis));
+            X.assign(Xnod * (separateGeometry ? bfs.back()->N : fe.basis(itgBasis)));
 
             // Evaluate the integrand and accumulate element contributions
             fe.detJxW *= dV*wg[0][i]*wg[1][j]*wg[2][k];
@@ -496,6 +499,8 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
 
   PROFILE2("ASMu3Dmx::integrate(B)");
 
+  const bool separateGeometry = this->getBasis(ASM::GEOMETRY_BASIS) != lrspline.get();
+
   // Get Gaussian quadrature points and weights
   int nGP = integrand.getBouIntegrationPoints(nGauss);
   const double* xg = GaussQuadrature::getCoord(nGP);
@@ -567,11 +572,12 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
     fe.v = gpar[1](1);
     fe.w = gpar[2](1);
 
-    std::vector<Matrix> dNxdu(m_basis.size());
+    std::vector<Matrix> dNxdu(m_basis.size() + separateGeometry);
     Matrix Xnod, Jac;
 
     double param[3] = { fe.u, fe.v, fe.w };
     Vec4   X(param,time.t);
+    Vector Ng;
     Vec3   normal;
     double dXidu[3];
 
@@ -630,9 +636,13 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
           fe.w = param[2] = gpar[2](k3+1);
         }
 
+        if (separateGeometry)
+          this->evaluateBasis(els.back()-1, fe.u, fe.v, fe.w,
+                              Ng, dNxdu.back(), ASM::GEOMETRY_BASIS);
+
         // Fetch basis function derivatives at current integration point
         for (size_t b = 1; b <= m_basis.size(); ++b)
-          this->evaluateBasis(iEl, fe, dNxdu[b-1], b);
+          this->evaluateBasis(els[b-1]-1, fe, dNxdu[b-1], b);
 
         // Compute basis function derivatives and the face normal
         if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
@@ -645,7 +655,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
           utl::getGmat(Jac,dXidu,fe.G);
 
         // Cartesian coordinates of current integration point
-        X.assign(Xnod * fe.basis(itgBasis));
+        X.assign(Xnod * (separateGeometry ? Ng : fe.basis(itgBasis)));
 
         // Evaluate the integrand and accumulate element contributions
         fe.detJxW *= dA*wg[i]*wg[j];
@@ -1012,6 +1022,9 @@ void ASMu3Dmx::getElementsAt (const RealArray& param,
     elms.push_back(1+iel);
     sizes.push_back(basis->getElement(iel)->nBasisFunctions());
   }
+  const LR::LRSplineVolume* geo = this->getBasis(ASM::GEOMETRY_BASIS);
+  if (geo != lrspline.get())
+    elms.push_back(1+geo->getElementContaining(param));
 }