Fixed: Clear the threadGroups before (re)generating them.

Added: Implementation of ASMunstruct::findClosestNode.
Changed: Some cosmetics in other ASMunstruct methods.

src/ASM/ASMunstruct.C

@@ -53,4 +53,10 @@ bool ASMunstruct::refine (const LR::RefineData&, Vectors&, const char*)
   return false;
 }
 
+
+std::pair<size_t,double> ASMunstruct::findClosestNode (const Vec3&) const
+{
+  return std::make_pair(0,-1.0);
+}
+
 #endif

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<size_t,double> findClosestNode(const Vec3& X) const;
+
 protected:
   //! \brief Refines the mesh adaptively.
   //! \param[in] prm Input data used to control the mesh refinement

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 <fstream>
-#include <numeric>
 
 #ifdef USE_OPENMP
 #include <omp.h>
@@ -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();
+  if (omp_get_max_threads() > 1)
-#endif
+  {
-  if (nt == 1) {
+    for (int i = 0; i < 2; i++)
-    threadGroups[0].resize(1);
+      threadGroups[i].clear();
-    threadGroups[0][0].resize(nElement);
+    IntMat& answer = threadGroups[0];
-    std::iota(threadGroups[0][0].begin(), threadGroups[0][0].end(), 0);
+
-  } else {
+    IntVec status(nElement,0); // status vector for elements:
-    std::vector<int> status(nElement, 0); // status vector for elements: -1 is unusable for current color, 0 is available, any other is assigned color
+    // -1 is unusable for current color, 0 is available,
-    std::vector<std::vector<int> > answer;
+    // any other value is the assigned color
+
     int fixedElements = 0;
-    int nColors       = 0;
+    for (int nColors = 0; fixedElements < nElement; nColors++)
-    while(fixedElements < nElement) {
+    {
       // reset un-assigned element tags
       for (int i=0; i<nElement; i++)
         if (status[i]<0)
           status[i] = 0;
-      std::vector<int> 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++;
     }
-
+    return;
-    threadGroups[0] = answer;
   }
+#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;
   }
-
+  else if (prm.errors.empty() && prm.elements.empty())
-  std::vector<int> 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
     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--) {
+    strcpy(fullFileName, "lrspline_");
-      sol[i].resize(nf[i]*geo->nBasisFunctions());
+    strcat(fullFileName, fName);
-      LR::contractControlPoints(geo,sol[i],nf[i]);
+    std::ofstream lrOut(fullFileName);
-    }
+    lrOut << *geo;
+    lrOut.close();
 
