Merge pull request #966 from su2code/contiguous_cpoint

Contiguous storage of CPoint(s)

Common/include/geometry/CGeometry.hpp

@@ -191,7 +191,7 @@ class CGeometry {
   CPrimalGrid** elem;                    /*!< \brief Element vector (primal grid information). */
   CPrimalGrid** face;                    /*!< \brief Face vector (primal grid information). */
   CPrimalGrid*** bound;	                 /*!< \brief Boundary vector (primal grid information). */
-  CPoint** node;                         /*!< \brief Node vector (dual grid information). */
+  CPoint* nodes;                         /*!< \brief Node vector (dual grid information). */
   CEdge* edges;                          /*!< \brief Edge vector (dual grid information). */
   CVertex*** vertex;                     /*!< \brief Boundary Vertex vector (dual grid information). */
   CTurboVertex**** turbovertex;          /*!< \brief Boundary Vertex vector ordered for turbomachinery calculation(dual grid information). */
@@ -585,9 +585,9 @@ class CGeometry {
   inline virtual void SetPositive_ZArea(CConfig *config) {}
 
   /*!
-   * \brief Setas connectivity between points.
+   * \brief Set connectivity between points.
    */
-  inline virtual void SetPoint_Connectivity(void) {}
+  inline virtual void SetPoint_Connectivity() {}
 
   /*!
    * \brief Orders the RCM.

Common/include/geometry/CPhysicalGeometry.hpp

@@ -391,7 +391,7 @@ class CPhysicalGeometry final : public CGeometry {
   /*!
    * \brief Set points which surround a point.
    */
-  void SetPoint_Connectivity(void) override;
+  void SetPoint_Connectivity() override;
 
   /*!
    * \brief Set a renumbering using a Reverse Cuthill-McKee Algorithm
@@ -800,7 +800,7 @@ class CPhysicalGeometry final : public CGeometry {
    */
   void SetWallDistance(su2double val) override {
     for (unsigned long iPoint = 0; iPoint < GetnPoint(); iPoint++){
-      node[iPoint]->SetWall_Distance(val);
+      nodes->SetWall_Distance(iPoint, val);
     }
   }
 

Common/include/omp_structure.hpp

@@ -156,7 +156,7 @@ inline size_t computeStaticChunkSize(size_t totalWork,
 template<class T, class U>
 void parallelCopy(size_t size, const T* src, U* dst)
 {
-  SU2_OMP_FOR_STAT(4196)
+  SU2_OMP_FOR_STAT(4096)
   for(size_t i=0; i<size; ++i) dst[i] = src[i];
 }
 
@@ -169,7 +169,7 @@ void parallelCopy(size_t size, const T* src, U* dst)
 template<class T, class U>
 void parallelSet(size_t size, T val, U* dst)
 {
-  SU2_OMP_FOR_STAT(4196)
+  SU2_OMP_FOR_STAT(4096)
   for(size_t i=0; i<size; ++i) dst[i] = val;
 }
 

Common/include/toolboxes/C2DContainer.hpp

@@ -485,6 +485,52 @@ class C2DContainer :
   }
 };
 
+/*!
+ * \brief Useful typedefs with default template parameters
+ */
+template<class T> using su2vector = C2DContainer<unsigned long, T, StorageType::ColumnMajor, 64, DynamicSize, 1>;
+template<class T> using su2matrix = C2DContainer<unsigned long, T, StorageType::RowMajor,    64, DynamicSize, DynamicSize>;
+
+using su2activevector = su2vector<su2double>;
+using su2activematrix = su2matrix<su2double>;
+
+using su2passivevector = su2vector<passivedouble>;
+using su2passivematrix = su2matrix<passivedouble>;
+
+/*!
+ * \class CVectorOfMatrix
+ * \brief This contrived container is used to store small matrices in a contiguous manner
+ *        but still present the "su2double**" interface to the outside world.
+ *        The "interface" part should be replaced by something more efficient, e.g. a "matrix view".
+ */
+struct CVectorOfMatrix {
+  su2activevector storage;
+  su2matrix<su2double*> interface;
+  unsigned long M, N;
+
+  CVectorOfMatrix() = default;
+
+  CVectorOfMatrix(unsigned long length, unsigned long rows, unsigned long cols, su2double value = 0.0) {
+    resize(length, rows, cols, value);
+  }
+
+  void resize(unsigned long length, unsigned long rows, unsigned long cols, su2double value = 0.0) {
+    M = rows;
+    N = cols;
+    storage.resize(length*rows*cols) = value;
+    interface.resize(length,rows);
+
+    for(unsigned long i=0; i<length; ++i)
+      for(unsigned long j=0; j<rows; ++j)
+        interface(i,j) = &(*this)(i,j,0);
+  }
+
+  su2double& operator() (unsigned long i, unsigned long j, unsigned long k) { return storage(i*M*N + j*N + k); }
+  const su2double& operator() (unsigned long i, unsigned long j, unsigned long k) const { return storage(i*M*N + j*N + k); }
+
+  su2double** operator[] (unsigned long i) { return interface[i]; }
+  const su2double* const* operator[] (unsigned long i) const { return interface[i]; }
+};
 
 /*!
  * \class C2DDummyLastView
@@ -549,15 +595,3 @@ struct C3DDummyMiddleView
     return data(i,k);
   }
 };
-
-/*!
- * \brief Useful typedefs with default template parameters
- */
-template<class T> using su2vector = C2DContainer<unsigned long, T, StorageType::ColumnMajor, 64, DynamicSize, 1>;
-template<class T> using su2matrix = C2DContainer<unsigned long, T, StorageType::RowMajor,    64, DynamicSize, DynamicSize>;
-
-using su2activevector = su2vector<su2double>;
-using su2activematrix = su2matrix<su2double>;
-
-using su2passivevector = su2vector<passivedouble>;
-using su2passivematrix = su2matrix<passivedouble>;

Common/include/toolboxes/graph_toolbox.hpp

@@ -67,6 +67,24 @@ class CCompressedSparsePattern {
 
   CCompressedSparsePattern() = default;
 
+  /*!
+   * \brief Construct with default inner indices.
+   * \param[in] outerPtrBegin - Start of outer pointers.
+   * \param[in] outerPtrEnd - End of outer pointers.
+   * \param[in] defaultInnerIdx - Default value for inner indices.
+   */
+  template<class Iterator>
+  CCompressedSparsePattern(Iterator outerPtrBegin, Iterator outerPtrEnd, Index_t defaultInnerIdx)
+  {
+    const auto size = outerPtrEnd - outerPtrBegin;
+    m_outerPtr.resize(size);
+    Index_t k = 0;
+    for(auto it = outerPtrBegin; it != outerPtrEnd; ++it)
+      m_outerPtr(k++) = *it;
+
+    m_innerIdx.resize(m_outerPtr(size-1)) = defaultInnerIdx;
+  }
+
   /*!
    * \brief Construct from rvalue refs.
    * \note This is the most efficient constructor as no data copy occurs.
@@ -102,6 +120,25 @@ class CCompressedSparsePattern {
     assert(m_innerIdx.size() == m_outerPtr(m_outerPtr.size()-1));
   }
 
+  /*!
+   * \brief Build from a "list of lists" type object
+   * \param[in] lil - An object with operator [][] and method "size" e.g. vector<vector<?> >.
+   */
+  template<class T>
+  CCompressedSparsePattern(const T& lil)
+  {
+    m_outerPtr.resize(lil.size()+1);
+    m_outerPtr(0) = 0;
+    for(Index_t i=1; i < Index_t(m_outerPtr.size()); ++i)
+      m_outerPtr(i) = m_outerPtr(i-1) + lil[i-1].size();
+
+    m_innerIdx.resize(m_outerPtr(lil.size()));
+    Index_t k = 0;
+    for(Index_t i=0; i < Index_t(lil.size()); ++i)
+      for(Index_t j=0; j < Index_t(lil[i].size()); ++j)
+        m_innerIdx(k++) = lil[i][j];
+  }
+
   /*!
    * \brief Build a list of pointers to the diagonal entries of the pattern.
    */
@@ -172,6 +209,16 @@ class CCompressedSparsePattern {
     return m_innerIdx(m_outerPtr(iOuterIdx) + iNonZero);
   }
 
+  /*!
+   * \param[in] iOuterIdx - Outer index.
+   * \param[in] iNonZero - Relative position of the inner index.
+   * \return The index of the i'th inner index associated with the outer index.
+   */
+  inline Index_t& getInnerIdx(Index_t iOuterIdx, Index_t iNonZero) {
+    assert(iNonZero >= 0 && iNonZero < getNumNonZeros(iOuterIdx));
+    return m_innerIdx(m_outerPtr(iOuterIdx) + iNonZero);
+  }
+
   /*!
    * \param[in] iOuterIdx - Outer index (row/col).
    * \param[in] iInnerIdx - Inner index (col/row).
@@ -285,6 +332,7 @@ using CEdgeToNonZeroMap = C2DContainer<unsigned long, Index_t, StorageType::RowM
 
 
 using CCompressedSparsePatternUL = CCompressedSparsePattern<unsigned long>;
+using CCompressedSparsePatternL = CCompressedSparsePattern<long>;
 using CEdgeToNonZeroMapUL = CEdgeToNonZeroMap<unsigned long>;
 
 
@@ -335,14 +383,12 @@ CCompressedSparsePattern<Index_t> buildCSRPattern(Geometry_t& geometry,
        *    neighbors, not duplicating any existing neighbor. ---*/
       for(auto jPoint : addedNeighbors)
       {
-        auto point = geometry.node[jPoint];
-
         if(type == ConnectivityType::FiniteVolume)
         {
           /*--- For FVM we know the neighbors of point j directly. ---*/
-          for(unsigned short iNeigh = 0; iNeigh < point->GetnPoint(); ++iNeigh)
+          for(unsigned short iNeigh = 0; iNeigh < geometry.nodes->GetnPoint(jPoint); ++iNeigh)
           {
-            Index_t kPoint = point->GetPoint(iNeigh);
+            Index_t kPoint = geometry.nodes->GetPoint(jPoint, iNeigh);
 
             if(neighbors.count(kPoint) == 0) // no duplication
               newNeighbors.insert(kPoint);
@@ -351,9 +397,9 @@ CCompressedSparsePattern<Index_t> buildCSRPattern(Geometry_t& geometry,
         else // FiniteElement
         {
           /*--- For FEM we need the nodes of all elements that contain point j. ---*/
-          for(unsigned short iNeigh = 0; iNeigh < point->GetnElem(); ++iNeigh)
+          for(unsigned short iNeigh = 0; iNeigh < geometry.nodes->GetnElem(jPoint); ++iNeigh)
           {
-            auto elem = geometry.elem[point->GetElem(iNeigh)];
+            auto elem = geometry.elem[geometry.nodes->GetElem(jPoint, iNeigh)];
 
             for(unsigned short iNode = 0; iNode < elem->GetnNodes(); ++iNode)
             {

Common/src/CMultiGridQueue.cpp

@@ -183,9 +183,9 @@ void CMultiGridQueue::Update(unsigned long updatePoint, CGeometry *fineGrid) {
 
   RemoveCV(updatePoint);
 
-  for (auto iNode = 0u; iNode < fineGrid->node[updatePoint]->GetnPoint(); ++iNode) {
-    const auto jPoint = fineGrid->node[updatePoint]->GetPoint(iNode);
-    if (!fineGrid->node[jPoint]->GetAgglomerate())
+  for (auto iNode = 0u; iNode < fineGrid->nodes->GetnPoint(updatePoint); ++iNode) {
+    const auto jPoint = fineGrid->nodes->GetPoint(updatePoint,iNode);
+    if (!fineGrid->nodes->GetAgglomerate(jPoint))
       IncrPriorityCV(jPoint);
   }
 }

Common/src/fem_geometry_structure.cpp

@@ -352,7 +352,7 @@ CMeshFEM::CMeshFEM(CGeometry *geometry, CConfig *config) {
         of the points.            ---*/
   map<unsigned long,unsigned long> globalPointIDToLocalInd;
   for(unsigned long i=0; i<geometry->GetnPoint(); ++i)
-    globalPointIDToLocalInd[geometry->node[i]->GetGlobalIndex()] = i;
+    globalPointIDToLocalInd[geometry->nodes->GetGlobalIndex(i)] = i;
 
   /*----------------------------------------------------------------------------*/
   /*--- Step 1: Communicate the elements and the boundary elements to the    ---*/
@@ -468,7 +468,7 @@ CMeshFEM::CMeshFEM(CGeometry *geometry, CConfig *config) {
 
       unsigned long ind = LMI->second;
       for(unsigned short l=0; l<nDim; ++l)
-        doubleSendBuf[i].push_back(geometry->node[ind]->GetCoord(l));
+        doubleSendBuf[i].push_back(geometry->nodes->GetCoord(ind, l));
     }
   }
 
@@ -6950,7 +6950,7 @@ CDummyMeshFEM_DG::CDummyMeshFEM_DG(CConfig *config): CMeshFEM_DG() {
   elem                = NULL;
   face                = NULL;
   bound               = NULL;
-  node                = NULL;
+  nodes               = NULL;
   edges               = NULL;
   vertex              = NULL;
   nVertex             = NULL;

Common/src/geometry/CDummyGeometry.cpp

@@ -46,7 +46,7 @@ CDummyGeometry::CDummyGeometry(CConfig *config){
   elem                = NULL;
   face                = NULL;
   bound               = NULL;
-  node                = NULL;
+  nodes               = NULL;
   edges               = NULL;
   vertex              = NULL;
   nVertex             = NULL;