Feature/mass center

src/mesh/MeshData.hpp

@@ -155,7 +155,7 @@ class MeshData : public IMeshData
   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;
@@ -190,8 +190,53 @@ class MeshData : public IMeshData
           });
         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;
+        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;
       }
     }
   }