Feature/mass center

src/mesh/MeshData.hpp

@@ -38,7 +38,8 @@ class MeshData : public IMeshData
   NodeValuePerCell<const double> m_ljr;
   NodeValuePerCell<const Rd> m_njr;
   FaceValue<const Rd> m_xl;
-  CellValue<const Rd> m_xj;
+  CellValue<const Rd> m_cell_centroid;
+  CellValue<const Rd> m_cell_iso_barycenter;
   CellValue<const double> m_Vj;
   FaceValue<const double> m_ll;
 
@@ -91,7 +92,8 @@ class MeshData : public IMeshData
     }
   }
 
-  PUGS_INLINE void
+  PUGS_INLINE
+  void
   _computeFaceIsobarycenter()
   {   // Computes vertices isobarycenter
     if constexpr (Dimension == 1) {
@@ -114,28 +116,29 @@ class MeshData : public IMeshData
     }
   }
 
+  PUGS_INLINE
   void
-  _computeCellIsobarycenter()
+  _computeCellIsoBarycenter()
   {   // Computes vertices isobarycenter
     if constexpr (Dimension == 1) {
       const NodeValue<const Rd>& xr = m_mesh.xr();
 
       const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
 
-      CellValue<Rd> xj(m_mesh.connectivity());
+      CellValue<Rd> cell_iso_barycenter{m_mesh.connectivity()};
       parallel_for(
         m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
           const auto& cell_nodes = cell_to_node_matrix[j];
-          xj[j]                  = 0.5 * (xr[cell_nodes[0]] + xr[cell_nodes[1]]);
+          cell_iso_barycenter[j] = 0.5 * (xr[cell_nodes[0]] + xr[cell_nodes[1]]);
         });
-      m_xj = xj;
+      m_cell_iso_barycenter = cell_iso_barycenter;
     } else {
       const NodeValue<const Rd>& xr = m_mesh.xr();
 
       const CellValue<const double>& inv_cell_nb_nodes = m_mesh.connectivity().invCellNbNodes();
 
       const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
-      CellValue<Rd> xj(m_mesh.connectivity());
+      CellValue<Rd> cell_iso_barycenter{m_mesh.connectivity()};
       parallel_for(
         m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
           Rd X                   = zero;
@@ -143,9 +146,53 @@ class MeshData : public IMeshData
           for (size_t R = 0; R < cell_nodes.size(); ++R) {
             X += xr[cell_nodes[R]];
           }
-          xj[j] = inv_cell_nb_nodes[j] * X;
+          cell_iso_barycenter[j] = inv_cell_nb_nodes[j] * X;
         });
-      m_xj = xj;
+      m_cell_iso_barycenter = cell_iso_barycenter;
+    }
+  }
+
+  PUGS_INLINE
+  void
+  _computeCellCentroid()
+  {   // Computes vertices isobarycenter
+    const CellValue<const Rd> cell_iso_barycenter = this->cellIsoBarycenter();
+    if constexpr (Dimension == 1) {
+      m_cell_centroid = cell_iso_barycenter;
+    } else {
+      if constexpr (Dimension == 2) {
+        const CellValue<const double> Vj = this->Vj();
+        const NodeValue<const Rd>& xr    = m_mesh.xr();
+        const auto& cell_to_node_matrix  = m_mesh.connectivity().cellToNodeMatrix();
+
+        CellValue<Rd> cell_centroid{m_mesh.connectivity()};
+        parallel_for(
+          m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+            Rd sum = zero;
+
+            const auto& cell_nodes = cell_to_node_matrix[j];
+
+            for (size_t R = 0; R < cell_nodes.size(); ++R) {
+              const Rd& xr0 = xr[cell_nodes[R]];
+              const Rd& xr1 = xr[cell_nodes[(R + 1) % cell_nodes.size()]];
+
+              Rd xjxr0 = xr[cell_nodes[R]] - cell_iso_barycenter[j];
+              Rd xjxr1 = xr[cell_nodes[(R + 1) % cell_nodes.size()]] - cell_iso_barycenter[j];
+
+              const double Vjl = 0.5 * (xjxr0[0] * xjxr1[1] - xjxr0[1] * xjxr1[0]);
+
+              sum += Vjl * (xr0 + xr1 + cell_iso_barycenter[j]);
+            }
+
+            sum *= 1 / (3 * Vj[j]);
+
+            cell_centroid[j] = sum;
+          });
+        m_cell_centroid = cell_centroid;
+      } else {
+        std::cout << rang::fgB::yellow << "Centroids are not computed correctly in 3d" << rang::style::reset << '\n';
+        m_cell_centroid = cell_iso_barycenter;
+      }
     }
   }
 
@@ -440,14 +487,24 @@ class MeshData : public IMeshData
     return m_xl;
   }
 
+  PUGS_INLINE
+  CellValue<const Rd>
+  cellIsoBarycenter()
+  {
+    if (not m_cell_iso_barycenter.isBuilt()) {
+      this->_computeCellIsoBarycenter();
+    }
+    return m_cell_iso_barycenter;
+  }
+
   PUGS_INLINE
   CellValue<const Rd>
   xj()
   {
-    if (not m_xj.isBuilt()) {
-      this->_computeCellIsobarycenter();
+    if (not m_cell_centroid.isBuilt()) {
+      this->_computeCellCentroid();
     }
-    return m_xj;
+    return m_cell_centroid;
   }
 
   PUGS_INLINE