From e7095543aa145d19db932325addc7ebc0318f4d6 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Mon, 3 Aug 2020 22:07:15 +0200
Subject: [PATCH] Add mass center calculation in 2D

It is replace the iso-barycenter when using MeshData::xj()

Iso-barycenter is now accessed through MeshData::cellIsoBarycenter() function

note that the functions of MeshData will be eventually renamed
---
 src/mesh/MeshData.hpp | 81 ++++++++++++++++++++++++++++++++++++-------
 1 file changed, 69 insertions(+), 12 deletions(-)

diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index aeddc582f..df4246a0e 100644
--- a/src/mesh/MeshData.hpp
+++ b/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
-- 
GitLab