added: class for simple global node number establishment for LR splines

Apps/HDF5toVTx/HDF5toVTx.C

@@ -80,6 +80,7 @@ bool readBasis (std::vector<ASMbase*>& result, const std::string& name,
     std::string out;
     hdf.readString(geom.str(),out);
     ptype = out.substr(0,10) == "# LRSPLINE" ? ASM::LRSpline : ASM::Spline;
+    isLR = ptype == ASM::LRSpline;
     int gdim = 0;
     if (ptype == ASM::LRSpline)
       gdim = out.substr(11,7) == "SURFACE" ? 2 : 3;

CMakeLists.txt

@@ -104,6 +104,9 @@ endif()
 if(LRSPLINE_FOUND OR LRSpline_FOUND)
   file(GLOB LR_SRCS ${PROJECT_SOURCE_DIR}/src/ASM/LR/*.C)
   list(APPEND IFEM_SRCS ${LR_SRCS})
+else()
+  list(REMOVE_ITEM IFEM_SRCS
+                   ${PROJECT_SOURCE_DIR}/src/ASM/GlobalNodes.C)
 endif()
 list(APPEND CHECK_SOURCES ${IFEM_SRCS})
 add_library(IFEM ${IFEM_SRCS} ${TINYXML_SRCS})

src/ASM/ASMunstruct.h

@@ -46,6 +46,8 @@ 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
+    IntVec pMLGN; //!< Parallel MLGN mapping to use with MPI
 
     //! \brief Default constructor.
     explicit RefineData(bool rs = false) : refShare(rs) {}
@@ -217,6 +219,9 @@ class ASMunstruct : public ASMbase
   //! \brief Finds the node that is closest to the given point \b X.
   virtual std::pair<size_t,double> findClosestNode(const Vec3& X) const;
 
+  //! \brief Obtain the refinement basis.
+  virtual const LR::LRSpline* getRefinementBasis() const { return geo; }
+
 protected:
   //! \brief Refines the mesh adaptively.
   //! \param[in] prm Input data used to control the mesh refinement

src/ASM/DomainDecomposition.C

@@ -322,8 +322,10 @@ void DomainDecomposition::setupNodeNumbers(int basis, IntVec& lNodes,
     // need to expand to all bases for corners and edges
     for (size_t b = 1; b <= pch->getNoBasis(); ++b)
       cbasis.insert(b);
-  else // directly add nodes, cbasis remains empty
+  else { // directly add nodes, cbasis remains empty
+    std::cout << "we go here yo" << std::endl;
     pch->getBoundaryNodes(lidx, lNodes, 0, thick, orient, false);
+  }
 
   const ASM2D* pch2D = dynamic_cast<const ASM2D*>(pch);
   const ASM3D* pch3D = dynamic_cast<const ASM3D*>(pch);

src/ASM/GlobalNodes.C

@@ -0,0 +1,259 @@
+// $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 "ASMunstruct.h"
+#include "SIMbase.h"
+#include "DomainDecomposition.h"
+#include "Utilities.h"
+
+
+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;
+      ASMunstruct::Sort(u, v, orient, edgeFunctions);
+    } else {
+      int dir = (lidx-1)/4;
+      int u = dir == 0;
+      int v = 1 + (dir != 2);
+      ASMunstruct::Sort(u, v, orient, edgeFunctions);
+    }
+  } else if (dim == 2) {
+    int dir = (lidx-1)/2;
+    int u = dir == 0;
+    int v = 1 + (dir != 2);
+    ASMunstruct::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;
+}
+
+
+/*!
+  \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;
+  }
+};
+
+
+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 LR::LRSpline* pch : pchs) {
+    it->resize(pch->nBasisFunctions());
+    std::iota(it->begin(), it->end(), nNodes);
+    nNodes += pch->nBasisFunctions();
+    ++it;
+  }
+
+  // remap common nodes
+  InterfaceOrder ifOrder;
+  std::set<ASM::Interface, InterfaceOrder> ifset(ifOrder);
+  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;
+}
+
+
+GlobalNodes::IntVec
+GlobalNodes::calcDDMapping(const GlobalNodes::LRSplineVec& pchs,
+                           const std::vector<GlobalNodes::IntVec>& MLGN,
+                           const SIMbase& sim, int& nNodes)
+{
+  const ProcessAdm& adm = sim.getProcessAdm();
+
+  int minNode = 0;
+  if (adm.getProcId() > 0)
+    adm.receive(minNode, adm.getProcId()-1);
+
+  IntVec result(*std::max_element(MLGN.back().begin(), MLGN.back().end()) + 1);
+  std::iota(result.begin(), result.end(), minNode);
+  int maxNode = adm.getProcId() == 0 ? 0 : minNode;
+
+  std::map<int,int> old2new;
+  for (const auto& it : adm.dd.ghostConnections) {
+    int sidx = sim.getLocalPatchIndex(it.slave);
+    if (sidx < 1)
+      continue;
+
+    IntVec lNodes = GlobalNodes::getBoundaryNodes(*pchs[sidx-1], it.dim, it.sidx, 0);
+    for (int& i : lNodes)
+      i = result[MLGN[sidx-1][i]];
+
+    int nRecv;
+    adm.receive(nRecv, adm.dd.getPatchOwner(it.master));
+    if (nRecv =! lNodes.size()) {
+      std::cerr <<"\n *** GlobalNodes::calcDDMapping(): "
+                <<" Topology error, boundary size "
+                << nRecv << ", expected " << lNodes.size() << std::endl;
+      return IntVec();
+    }
+    IntVec glbNodes(lNodes.size());
+    adm.receive(glbNodes, adm.dd.getPatchOwner(it.master));
+    for (size_t i = 0; i < lNodes.size(); ++i)
+      old2new[lNodes[i]] = glbNodes[i];
+  }
+
+  // remap ghost nodes
+  for (auto& it : result)
+    utl::renumber(it, old2new, false);
+
+  // remap rest of our nodes
+  for (int i = 0; i < (int)result.size(); ++i)
+    if (old2new.find(i + minNode) == old2new.end()) {
+      std::map<int,int> old2new2;
+      old2new2[i + minNode] = maxNode++;
+      for (auto& it : result)
+        utl::renumber(it, old2new2, false);
+    }
+
+  if (adm.getProcId() < adm.getNoProcs()-1)
+    adm.send(maxNode, adm.getProcId()+1);
+
+  for (const auto& it : adm.dd.ghostConnections) {
+    int midx = sim.getLocalPatchIndex(it.master);
+    if (midx < 1)
+      continue;
+
+    IntVec glbNodes = GlobalNodes::getBoundaryNodes(*pchs[midx-1], it.dim,
+                                                    it.midx, it.orient);
+    for (size_t i = 0; i < glbNodes.size(); ++i)
+      glbNodes[i] = result[MLGN[midx-1][glbNodes[i]]];
+
+    adm.send(int(glbNodes.size()), adm.dd.getPatchOwner(it.slave));
+    adm.send(glbNodes, adm.dd.getPatchOwner(it.slave));
+  }
+
+#ifdef HAVE_MPI
+  nNodes = adm.allReduce(adm.getProcId() == adm.getNoProcs()-1 ? maxNode : 0, MPI_SUM);
+#endif
+
+  return result;
+}

src/ASM/GlobalNodes.h

@@ -0,0 +1,64 @@
+// $Id$
+//==============================================================================
+//!
+//! \file GlobalNodes.h
+//!
+//! \date Mar 13 2018
+//!
+//! \author Arne Morten Kvarving / SINTEF
+//!
+//! \brief Simple global node establishment for unstructured FE models.
+//!
+//==============================================================================
+
+#ifndef _GLOBAL_NODES_H_
+#define _GLOBAL_NODES_H_
+
+#include "Interface.h"
+#include <LRSpline/LRSplineSurface.h>
+
+#include <vector>
+
+class DomainDecomposition;
+class ProcessAdm;
+class SIMbase;
+
+
+/*!
+  \brief Class establishing global node numbers for unstructed FE models.
+*/
+
+class GlobalNodes
+{
+public:
+  typedef std::vector<int> IntVec; //!< Convenience typedef
+  typedef std::vector<const LR::LRSpline*> LRSplineVec; //!< Convenience typedef
+  typedef std::vector<ASM::Interface> InterfaceVec; //!< Convenience typedef
+
+  //! \brief Extract local boundary nodes for a LR spline.
+  //! \param lr The LR spline to extract boundary nodes for
+  //! \param dim The dimension of the boundary to extract
+  //! \param lidx The local index of the boundary to extract
+  //! \param orient Orientation of nodes on boundary
+  static IntVec getBoundaryNodes(const LR::LRSpline& lr,
+                                 int dim, int lidx, int orient);
+
+  //! \brief Calculate global node numbers for a FE model.
+   //! \param pchs The spline patches in the model
+  //! \param interfaces The topological connections for the spline patches
+  static std::vector<IntVec> calcGlobalNodes(const LRSplineVec& pchs,
+                                             const InterfaceVec& interfaces);
+
+  //! \brief Calculate parallel global node numbers for a FE model.
+  //! \param pchs The spline patches in the model
+  //! \param MLGN Process-global node numbers
+  //! \param sim The simulator holding patch owner information
+  //! \param adm The parallel process administrator
+  //! \param dd The domain decomposition holding ghost connections
+  //! \param[out] nNodes Total number of nodes in the model
+  static IntVec calcDDMapping(const LRSplineVec& pchs,
+                              const std::vector<GlobalNodes::IntVec>& MLGN,
+                              const SIMbase& sim, int& nNodes);
+};
+
+#endif