diff --git a/src/language/modules/MeshModule.hpp b/src/language/modules/MeshModule.hpp
index 64e6194946b6f6f9ceed019f513730d6fd3de8c0..f2c3d3d6360e5fd378d128d3b6b57e0de5a04509 100644
--- a/src/language/modules/MeshModule.hpp
+++ b/src/language/modules/MeshModule.hpp
@@ -5,7 +5,7 @@
#include <language/utils/ASTNodeDataTypeTraits.hpp>
#include <utils/PugsMacros.hpp>
-struct IMesh;
+class IMesh;
template <>
inline ASTNodeDataType ast_node_data_type_from<std::shared_ptr<const IMesh>> = {ASTNodeDataType::type_id_t, "mesh"};
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index f8cd50f30bed85e3f905aa6efc877b2612e669e0..01eccadce6fcedb1cfaea71d87e959781190cc25 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -22,6 +22,7 @@ class MeshData : public IMeshData
{
public:
static_assert(Dimension > 0, "dimension must be strictly positive");
+ static_assert((Dimension <= 3), "only 1d, 2d and 3d are implemented");
using MeshType = Mesh<Connectivity<Dimension>>;
@@ -36,10 +37,11 @@ class MeshData : public IMeshData
std::shared_ptr<NodeValuePerCell<const Rd>> m_njr;
std::shared_ptr<CellValue<const Rd>> m_xj;
std::shared_ptr<CellValue<const double>> m_Vj;
+ std::shared_ptr<NodeValuePerFace<const Rd>> m_Nlr;
PUGS_INLINE
void
- _updateCenter()
+ _computeIsobarycenter()
{ // Computes vertices isobarycenter
if constexpr (Dimension == 1) {
const NodeValue<const Rd>& xr = m_mesh.xr();
@@ -75,11 +77,13 @@ class MeshData : public IMeshData
PUGS_INLINE
void
- _updateVolume()
+ _computeCellVolume()
{
const NodeValue<const Rd>& xr = m_mesh.xr();
const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+ auto Cjr = this->Cjr();
+
CellValue<double> Vj(m_mesh.connectivity());
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
@@ -87,7 +91,7 @@ class MeshData : public IMeshData
const auto& cell_nodes = cell_to_node_matrix[j];
for (size_t R = 0; R < cell_nodes.size(); ++R) {
- sum_cjr_xr += (xr[cell_nodes[R]], (*m_Cjr)(j, R));
+ sum_cjr_xr += (xr[cell_nodes[R]], Cjr(j, R));
}
Vj[j] = inv_Dimension * sum_cjr_xr;
});
@@ -96,43 +100,13 @@ class MeshData : public IMeshData
PUGS_INLINE
void
- _updateCjr()
+ _computeNlr()
{
if constexpr (Dimension == 1) {
- // Cjr/njr/ljr are constant overtime
+ static_assert(Dimension != 1, "Nlr does not make sense in 1d");
} else if constexpr (Dimension == 2) {
- const NodeValue<const Rd>& xr = m_mesh.xr();
- const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
-
- {
- NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
- for (size_t R = 0; R < cell_nodes.size(); ++R) {
- int Rp1 = (R + 1) % cell_nodes.size();
- int Rm1 = (R + cell_nodes.size() - 1) % cell_nodes.size();
- Rd half_xrp_xrm = 0.5 * (xr[cell_nodes[Rp1]] - xr[cell_nodes[Rm1]]);
- Cjr(j, R) = Rd{half_xrp_xrm[1], -half_xrp_xrm[0]};
- }
- });
- m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
- }
-
- {
- NodeValuePerCell<double> ljr(m_mesh.connectivity());
- parallel_for(
- m_Cjr->numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm((*m_Cjr)[jr]); });
- m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
- }
-
- {
- NodeValuePerCell<Rd> njr(m_mesh.connectivity());
- parallel_for(
- m_Cjr->numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / (*m_ljr)[jr]) * (*m_Cjr)[jr]; });
- m_njr = std::make_shared<NodeValuePerCell<const Rd>>(njr);
- }
- } else if (Dimension == 3) {
+ throw NotImplementedError("Nlr are not yet computed in 2d");
+ } else {
const NodeValue<const Rd>& xr = m_mesh.xr();
NodeValuePerFace<Rd> Nlr(m_mesh.connectivity());
@@ -158,78 +132,133 @@ class MeshData : public IMeshData
Nlr(l, r) = Nr;
}
});
+ m_Nlr = std::make_shared<NodeValuePerFace<const Rd>>(Nlr);
+ }
+ }
+ PUGS_INLINE
+ void
+ _computeCjr()
+ {
+ if constexpr (Dimension == 1) {
+ NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ Cjr(j, 0) = -1;
+ Cjr(j, 1) = 1;
+ });
+ m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
+ } else if constexpr (Dimension == 2) {
+ const NodeValue<const Rd>& xr = m_mesh.xr();
const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
- const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
+ NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ int Rp1 = (R + 1) % cell_nodes.size();
+ int Rm1 = (R + cell_nodes.size() - 1) % cell_nodes.size();
+ Rd half_xrp_xrm = 0.5 * (xr[cell_nodes[Rp1]] - xr[cell_nodes[Rm1]]);
+ Cjr(j, R) = Rd{half_xrp_xrm[1], -half_xrp_xrm[0]};
+ }
+ });
+ m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
+ } else if (Dimension == 3) {
+ auto Nlr = this->Nlr();
+ const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+ const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
- {
- NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
- parallel_for(
- Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Cjr[jr] = zero; });
+ NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
+ parallel_for(
+ Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Cjr[jr] = zero; });
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
- const auto& cell_faces = cell_to_face_matrix[j];
- const auto& face_is_reversed = cell_face_is_reversed.itemValues(j);
+ const auto& cell_faces = cell_to_face_matrix[j];
+ const auto& face_is_reversed = cell_face_is_reversed.itemValues(j);
- for (size_t L = 0; L < cell_faces.size(); ++L) {
- const FaceId& l = cell_faces[L];
- const auto& face_nodes = face_to_node_matrix[l];
+ for (size_t L = 0; L < cell_faces.size(); ++L) {
+ const FaceId& l = cell_faces[L];
+ const auto& face_nodes = face_to_node_matrix[l];
- auto local_node_number_in_cell = [&](NodeId node_number) {
- for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
- if (node_number == cell_nodes[i_node]) {
- return i_node;
- }
+ auto local_node_number_in_cell = [&](NodeId node_number) {
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ if (node_number == cell_nodes[i_node]) {
+ return i_node;
}
- return std::numeric_limits<size_t>::max();
- };
+ }
+ return std::numeric_limits<size_t>::max();
+ };
- if (face_is_reversed[L]) {
- for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
- const size_t R = local_node_number_in_cell(face_nodes[rl]);
- Cjr(j, R) -= Nlr(l, rl);
- }
- } else {
- for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
- const size_t R = local_node_number_in_cell(face_nodes[rl]);
- Cjr(j, R) += Nlr(l, rl);
- }
+ if (face_is_reversed[L]) {
+ for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
+ Cjr(j, R) -= Nlr(l, rl);
+ }
+ } else {
+ for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
+ Cjr(j, R) += Nlr(l, rl);
}
}
- });
+ }
+ });
- m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
- }
+ m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
+ }
+ }
- {
- NodeValuePerCell<double> ljr(m_mesh.connectivity());
- parallel_for(
- m_Cjr->numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm((*m_Cjr)[jr]); });
- m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
- }
+ PUGS_INLINE
+ void
+ _compute_ljr()
+ {
+ auto Cjr = this->Cjr();
+ if constexpr (Dimension == 1) {
+ NodeValuePerCell<double> ljr(m_mesh.connectivity());
+ parallel_for(
+ ljr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = 1; });
+ m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
- {
- NodeValuePerCell<Rd> njr(m_mesh.connectivity());
- parallel_for(
- m_Cjr->numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / (*m_ljr)[jr]) * (*m_Cjr)[jr]; });
- m_njr = std::make_shared<NodeValuePerCell<const Rd>>(njr);
- }
+ } else {
+ NodeValuePerCell<double> ljr(m_mesh.connectivity());
+ parallel_for(
+ Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm(Cjr[jr]); });
+ m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
}
- static_assert((Dimension <= 3), "only 1d, 2d and 3d are implemented");
}
+ PUGS_INLINE
void
- _checkCellVolume() const
+ _compute_njr()
{
+ auto Cjr = this->Cjr();
+ if constexpr (Dimension == 1) {
+ // in 1d njr=Cjr (here performs a shallow copy)
+ m_njr = m_Cjr;
+ } else {
+ auto ljr = this->ljr();
+
+ NodeValuePerCell<Rd> njr(m_mesh.connectivity());
+ parallel_for(
+ Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / ljr[jr]) * Cjr[jr]; });
+ m_njr = std::make_shared<NodeValuePerCell<const Rd>>(njr);
+ }
+ }
+
+ void
+ _checkCellVolume()
+ {
+ auto Vj = this->Vj();
+
bool is_valid = [&] {
for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
- if ((*m_Vj)[j] <= 0) {
+ if (Vj[j] <= 0) {
return false;
}
}
@@ -250,76 +279,73 @@ class MeshData : public IMeshData
}
PUGS_INLINE
- const NodeValuePerCell<const Rd>
- Cjr() const
+ NodeValuePerFace<const Rd>
+ Nlr()
+ {
+ if (not m_Nlr) {
+ this->_computeNlr();
+ }
+ return *m_Nlr;
+ }
+
+ PUGS_INLINE
+ NodeValuePerCell<const Rd>
+ Cjr()
{
+ if (not m_Cjr) {
+ this->_computeCjr();
+ }
return *m_Cjr;
}
PUGS_INLINE
- const NodeValuePerCell<const double>
- ljr() const
+ NodeValuePerCell<const double>
+ ljr()
{
+ if (not m_ljr) {
+ this->_compute_ljr();
+ }
return *m_ljr;
}
PUGS_INLINE
- const NodeValuePerCell<const Rd>
- njr() const
+ NodeValuePerCell<const Rd>
+ njr()
{
+ if (not m_njr) {
+ this->_compute_njr();
+ }
return *m_njr;
}
PUGS_INLINE
- const CellValue<const Rd>
+ CellValue<const Rd>
xj()
{
if (not m_xj) {
- this->_updateCenter();
+ this->_computeIsobarycenter();
}
return *m_xj;
}
PUGS_INLINE
- const CellValue<const double>
- Vj() const
+ CellValue<const double>
+ Vj()
{
+ if (not m_Vj) {
+ this->_computeCellVolume();
+ this->_checkCellVolume();
+ }
return *m_Vj;
}
void
updateAllData()
{
- this->_updateCjr();
- // this->_updateCenter();
- this->_updateVolume();
- this->_checkCellVolume();
+ ;
}
- MeshData(const MeshType& mesh) : m_mesh(mesh)
- {
- if constexpr (Dimension == 1) {
- // in 1d Cjr are computed once for all
- {
- NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- Cjr(j, 0) = -1;
- Cjr(j, 1) = 1;
- });
- m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
- }
- // in 1d njr=Cjr (here performs a shallow copy)
- m_njr = m_Cjr;
- {
- NodeValuePerCell<double> ljr(m_mesh.connectivity());
- parallel_for(
- ljr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = 1; });
- m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
- }
- }
- this->updateAllData();
- }
+ MeshData(const MeshType& mesh) : m_mesh(mesh) {}
MeshData() = delete;
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index c531154bd5db52a3f102dfc8e0dd455157e01351..971c8fd2458c35d55013a084a4919d7dce026ad6 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -279,7 +279,7 @@ class AcousticSolver
CellValue<double>& pj = unknowns.pj();
CellValue<double>& cj = unknowns.cj();
- const MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+ MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
const CellValue<const double> Vj = mesh_data.Vj();
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();