enable adaptive refinement for mixed multi-patch models

in particular, the refinements basis differ from the finite element
bases for some formulations, so a simple node numbering scheme for
the refinement basis is required to perform the multi-patch refinements.
use the GlobalNodes class for this

src/ASM/ASMunstruct.h

@@ -46,6 +46,7 @@ namespace LR //! Utilities for LR-splines.
     IntVec    options;  //!< Parameters used to control the refinement
     IntVec    elements; //!< 0-based indices of the elements to refine
     RealArray errors;   //!< List of error indicators for the elements
+    std::vector<IntVec> MLGN; //!< MLGN mapping to use for multipatch
 
     //! \brief Default constructor.
     explicit RefineData(bool rs = false) : refShare(rs) {}

src/ASM/LR/ASMLRSpline.C

@@ -16,6 +16,7 @@
 #include "LRSpline/LRSplineVolume.h"
 
 #include "ASMLRSpline.h"
+#include "GlobalNodes.h"
 #include "Vec3.h"
 #include "Vec3Oper.h"
 #include "ThreadGroups.h"
@@ -387,12 +388,12 @@ IntVec ASMLRSpline::getBoundaryCovered (const IntSet& nodes) const
   int numbEdges = (this->getNoParamDim() == 2) ? 4 : 6;
   for (int edge = 1; edge <= numbEdges; edge++)
   {
-    IntVec oneBoundary; // 1-based list of boundary nodes
-    this->getBoundaryNodes(edge,oneBoundary,1,1,0,true);
-    for (int i : nodes)
-      for (int j : oneBoundary)
-        if (geomB->getBasisfunction(i)->contains(*geomB->getBasisfunction(j-1)))
-          result.insert(j-1);
+    IntVec bnd = GlobalNodes::getBoundaryNodes(*refB,
+                                               this->getNoParamDim()-1, edge, 0);
+    for (const int i : nodes)
+      for (const int j : bnd)
+        if (refB->getBasisfunction(i)->contains(*refB->getBasisfunction(j)))
+          result.insert(j);
   }
 
   return IntVec(result.begin(), result.end());
