Skip to content
Snippets Groups Projects
Commit 3c8c4274 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Merge branch 'feature/mass-center' into 'develop'

Feature/mass center

See merge request !47
parents b7527d95 8285983f
No related branches found
No related tags found
1 merge request!47Feature/mass center
...@@ -38,7 +38,8 @@ class MeshData : public IMeshData ...@@ -38,7 +38,8 @@ class MeshData : public IMeshData
NodeValuePerCell<const double> m_ljr; NodeValuePerCell<const double> m_ljr;
NodeValuePerCell<const Rd> m_njr; NodeValuePerCell<const Rd> m_njr;
FaceValue<const Rd> m_xl; 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; CellValue<const double> m_Vj;
FaceValue<const double> m_ll; FaceValue<const double> m_ll;
...@@ -91,7 +92,8 @@ class MeshData : public IMeshData ...@@ -91,7 +92,8 @@ class MeshData : public IMeshData
} }
} }
PUGS_INLINE void PUGS_INLINE
void
_computeFaceIsobarycenter() _computeFaceIsobarycenter()
{ // Computes vertices isobarycenter { // Computes vertices isobarycenter
if constexpr (Dimension == 1) { if constexpr (Dimension == 1) {
...@@ -114,28 +116,29 @@ class MeshData : public IMeshData ...@@ -114,28 +116,29 @@ class MeshData : public IMeshData
} }
} }
PUGS_INLINE
void void
_computeCellIsobarycenter() _computeCellIsoBarycenter()
{ // Computes vertices isobarycenter { // Computes vertices isobarycenter
if constexpr (Dimension == 1) { if constexpr (Dimension == 1) {
const NodeValue<const Rd>& xr = m_mesh.xr(); const NodeValue<const Rd>& xr = m_mesh.xr();
const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix(); 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( parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[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 { } else {
const NodeValue<const Rd>& xr = m_mesh.xr(); const NodeValue<const Rd>& xr = m_mesh.xr();
const CellValue<const double>& inv_cell_nb_nodes = m_mesh.connectivity().invCellNbNodes(); const CellValue<const double>& inv_cell_nb_nodes = m_mesh.connectivity().invCellNbNodes();
const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix(); 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( parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
Rd X = zero; Rd X = zero;
...@@ -143,9 +146,98 @@ class MeshData : public IMeshData ...@@ -143,9 +146,98 @@ class MeshData : public IMeshData
for (size_t R = 0; R < cell_nodes.size(); ++R) { for (size_t R = 0; R < cell_nodes.size(); ++R) {
X += xr[cell_nodes[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()
{
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 {
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;
}
} }
} }
...@@ -440,14 +532,24 @@ class MeshData : public IMeshData ...@@ -440,14 +532,24 @@ class MeshData : public IMeshData
return m_xl; 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 PUGS_INLINE
CellValue<const Rd> CellValue<const Rd>
xj() xj()
{ {
if (not m_xj.isBuilt()) { if (not m_cell_centroid.isBuilt()) {
this->_computeCellIsobarycenter(); this->_computeCellCentroid();
} }
return m_xj; return m_cell_centroid;
} }
PUGS_INLINE PUGS_INLINE
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment