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,98 @@ 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()
+  {
+    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 {
+        const auto& face_center           = this->xl();
+        const CellValue<const double> Vj  = this->Vj();
+        const NodeValue<const Rd>& xr     = m_mesh.xr();
+        const auto& cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
+        const auto& face_to_node_matrix   = m_mesh.connectivity().faceToNodeMatrix();
+        const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
+
+        CellValue<Rd> cell_centroid{m_mesh.connectivity()};
+        parallel_for(
+          m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+            const Rd xj = m_cell_iso_barycenter[j];
+
+            const auto& cell_faces = cell_to_face_matrix[j];
+
+            Rd sum = zero;
+            for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
+              const FaceId face_id = cell_faces[i_face];
+
+              const Rd xl = face_center[face_id];
+
+              const Rd xjxl = xl - xj;
+
+              const auto& face_nodes = face_to_node_matrix[face_id];
+
+              const Rd xl_plus_xj = xl + xj;
+
+              double sign = (cell_face_is_reversed(j, i_face)) ? -1 : 1;
+
+              for (size_t i_face_node = 0; i_face_node < face_nodes.size(); ++i_face_node) {
+                const Rd& xr0 = xr[face_nodes[i_face_node]];
+                const Rd& xr1 = xr[face_nodes[(i_face_node + 1) % face_nodes.size()]];
+
+                const Rd xjxr0 = xr0 - xj;
+                const Rd xjxr1 = xr1 - xj;
+
+                const double Vjlr = (crossProduct(xjxr0, xjxr1), xjxl);
+
+                sum += (sign * Vjlr) * (xl_plus_xj + xr0 + xr1);
+              }
+            }
+
+            sum *= 1 / (24 * Vj[j]);
+
+            cell_centroid[j] = sum;
+          });
+
+        m_cell_centroid = cell_centroid;
+      }
     }
   }
 
@@ -440,14 +532,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