From c97035a74946ad30e186712c479d4b07cc6803d5 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 31 Aug 2023 09:43:26 +0200 Subject: [PATCH] ASMu3Dmx::integrate: add support for a separate geometry basis --- src/ASM/LR/ASMu3Dmx.C | 41 +++++++++++++++++++++++++++-------------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/src/ASM/LR/ASMu3Dmx.C b/src/ASM/LR/ASMu3Dmx.C index c42dac7c..b670829c 100644 --- a/src/ASM/LR/ASMu3Dmx.C +++ b/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(*this, cachePolicy, 1, ASMmxBase::Type != SUBGRID)); - for (size_t b = 2; b <= this->getNoBasis(); ++b) { - const BasisFunctionCache& c = static_cast(*myCache.front()); + const BasisFunctionCache& c = static_cast(*myCache.front()); + for (size_t b = 2; b <= this->getNoBasis(); ++b) myCache.emplace_back(std::make_unique(c,b)); - } + if (separateGeometry) + myCache.emplace_back(std::make_unique(c,ASM::GEOMETRY_BASIS)); } for (std::unique_ptr& 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 els; std::vector 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,10 +444,11 @@ 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 bfs(this->getNoBasis()); - for (size_t b = 0; b < m_basis.size(); ++b) { + std::vector bfs(myCache.size()); + for (size_t b = 0; b < myCache.size(); ++b) { bfs[b] = &myCache[b]->getVals(iEl,ig); - fe.basis(b+1) = bfs[b]->N; + if (b < m_basis.size()) + fe.basis(b+1) = bfs[b]->N; } // Compute Jacobian inverse of the coordinate mapping and @@ -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 dNxdu(m_basis.size()); + std::vector 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)); }