diff --git a/src/ASM/LR/ASMLRSpline.C b/src/ASM/LR/ASMLRSpline.C index ef1cb64c..077158cd 100644 --- a/src/ASM/LR/ASMLRSpline.C +++ b/src/ASM/LR/ASMLRSpline.C @@ -187,20 +187,20 @@ ASMLRSpline::ASMLRSpline (unsigned char n_p, unsigned char n_s, unsigned char n_f) : ASMbase(n_p,n_s,n_f) { - geo = nullptr; + geomB = nullptr; } ASMLRSpline::ASMLRSpline (const ASMLRSpline& patch, unsigned char n_f) : ASMbase(patch,n_f) { - geo = patch.geo; + geomB = patch.geomB; } bool ASMLRSpline::refine (const LR::RefineData& prm, Vectors& sol) { - if (!geo) + if (!geomB) return false; if (shareFE && !prm.refShare) @@ -208,8 +208,8 @@ bool ASMLRSpline::refine (const LR::RefineData& prm, Vectors& sol) // This patch shares spline object with another patch // (in another simulator on the same mesh), // and is assumed to have been refined already - nnod = geo->nBasisFunctions(); - nel = geo->nElements(); + nnod = geomB->nBasisFunctions(); + nel = geomB->nElements(); return true; } @@ -220,32 +220,32 @@ bool ASMLRSpline::refine (const LR::RefineData& prm, Vectors& sol) IntVec nf(sol.size()); for (size_t j = 0; j < sol.size(); j++) - if ((nf[j] = LR::extendControlPoints(geo,sol[j],this->getNoFields(1))) < 0) + if ((nf[j] = LR::extendControlPoints(geomB.get(),sol[j],this->getNoFields(1))) < 0) return false; - if (!this->doRefine(prm,geo)) + if (!this->doRefine(prm,geomB.get())) return false; - nnod = geo->nBasisFunctions(); - nel = geo->nElements(); + nnod = geomB->nBasisFunctions(); + nel = geomB->nElements(); IFEM::cout <<"Refined mesh: "<< nel <<" elements "<< nnod <<" nodes."<< std::endl; for (int i = sol.size()-1; i >= 0; i--) if (nf[i] > 0) { sol[i].resize(nf[i]*nnod); - LR::contractControlPoints(geo,sol[i],nf[i]); + LR::contractControlPoints(geomB.get(),sol[i],nf[i]); } bool linIndepTest = prm.options.size() > 3 ? prm.options[3] != 0 : false; if (linIndepTest) { std::cout <<"Testing for linear independence by overloading"<< std::endl; - bool isLinIndep = geo->isLinearIndepByOverloading(false); + bool isLinIndep = geomB->isLinearIndepByOverloading(false); if (!isLinIndep) { std::cout <<"Inconclusive..."<< std::endl; #ifdef HAS_BOOST std::cout <<"Testing for linear independence by full tensor expansion"<< std::endl; - isLinIndep = geo->isLinearIndepByMappingMatrix(false); + isLinIndep = geomB->isLinearIndepByMappingMatrix(false); #endif } if (isLinIndep) @@ -337,12 +337,12 @@ void ASMLRSpline::getFunctionsForElements (IntSet& functions, const IntVec& elements, bool globalId) const { - geo->generateIDs(); + geomB->generateIDs(); for (int elmId : elements) { int iel = globalId ? utl::findIndex(MLGE,1+elmId) : elmId; - if (iel >= 0 && iel < geo->nElements()) - for (LR::Basisfunction* b : geo->getElement(iel)->support()) + if (iel >= 0 && iel < geomB->nElements()) + for (LR::Basisfunction* b : geomB->getElement(iel)->support()) functions.insert(globalId ? this->getNodeID(b->getId()+1)-1:b->getId()); } } @@ -358,7 +358,7 @@ IntVec ASMLRSpline::getBoundaryCovered (const IntSet& nodes) const this->getBoundaryNodes(edge,oneBoundary,1,1,0,true); for (int i : nodes) for (int j : oneBoundary) - if (geo->getBasisfunction(i)->contains(*geo->getBasisfunction(j-1))) + if (geomB->getBasisfunction(i)->contains(*geomB->getBasisfunction(j-1))) result.insert(j-1); } @@ -371,7 +371,7 @@ IntVec ASMLRSpline::getOverlappingNodes (const IntSet& nodes, int dir) const IntSet result; for (int i : nodes) { - LR::Basisfunction* b = geo->getBasisfunction(i); + LR::Basisfunction* b = geomB->getBasisfunction(i); for (LR::Element* el : b->support()) for (LR::Basisfunction* basis : el->support()) { @@ -432,7 +432,7 @@ std::pair ASMLRSpline::findClosestNode (const Vec3& X) const { double distance = 0.0; size_t inod = 0, iclose = 0; - for (LR::Basisfunction* b : geo->getAllBasisfunctions()) + for (LR::Basisfunction* b : geomB->getAllBasisfunctions()) { double d = (X-Vec3(&(*b->cp()),nsd)).length(); if (++inod == 1 || d < distance) @@ -449,16 +449,16 @@ std::pair ASMLRSpline::findClosestNode (const Vec3& X) const Vec3 ASMLRSpline::getElementCenter (int iel) const { #ifdef INDEX_CHECK - if (iel < 1 || iel > geo->nElements()) + if (iel < 1 || iel > geomB->nElements()) { std::cerr <<" *** ASMLRSpline::getElementCenter: Element index "<< iel - <<" out of range [1,"<< geo->nElements() <<"]."<< std::endl; + <<" out of range [1,"<< geomB->nElements() <<"]."<< std::endl; return Vec3(); } #endif double u0[3] = { 0.0, 0.0, 0.0 }; - LR::Element* elm = geo->getElement(iel-1); + LR::Element* elm = geomB->getElement(iel-1); for (unsigned char d = 0; d < ndim; d++) u0[d] = 0.5*(elm->getParmin(d) + elm->getParmax(d)); @@ -524,23 +524,23 @@ void ASMLRSpline::analyzeThreadGroups (const IntMat& groups) bool ASMLRSpline::getParameterDomain (Real2DMat& u, IntVec* corners) const { - u.resize(geo->nVariate(),RealArray(2)); - for (int i = 0; i < geo->nVariate(); ++i) { - u[i][0] = geo->startparam(i); - u[i][1] = geo->endparam(i); + u.resize(geomB->nVariate(),RealArray(2)); + for (int i = 0; i < geomB->nVariate(); ++i) { + u[i][0] = geomB->startparam(i); + u[i][1] = geomB->endparam(i); } if (corners) { - if (geo->nVariate() == 2) { + if (geomB->nVariate() == 2) { for (int J = 0; J < 2; ++J) for (int I = 0; I < 2; ++I) { std::vector funcs; int dir = (I > 0 ? LR::EAST : LR::WEST) | (J > 0 ? LR::NORTH : LR::SOUTH); - geo->getEdgeFunctions(funcs, static_cast(dir)); + geomB->getEdgeFunctions(funcs, static_cast(dir)); corners->push_back(funcs.front()->getId()+1); } - } else if (geo->nVariate() == 3) { + } else if (geomB->nVariate() == 3) { for (int K = 0; K < 2; ++K) for (int J = 0; J < 2; ++J) for (int I = 0; I < 2; ++I) { @@ -548,7 +548,7 @@ bool ASMLRSpline::getParameterDomain (Real2DMat& u, IntVec* corners) const int dir = (I > 0 ? LR::EAST : LR::WEST) | (J > 0 ? LR::NORTH : LR::SOUTH) | (K > 0 ? LR::TOP : LR::BOTTOM); - geo->getEdgeFunctions(funcs, static_cast(dir)); + geomB->getEdgeFunctions(funcs, static_cast(dir)); corners->push_back(funcs.front()->getId()+1); } } diff --git a/src/ASM/LR/ASMLRSpline.h b/src/ASM/LR/ASMLRSpline.h index 257ce3b0..90372a71 100644 --- a/src/ASM/LR/ASMLRSpline.h +++ b/src/ASM/LR/ASMLRSpline.h @@ -85,7 +85,7 @@ public: virtual ~ASMLRSpline() {} //! \brief Checks if the patch is empty. - virtual bool empty() const { return geo == nullptr; } + virtual bool empty() const { return geomB == nullptr; } //! \brief Returns parameter values and node numbers of the domain corners. virtual bool getParameterDomain(Real2DMat&, IntVec*) const; @@ -203,7 +203,7 @@ protected: //! \param[in] groups The generated thread groups static void analyzeThreadGroups(const IntMat& groups); - LR::LRSpline* geo; //!< Pointer to the actual spline geometry object + std::shared_ptr geomB; //!< Pointer to the actual spline geometry object }; #endif diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index a813076b..1c3d3e62 100644 --- a/src/ASM/LR/ASMu2D.C +++ b/src/ASM/LR/ASMu2D.C @@ -116,7 +116,7 @@ bool ASMu2D::read (std::istream& is) nsd = readDim; } - geo = lrspline.get(); + geomB = lrspline; return true; } @@ -140,7 +140,7 @@ void ASMu2D::clear (bool retainGeometry) delete tensorspline; delete tensorPrjBas; } - geo = nullptr; + geomB = nullptr; tensorspline = tensorPrjBas = nullptr; } @@ -402,7 +402,7 @@ bool ASMu2D::createProjectionBasis (bool init) std::swap(tensorspline,tensorPrjBas); std::swap(lrspline,projBasis); - geo = lrspline.get(); + geomB = lrspline; return true; } @@ -507,7 +507,7 @@ LR::LRSplineSurface* ASMu2D::createLRNurbs (const Go::SplineSurface& srf) } -LR::LRSplineSurface* ASMu2D::createLRfromTensor () +std::shared_ptr ASMu2D::createLRfromTensor () { if (tensorspline) { @@ -541,13 +541,13 @@ LR::LRSplineSurface* ASMu2D::createLRfromTensor () tensorspline = nullptr; } - return lrspline.get(); + return lrspline; } bool ASMu2D::generateFEMTopology () { - geo = this->createLRfromTensor(); + geomB = this->createLRfromTensor(); if (tensorPrjBas) { @@ -1687,7 +1687,7 @@ bool ASMu2D::integrate (Integrand& integrand, const TimeDomain& time, const ASM::InterfaceChecker& iChkgen) { - if (!geo) return true; // silently ignore empty patches + if (!geomB) return true; // silently ignore empty patches if (!(integrand.getIntegrandType() & Integrand::INTERFACE_TERMS)) return true; PROFILE2("ASMu2D::integrate(J)"); @@ -2322,8 +2322,8 @@ void ASMu2D::getBoundaryNodes (int lIndex, IntVec& nodes, int basis, bool ASMu2D::getOrder (int& p1, int& p2, int& p3) const { - p1 = geo->order(0); - p2 = geo->order(1); + p1 = geomB->order(0); + p2 = geomB->order(1); p3 = 0; return true; @@ -2700,13 +2700,13 @@ ASMu2D::InterfaceChecker::InterfaceChecker (const ASMu2D& pch) : myPatch(pch) short int ASMu2D::InterfaceChecker::hasContribution (int e, int, int, int) const { - const LR::Element* elm = myPatch.geo->getElement(e-1); + const LR::Element* elm = myPatch.geomB->getElement(e-1); bool neighbor[4]; - neighbor[0] = elm->getParmin(0) != myPatch.geo->startparam(0); // West - neighbor[1] = elm->getParmax(0) != myPatch.geo->endparam(0); // East - neighbor[2] = elm->getParmin(1) != myPatch.geo->startparam(1); // South - neighbor[3] = elm->getParmax(1) != myPatch.geo->endparam(1); // North + neighbor[0] = elm->getParmin(0) != myPatch.geomB->startparam(0); // West + neighbor[1] = elm->getParmax(0) != myPatch.geomB->endparam(0); // East + neighbor[2] = elm->getParmin(1) != myPatch.geomB->startparam(1); // South + neighbor[3] = elm->getParmax(1) != myPatch.geomB->endparam(1); // North // Check for existing neighbors short int status = 0, s = 1; @@ -2763,8 +2763,8 @@ bool ASMu2D::refine (const LR::RefineData& prm, Vectors& sol) void ASMu2D::generateBezierBasis () { - bezier_u = this->getBezierBasis(geo->order(0)); - bezier_v = this->getBezierBasis(geo->order(1)); + bezier_u = this->getBezierBasis(geomB->order(0)); + bezier_v = this->getBezierBasis(geomB->order(1)); } @@ -2772,16 +2772,16 @@ void ASMu2D::generateBezierExtraction () { PROFILE2("Bezier extraction"); - const int p1 = geo->order(0); - const int p2 = geo->order(1); + const int p1 = geomB->order(0); + const int p2 = geomB->order(1); - myBezierExtract.resize(geo->nElements()); + myBezierExtract.resize(geomB->nElements()); RealArray extrMat; int iel = 0; - for (const LR::Element* elm : geo->getAllElements()) + for (const LR::Element* elm : geomB->getAllElements()) { // Get bezier extraction matrix - geo->getBezierExtraction(iel,extrMat); + geomB->getBezierExtraction(iel,extrMat); myBezierExtract[iel].resize(elm->nBasisFunctions(),p1*p2); myBezierExtract[iel++].fill(extrMat.data(),extrMat.size()); } diff --git a/src/ASM/LR/ASMu2D.h b/src/ASM/LR/ASMu2D.h index 01cbd4e4..9b6ab9db 100644 --- a/src/ASM/LR/ASMu2D.h +++ b/src/ASM/LR/ASMu2D.h @@ -131,10 +131,10 @@ public: //! \brief Copy constructor. ASMu2D(const ASMu2D& patch, unsigned char n_f = 0); //! \brief Empty destructor. - virtual ~ASMu2D() { geo = nullptr; } + virtual ~ASMu2D() { geomB = nullptr; } //! \brief Returns the spline surface representing the geometry of this patch. - LR::LRSplineSurface* getSurface() { return this->createLRfromTensor(); } + LR::LRSplineSurface* getSurface() { return this->createLRfromTensor().get(); } //! \brief Returns the spline surface representing the geometry of this patch. const LR::LRSplineSurface* getSurface() const { return lrspline.get(); } @@ -687,7 +687,7 @@ protected: const IntSet& neighborIndices) const; //! \brief Converts current tensor spline object to LR-spline. - LR::LRSplineSurface* createLRfromTensor(); + std::shared_ptr createLRfromTensor(); //! \brief Converts a rational spline surface to a LR NURBS surface. static LR::LRSplineSurface* createLRNurbs(const Go::SplineSurface& srf); diff --git a/src/ASM/LR/ASMu2DC1.C b/src/ASM/LR/ASMu2DC1.C index 919fd5fd..e85b4715 100644 --- a/src/ASM/LR/ASMu2DC1.C +++ b/src/ASM/LR/ASMu2DC1.C @@ -18,13 +18,14 @@ bool ASMu2DC1::generateFEMTopology () { - if (!(geo = this->createLRfromTensor())) + geomB = this->createLRfromTensor(); + if (!geomB) return false; - else if (geo->order(0) > 2 || geo->order(1) > 2) + else if (geomB->order(0) > 2 || geomB->order(1) > 2) return this->ASMu2D::generateFEMTopology(); std::cerr <<" *** ASMu2DC1::generateFEMTopology:" - <<" The polynomial order "<< geo->order(0) <<"x"<< geo->order(1) + <<" The polynomial order "<< geomB->order(0) <<"x"<< geomB->order(1) <<" is too low.\n C1-continuity requires" <<" at least quadratic order."<< std::endl; return false; diff --git a/src/ASM/LR/ASMu2Dmx.C b/src/ASM/LR/ASMu2Dmx.C index b29b037e..c8da7e6d 100644 --- a/src/ASM/LR/ASMu2Dmx.C +++ b/src/ASM/LR/ASMu2Dmx.C @@ -302,7 +302,7 @@ bool ASMu2Dmx::generateFEMTopology () std::cout <<"NEL = "<< nel <<" NNOD = "<< nnod << std::endl; #endif - geo = m_basis[geoBasis-1].get(); + geomB = m_basis[geoBasis-1]; this->generateBezierBasis(); this->generateBezierExtraction(); @@ -397,8 +397,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand, if (integrand.getIntegrandType() & Integrand::G_MATRIX) { // Element size in parametric space - dXidu[0] = geo->getElement(geoEl-1)->umax()-geo->getElement(geoEl-1)->umin(); - dXidu[1] = geo->getElement(geoEl-1)->vmax()-geo->getElement(geoEl-1)->vmin(); + dXidu[0] = geomB->getElement(geoEl-1)->umax() - geomB->getElement(geoEl-1)->umin(); + dXidu[1] = geomB->getElement(geoEl-1)->vmax() - geomB->getElement(geoEl-1)->vmin(); } // Initialize element quantities @@ -1220,7 +1220,7 @@ void ASMu2Dmx::swapProjectionBasis () std::swap(projBasis, altProjBasis); std::swap(projThreadGroups, altProjThreadGroups); lrspline = m_basis[ASMmxBase::geoBasis-1]; - geo = lrspline.get(); + geomB = lrspline; this->generateBezierBasis(); this->generateBezierExtraction(); } diff --git a/src/ASM/LR/ASMu2Dmx.h b/src/ASM/LR/ASMu2Dmx.h index baa4af2b..672f28c2 100644 --- a/src/ASM/LR/ASMu2Dmx.h +++ b/src/ASM/LR/ASMu2Dmx.h @@ -63,7 +63,7 @@ public: //! \brief Copy constructor. ASMu2Dmx(const ASMu2Dmx& patch, const CharVec& n_f = CharVec(2,0)); //! \brief Empty destructor. - virtual ~ASMu2Dmx() { lrspline = nullptr; geo = nullptr; } + virtual ~ASMu2Dmx() {} //! \brief Returns the spline surface representing the basis of this patch. virtual LR::LRSplineSurface* getBasis(int basis = 1); diff --git a/src/ASM/LR/ASMu2Dnurbs.C b/src/ASM/LR/ASMu2Dnurbs.C index b1fc7d54..13b86ca8 100644 --- a/src/ASM/LR/ASMu2Dnurbs.C +++ b/src/ASM/LR/ASMu2Dnurbs.C @@ -298,7 +298,7 @@ void ASMu2D::writePostscriptElementsNurbs (std::shared_ptr if (mesh != lrspline) { mesh.swap(lrspline); - geo = lrspline.get(); + geomB = lrspline; } // This is done unconditionally as it will not have been performed @@ -383,7 +383,7 @@ void ASMu2D::writePostscriptElementsNurbs (std::shared_ptr if (mesh != lrspline) { mesh.swap(lrspline); - geo = lrspline.get(); + geomB = lrspline; this->generateBezierBasis(); this->generateBezierExtraction(); } diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C index c9ec77ab..1d2aa6bc 100644 --- a/src/ASM/LR/ASMu3D.C +++ b/src/ASM/LR/ASMu3D.C @@ -105,7 +105,7 @@ bool ASMu3D::read (std::istream& is) return false; } - geo = lrspline.get(); + geomB = lrspline; return true; } @@ -129,7 +129,7 @@ void ASMu3D::clear (bool retainGeometry) delete tensorspline; delete tensorPrjBas; } - geo = nullptr; + geomB = nullptr; tensorspline = tensorPrjBas = nullptr; } @@ -231,12 +231,12 @@ bool ASMu3D::createProjectionBasis (bool init) std::swap(tensorspline,tensorPrjBas); std::swap(lrspline,projBasis); - geo = lrspline.get(); + geomB = lrspline; return true; } -LR::LRSplineVolume* ASMu3D::createLRfromTensor () +std::shared_ptr ASMu3D::createLRfromTensor () { if (tensorspline) { @@ -245,13 +245,13 @@ LR::LRSplineVolume* ASMu3D::createLRfromTensor () tensorspline = nullptr; } - return lrspline.get(); + return lrspline; } bool ASMu3D::generateFEMTopology () { - geo = this->createLRfromTensor(); + geomB = this->createLRfromTensor(); if (tensorPrjBas) { @@ -1836,9 +1836,9 @@ void ASMu3D::getBoundaryNodes (int lIndex, IntVec& nodes, int basis, bool ASMu3D::getOrder (int& p1, int& p2, int& p3) const { - p1 = geo->order(0); - p2 = geo->order(1); - p3 = geo->order(2); + p1 = geomB->order(0); + p2 = geomB->order(1); + p3 = geomB->order(2); return true; } diff --git a/src/ASM/LR/ASMu3D.h b/src/ASM/LR/ASMu3D.h index 3abce894..8717aa0f 100644 --- a/src/ASM/LR/ASMu3D.h +++ b/src/ASM/LR/ASMu3D.h @@ -117,7 +117,7 @@ public: virtual ~ASMu3D() {} //! \brief Returns the spline volume representing the geometry of this patch. - LR::LRSplineVolume* getVolume() { return this->createLRfromTensor(); } + LR::LRSplineVolume* getVolume() { return this->createLRfromTensor().get(); } //! \brief Returns the spline volume representing the geometry of this patch. const LR::LRSplineVolume* getVolume() const { return lrspline.get(); } @@ -685,7 +685,7 @@ protected: const IntSet& neighborIndices) const; //! \brief Converts current tensor spline object to LR-spline. - LR::LRSplineVolume* createLRfromTensor(); + std::shared_ptr createLRfromTensor(); public: //! \brief Returns the number of elements on a boundary. diff --git a/src/ASM/LR/ASMu3Dmx.C b/src/ASM/LR/ASMu3Dmx.C index 7cfd1def..2939b157 100644 --- a/src/ASM/LR/ASMu3Dmx.C +++ b/src/ASM/LR/ASMu3Dmx.C @@ -311,7 +311,7 @@ bool ASMu3Dmx::generateFEMTopology () std::cout <<"NEL = "<< nel <<" NNOD = "<< nnod << std::endl; #endif - geo = m_basis[geoBasis-1].get(); + geomB = m_basis[geoBasis-1]; return true; } diff --git a/src/ASM/LR/ASMu3Dmx.h b/src/ASM/LR/ASMu3Dmx.h index ae289cbf..b3f89ed9 100644 --- a/src/ASM/LR/ASMu3Dmx.h +++ b/src/ASM/LR/ASMu3Dmx.h @@ -63,7 +63,7 @@ public: //! \brief Copy constructor. ASMu3Dmx(const ASMu3Dmx& patch, const CharVec& n_f = CharVec(3,0)); //! \brief Empty destructor. - virtual ~ASMu3Dmx() { lrspline = nullptr; geo = nullptr; } + virtual ~ASMu3Dmx() {} //! \brief Returns the spline surface representing the basis of this patch. virtual LR::LRSplineVolume* getBasis(int basis = 1);