added: support for mixed multipatch adaptive simulations

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/ASMLRSpline.h

@@ -108,6 +108,9 @@ class ASMLRSpline : public ASMbase, public ASMunstruct
   //! \brief Returns parameter values and node numbers of the domain corners.
   virtual bool getParameterDomain(Real2DMat&, IntVec*) const;
 
+  //! \brief Returns a const pointer to refinement basis.
+  const LR::LRSpline* getRefinementBasis() const { return refB.get(); }
+
   //! \brief Refines the mesh adaptively.
   //! \param[in] prm Input data used to control the mesh refinement
   //! \param sol Control point results values that are transferred to new mesh

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());
@@ -2464,9 +2467,6 @@ Field* ASMu2D::getProjectedField (const Vector& coefs) const
 
 Fields* ASMu2D::getProjectedFields (const Vector& coefs, size_t) const
 {
-  if (!this->separateProjectionBasis())
-    return nullptr;
-
   size_t ncmp = coefs.size() / this->getNoProjectionNodes();
   if (ncmp*this->getNoProjectionNodes() == coefs.size())
     return new LRSplineFields2D(this->getBasis(ASM::PROJECTION_BASIS),

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"
@@ -2004,9 +2005,6 @@ Field* ASMu3D::getProjectedField (const Vector& coefs) const
 
 Fields* ASMu3D::getProjectedFields (const Vector& coefs, size_t) const
 {
-  if (!this->separateProjectionBasis())
-    return nullptr;
-
   size_t ncmp = coefs.size() / this->getNoProjectionNodes();
   if (ncmp*this->getNoProjectionNodes() == coefs.size())
     return new LRSplineFields3D(this->getBasis(ASM::PROJECTION_BASIS),coefs,ncmp);
@@ -2254,18 +2252,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)
@@ -2275,7 +2273,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());
@@ -2288,7 +2286,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
@@ -2306,7 +2304,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/ASM/LR/GlobalNodes.C

@@ -0,0 +1,186 @@
+// $Id$
+//==============================================================================
+//!
+//! \file GlobalNodes.C
+//!
+//! \date Mar 13 2018
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Simple global node establishment for unstructured FE models.
+//!
+//==============================================================================
+
+#include "GlobalNodes.h"
+
+#include "ASMLRSpline.h"
+#include "Utilities.h"
+
+#include <numeric>
+
+namespace {
+
+/*!
+  \brief Class for ordering interfaces in processing order.
+*/
+
+class InterfaceOrder {
+public:
+  //! \brief Comparison operator for interfaces
+  //! \param A First interface
+  //! \param B Second interface
+  bool operator()(const ASM::Interface& A, const ASM::Interface& B) const
+  {
+    if (A.master != B.master)
+      return A.master < B.master;
+
+    if (A.slave != B.slave)
+      return A.slave < B.slave;
+
+    if (A.dim != B.dim)
+      return A.dim < B.dim;
+
+    return A.midx < B.midx;
+  }
+};
+
+}
+
+
+GlobalNodes::IntVec
+GlobalNodes::getBoundaryNodes (const LR::LRSpline& lr,
+                               int dim, int lidx, int orient)
+{
+  LR::parameterEdge edge{};
+  if (dim == 0) {
+    if (lr.nVariate() == 2) {
+      switch (lidx) {
+        case 1: edge = LR::WEST | LR::SOUTH; break;
+        case 2: edge = LR::EAST | LR::SOUTH; break;
+        case 3: edge = LR::WEST | LR::NORTH; break;
+        case 4: edge = LR::EAST | LR::NORTH; break;
+      }
+    } else {
+      switch (lidx) {
+        case 1: edge = LR::WEST | LR::SOUTH | LR::BOTTOM; break;
+        case 2: edge = LR::EAST | LR::SOUTH | LR::BOTTOM; break;
+        case 3: edge = LR::WEST | LR::NORTH | LR::BOTTOM; break;
+        case 4: edge = LR::EAST | LR::NORTH | LR::BOTTOM; break;
+        case 5: edge = LR::WEST | LR::SOUTH | LR::TOP; break;
+        case 6: edge = LR::EAST | LR::SOUTH | LR::TOP; break;
+        case 7: edge = LR::WEST | LR::NORTH | LR::TOP; break;
+        case 8: edge = LR::EAST | LR::NORTH | LR::TOP; break;
+      }
+    }
+  } else if (dim == 1) {
+    if (lr.nVariate() == 2) {
+      switch (lidx) {
+        case 1: edge = LR::WEST; break;
+        case 2: edge = LR::EAST; break;
+        case 3: edge = LR::SOUTH; break;
+        case 4: edge = LR::NORTH; break;
+        default: break;
+      }
+    } else {
+      switch (lidx) {
+        case 1: edge = LR::BOTTOM | LR::SOUTH; break;
+        case 2: edge = LR::BOTTOM | LR::NORTH; break;
+        case 3: edge = LR::TOP | LR::SOUTH; break;
+        case 4: edge = LR::TOP | LR::NORTH; break;
+        case 5: edge = LR::BOTTOM | LR::WEST; break;
+        case 6: edge = LR::BOTTOM | LR::EAST; break;
+        case 7: edge = LR::TOP | LR::WEST; break;
+        case 8: edge = LR::TOP | LR::WEST; break;
+        case 9: edge = LR::SOUTH | LR::WEST; break;
+        case 10: edge = LR::SOUTH | LR::EAST; break;
+        case 11: edge = LR::NORTH | LR::WEST; break;
+        case 12: edge = LR::NORTH | LR::EAST; break;
+      }
+    }
+  } else if (dim == 2) {
+    switch (lidx) {
+      case 1: edge = LR::WEST; break;
+      case 2: edge = LR::EAST; break;
+      case 3: edge = LR::SOUTH; break;
+      case 4: edge = LR::NORTH; break;
+      case 5: edge = LR::BOTTOM; break;
+      case 6: edge = LR::TOP; break;
+    }
+  }
+
+  std::vector<LR::Basisfunction*> edgeFunctions;
+  lr.getEdgeFunctions(edgeFunctions, edge);
+
+  if (dim == 1) {
+    if (lr.nVariate() == 2) {
+      int v = (lidx == 1 || lidx == 2) ? 0 : 1;
+      int u = 1-v;
+      ASMLRSpline::Sort(u, v, orient, edgeFunctions);
+    } else {
+      int dir = (lidx-1)/4;
+      int u = dir == 0;
+      int v = 1 + (dir != 2);
+      ASMLRSpline::Sort(u, v, orient, edgeFunctions);
+    }
+  } else if (dim == 2) {
+    int dir = (lidx-1)/2;
+    int u = dir == 0;
+    int v = 1 + (dir != 2);
+    ASMLRSpline::Sort(u, v, orient, edgeFunctions);
+  }
+
+  GlobalNodes::IntVec lNodes;
+  lNodes.reserve(edgeFunctions.size());
+  for (const LR::Basisfunction* func : edgeFunctions)
+    lNodes.push_back(func->getId());
+
+  return lNodes;
+}
+
+
+std::vector<GlobalNodes::IntVec>
+GlobalNodes::calcGlobalNodes (const GlobalNodes::LRSplineVec& pchs,
+                              const GlobalNodes::InterfaceVec& interfaces)
+{
+  // count total number of nodes
+  size_t nNodes = 0;
+  std::vector<GlobalNodes::IntVec> result(pchs.size());
+  auto it = result.begin();
+  for (const auto& pch : pchs) {
+    it->resize(pch->nBasisFunctions());
+    std::iota(it->begin(), it->end(), nNodes);
+    nNodes += pch->nBasisFunctions();
+    ++it;
+  }
+
+  // remap common nodes
+  std::set<ASM::Interface, InterfaceOrder> ifset;
+  for (const ASM::Interface& it : interfaces)
+    ifset.insert(it);
+  for (size_t i = 0; i < pchs.size(); ++i) {
+    std::map<int,int> old2new;
+    for (const ASM::Interface& it : ifset) {
+      if (it.master != (int)i+1)
+        continue;
+
+      IntVec mNodes = getBoundaryNodes(*pchs[i], it.dim, it.midx, 0);
+      IntVec sNodes = getBoundaryNodes(*pchs[it.slave-1], it.dim, it.sidx, it.orient);
+      for (size_t n = 0; n < mNodes.size(); ++n)
+        old2new[result[it.slave-1][sNodes[n]]] = result[i][mNodes[n]];
+    }
+
+    // renumber
+    for (size_t j = i; j < pchs.size(); ++j)
+      for (int& it : result[j])
+        utl::renumber(it, old2new, false);
+
+    // compress
+    int maxNode = *std::max_element(result[i].begin(), result[i].end());
+    for (size_t j = i+1; j < pchs.size(); ++j)
+      for (int& n : result[j])
+        if (n > maxNode)
+          n = ++maxNode;
+  }
+
+  return result;
+}