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

Add calculation of face area/lenght in mesh data 

Also replace use of unnecessary `std::shared_ptr` which was duplicate
reference counters of `(Sub)ItemValue` internal ones. Now check if the
variable is built which is even cleaner.
parent cd9a2dbd
No related branches found
No related tags found
1 merge request!43Feature/glace boundary conditions
Show whitespace changes
Inline Side-by-side
Showing
with 69 additions and 32 deletions
src/mesh/MeshData.hpp
+ 69
− 32
View file @ 468b5bc8
...@@ -32,16 +32,44 @@ class MeshData : public IMeshData ...@@ -32,16 +32,44 @@ class MeshData : public IMeshData
private: private:
const MeshType& m_mesh; const MeshType& m_mesh;
std::shared_ptr<NodeValuePerCell<const Rd>> m_Cjr; NodeValuePerCell<const Rd> m_Cjr;
std::shared_ptr<NodeValuePerCell<const double>> m_ljr; NodeValuePerCell<const double> m_ljr;
std::shared_ptr<NodeValuePerCell<const Rd>> m_njr; NodeValuePerCell<const Rd> m_njr;
FaceValue<const Rd> m_xl; FaceValue<const Rd> m_xl;
std::shared_ptr<CellValue<const Rd>> m_xj; CellValue<const Rd> m_xj;
std::shared_ptr<CellValue<const double>> m_Vj; CellValue<const double> m_Vj;
std::shared_ptr<NodeValuePerFace<const Rd>> m_Nlr; FaceValue<const double> m_ll;
NodeValuePerFace<const Rd> m_Nlr;
PUGS_INLINE PUGS_INLINE
void void
_compute_ll()
{
if constexpr (Dimension == 1) {
static_assert(Dimension != 1, "ll does not make sense in 1d");
} else {
const auto& nlr = this->Nlr();
FaceValue<double> ll{m_mesh.connectivity()};
const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
parallel_for(
m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
const auto& face_nodes = face_to_node_matrix[face_id];
double lenght = 0;
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
lenght += l2Norm(nlr(face_id, i_node));
}
ll[face_id] = lenght;
});
m_ll = ll;
}
}
PUGS_INLINE void
_computeFaceIsobarycenter() _computeFaceIsobarycenter()
{ // Computes vertices isobarycenter { // Computes vertices isobarycenter
if constexpr (Dimension == 1) { if constexpr (Dimension == 1) {
...@@ -78,7 +106,7 @@ class MeshData : public IMeshData ...@@ -78,7 +106,7 @@ class MeshData : public IMeshData
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]]); xj[j] = 0.5 * (xr[cell_nodes[0]] + xr[cell_nodes[1]]);
}); });
m_xj = std::make_shared<CellValue<const Rd>>(xj); m_xj = xj;
} else { } else {
const NodeValue<const Rd>& xr = m_mesh.xr(); const NodeValue<const Rd>& xr = m_mesh.xr();
...@@ -95,7 +123,7 @@ class MeshData : public IMeshData ...@@ -95,7 +123,7 @@ class MeshData : public IMeshData
} }
xj[j] = inv_cell_nb_nodes[j] * X; xj[j] = inv_cell_nb_nodes[j] * X;
}); });
m_xj = std::make_shared<CellValue<const Rd>>(xj); m_xj = xj;
} }
} }
...@@ -119,7 +147,7 @@ class MeshData : public IMeshData ...@@ -119,7 +147,7 @@ class MeshData : public IMeshData
} }
Vj[j] = inv_Dimension * sum_cjr_xr; Vj[j] = inv_Dimension * sum_cjr_xr;
}); });
m_Vj = std::make_shared<CellValue<const double>>(Vj); m_Vj = Vj;
} }
PUGS_INLINE PUGS_INLINE
...@@ -147,7 +175,7 @@ class MeshData : public IMeshData ...@@ -147,7 +175,7 @@ class MeshData : public IMeshData
Nlr(l, 0) = Nr; Nlr(l, 0) = Nr;
Nlr(l, 1) = Nr; Nlr(l, 1) = Nr;
}); });
m_Nlr = std::make_shared<NodeValuePerFace<const Rd>>(Nlr); m_Nlr = Nlr;
} else { } else {
const NodeValue<const Rd>& xr = m_mesh.xr(); const NodeValue<const Rd>& xr = m_mesh.xr();
...@@ -174,7 +202,7 @@ class MeshData : public IMeshData ...@@ -174,7 +202,7 @@ class MeshData : public IMeshData
Nlr(l, r) = Nr; Nlr(l, r) = Nr;
} }
}); });
m_Nlr = std::make_shared<NodeValuePerFace<const Rd>>(Nlr); m_Nlr = Nlr;
} }
} }
...@@ -189,7 +217,7 @@ class MeshData : public IMeshData ...@@ -189,7 +217,7 @@ class MeshData : public IMeshData
Cjr(j, 0) = -1; Cjr(j, 0) = -1;
Cjr(j, 1) = 1; Cjr(j, 1) = 1;
}); });
m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr); m_Cjr = Cjr;
} else if constexpr (Dimension == 2) { } else if constexpr (Dimension == 2) {
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();
...@@ -205,7 +233,7 @@ class MeshData : public IMeshData ...@@ -205,7 +233,7 @@ class MeshData : public IMeshData
Cjr(j, R) = Rd{half_xrp_xrm[1], -half_xrp_xrm[0]}; Cjr(j, R) = Rd{half_xrp_xrm[1], -half_xrp_xrm[0]};
} }
}); });
m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr); m_Cjr = Cjr;
} else if (Dimension == 3) { } else if (Dimension == 3) {
auto Nlr = this->Nlr(); auto Nlr = this->Nlr();
...@@ -252,7 +280,7 @@ class MeshData : public IMeshData ...@@ -252,7 +280,7 @@ class MeshData : public IMeshData
} }
}); });
m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr); m_Cjr = Cjr;
} }
} }
...@@ -265,13 +293,13 @@ class MeshData : public IMeshData ...@@ -265,13 +293,13 @@ class MeshData : public IMeshData
NodeValuePerCell<double> ljr(m_mesh.connectivity()); NodeValuePerCell<double> ljr(m_mesh.connectivity());
parallel_for( parallel_for(
ljr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = 1; }); ljr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = 1; });
m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr); m_ljr = ljr;
} else { } else {
NodeValuePerCell<double> ljr(m_mesh.connectivity()); NodeValuePerCell<double> ljr(m_mesh.connectivity());
parallel_for( parallel_for(
Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm(Cjr[jr]); }); Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm(Cjr[jr]); });
m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr); m_ljr = ljr;
} }
} }
...@@ -289,19 +317,18 @@ class MeshData : public IMeshData ...@@ -289,19 +317,18 @@ class MeshData : public IMeshData
NodeValuePerCell<Rd> njr(m_mesh.connectivity()); NodeValuePerCell<Rd> njr(m_mesh.connectivity());
parallel_for( parallel_for(
Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / ljr[jr]) * Cjr[jr]; }); Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / ljr[jr]) * Cjr[jr]; });
m_njr = std::make_shared<NodeValuePerCell<const Rd>>(njr); m_njr = njr;
} }
} }
void void
_checkCellVolume() _checkCellVolume()
{ {
Assert(m_Vj); Assert(m_Vj.isBuilt());
auto& Vj = *m_Vj;
bool is_valid = [&] { bool is_valid = [&] {
for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) { for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
if (Vj[j] <= 0) { if (m_Vj[j] <= 0) {
return false; return false;
} }
} }
...@@ -321,44 +348,54 @@ class MeshData : public IMeshData ...@@ -321,44 +348,54 @@ class MeshData : public IMeshData
return m_mesh; return m_mesh;
} }
PUGS_INLINE
FaceValue<const double>
ll()
{
if (not m_ll.isBuilt()) {
this->_compute_ll();
}
return m_ll;
}
PUGS_INLINE PUGS_INLINE
NodeValuePerFace<const Rd> NodeValuePerFace<const Rd>
Nlr() Nlr()
{ {
if (not m_Nlr) { if (not m_Nlr.isBuilt()) {
this->_computeNlr(); this->_computeNlr();
} }
return *m_Nlr; return m_Nlr;
} }
PUGS_INLINE PUGS_INLINE
NodeValuePerCell<const Rd> NodeValuePerCell<const Rd>
Cjr() Cjr()
{ {
if (not m_Cjr) { if (not m_Cjr.isBuilt()) {
this->_computeCjr(); this->_computeCjr();
} }
return *m_Cjr; return m_Cjr;
} }
PUGS_INLINE PUGS_INLINE
NodeValuePerCell<const double> NodeValuePerCell<const double>
ljr() ljr()
{ {
if (not m_ljr) { if (not m_ljr.isBuilt()) {
this->_compute_ljr(); this->_compute_ljr();
} }
return *m_ljr; return m_ljr;
} }
PUGS_INLINE PUGS_INLINE
NodeValuePerCell<const Rd> NodeValuePerCell<const Rd>
njr() njr()
{ {
if (not m_njr) { if (not m_njr.isBuilt()) {
this->_compute_njr(); this->_compute_njr();
} }
return *m_njr; return m_njr;
} }
PUGS_INLINE PUGS_INLINE
...@@ -375,21 +412,21 @@ class MeshData : public IMeshData ...@@ -375,21 +412,21 @@ class MeshData : public IMeshData
CellValue<const Rd> CellValue<const Rd>
xj() xj()
{ {
if (not m_xj) { if (not m_xj.isBuilt()) {
this->_computeCellIsobarycenter(); this->_computeCellIsobarycenter();
} }
return *m_xj; return m_xj;
} }
PUGS_INLINE PUGS_INLINE
CellValue<const double> CellValue<const double>
Vj() Vj()
{ {
if (not m_Vj) { if (not m_Vj.isBuilt()) {
this->_computeCellVolume(); this->_computeCellVolume();
this->_checkCellVolume(); this->_checkCellVolume();
} }
return *m_Vj; return m_Vj;
} }
private: private:
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment