From 8285983f65d93e1a4d13cf3018c489ce327b30a0 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 4 Aug 2020 23:30:13 +0200
Subject: [PATCH] Add centroid calculation in 3d
---
src/mesh/MeshData.hpp | 51 ++++++++++++++++++++++++++++++++++++++++---
1 file changed, 48 insertions(+), 3 deletions(-)
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index df4246a0e..c82e9ce06 100644
--- a/src/mesh/MeshData.hpp
+++ b/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;
}
}
}
--
GitLab