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

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
parent ec8e0e2d
No related branches found
No related tags found
1 merge request!47Feature/mass center
......@@ -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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment