diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp index df4246a0ee36f2277e50d34398bf6de7a8d5ad3f..c82e9ce06ec8b7d43ce1c6cc2fec282f89dc1fee 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; } } }