added: BasisFunctionCache in ASMs3DmxLag

This commit is contained in:
Arne Morten Kvarving 2022-06-08 13:29:06 +02:00
parent e595bd3556
commit bffcfa640f
2 changed files with 101 additions and 31 deletions

View File

@ -266,20 +266,28 @@ bool ASMs3DmxLag::integrate (Integrand& integrand,
{
if (this->empty()) return true; // silently ignore empty patches
// Get Gaussian quadrature points and weights
const double* x = GaussQuadrature::getCoord(nGauss);
const double* w = GaussQuadrature::getWeight(nGauss);
if (!x || !w) return false;
if (myCache.empty()) {
myCache.emplace_back(std::make_unique<BasisFunctionCache>(*this, cachePolicy, 1));
const BasisFunctionCache& front = static_cast<const BasisFunctionCache&>(*myCache.front());
for (size_t b = 2; b <= this->getNoBasis(); ++b)
myCache.emplace_back(std::make_unique<BasisFunctionCache>(front, b));
}
// Get parametric coordinates of the elements
RealArray upar, vpar, wpar;
this->getGridParameters(upar,0,1);
this->getGridParameters(vpar,1,1);
this->getGridParameters(wpar,2,1);
for (std::unique_ptr<ASMs3D::BasisFunctionCache>& cache : myCache) {
cache->setIntegrand(&integrand);
cache->init(1);
}
BasisFunctionCache& cache = static_cast<BasisFunctionCache&>(*myCache.front());
// Get Gaussian quadrature points and weights
const std::array<int,3>& ng = cache.nGauss();
const std::array<const double*,3>& xg = cache.coord();
const std::array<const double*,3>& wg = cache.weight();
// Number of elements in each direction
const int nelx = upar.size() - 1;
const int nely = vpar.size() - 1;
const int nelx = cache.noElms()[0];
const int nely = cache.noElms()[1];
// === Assembly loop over all elements in the patch ==========================
@ -321,41 +329,40 @@ bool ASMs3DmxLag::integrate (Integrand& integrand,
// --- Integration loop over all Gauss points in each direction --------
int jp = ((i3*nely + i2)*nelx + i1)*nGauss*nGauss*nGauss;
size_t ip = 0;
int jp = ((i3*nely + i2)*nelx + i1)*ng[0]*ng[1]*ng[2];
fe.iGP = firstIp + jp; // Global integration point counter
for (int k = 0; k < nGauss; k++)
for (int j = 0; j < nGauss; j++)
for (int i = 0; i < nGauss; i++, 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++, fe.iGP++)
{
// Parameter value of current integration point
fe.u = 0.5*(upar[i1]*(1.0-x[i]) + upar[i1+1]*(1.0+x[i]));
fe.v = 0.5*(vpar[i2]*(1.0-x[j]) + vpar[i2+1]*(1.0+x[j]));
fe.w = 0.5*(wpar[i3]*(1.0-x[k]) + wpar[i3+1]*(1.0+x[k]));
fe.u = cache.getParam(0,i1,i);
fe.v = cache.getParam(1,i2,j);
fe.w = cache.getParam(2,i3,k);
// Local coordinate of current integration point
fe.xi = x[i];
fe.eta = x[j];
fe.zeta = x[k];
fe.xi = xg[0][i];
fe.eta = xg[1][j];
fe.zeta = xg[2][k];
// 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],x[i],
elem_sizes[b][1],x[j],
elem_sizes[b][2],x[k]))
ok = false;
std::vector<const BasisFunctionVals*> 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 *= w[i]*w[j]*w[k];
fe.detJxW *= wg[0][i]*wg[1][j]*wg[2][k];
if (!integrand.evalIntMx(*A,fe,time,X))
ok = false;
}
@ -596,7 +603,7 @@ bool ASMs3DmxLag::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
@ -618,3 +625,38 @@ bool ASMs3DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
return true;
}
ASMs3DmxLag::BasisFunctionCache::BasisFunctionCache (const ASMs3DLag& pch,
ASM::CachePolicy plcy,
int b) :
ASMs3DLag::BasisFunctionCache(pch,plcy,b)
{
}
ASMs3DmxLag::BasisFunctionCache::BasisFunctionCache (const BasisFunctionCache& cache,
int b) :
ASMs3DLag::BasisFunctionCache(cache,b)
{
}
BasisFunctionVals ASMs3DmxLag::BasisFunctionCache::calculatePt (size_t el,
size_t gp,
bool reduced) const
{
std::array<size_t,3> gpIdx = this->gpIndex(gp,reduced);
const Quadrature& q = reduced ? *reducedQ : *mainQ;
const ASMs3DmxLag& pch = static_cast<const ASMs3DmxLag&>(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]],
pch.elem_sizes[basis-1][2],q.xg[2][gpIdx[2]]);
return result;
}

View File

@ -28,6 +28,32 @@
class ASMs3DmxLag : public ASMs3DLag, private ASMmxBase
{
//! \brief Implementation of basis function cache.
class BasisFunctionCache : public ASMs3DLag::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 ASMs3DLag& 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.
explicit ASMs3DmxLag(const CharVec& n_f);
@ -50,6 +76,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.