From 982eb0ce3eec99e908f423ea429b4d4833fb2700 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 30 Nov 2017 10:29:55 +0100 Subject: [PATCH] changed: virtualize transferGaussPtVars / transferCtrlPtVars and implement for 3D --- src/ASM/ASMunstruct.h | 31 +++++++ src/ASM/LR/ASMLRSpline.C | 10 ++ src/ASM/LR/ASMu2D.C | 19 ++-- src/ASM/LR/ASMu2D.h | 32 +++---- src/ASM/LR/ASMu3D.C | 175 +++++++++++++++++++++++++++++++++++ src/ASM/LR/ASMu3D.h | 24 +++++ src/ASM/LR/Test/TestASMu3D.C | 71 ++++++++++++++ 7 files changed, 330 insertions(+), 32 deletions(-) diff --git a/src/ASM/ASMunstruct.h b/src/ASM/ASMunstruct.h index c57b330c..72cb8e00 100644 --- a/src/ASM/ASMunstruct.h +++ b/src/ASM/ASMunstruct.h @@ -172,6 +172,37 @@ public: virtual void remapErrors(RealArray& errors, const RealArray& orig, bool elemErrors = false) const = 0; + //! \brief Transfers Gauss point variables from old basis to this patch. + //! \param[in] oldBasis The LR-spline basis to transfer from + //! \param[in] oldVars Gauss point variables associated with \a oldBasis + //! \param[out] newVars Gauss point variables associated with this patch. + //! \param[in] nGauss Number of Gauss points along a knot-span + virtual bool transferGaussPtVars(const LR::LRSpline* oldBasis, + const RealArray& oldVars, RealArray& newVars, + int nGauss) const = 0; + //! \brief Transfers Gauss point variables from old basis to this patch. + //! \param[in] oldBasis The LR-spline basis to transfer from + //! \param[in] oldVars Gauss point variables associated with \a oldBasis + //! \param[out] newVars Gauss point variables associated with this patch. + //! \param[in] nGauss Number of Gauss points along a knot-span + virtual bool transferGaussPtVarsN(const LR::LRSpline* oldBasis, + const RealArray& oldVars, RealArray& newVars, + int nGauss) const = 0; + //! \brief Transfers control point variables from old basis to this patch. + //! \param[in] oldBasis The LR-spline basis to transfer from + //! \param[out] newVars Gauss point variables associated with this patch. + //! \param[in] nGauss Number of Gauss points along a knot-span + virtual bool transferCntrlPtVars(const LR::LRSpline* oldBasis, + RealArray& newVars, int nGauss) const = 0; + //! \brief Transfers control point variables from old basis to this patch. + //! \param[in] oldBasis The LR-spline basis to transfer from + //! \param[in] oldVars Control point variables associated with \a oldBasis + //! \param[out] newVars Gauss point variables associated with this patch. + //! \param[in] nGauss Number of Gauss points along a knot-span + bool transferCntrlPtVars(LR::LRSpline* oldBasis, + const RealArray& oldVars, RealArray& newVars, + int nGauss, int nf = 1) const; + protected: LR::LRSpline* geo; //!< Pointer to the actual spline geometry object diff --git a/src/ASM/LR/ASMLRSpline.C b/src/ASM/LR/ASMLRSpline.C index 8d44a07d..da46c89a 100644 --- a/src/ASM/LR/ASMLRSpline.C +++ b/src/ASM/LR/ASMLRSpline.C @@ -384,3 +384,13 @@ void ASMunstruct::Sort(int u, int v, int orient, return false; }); } + + +bool ASMunstruct::transferCntrlPtVars (LR::LRSpline* oldBasis, + const RealArray& oldVars, RealArray& newVars, + int nGauss, int nf) const +{ + oldBasis->rebuildDimension(nf); + oldBasis->setControlPoints(const_cast(oldVars)); + return this->transferCntrlPtVars(oldBasis,newVars,nGauss); +} diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index fb5b2ea9..a02934ee 100644 --- a/src/ASM/LR/ASMu2D.C +++ b/src/ASM/LR/ASMu2D.C @@ -1992,11 +1992,12 @@ size_t ASMu2D::getNoNodes (int) const } -bool ASMu2D::transferGaussPtVars (const LR::LRSplineSurface* oldBasis, +bool ASMu2D::transferGaussPtVars (const LR::LRSpline* old_basis, const RealArray& oldVars, RealArray& newVars, int nGauss) const { const LR::LRSplineSurface* newBasis = this->getBasis(); + const LR::LRSplineSurface* oldBasis = static_cast(old_basis); size_t nGp = nGauss*nGauss; newVars.resize(newBasis->nElements()*nGp); @@ -2044,11 +2045,12 @@ bool ASMu2D::transferGaussPtVars (const LR::LRSplineSurface* oldBasis, } -bool ASMu2D::transferGaussPtVarsN (const LR::LRSplineSurface* oldBasis, +bool ASMu2D::transferGaussPtVarsN (const LR::LRSpline* old_basis, const RealArray& oldVars, RealArray& newVars, int nGauss) const { const LR::LRSplineSurface* newBasis = this->getBasis(); + const LR::LRSplineSurface* oldBasis = static_cast(old_basis); size_t nGP = nGauss*nGauss; newVars.clear(); @@ -2110,20 +2112,11 @@ bool ASMu2D::transferGaussPtVarsN (const LR::LRSplineSurface* oldBasis, } -bool ASMu2D::transferCntrlPtVars (LR::LRSplineSurface* oldBasis, - const RealArray& oldVars, RealArray& newVars, - int nGauss) const -{ - oldBasis->rebuildDimension(1); - oldBasis->setControlPoints(const_cast(oldVars)); - return this->transferCntrlPtVars(oldBasis,newVars,nGauss); -} - - -bool ASMu2D::transferCntrlPtVars (const LR::LRSplineSurface* oldBasis, +bool ASMu2D::transferCntrlPtVars (const LR::LRSpline* old_basis, RealArray& newVars, int nGauss) const { const LR::LRSplineSurface* newBasis = this->getBasis(); + const LR::LRSplineSurface* oldBasis = static_cast(old_basis); newVars.clear(); newVars.reserve(newBasis->nElements()*nGauss*nGauss*oldBasis->dimension()); diff --git a/src/ASM/LR/ASMu2D.h b/src/ASM/LR/ASMu2D.h index 6ee780df..28381c07 100644 --- a/src/ASM/LR/ASMu2D.h +++ b/src/ASM/LR/ASMu2D.h @@ -424,35 +424,29 @@ public: double time) const; //! \brief Transfers Gauss point variables from old basis to this patch. - //! \param[in] oldBasis The LR-spline basis to transfer from + //! \param[in] old_basis The LR-spline basis to transfer from //! \param[in] oldVars Gauss point variables associated with \a oldBasis //! \param[out] newVars Gauss point variables associated with this patch. //! \param[in] nGauss Number of Gauss points along a knot-span - bool transferGaussPtVars(const LR::LRSplineSurface* oldBasis, - const RealArray& oldVars, RealArray& newVars, - int nGauss) const; + virtual bool transferGaussPtVars(const LR::LRSpline* old_basis, + const RealArray& oldVars, RealArray& newVars, + int nGauss) const; //! \brief Transfers Gauss point variables from old basis to this patch. - //! \param[in] oldBasis The LR-spline basis to transfer from + //! \param[in] old_basis The LR-spline basis to transfer from //! \param[in] oldVars Gauss point variables associated with \a oldBasis //! \param[out] newVars Gauss point variables associated with this patch. //! \param[in] nGauss Number of Gauss points along a knot-span - bool transferGaussPtVarsN(const LR::LRSplineSurface* oldBasis, - const RealArray& oldVars, RealArray& newVars, - int nGauss) const; + virtual bool transferGaussPtVarsN(const LR::LRSpline* old_basis, + const RealArray& oldVars, RealArray& newVars, + int nGauss) const; + + using ASMunstruct::transferCntrlPtVars; //! \brief Transfers control point variables from old basis to this patch. - //! \param[in] oldBasis The LR-spline basis to transfer from - //! \param[in] oldVars Control point variables associated with \a oldBasis + //! \param[in] old_basis The LR-spline basis to transfer from //! \param[out] newVars Gauss point variables associated with this patch. //! \param[in] nGauss Number of Gauss points along a knot-span - bool transferCntrlPtVars(LR::LRSplineSurface* oldBasis, - const RealArray& oldVars, RealArray& newVars, - int nGauss) const; - //! \brief Transfers control point variables from old basis to this patch. - //! \param[in] oldBasis The LR-spline basis to transfer from - //! \param[out] newVars Gauss point variables associated with this patch. - //! \param[in] nGauss Number of Gauss points along a knot-span - bool transferCntrlPtVars(const LR::LRSplineSurface* oldBasis, - RealArray& newVars, int nGauss) const; + virtual bool transferCntrlPtVars(const LR::LRSpline* old_basis, + RealArray& newVars, int nGauss) const; protected: diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C index e35eab73..79d492a9 100644 --- a/src/ASM/LR/ASMu3D.C +++ b/src/ASM/LR/ASMu3D.C @@ -22,6 +22,7 @@ #include "TimeDomain.h" #include "FiniteElement.h" #include "GlobalIntegral.h" +#include "LagrangeInterpolator.h" #include "LocalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -2234,3 +2235,177 @@ bool ASMu3D::evaluate (const FunctionBase* func, RealArray& vec, vec = ctrlPvals; return ok; } + + +bool ASMu3D::transferGaussPtVars (const LR::LRSpline* old_basis, + const RealArray& oldVars, RealArray& newVars, + int nGauss) const +{ + const LR::LRSplineVolume* newBasis = this->getBasis(); + const LR::LRSplineVolume* oldBasis = static_cast(old_basis); + + size_t nGp = nGauss*nGauss*nGauss; + newVars.resize(newBasis->nElements()*nGp); + + const double* xi = GaussQuadrature::getCoord(nGauss); + LagrangeInterpolator interp(Vector(xi,nGauss)); + + for (int iEl = 0; iEl < newBasis->nElements(); iEl++) + { + const LR::Element* newEl = newBasis->getElement(iEl); + double u_center = 0.5*(newEl->umin() + newEl->umax()); + double v_center = 0.5*(newEl->vmin() + newEl->vmax()); + double w_center = 0.5*(newEl->wmin() + newEl->wmax()); + int iOld = oldBasis->getElementContaining(u_center,v_center,w_center); + if (iOld < 0) + { + std::cerr <<" *** ASMu3D: Failed to locate element "<< iEl + <<" of the old mesh in the new mesh."<< std::endl; + return false; + } + const LR::Element* oldEl = oldBasis->getElement(iOld); + + std::array I; + for (int i = 0; i < 3; i++) + { + RealArray UGP; + LR::getGaussPointParameters(newBasis, UGP, i, nGauss, iEl+1, xi); + double pmin = oldEl->getParmin(i); + double pmax = oldEl->getParmax(i); + for (size_t j = 0; j < UGP.size(); j++) + UGP[j] = -1.0 + 2.0*(UGP[j]-pmin)/(pmax-pmin); + + // lagrangian interpolation + I[i] = interp.get(UGP); + } + + Matrix data(nGauss,nGauss*nGauss); + data.fill(&oldVars[nGp*iOld],nGp); + + Matrix I1 = I[0]*data; + for (int i = 0; i < nGauss; ++i) + for (int j = 0; j < nGauss; ++j) + for (int k = 0; k < nGauss; ++k) + data(j + 1, i + k*(nGauss) + 1) = I1(i+1, 1 + j + k*(nGauss)); + I1 = I[1]*data; + for (int i = 0; i < nGauss; ++i) + for (int j = 0; j < nGauss; ++j) + for (int k = 0; k < nGauss; ++k) + data(k + 1, i + j*(nGauss) + 1) = I1(j+1, 1 + i + k*(nGauss)); + + I1 = I[2]*data; + for (int i = 0; i < nGauss; ++i) + for (int j = 0; j < nGauss; ++j) + for (int k = 0; k < nGauss; ++k) + newVars[iEl*nGp+i+(j + k*nGauss)*nGauss] = I1(k+1, 1+i + j*nGauss); + } + + return true; +} + + +bool ASMu3D::transferGaussPtVarsN (const LR::LRSpline* old_basis, + const RealArray& oldVars, RealArray& newVars, + int nGauss) const +{ + const LR::LRSplineVolume* newBasis = this->getBasis(); + const LR::LRSplineVolume* oldBasis = static_cast(old_basis); + + size_t nGP = nGauss*nGauss*nGauss; + newVars.clear(); + newVars.reserve(nGP*newBasis->nElements()); + + struct Param{ double u,v,w; }; + std::vector oGP(nGP); + + const double* xi = GaussQuadrature::getCoord(nGauss); + for (int iEl = 0; iEl < newBasis->nElements(); iEl++) + { + const LR::Element* newEl = newBasis->getElement(iEl); + double u_center = 0.5*(newEl->umin() + newEl->umax()); + double v_center = 0.5*(newEl->vmin() + newEl->vmax()); + double w_center = 0.5*(newEl->wmin() + newEl->wmax()); + int iOld = oldBasis->getElementContaining(u_center,v_center,w_center); + if (iOld < 0) + { + std::cerr <<" *** ASMu3D: Failed to locate element "<< iEl + <<" of the old mesh in the new mesh."<< std::endl; + return false; + } + const LR::Element* oldEl = oldBasis->getElement(iOld); + + // find parameters of old gauss points + double umin = oldEl->umin(); + double vmin = oldEl->vmin(); + double wmin = oldEl->wmin(); + double du = 0.5*(oldEl->umax() - umin); + double dv = 0.5*(oldEl->vmax() - vmin); + double dw = 0.5*(oldEl->wmax() - wmin); + size_t l = 0; + for (int k = 0; k < nGauss; ++k) + for (int j = 0; j < nGauss; ++j) + for (int i = 0; i < nGauss; ++i, ++l) { + oGP[l].u = umin + du * (xi[i] + 1.0); + oGP[l].v = vmin + dv * (xi[j] + 1.0); + oGP[l].w = wmin + dw * (xi[k] + 1.0); + } + + // parameters of new gauss points + umin = newEl->umin(); + vmin = newEl->vmin(); + wmin = newEl->wmin(); + du = 0.5*(newEl->umax() - umin); + dv = 0.5*(newEl->vmax() - vmin); + dw = 0.5*(newEl->wmax() - wmin); + for (int k = 0; k < nGauss; ++k) + for (int j = 0; j < nGauss; ++j) + for (int i = 0; i < nGauss; ++i) { + double u = umin + du * (xi[i] + 1.0); + double v = vmin + dv * (xi[j] + 1.0); + double w = wmin + dw * (xi[k] + 1.0); + double dist = 1.0e16; + size_t near = 0; + for (size_t l = 0; l < nGP; ++l) { + Vec3 d(oGP[l].u-u, oGP[l].v-v, oGP[l].w-w); + double nd = d.length(); + if (nd < dist) { + near = l; + dist = nd; + } + } + newVars.push_back(oldVars[iOld*nGP+near]); + } + } + + return true; +} + + +bool ASMu3D::transferCntrlPtVars (const LR::LRSpline* old_basis, + RealArray& newVars, int nGauss) const +{ + const LR::LRSplineVolume* newBasis = this->getBasis(); + const LR::LRSplineVolume* oldBasis = static_cast(old_basis); + + newVars.clear(); + newVars.reserve(newBasis->nElements()*nGauss*nGauss*nGauss*oldBasis->dimension()); + const double* xi = GaussQuadrature::getCoord(nGauss); + + for (int iEl = 0; iEl < newBasis->nElements(); iEl++) + { + RealArray U, V, W, ptVar; + LR::getGaussPointParameters(newBasis, U, 0, nGauss, iEl+1, xi); + LR::getGaussPointParameters(newBasis, V, 1, nGauss, iEl+1, xi); + LR::getGaussPointParameters(newBasis, W, 2, nGauss, iEl+1, xi); + for (int k = 0; k < nGauss; k++) + for (int j = 0; j < nGauss; j++) + for (int i = 0; i < nGauss; i++) + { + oldBasis->point(ptVar,U[i],V[j],W[k]); + for (size_t l = 0; l < ptVar.size(); l++) + newVars.push_back(ptVar[l]); + } + } + + return true; +} diff --git a/src/ASM/LR/ASMu3D.h b/src/ASM/LR/ASMu3D.h index b404394c..cc129ec9 100644 --- a/src/ASM/LR/ASMu3D.h +++ b/src/ASM/LR/ASMu3D.h @@ -339,6 +339,30 @@ public: virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand, const int* npe, char project = '\0') const; + //! \brief Transfers Gauss point variables from old basis to this patch. + //! \param[in] old_basis The LR-spline basis to transfer from + //! \param[in] oldVars Gauss point variables associated with \a oldBasis + //! \param[out] newVars Gauss point variables associated with this patch. + //! \param[in] nGauss Number of Gauss points along a knot-span + virtual bool transferGaussPtVars(const LR::LRSpline* old_basis, + const RealArray& oldVars, RealArray& newVars, + int nGauss) const; + //! \brief Transfers Gauss point variables from old basis to this patch. + //! \param[in] old_basis The LR-spline basis to transfer from + //! \param[in] oldVars Gauss point variables associated with \a oldBasis + //! \param[out] newVars Gauss point variables associated with this patch. + //! \param[in] nGauss Number of Gauss points along a knot-span + virtual bool transferGaussPtVarsN(const LR::LRSpline* old_basis, + const RealArray& oldVars, RealArray& newVars, + int nGauss) const; + using ASMunstruct::transferCntrlPtVars; + //! \brief Transfers control point variables from old basis to this patch. + //! \param[in] old_basis The LR-spline basis to transfer from + //! \param[out] newVars Gauss point variables associated with this patch. + //! \param[in] nGauss Number of Gauss points along a knot-span + virtual bool transferCntrlPtVars(const LR::LRSpline* old_basis, + RealArray& newVars, int nGauss) const; + private: //! \brief Struct representing an inhomogeneous Dirichlet boundary condition. struct DirichletFace diff --git a/src/ASM/LR/Test/TestASMu3D.C b/src/ASM/LR/Test/TestASMu3D.C index d87d4782..66144480 100644 --- a/src/ASM/LR/Test/TestASMu3D.C +++ b/src/ASM/LR/Test/TestASMu3D.C @@ -11,7 +11,9 @@ //============================================================================== #include "ASMu3D.h" +#include "GaussQuadrature.h" #include "SIM3D.h" +#include "LRSpline/LRSplineVolume.h" #include "gtest/gtest.h" @@ -59,3 +61,72 @@ const std::vector tests = {1,2,3,4,5,6}; INSTANTIATE_TEST_CASE_P(TestASMu3D, TestASMu3D, testing::ValuesIn(tests)); + + +TEST(TestASMu3D, TransferGaussPtVars) +{ + SIM3D sim(1); + sim.opt.discretization = ASM::LRSpline; + sim.createDefaultModel(); + ASMu3D* pch = static_cast(sim.getPatch(1)); + size_t id[3]; + const double* xi = GaussQuadrature::getCoord(3); + RealArray oldAr(3*3*3), newAr; + for (size_t idx = 0; idx < 3; ++idx) { + SIM3D simNew(1); + simNew.opt.discretization = ASM::LRSpline; + simNew.createDefaultModel(); + ASMu3D* pchNew = static_cast(simNew.getPatch(1)); + pchNew->uniformRefine(idx, 1); + + for (id[2] = 0; id[2] < 3; ++id[2]) + for (id[1] = 0; id[1] < 3; ++id[1]) + for (id[0] = 0; id[0] < 3; ++id[0]) + oldAr[id[0]+(id[1]+id[2]*3)*3] = (1.0 + xi[id[idx]]) / 2.0; + + pchNew->transferGaussPtVars(pch->getVolume(), oldAr, newAr, 3); + size_t k = 0; + for (size_t iEl = 0; iEl < 2; ++iEl) + for (id[2] = 0; id[2] < 3; ++id[2]) + for (id[1] = 0; id[1] < 3; ++id[1]) + for (id[0] = 0; id[0] < 3; ++id[0], ++k) + EXPECT_FLOAT_EQ(newAr[k], 0.5*iEl + 0.5*(xi[id[idx]] + 1.0) / 2.0); + } +} + + +TEST(TestASMu3D, TransferGaussPtVarsN) +{ + SIM3D sim(1), sim2(1); + sim.opt.discretization = sim2.opt.discretization = ASM::LRSpline; + sim.createDefaultModel(); + sim2.createDefaultModel(); + ASMu3D* pch = static_cast(sim.getPatch(1)); + ASMu3D* pchNew = static_cast(sim2.getPatch(1)); + pchNew->uniformRefine(0, 1); + RealArray oldAr(3*3*3), newAr; + std::iota(oldAr.begin(), oldAr.end(), 1); + + pchNew->transferGaussPtVarsN(pch->getVolume(), oldAr, newAr, 3); + static RealArray refAr = {{ 1.0, 1.0, 2.0, + 4.0, 4.0, 5.0, + 7.0, 7.0, 8.0, + 10.0, 10.0, 11.0, + 13.0, 13.0, 14.0, + 16.0, 16.0, 17.0, + 19.0, 19.0, 20.0, + 22.0, 22.0, 23.0, + 25.0, 25.0, 26.0, + 2.0, 3.0, 3.0, + 5.0, 6.0, 6.0, + 8.0, 9.0, 9.0, + 11.0, 12.0, 12.0, + 14.0, 15.0, 15.0, + 17.0, 18.0, 18.0, + 20.0, 21.0, 21.0, + 23.0, 24.0, 24.0, + 26.0, 27.0, 27.0}}; + EXPECT_EQ(refAr.size(), newAr.size()); + for (size_t i = 0; i < refAr.size(); ++i) + EXPECT_FLOAT_EQ(refAr[i], newAr[i]); +}