@@ -404,7 +405,7 @@ IntVec ASMLRSpline::getOverlappingNodes (const IntSet& nodes, int dir) const
   IntSet result;
   for (int i : nodes)
   {
-    LR::Basisfunction* b = geomB->getBasisfunction(i);
+    const LR::Basisfunction* b = refB->getBasisfunction(i);
     for (LR::Element* el : b->support())
       for (LR::Basisfunction* basis : el->support())
       {

src/ASM/LR/ASMu2D.C

@@ -21,6 +21,7 @@
 #include "TimeDomain.h"
 #include "FiniteElement.h"
 #include "GlobalIntegral.h"
+#include "GlobalNodes.h"
 #include "LocalIntegral.h"
 #include "IntegrandBase.h"
 #include "CoordinateMapping.h"
@@ -363,10 +364,12 @@ void ASMu2D::extendRefinementDomain (IntSet& refineIndices,
 
   IntVec              bndry0;
   std::vector<IntVec> bndry1(4);
-  for (int i = 0; i < 4; i++)
+  for (int i = 1; i <= 4; i++)
   {
-    bndry0.push_back(this->getCorner((i%2)*2-1, (i/2)*2-1, 1));
-    this->getBoundaryNodes(1+i, bndry1[i], 1, 1, 0, true);
+    bndry0.push_back(GlobalNodes::getBoundaryNodes(*this->getBasis(ASM::REFINEMENT_BASIS),
+                                                   0, i, 0).front());
+    bndry1[i-1] = GlobalNodes::getBoundaryNodes(*this->getBasis(ASM::REFINEMENT_BASIS),
+                                                  1, i, 0);
   }
 
   // Add refinement from neighbors
@@ -377,7 +380,7 @@ void ASMu2D::extendRefinementDomain (IntSet& refineIndices,
     // Check if node is a corner node,
     // compute large extended domain (all directions)
     for (int edgeNode : bndry0)
-      if (edgeNode-1 == j)
+      if (edgeNode == j)
       {
         IntVec secondary = this->getOverlappingNodes(j);
         refineIndices.insert(secondary.begin(),secondary.end());
@@ -389,7 +392,7 @@ void ASMu2D::extendRefinementDomain (IntSet& refineIndices,
     // compute small extended domain (one direction)
     for (int edge = 0; edge < 4 && !done_with_this_node; edge++)
       for (int edgeNode : bndry1[edge])
-        if (edgeNode-1 == j)
+        if (edgeNode == j)
         {
           IntVec secondary = this->getOverlappingNodes(j, edge/2+1);
           refineIndices.insert(secondary.begin(),secondary.end());

src/ASM/LR/ASMu3D.C

@@ -21,6 +21,7 @@
 #include "TimeDomain.h"
 #include "FiniteElement.h"
 #include "GlobalIntegral.h"
+#include "GlobalNodes.h"
 #include "LocalIntegral.h"
 #include "IntegrandBase.h"
 #include "CoordinateMapping.h"
@@ -2253,18 +2254,18 @@ void ASMu3D::extendRefinementDomain (IntSet& refineIndices,
   const int nface =  6;
 
   IntVec bndry0;
-  for (int K = -1; K < 2; K += 2)
-    for (int J = -1; J < 2; J += 2)
-      for (int I = -1; I < 2; I += 2)
-        bndry0.push_back(this->getCorner(I,J,K,1));
+  const LR::LRSplineVolume* ref = this->getBasis(ASM::REFINEMENT_BASIS);
+  for (size_t c = 1; c <= 8; ++c)
+      bndry0.push_back(GlobalNodes::getBoundaryNodes(*ref,
+                                                     0, c, 0).front());
 
   std::vector<IntVec> bndry1(nedge);
   for (int j = 1; j <= nedge; j++)
-    this->getBoundary1Nodes(j , bndry1[j-1], 1, 0, true);
+    bndry1[j-1] = GlobalNodes::getBoundaryNodes(*ref, 1, j, 0);
 
   std::vector<IntVec> bndry2(nface);
   for (int j = 1; j <= nface; j++)
-    this->getBoundaryNodes(j, bndry2[j-1], 1, 1, 0, true);
+    bndry2[j-1] = GlobalNodes::getBoundaryNodes(*ref, 2, j, 0);
 
   // Add refinement from neighbors
   for (int j : neighborIndices)
@@ -2274,7 +2275,7 @@ void ASMu3D::extendRefinementDomain (IntSet& refineIndices,
     // Check if node is a corner node,
     // compute large extended domain (all directions)
     for (int edgeNode : bndry0)
-      if (edgeNode-1 == j)
+      if (edgeNode == j)
       {
         IntVec secondary = this->getOverlappingNodes(j);
         refineIndices.insert(secondary.begin(),secondary.end());
@@ -2287,7 +2288,7 @@ void ASMu3D::extendRefinementDomain (IntSet& refineIndices,
     int allowedDir;
     for (int edge = 0; edge < nedge && !done_with_this_node; edge++)
       for (int edgeNode : bndry1[edge])
-        if (edgeNode-1 == j)
+        if (edgeNode == j)
         {
           if (edge < 4)
             allowedDir = 6; // bin(110), allowed to grow in v- and w-direction
@@ -2305,7 +2306,7 @@ void ASMu3D::extendRefinementDomain (IntSet& refineIndices,
     // compute small extended domain (1 direction)
     for (int face = 0; face < nface && !done_with_this_node; face++)
       for (int edgeNode : bndry2[face])
-        if (edgeNode-1 == j)
+        if (edgeNode == j)
         {
           allowedDir = 1 << face/2;
           IntVec secondary = this->getOverlappingNodes(j,allowedDir);

src/SIM/AdaptiveSetup.C

@@ -21,7 +21,10 @@
 #include "IFEM.h"
 #include "tinyxml2.h"
 #ifdef HAS_LRSPLINE
+#include "ASMLRSpline.h"
+#include "GlobalNodes.h"
 #include <LRSpline/Basisfunction.h>
+#include <LRSpline/LRSplineSurface.h>
 #endif
 #include <fstream>
 #include <sstream>
@@ -331,16 +334,21 @@ int AdaptiveSetup::calcRefinement (LR::RefineData& prm, int iStep,
 
   if (scheme == ISOTROPIC_FUNCTION) // use errors per function
   {
-    if (thePatch->getNoRefineNodes() == thePatch->getNoNodes(1))
-      error.resize(model.getNoNodes(1),DblIdx(0.0,0));
-    else if (model.getNoPatches() == 1)
-      error.resize(thePatch->getNoRefineNodes(),DblIdx(0.0,0));
-    else
-    {
-      std::cerr <<" *** AdaptiveSetup::calcRefinement: Multi-patch refinement"
-                <<" is not available for mixed models."<< std::endl;
-      return -3;
-    }
+      size_t nNodes = 0;
+      if (model.getNoPatches() > 1) {
+#ifdef HAS_LRSPLINE
+        std::vector<const LR::LRSpline*> refBasis;
+        for (const ASMbase* pch : model.getFEModel())
+          refBasis.push_back(dynamic_cast<const ASMLRSpline*>(pch)->getRefinementBasis());
+
+        prm.MLGN = GlobalNodes::calcGlobalNodes(refBasis, model.getInterfaces());
+        nNodes = *std::max_element(prm.MLGN.back().begin(), prm.MLGN.back().end()) + 1;
+#endif
+      } else
+        nNodes = model.getPatch(1)->getNoRefineNodes();
+      error.reserve(nNodes);
+      for (i = 0; i < nNodes; i++)
+        error.push_back(DblIdx(0.0,i));
 
     for (i = 0; i < error.size(); i++)
       error[i].second = i;
@@ -359,7 +367,7 @@ int AdaptiveSetup::calcRefinement (LR::RefineData& prm, int iStep,
       // Insert into global error array
       for (i = 0; i < locErr.size(); i++)
         if (model.getNoPatches() > 1)
-          error[patch->getNodeID(i+1)-1].first += locErr[i];
+          error[prm.MLGN[patch->idx][i]].first += locErr[i];
         else
           error[i].first += locErr[i];
     }

src/SIM/SIMinput.C

@@ -1604,9 +1604,16 @@ bool SIMinput::refine (const LR::RefineData& prm, Vectors& sol)
       // Extract local indices from the vector of global indices
       int locId;
       for (int k : prm.elements)
-        if ((locId = myModel[i]->getNodeIndex(k+1)) > 0)
-          if (refineIndices[i].insert(locId-1).second)
+        if (!prm.MLGN.empty()) {
+          const auto it = std::find(prm.MLGN[i].begin(), prm.MLGN[i].end(), k);
+          if (it != prm.MLGN[i].end())
+            if (refineIndices[i].insert(std::distance(prm.MLGN[i].begin(), it)).second)
+              changed = true;
+        } else {
+         if ((locId = myModel[i]->getNodeIndex(k+1)) > 0)
+           if (refineIndices[i].insert(locId-1).second)
             changed = true;
+        }
 
       // Fetch boundary nodes covered (may need to pass this to other patches)
       IntVec bndry_nodes = pch[i]->getBoundaryCovered(refineIndices[i]);
@@ -1627,13 +1634,27 @@ bool SIMinput::refine (const LR::RefineData& prm, Vectors& sol)
       for (int k : bndry_nodes)
       {
         // Check if this boundary node appears on other patches
-        int globId = myModel[i]->getNodeID(k+1);
+        int globId;
+        if (prm.MLGN.empty())
+          globId = myModel[i]->getNodeID(k+1);
+        else
+          globId = prm.MLGN[i][k];
         for (size_t j = 0; j < myModel.size(); j++)
-          if (j != i && (locId = myModel[j]->getNodeIndex(globId)) > 0)
-            if (conformingIndices[j].insert(locId-1).second) {
-              changed = true;
-              conformingIndices[i].insert(k);
+          if (j != i) {
+            if (prm.MLGN.empty()) {
+              if ((locId = myModel[j]->getNodeIndex(globId)) > 0) {
+                conformingIndices[j].insert(locId-1);
+                conformingIndices[i].insert(k);
+              }
+            } else {
+              const auto it = std::find(prm.MLGN[j].begin(),
+                                        prm.MLGN[j].end(), globId);
+              if (it != prm.MLGN[j].end()) {
+                conformingIndices[j].insert(std::distance(prm.MLGN[j].begin(), it));
+                conformingIndices[i].insert(k);
+              }
             }
+          }
       }
     }
 

src/SIM/SIMinput.h

@@ -208,6 +208,10 @@ class SIMinput : public SIMbase
                                 coordCheck);
   }
 
+  //! \brief Obtain a const reference to model topology.
+  const std::vector<ASM::Interface>& getInterfaces() const
+  { return myInterfaces; }
+
 protected:
   //! \brief Helper method returning a stream for patch geometry input.
   //! \param[in] tag The XML-tag containing the patch geometry definition