From 972ab2f158f193ecdc14247c01ea434e6da66e6e Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 7 Jul 2020 00:41:16 +0200
Subject: [PATCH] Compute mesh data attributes on demand
Separate the notion of Nlr (Eucclhyd normals) of the calculation of
Cjr in 3D. It remain to implement its calculation in 2D to finish this
part.
This should allow a trivial implementation of the Eucclhyd scheme and
simplify coding off pressure boundary conditions.
---
src/language/modules/MeshModule.hpp | 2 +-
src/mesh/MeshData.hpp | 272 +++++++++++++++-------------
src/scheme/AcousticSolver.hpp | 2 +-
3 files changed, 151 insertions(+), 125 deletions(-)
diff --git a/src/language/modules/MeshModule.hpp b/src/language/modules/MeshModule.hpp
index 64e619494..f2c3d3d63 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 f8cd50f30..01eccadce 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 c531154bd..971c8fd24 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();
--
GitLab