diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp index c15bb1a66ddacc34fe392a08d133e6f683b0e757..b312e4a8e3296fceb7ca0b747bdf5841200b8274 100644 --- a/src/mesh/MeshData.hpp +++ b/src/mesh/MeshData.hpp @@ -199,85 +199,83 @@ class MeshData : public IMeshData 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(); + } 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; + 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]; + 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()]]; + 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]; + 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]); + const double Vjl = 0.5 * (xjxr0[0] * xjxr1[1] - xjxr0[1] * xjxr1[0]); - sum += Vjl * (xr0 + xr1 + cell_iso_barycenter[j]); - } + sum += Vjl * (xr0 + xr1 + cell_iso_barycenter[j]); + } - sum *= 1 / (3 * Vj[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(); + 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]; + 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]; + 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]; + 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 xl = face_center[face_id]; - const Rd xjxl = xl - xj; + const Rd xjxl = xl - xj; - const auto& face_nodes = face_to_node_matrix[face_id]; + const auto& face_nodes = face_to_node_matrix[face_id]; - const Rd xl_plus_xj = xl + xj; + const Rd xl_plus_xj = xl + xj; - double sign = (cell_face_is_reversed(j, i_face)) ? -1 : 1; + 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()]]; + 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 Rd xjxr0 = xr0 - xj; + const Rd xjxr1 = xr1 - xj; - const double Vjlr = dot(crossProduct(xjxr0, xjxr1), xjxl); + const double Vjlr = dot(crossProduct(xjxr0, xjxr1), xjxl); - sum += (sign * Vjlr) * (xl_plus_xj + xr0 + xr1); - } + sum += (sign * Vjlr) * (xl_plus_xj + xr0 + xr1); } + } - sum *= 1 / (24 * Vj[j]); + sum *= 1 / (24 * Vj[j]); - cell_centroid[j] = sum; - }); + cell_centroid[j] = sum; + }); - m_cell_centroid = cell_centroid; - } + m_cell_centroid = cell_centroid; } }