-    if (fName)
+    LR::LRSplineSurface* lr = dynamic_cast<LR::LRSplineSurface*>(geo);
-    {
+    if (lr) {
-      char fullFileName[256];
+      // open files for writing
-
+      strcpy(fullFileName, "param_");
-      strcpy(fullFileName, "lrspline_");
       strcat(fullFileName, fName);
-      std::ofstream lrOut(fullFileName);
+      std::ofstream paramMeshFile(fullFileName);
-      lrOut << *geo;
-      lrOut.close();
 
-      LR::LRSplineSurface* lr = dynamic_cast<LR::LRSplineSurface*>(geo);
+      strcpy(fullFileName, "physical_");
-      if (lr) {
+      strcat(fullFileName, fName);
-        // open files for writing
+      std::ofstream physicalMeshFile(fullFileName);
-        strcpy(fullFileName, "param_");
-        strcat(fullFileName, fName);
-        std::ofstream paramMeshFile(fullFileName);
 
-        strcpy(fullFileName, "physical_");
+      strcpy(fullFileName, "param_dot_");
-        strcat(fullFileName, fName);
+      strcat(fullFileName, fName);
-        std::ofstream physicalMeshFile(fullFileName);
+      std::ofstream paramDotMeshFile(fullFileName);
 
-        strcpy(fullFileName, "param_dot_");
+      strcpy(fullFileName, "physical_dot_");
-        strcat(fullFileName, fName);
+      strcat(fullFileName, fName);
-        std::ofstream paramDotMeshFile(fullFileName);
+      std::ofstream physicalDotMeshFile(fullFileName);
 
-        strcpy(fullFileName, "physical_dot_");
+      lr->writePostscriptMesh(paramMeshFile);
-        strcat(fullFileName, fName);
+      lr->writePostscriptElements(physicalMeshFile);
-        std::ofstream physicalDotMeshFile(fullFileName);
+      lr->writePostscriptFunctionSpace(paramDotMeshFile);
+      lr->writePostscriptMeshWithControlPoints(physicalDotMeshFile);
 
-        lr->writePostscriptMesh(paramMeshFile);
+      // close all files
-        lr->writePostscriptElements(physicalMeshFile);
+      paramMeshFile.close();
-        lr->writePostscriptFunctionSpace(paramDotMeshFile);
+      physicalMeshFile.close();
-        lr->writePostscriptMeshWithControlPoints(physicalDotMeshFile);
+      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;
+  // set refinement parameters
-  if (!prm.errors.empty())
+  if (maxTjoints > 0)
-    doRefine = 'E'; // Refine based on error indicators
+    lrspline->setMaxTjoints(maxTjoints);
-  else if (!prm.elements.empty())
+  if (maxAspectRatio > 0)
-    doRefine = 'I'; // Refine the specified elements
+    lrspline->setMaxAspectRatio((double)maxAspectRatio);
+  lrspline->setCloseGaps(closeGaps);
+  lrspline->setRefMultiplicity(multiplicity);
+  lrspline->setRefStrat(strat);
 
-  if (doRefine) {
+  // do actual refinement
-    // set refinement parameters
+  if (doRefine == 'E')
-    if (maxTjoints > 0)
+    lrspline->refineByDimensionIncrease(prm.errors,beta);
-      lrspline->setMaxTjoints(maxTjoints);
+  else if (strat == LR_STRUCTURED_MESH)
-    if (maxAspectRatio > 0)
+    lrspline->refineBasisFunction(prm.elements);
-      lrspline->setMaxAspectRatio((double)maxAspectRatio);
+  else
-    lrspline->setCloseGaps(closeGaps);
+    lrspline->refineElement(prm.elements);
-    lrspline->setRefMultiplicity(multiplicity);
-    lrspline->setRefStrat(strat);
 
-    // do actual refinement
+  lrspline->generateIDs();
-    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();
-  }
 
   return doRefine;
 }
@@ -375,26 +373,25 @@ IntVec ASMunstruct::getOverlappingNodes (const IntSet& nodes, int dir) const
 }
 
 
-void ASMunstruct::Sort(int u, int v, int orient,
+void ASMunstruct::Sort (int u, int v, int orient,
-                       std::vector<LR::Basisfunction*>& functions)
+                        std::vector<LR::Basisfunction*>& functions)
 {
   std::sort(functions.begin(), functions.end(),
             [u,v,orient](const LR::Basisfunction* a, const LR::Basisfunction* b)
             {
-              int p1 = a->getOrder(u);
+              int i,p = a->getOrder(orient < 4 ? v : u);
-              int p2 = a->getOrder(v);
+              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;
+              p   = a->getOrder(orient < 4 ? u : v);
-              int idx2 = orient & 4 ? u : v;
+              idx = orient & 4 ? u : v;
-
+              for (i = 0; i <= p; i++)
-              for (int i = 0; i < 1 + (orient < 4 ? p2 : p1); ++i)
+                if ((*a)[idx][i] != (*b)[idx][i])
-                if ((*a)[idx1][i] != (*b)[idx1][i])
+                  return orient & 1 ? (*a)[idx][i] > (*b)[idx][i]
-                  return orient & 2 ? (*a)[idx1][i] > (*b)[idx1][i]
+                                    : (*a)[idx][i] < (*b)[idx][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];
 
               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<RealArray&>(oldVars));
   return this->transferCntrlPtVars(oldBasis,newVars,nGauss);
 }
+
+
+std::pair<size_t,double> 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);
+}