diff --git a/src/ASM/ASMunstruct.C b/src/ASM/ASMunstruct.C index a1efa2e0..99c59042 100644 --- a/src/ASM/ASMunstruct.C +++ b/src/ASM/ASMunstruct.C @@ -53,4 +53,10 @@ bool ASMunstruct::refine (const LR::RefineData&, Vectors&, const char*) return false; } + +std::pair ASMunstruct::findClosestNode (const Vec3&) const +{ + return std::make_pair(0,-1.0); +} + #endif diff --git a/src/ASM/ASMunstruct.h b/src/ASM/ASMunstruct.h index 5ef29e86..34a2ce35 100644 --- a/src/ASM/ASMunstruct.h +++ b/src/ASM/ASMunstruct.h @@ -214,6 +214,9 @@ public: //! \param[in] refTol Mesh refinement threshold virtual bool refine(const RealFunc& refC, double refTol) = 0; + //! \brief Finds the node that is closest to the given point \b X. + virtual std::pair findClosestNode(const Vec3& X) const; + protected: //! \brief Refines the mesh adaptively. //! \param[in] prm Input data used to control the mesh refinement diff --git a/src/ASM/LR/ASMLRSpline.C b/src/ASM/LR/ASMLRSpline.C index f9f3676c..4651d0b2 100644 --- a/src/ASM/LR/ASMLRSpline.C +++ b/src/ASM/LR/ASMLRSpline.C @@ -15,10 +15,11 @@ #include "IFEM.h" #include "LRSpline/LRSplineSurface.h" #include "LRSpline/Basisfunction.h" -#include "Profiler.h" +#include "Vec3.h" +#include "Vec3Oper.h" #include "ThreadGroups.h" +#include "Profiler.h" #include -#include #ifdef USE_OPENMP #include @@ -88,30 +89,31 @@ void LR::getGaussPointParameters (const LR::LRSpline* lrspline, RealArray& uGP, } -void LR::generateThreadGroups (ThreadGroups& threadGroups, const LR::LRSpline* lr) +void LR::generateThreadGroups (ThreadGroups& threadGroups, + const LR::LRSpline* lr) { int nElement = lr->nElements(); - int nt = 1; #ifdef USE_OPENMP - nt = omp_get_max_threads(); -#endif - if (nt == 1) { - threadGroups[0].resize(1); - threadGroups[0][0].resize(nElement); - std::iota(threadGroups[0][0].begin(), threadGroups[0][0].end(), 0); - } else { - std::vector status(nElement, 0); // status vector for elements: -1 is unusable for current color, 0 is available, any other is assigned color - std::vector > answer; + if (omp_get_max_threads() > 1) + { + for (int i = 0; i < 2; i++) + threadGroups[i].clear(); + IntMat& answer = threadGroups[0]; + + IntVec status(nElement,0); // status vector for elements: + // -1 is unusable for current color, 0 is available, + // any other value is the assigned color + int fixedElements = 0; - int nColors = 0; - while(fixedElements < nElement) { + for (int nColors = 0; fixedElements < nElement; nColors++) + { // reset un-assigned element tags for (int i=0; i thisColor; // look for available elements + IntVec thisColor; for (auto e : lr->getAllElements() ) { int i = e->getId(); if (status[i] == 0) { @@ -119,7 +121,7 @@ void LR::generateThreadGroups (ThreadGroups& threadGroups, const LR::LRSpline* l thisColor.push_back(i); fixedElements++; for (auto b : e->support()) // for all basisfunctions with support here - for (auto el2 : b->support()) {// for all element this function supports + for (auto el2 : b->support()) {// for all elements this function supports int j = el2->getId(); if (status[j] == 0) // if not assigned a color yet status[j] = -1; // set as unavailable (with current color) @@ -127,11 +129,12 @@ void LR::generateThreadGroups (ThreadGroups& threadGroups, const LR::LRSpline* l } } answer.push_back(thisColor); - nColors++; } - - threadGroups[0] = answer; + return; } +#endif + + threadGroups.oneGroup(nElement); // No threading, all elements in one group } @@ -147,102 +150,105 @@ bool ASMunstruct::refine (const LR::RefineData& prm, nnod = geo->nBasisFunctions(); return true; } - - std::vector nf(sol.size()); - if (!prm.errors.empty() || !prm.elements.empty()) { - for (size_t j = 0; j < sol.size(); j++) - if (!(nf[j] = LR::extendControlPoints(geo,sol[j],this->getNoFields(1)))) - return false; - } - else + else if (prm.errors.empty() && prm.elements.empty()) return true; - if (doRefine(prm, geo)) + IntVec nf(sol.size()); + for (size_t j = 0; j < sol.size(); j++) + if (!(nf[j] = LR::extendControlPoints(geo,sol[j],this->getNoFields(1)))) + return false; + + if (!this->doRefine(prm,geo)) + return false; + + nnod = geo->nBasisFunctions(); + for (int i = sol.size()-1; i >= 0; i--) { + sol[i].resize(nf[i]*geo->nBasisFunctions()); + LR::contractControlPoints(geo,sol[i],nf[i]); + } + + if (fName) { - nnod = geo->nBasisFunctions(); + char fullFileName[256]; - for (int i = sol.size()-1; i >= 0; i--) { - sol[i].resize(nf[i]*geo->nBasisFunctions()); - LR::contractControlPoints(geo,sol[i],nf[i]); - } + strcpy(fullFileName, "lrspline_"); + strcat(fullFileName, fName); + std::ofstream lrOut(fullFileName); + lrOut << *geo; + lrOut.close(); - if (fName) - { - char fullFileName[256]; - - strcpy(fullFileName, "lrspline_"); + LR::LRSplineSurface* lr = dynamic_cast(geo); + if (lr) { + // open files for writing + strcpy(fullFileName, "param_"); strcat(fullFileName, fName); - std::ofstream lrOut(fullFileName); - lrOut << *geo; - lrOut.close(); + std::ofstream paramMeshFile(fullFileName); - LR::LRSplineSurface* lr = dynamic_cast(geo); - if (lr) { - // open files for writing - strcpy(fullFileName, "param_"); - strcat(fullFileName, fName); - std::ofstream paramMeshFile(fullFileName); + strcpy(fullFileName, "physical_"); + strcat(fullFileName, fName); + std::ofstream physicalMeshFile(fullFileName); - strcpy(fullFileName, "physical_"); - strcat(fullFileName, fName); - std::ofstream physicalMeshFile(fullFileName); + strcpy(fullFileName, "param_dot_"); + strcat(fullFileName, fName); + std::ofstream paramDotMeshFile(fullFileName); - strcpy(fullFileName, "param_dot_"); - strcat(fullFileName, fName); - std::ofstream paramDotMeshFile(fullFileName); + strcpy(fullFileName, "physical_dot_"); + strcat(fullFileName, fName); + std::ofstream physicalDotMeshFile(fullFileName); - strcpy(fullFileName, "physical_dot_"); - strcat(fullFileName, fName); - std::ofstream physicalDotMeshFile(fullFileName); + lr->writePostscriptMesh(paramMeshFile); + lr->writePostscriptElements(physicalMeshFile); + lr->writePostscriptFunctionSpace(paramDotMeshFile); + lr->writePostscriptMeshWithControlPoints(physicalDotMeshFile); - lr->writePostscriptMesh(paramMeshFile); - lr->writePostscriptElements(physicalMeshFile); - lr->writePostscriptFunctionSpace(paramDotMeshFile); - lr->writePostscriptMeshWithControlPoints(physicalDotMeshFile); - - // close all files - paramMeshFile.close(); - physicalMeshFile.close(); - paramDotMeshFile.close(); - physicalDotMeshFile.close(); - } + // close all files + paramMeshFile.close(); + physicalMeshFile.close(); + paramDotMeshFile.close(); + physicalDotMeshFile.close(); } - - IFEM::cout <<"Refined mesh: "<< geo->nElements() <<" elements " - << geo->nBasisFunctions() <<" nodes."<< std::endl; - - 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); - 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); - #endif - } - if (isLinIndep) - std::cout <<"...Passed."<< std::endl; - else { - std::cout <<"FAILED!!!"<< std::endl; - exit(228); - } - } - - return true; } - return false; + IFEM::cout <<"Refined mesh: "<< geo->nElements() <<" elements " + << geo->nBasisFunctions() <<" nodes."<< std::endl; + + 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); + 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); +#endif + } + if (isLinIndep) + std::cout <<"...Passed."<< std::endl; + else { + std::cout <<"FAILED!!!"<< std::endl; + exit(228); + } + } + + return true; } -bool ASMunstruct::doRefine(const LR::RefineData& prm, - LR::LRSpline* lrspline) + +bool ASMunstruct::doRefine (const LR::RefineData& prm, LR::LRSpline* lrspline) { // to pick up if LR splines get stuck while doing refinement, // print entry and exit point of this function + char doRefine = 0; + if (!prm.errors.empty()) + doRefine = 'E'; // Refine based on error indicators + else if (!prm.elements.empty()) + doRefine = 'I'; // Refine the specified elements + else + return doRefine; + double beta = prm.options.size() > 0 ? prm.options[0]/100.0 : 0.10; int multiplicity = prm.options.size() > 1 ? prm.options[1] : 1; if (multiplicity > 1) @@ -262,32 +268,24 @@ bool ASMunstruct::doRefine(const LR::RefineData& prm, int maxAspectRatio = prm.options.size() > 5 ? prm.options[5] : -1; bool closeGaps = prm.options.size() > 6 ? prm.options[6] != 0 : false; - char doRefine = 0; - if (!prm.errors.empty()) - doRefine = 'E'; // Refine based on error indicators - else if (!prm.elements.empty()) - doRefine = 'I'; // Refine the specified elements + // set refinement parameters + if (maxTjoints > 0) + lrspline->setMaxTjoints(maxTjoints); + if (maxAspectRatio > 0) + lrspline->setMaxAspectRatio((double)maxAspectRatio); + lrspline->setCloseGaps(closeGaps); + lrspline->setRefMultiplicity(multiplicity); + lrspline->setRefStrat(strat); - if (doRefine) { - // set refinement parameters - if (maxTjoints > 0) - lrspline->setMaxTjoints(maxTjoints); - if (maxAspectRatio > 0) - lrspline->setMaxAspectRatio((double)maxAspectRatio); - lrspline->setCloseGaps(closeGaps); - lrspline->setRefMultiplicity(multiplicity); - lrspline->setRefStrat(strat); + // do actual refinement + if (doRefine == 'E') + lrspline->refineByDimensionIncrease(prm.errors,beta); + else if (strat == LR_STRUCTURED_MESH) + lrspline->refineBasisFunction(prm.elements); + else + lrspline->refineElement(prm.elements); - // do actual refinement - if (doRefine == 'E') - lrspline->refineByDimensionIncrease(prm.errors,beta); - else if (strat == LR_STRUCTURED_MESH) - lrspline->refineBasisFunction(prm.elements); - else - lrspline->refineElement(prm.elements); - - lrspline->generateIDs(); - } + lrspline->generateIDs(); return doRefine; } @@ -375,26 +373,25 @@ IntVec ASMunstruct::getOverlappingNodes (const IntSet& nodes, int dir) const } -void ASMunstruct::Sort(int u, int v, int orient, - std::vector& functions) +void ASMunstruct::Sort (int u, int v, int orient, + std::vector& functions) { std::sort(functions.begin(), functions.end(), [u,v,orient](const LR::Basisfunction* a, const LR::Basisfunction* b) { - int p1 = a->getOrder(u); - int p2 = a->getOrder(v); + int i,p = a->getOrder(orient < 4 ? v : u); + int idx = orient & 4 ? v : u; + for (i = 0; i <= p; i++) + if ((*a)[idx][i] != (*b)[idx][i]) + return orient & 2 ? (*a)[idx][i] > (*b)[idx][i] + : (*a)[idx][i] < (*b)[idx][i]; - int idx1 = orient & 4 ? v : u; - int idx2 = orient & 4 ? u : v; - - for (int i = 0; i < 1 + (orient < 4 ? p2 : p1); ++i) - if ((*a)[idx1][i] != (*b)[idx1][i]) - return orient & 2 ? (*a)[idx1][i] > (*b)[idx1][i] - : (*a)[idx1][i] < (*b)[idx1][i]; - for (int i = 0; i < 1 + (orient < 4 ? p1 : p2); ++i) - if ((*a)[idx2][i] != (*b)[idx2][i]) - return orient & 1 ? (*a)[idx2][i] > (*b)[idx2][i] - : (*a)[idx2][i] < (*b)[idx2][i]; + p = a->getOrder(orient < 4 ? u : v); + idx = orient & 4 ? u : v; + for (i = 0; i <= p; i++) + if ((*a)[idx][i] != (*b)[idx][i]) + return orient & 1 ? (*a)[idx][i] > (*b)[idx][i] + : (*a)[idx][i] < (*b)[idx][i]; return false; }); @@ -402,10 +399,29 @@ void ASMunstruct::Sort(int u, int v, int orient, bool ASMunstruct::transferCntrlPtVars (LR::LRSpline* oldBasis, - const RealArray& oldVars, RealArray& newVars, + const RealArray& oldVars, + RealArray& newVars, int nGauss, int nf) const { oldBasis->rebuildDimension(nf); oldBasis->setControlPoints(const_cast(oldVars)); return this->transferCntrlPtVars(oldBasis,newVars,nGauss); } + + +std::pair ASMunstruct::findClosestNode (const Vec3& X) const +{ + double distance = 0.0; + size_t inod = 0, iclose = 0; + for (LR::Basisfunction* b : geo->getAllBasisfunctions()) + { + double d = (X-Vec3(&(*b->cp()),nsd)).length(); + if (++inod == 1 || d < distance) + { + iclose = inod; + distance = d; + } + } + + return std::make_pair(iclose,distance); +}