From 6827453eb9bdf9b7101e044a07adff1c6537f9a0 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 19 Mar 2021 09:07:50 +0100 Subject: [PATCH] add evaluate for ASMsxDLag also expand evalSolution to respect evaluation points --- src/ASM/ASMs2DLag.C | 94 ++++++++++++++++++++++++++++++++++-------- src/ASM/ASMs2DLag.h | 8 ++++ src/ASM/ASMs3DLag.C | 99 ++++++++++++++++++++++++++++++++++++--------- src/ASM/ASMs3DLag.h | 8 ++++ 4 files changed, 172 insertions(+), 37 deletions(-) diff --git a/src/ASM/ASMs2DLag.C b/src/ASM/ASMs2DLag.C index 4f46f9ce..c7c9ddf9 100644 --- a/src/ASM/ASMs2DLag.C +++ b/src/ASM/ASMs2DLag.C @@ -682,18 +682,57 @@ bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol, bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol, - const RealArray*, bool, int, int) const + const RealArray* gpar, bool, int, int) const { - size_t nPoints = coord.size(); - size_t nNodes = this->getNoNodes(-1); - size_t nComp = locSol.size() / nNodes; - if (nNodes < nPoints || nComp*nNodes != locSol.size()) - return false; + if (!gpar) { + size_t nPoints = coord.size(); + size_t nNodes = this->getNoNodes(-1); + size_t nComp = locSol.size() / nNodes; + if (nNodes < nPoints || nComp*nNodes != locSol.size()) + return false; - sField.resize(nComp,nPoints); - const double* u = locSol.ptr(); - for (size_t n = 1; n <= nPoints; n++, u += nComp) - sField.fillColumn(n,u); + sField.resize(nComp,nPoints); + const double* u = locSol.ptr(); + for (size_t n = 1; n <= nPoints; n++, u += nComp) + sField.fillColumn(n,u); + + return true; + } + + size_t nNodes = gpar[0].size()*gpar[1].size(); + size_t nComp = locSol.size() / this->getNoNodes(); + + sField.resize(nComp, nNodes); + + RealArray N(p1*p2); + double xi, eta; + + size_t n = 1; + for (size_t j = 0; j < gpar[1].size(); ++j) + for (size_t i = 0; i < gpar[0].size(); ++i, ++n) { + int iel = this->findElement(gpar[0][i], gpar[1][j], &xi, &eta); + + if (iel < 1 || iel > int(MNPC.size())) + return false; + + if (!Lagrange::computeBasis(N,p1,xi,p2,eta)) + return false; + + Matrix elmSol(nComp, p1*p2); + const IntVec& mnpc = MNPC[iel-1]; + + size_t idx = 1; + for (const int& m : mnpc) { + for (size_t c = 1; c <= nComp; ++c) + elmSol(c,idx) = locSol(m*nComp+c); + ++idx; + } + + Vector val; + elmSol.multiply(N, val); + for (size_t c = 1; c <= nComp; ++c) + sField(c,n) = val(c); + } return true; } @@ -778,16 +817,37 @@ bool ASMs2DLag::write(std::ostream& os, int) const int ASMs2DLag::findElement(double u, double v, double* xi, double* eta) const { - double du = 1.0 / (nx-1); - double dv = 1.0 / (ny-1); - - int elmx = std::min(nx-2.0, floor(u / du)); - int elmy = std::min(ny-2.0, floor(v / dv)); + int elmx = std::min(nx-2.0, floor(u*(nx-1))); + int elmy = std::min(ny-2.0, floor(v*(ny-1))); if (xi) - *xi = -1.0 + (u - elmx*du)*2.0 / du; + *xi = -1.0 + (u*(nx-1) - elmx)*2.0; if (eta) - *eta = -1.0 + (v - elmy*dv)*2.0 / dv; + *eta = -1.0 + (v*(ny-1) - elmy)*2.0; return 1 + elmx + elmy*(nx-1); } + + +bool ASMs2DLag::evaluate (const ASMbase* basis, const Vector& locVec, + RealArray& vec, int basisNum) const +{ + // Compute parameter values of the result sampling points + std::array gpar; + for (int dir = 0; dir < 2; dir++) { + int nel1 = dir == 0 ? (nx-1)/(p1-1) : (ny-1)/(p2-1); + gpar[dir].resize(nel1+1); + double du = 1.0 / nel1; + for (int i = 0; i <= nel1; ++i) + gpar[dir][i] = i*du; + } + + // Evaluate the result field at all sampling points. + // Note: it is here assumed that *basis and *this have spline bases + // defined over the same parameter domain. + Matrix sValues; + if (!basis->evalSolution(sValues,locVec,gpar.data())) + return false; + vec = sValues; + return true; +} diff --git a/src/ASM/ASMs2DLag.h b/src/ASM/ASMs2DLag.h index d111bf35..8d1b813a 100644 --- a/src/ASM/ASMs2DLag.h +++ b/src/ASM/ASMs2DLag.h @@ -158,6 +158,14 @@ public: virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand, const RealArray*, bool = false) const; + //! \brief Evaluates and interpolates a field over a given geometry. + //! \param[in] basis The basis of the field to evaluate + //! \param[in] locVec The coefficients of the field to evaluate + //! \param[out] vec The obtained coefficients after interpolation + //! \param[in] basisNum The basis to evaluate for (mixed) + virtual bool evaluate (const ASMbase* basis, const Vector& locVec, + RealArray& vec, int basisNum) const; + using ASMs2D::getSize; //! \brief Returns the number of nodal points in each parameter direction. //! \param[out] n1 Number of nodes in first (u) direction diff --git a/src/ASM/ASMs3DLag.C b/src/ASM/ASMs3DLag.C index 7087e0c4..1b47e80d 100644 --- a/src/ASM/ASMs3DLag.C +++ b/src/ASM/ASMs3DLag.C @@ -957,38 +957,73 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const Vector& locSol, int ASMs3DLag::findElement(double u, double v, double w, double* xi, double* eta, double* zeta) const { - double du = 1.0 / (nx-1); - double dv = 1.0 / (ny-1); - double dw = 1.0 / (nz-1); - - int elmx = std::min(nx-2.0, floor(u / du)); - int elmy = std::min(ny-2.0, floor(v / dv)); - int elmz = std::min(nz-2.0, floor(w / dw)); + int elmx = std::min(nx-2.0, floor(u*(nx-1))); + int elmy = std::min(ny-2.0, floor(v*(ny-1))); + int elmz = std::min(nz-2.0, floor(w*(nz-1))); if (xi) - *xi = -1.0 + (u - elmx*du)*2.0 / du; + *xi = -1.0 + (u*(nx-1) - elmx)*2.0; if (eta) - *eta = -1.0 + (v - elmy*dv)*2.0 / dv; + *eta = -1.0 + (v*(ny-1) - elmy)*2.0; if (zeta) - *zeta = -1.0 + (w - elmz*dw)*2.0 / dw; + *zeta = -1.0 + (w*(nz-1) - elmz)*2.0; return 1 + elmx + elmy*(nx-1) + elmz*(ny-1)*(nx-1); } bool ASMs3DLag::evalSolution (Matrix& sField, const Vector& locSol, - const RealArray*, bool, int, int) const + const RealArray* gpar, bool, int, int) const { - size_t nPoints = coord.size(); - size_t nNodes = this->getNoNodes(-1); - size_t nComp = locSol.size() / nNodes; - if (nNodes < nPoints || nComp*nNodes != locSol.size()) - return false; + if (!gpar) { + size_t nPoints = coord.size(); + size_t nNodes = this->getNoNodes(-1); + size_t nComp = locSol.size() / nNodes; + if (nNodes < nPoints || nComp*nNodes != locSol.size()) + return false; + sField.resize(nComp,nPoints); + const double* u = locSol.ptr(); + for (size_t n = 1; n <= nPoints; n++, u += nComp) + sField.fillColumn(n,u); + return true; + } - sField.resize(nComp,nPoints); - const double* u = locSol.ptr(); - for (size_t n = 1; n <= nPoints; n++, u += nComp) - sField.fillColumn(n,u); + size_t nNodes = gpar[0].size()*gpar[1].size()*gpar[2].size(); + size_t nComp = locSol.size() / this->getNoNodes(); + + sField.resize(nComp, nNodes); + + RealArray N(p1*p2*p3); + double xi, eta, zeta; + + size_t n = 1; + for (size_t k = 0; k < gpar[2].size(); ++k) + for (size_t j = 0; j < gpar[1].size(); ++j) + for (size_t i = 0; i < gpar[0].size(); ++i, ++n) { + int iel = this->findElement(gpar[0][i], gpar[1][j], gpar[2][k], + &xi, &eta, &zeta); + + if (iel < 1 || iel > int(MNPC.size())) + return false; + + if (!Lagrange::computeBasis(N,p1,xi,p2,eta,p3,zeta)) + return false; + + Matrix elmSol(nComp, p1*p2*p3); + const IntVec& mnpc = MNPC[iel-1]; + + size_t idx = 1; + for (const int& m : mnpc) { + for (size_t c = 1; c <= nComp; ++c) + elmSol(c,idx) = locSol(m*nComp+c); + ++idx; + } + + Vector val; + elmSol.multiply(N, val); + for (size_t c = 1; c <= nComp; ++c) + sField(c,n) = val(c); + } return true; } @@ -1152,3 +1187,27 @@ bool ASMs3DLag::write (std::ostream& os, int) const { return this->writeLagBasis(os,"hexahedron"); } + + +bool ASMs3DLag::evaluate (const ASMbase* basis, const Vector& locVec, + RealArray& vec, int basisNum) const +{ + // Compute parameter values of the result sampling points + std::array gpar; + for (int dir = 0; dir < 3; dir++) { + int nel1 = dir == 0 ? (nx-1)/(p1-1) : (dir == 1 ? (ny-1)/(p2-1) : (nz-1)/(p3-1)); + gpar[dir].resize(nel1+1); + double du = 1.0 / nel1; + for (int i = 0; i <= nel1; ++i) + gpar[dir][i] = i*du; + } + + // Evaluate the result field at all sampling points. + // Note: it is here assumed that *basis and *this have spline bases + // defined over the same parameter domain. + Matrix sValues; + if (!basis->evalSolution(sValues,locVec,gpar.data())) + return false; + vec = sValues; + return true; +} diff --git a/src/ASM/ASMs3DLag.h b/src/ASM/ASMs3DLag.h index a079efe2..0bbfd208 100644 --- a/src/ASM/ASMs3DLag.h +++ b/src/ASM/ASMs3DLag.h @@ -172,6 +172,14 @@ public: virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand, const RealArray*, bool = false) const; + //! \brief Evaluates and interpolates a field over a given geometry. + //! \param[in] basis The basis of the field to evaluate + //! \param[in] locVec The coefficients of the field to evaluate + //! \param[out] vec The obtained coefficients after interpolation + //! \param[in] basisNum The basis to evaluate for (mixed) + virtual bool evaluate (const ASMbase* basis, const Vector& locVec, + RealArray& vec, int basisNum) const; + //! \brief Returns the number of nodal points in each parameter direction. //! \param[out] n1 Number of nodes in first (u) direction //! \param[out] n2 Number of nodes in second (v) direction