diff --git a/CMakeLists.txt b/CMakeLists.txt
index 127b72af3dbdccab17393fb27568ba7b3bdc23db..73cda5739485b07c04c44c93ae6f4e9533d0f055 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -425,6 +425,7 @@ target_link_libraries(
PugsLanguage
PugsLanguageAST
PugsLanguageModules
+ PugsLanguageAlgorithms
PugsMesh
PugsUtils
PugsLanguageUtils
diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp
index c931385b214c91ca2807db4f07fca027ebd79b25..bc822668a7f63c4d5d7cc4c45f4d1b1d70084fc0 100644
--- a/src/algebra/BiCGStab.hpp
+++ b/src/algebra/BiCGStab.hpp
@@ -5,14 +5,19 @@
#include <iomanip>
#include <iostream>
+#include <rang.hpp>
+
+template <bool verbose = true>
struct BiCGStab
{
template <typename VectorType, typename MatrixType>
BiCGStab(const VectorType& b, const MatrixType& A, VectorType& x, const size_t max_iter, const double epsilon = 1e-6)
{
- std::cout << "- bi-conjugate gradient stabilized\n";
- std::cout << " epsilon = " << epsilon << '\n';
- std::cout << " maximum number of iterations: " << max_iter << '\n';
+ if constexpr (verbose) {
+ std::cout << "- bi-conjugate gradient stabilized\n";
+ std::cout << " epsilon = " << epsilon << '\n';
+ std::cout << " maximum number of iterations: " << max_iter << '\n';
+ }
VectorType r_k_1{b.size()};
@@ -33,10 +38,14 @@ struct BiCGStab
VectorType r_k{x.size()};
- std::cout << " initial residu: " << resid0 << '\n';
+ if constexpr (verbose) {
+ std::cout << " initial residu: " << resid0 << '\n';
+ }
for (size_t i = 1; i <= max_iter; ++i) {
- std::cout << " - iteration: " << std::setw(6) << i << "\tresidu: " << residu / resid0
- << "\tabsolute: " << residu << '\n';
+ if constexpr (verbose) {
+ std::cout << " - iteration: " << std::setw(6) << i << "\tresidu: " << residu / resid0
+ << "\tabsolute: " << residu << '\n';
+ }
Ap_k = A * p_k;
@@ -64,14 +73,14 @@ struct BiCGStab
}
if (residu / resid0 > epsilon) {
- std::cout << " bi_conjugate gradient stabilized: *NOT CONVERGED*\n";
+ std::cout << " bi-conjugate gradient stabilized: " << rang::fgB::red << "*NOT CONVERGED*" << rang::style::reset
+ << '\n';
+ ;
std::cout << " - epsilon: " << epsilon << '\n';
std::cout << " - relative residu : " << residu / resid0 << '\n';
std::cout << " - absolute residu : " << residu << '\n';
}
}
-
- std::cout << "- bi-conjugate gradient stabilized: done\n";
}
};
diff --git a/src/algebra/PCG.hpp b/src/algebra/PCG.hpp
index 74c17ee40971b96c512c196d630f1edfadfa6ddc..36a4635b313d2bba13d97e1b7a4ea1882e2a66a2 100644
--- a/src/algebra/PCG.hpp
+++ b/src/algebra/PCG.hpp
@@ -4,6 +4,9 @@
#include <iomanip>
#include <iostream>
+#include <rang.hpp>
+
+template <bool verbose = true>
struct PCG
{
template <typename VectorType, typename MatrixType, typename PreconditionerType>
@@ -14,14 +17,16 @@ struct PCG
const size_t maxiter,
const double epsilon = 1e-6)
{
- std::cout << "- conjugate gradient\n";
- std::cout << " epsilon = " << epsilon << '\n';
- std::cout << " maximum number of iterations: " << maxiter << '\n';
+ if constexpr (verbose) {
+ std::cout << "- conjugate gradient\n";
+ std::cout << " epsilon = " << epsilon << '\n';
+ std::cout << " maximum number of iterations: " << maxiter << '\n';
+ }
VectorType h{f.size()};
VectorType b = copy(f);
- if (true) {
+ if constexpr (verbose) {
h = A * x;
h -= f;
std::cout << "- initial *real* residu : " << (h, h) << '\n';
@@ -69,10 +74,14 @@ struct PCG
if ((i == 1) && (gcg != 0)) {
relativeEpsilon = epsilon * gcg;
gcg0 = gcg;
- std::cout << " initial residu: " << gcg << '\n';
+ if constexpr (verbose) {
+ std::cout << " initial residu: " << gcg << '\n';
+ }
+ }
+ if constexpr (verbose) {
+ std::cout << " - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0;
+ std::cout << "\tabsolute: " << gcg << '\n';
}
- std::cout << " - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0;
- std::cout << "\tabsolute: " << gcg << '\n';
if (gcg < relativeEpsilon) {
break;
@@ -85,7 +94,7 @@ struct PCG
}
if (gcg > relativeEpsilon) {
- std::cout << " conjugate gradient: *NOT CONVERGED*\n";
+ std::cout << " conjugate gradient: " << rang::fgB::red << "*NOT CONVERGED*" << rang::style::reset << '\n';
std::cout << " - epsilon: " << epsilon << '\n';
std::cout << " - relative residu : " << gcg / gcg0 << '\n';
std::cout << " - absolute residu : " << gcg << '\n';
diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp
index 3533a4017011b0c5c1f5d7f04ae2b893fa8a7d0b..100000366d3bd0e5a072084864e82cdb4b73ba29 100644
--- a/src/algebra/SparseMatrixDescriptor.hpp
+++ b/src/algebra/SparseMatrixDescriptor.hpp
@@ -28,14 +28,16 @@ class SparseMatrixDescriptor
return m_id_value_map.size();
}
- const DataType& operator[](const IndexType& j) const
+ const DataType&
+ operator[](const IndexType& j) const
{
auto i_value = m_id_value_map.find(j);
Assert(i_value != m_id_value_map.end());
return i_value->second;
}
- DataType& operator[](const IndexType& j)
+ DataType&
+ operator[](const IndexType& j)
{
auto i_value = m_id_value_map.find(j);
if (i_value != m_id_value_map.end()) {
@@ -51,7 +53,7 @@ class SparseMatrixDescriptor
operator<<(std::ostream& os, const SparseRowDescriptor& row)
{
for (auto [j, value] : row.m_id_value_map) {
- os << ' ' << j << ':' << value;
+ os << ' ' << static_cast<size_t>(j) << ':' << value;
}
return os;
}
@@ -73,24 +75,30 @@ class SparseMatrixDescriptor
SparseRowDescriptor&
row(const IndexType i)
{
+ Assert(i < m_row_array.size());
return m_row_array[i];
}
const SparseRowDescriptor&
row(const IndexType i) const
{
+ Assert(i < m_row_array.size());
return m_row_array[i];
}
DataType&
operator()(const IndexType& i, const IndexType& j)
{
+ Assert(i < m_row_array.size());
+ Assert(j < m_row_array.size());
return m_row_array[i][j];
}
const DataType&
operator()(const IndexType& i, const IndexType& j) const
{
+ Assert(i < m_row_array.size());
+ Assert(j < m_row_array.size());
const auto& r = m_row_array[i]; // split to ensure const-ness of call
return r[j];
}
@@ -99,7 +107,7 @@ class SparseMatrixDescriptor
operator<<(std::ostream& os, const SparseMatrixDescriptor& M)
{
for (IndexType i = 0; i < M.m_row_array.size(); ++i) {
- os << i << " |" << M.m_row_array[i] << '\n';
+ os << static_cast<size_t>(i) << " |" << M.m_row_array[i] << '\n';
}
return os;
}
diff --git a/src/language/CMakeLists.txt b/src/language/CMakeLists.txt
index b9c6094e4403cb1d2289d8cfcae58de9e2e3c4e5..3565082161700da019088f7bfca70b0d9dfe6fd8 100644
--- a/src/language/CMakeLists.txt
+++ b/src/language/CMakeLists.txt
@@ -1,5 +1,6 @@
# ------------------- Source files --------------------
+add_subdirectory(algorithms)
add_subdirectory(ast)
add_subdirectory(modules)
add_subdirectory(node_processor)
@@ -12,6 +13,7 @@ add_library(
# Additional dependencies
add_dependencies(PugsLanguage
PugsLanguage
+ PugsLanguageAlgorithms
PugsLanguageAST
PugsLanguageModules
PugsLanguageUtils
diff --git a/src/language/algorithms/AcousticSolverAlgorithm.cpp b/src/language/algorithms/AcousticSolverAlgorithm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9811c90f0b0f2c52d5f72369efea4a4398f007e3
--- /dev/null
+++ b/src/language/algorithms/AcousticSolverAlgorithm.cpp
@@ -0,0 +1,245 @@
+#include <language/algorithms/AcousticSolverAlgorithm.hpp>
+
+#include <language/utils/InterpolateItemValue.hpp>
+#include <output/VTKWriter.hpp>
+#include <scheme/AcousticSolver.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+
+template <size_t Dimension>
+AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
+ const AcousticSolverType& solver_type,
+ std::shared_ptr<const IMesh> i_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& rho_id,
+ const FunctionSymbolId& u_id,
+ const FunctionSymbolId& p_id)
+{
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<Dimension>;
+ using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>;
+
+ std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+ typename AcousticSolver<MeshType>::BoundaryConditionList bc_list;
+ {
+ constexpr ItemType FaceType = [] {
+ if constexpr (Dimension > 1) {
+ return ItemType::face;
+ } else {
+ return ItemType::node;
+ }
+ }();
+
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ using SymmetryBoundaryCondition = typename AcousticSolver<MeshType>::SymmetryBoundaryCondition;
+
+ const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+ dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+ bc_list.push_back(
+ SymmetryBoundaryCondition{MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list)});
+ }
+ }
+ is_valid_boundary_condition = true;
+ break;
+ }
+
+ case IBoundaryConditionDescriptor::Type::dirichlet: {
+ const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+ dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+ if (dirichlet_bc_descriptor.name() == "velocity") {
+ using VelocityBoundaryCondition = typename AcousticSolver<MeshType>::VelocityBoundaryCondition;
+
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+ MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
+
+ const FunctionSymbolId velocity_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ const auto& node_list = mesh_node_boundary.nodeList();
+
+ Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, mesh->xr(), node_list);
+
+ bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
+ }
+ }
+ } else if (dirichlet_bc_descriptor.name() == "pressure") {
+ using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition;
+
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+ const auto& face_list = ref_face_list.list();
+
+ const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ Array<const double> face_values = [&] {
+ if constexpr (Dimension == 1) {
+ return InterpolateItemValue<double(
+ TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh->xr(), face_list);
+ } else {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ return InterpolateItemValue<double(
+ TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(), face_list);
+ }
+ }();
+
+ bc_list.push_back(PressureBoundaryCondition{face_list, face_values});
+ }
+ }
+ } else {
+ is_valid_boundary_condition = false;
+ }
+ break;
+ }
+ default: {
+ is_valid_boundary_condition = false;
+ }
+ }
+ if (not is_valid_boundary_condition) {
+ std::ostringstream error_msg;
+ error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver";
+ throw NormalError(error_msg.str());
+ }
+ }
+ }
+
+ UnknownsType unknowns(*mesh);
+
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ unknowns.rhoj() =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id, mesh_data.xj());
+
+ unknowns.pj() =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj());
+
+ unknowns.uj() =
+ InterpolateItemValue<TinyVector<Dimension>(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id,
+ mesh_data
+ .xj());
+ }
+ unknowns.gammaj().fill(1.4);
+
+ AcousticSolver acoustic_solver(mesh, bc_list, solver_type);
+
+ const double tmax = 0.2;
+ double t = 0;
+
+ int itermax = std::numeric_limits<int>::max();
+ int iteration = 0;
+
+ CellValue<double>& rhoj = unknowns.rhoj();
+ CellValue<double>& ej = unknowns.ej();
+ CellValue<double>& pj = unknowns.pj();
+ CellValue<double>& gammaj = unknowns.gammaj();
+ CellValue<double>& cj = unknowns.cj();
+ CellValue<TinyVector<Dimension>>& uj = unknowns.uj();
+ CellValue<double>& Ej = unknowns.Ej();
+ CellValue<double>& mj = unknowns.mj();
+ CellValue<double>& inv_mj = unknowns.invMj();
+
+ BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+ block_eos.updateEandCFromRhoP();
+
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ const CellValue<const double> Vj = mesh_data.Vj();
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; });
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
+ }
+
+ VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
+
+ while ((t < tmax) and (iteration < itermax)) {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ vtk_writer.write(mesh,
+ {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+ NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+ NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+ NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
+ t);
+
+ double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
+ if (t + dt > tmax) {
+ dt = tmax - t;
+ }
+
+ std::cout.setf(std::cout.scientific);
+ std::cout << "iteration " << rang::fg::cyan << std::setw(4) << iteration << rang::style::reset
+ << " time=" << rang::fg::green << t << rang::style::reset << " dt=" << rang::fgB::blue << dt
+ << rang::style::reset << '\n';
+
+ mesh = acoustic_solver.computeNextStep(dt, unknowns);
+
+ block_eos.updatePandCFromRhoE();
+
+ t += dt;
+ ++iteration;
+ }
+ std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after "
+ << rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations"
+ << rang::style::reset << '\n';
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ vtk_writer.write(mesh,
+ {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+ NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+ NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+ NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
+ t, true); // forces last output
+ }
+}
+
+template AcousticSolverAlgorithm<1>::AcousticSolverAlgorithm(
+ const AcousticSolverType&,
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&);
+
+template AcousticSolverAlgorithm<2>::AcousticSolverAlgorithm(
+ const AcousticSolverType&,
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&);
+
+template AcousticSolverAlgorithm<3>::AcousticSolverAlgorithm(
+ const AcousticSolverType&,
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&);
diff --git a/src/language/algorithms/AcousticSolverAlgorithm.hpp b/src/language/algorithms/AcousticSolverAlgorithm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a6973480c4cd913c6463d39012f6eea906540003
--- /dev/null
+++ b/src/language/algorithms/AcousticSolverAlgorithm.hpp
@@ -0,0 +1,20 @@
+#ifndef ACOUSTIC_SOLVER_ALGORITHM_HPP
+#define ACOUSTIC_SOLVER_ALGORITHM_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+#include <scheme/AcousticSolverType.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+template <size_t Dimension>
+struct AcousticSolverAlgorithm
+{
+ AcousticSolverAlgorithm(const AcousticSolverType& solver_type,
+ std::shared_ptr<const IMesh> i_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& rho_id,
+ const FunctionSymbolId& u_id,
+ const FunctionSymbolId& p_id);
+};
+
+#endif // ACOUSTIC_SOLVER_ALGORITHM_HPP
diff --git a/src/language/algorithms/CMakeLists.txt b/src/language/algorithms/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a8599f575058185b5ceeed0814b242ae6ec70b0a
--- /dev/null
+++ b/src/language/algorithms/CMakeLists.txt
@@ -0,0 +1,10 @@
+# ------------------- Source files --------------------
+
+add_library(PugsLanguageAlgorithms
+ AcousticSolverAlgorithm.cpp
+)
+
+
+add_dependencies(PugsLanguageAlgorithms
+ PugsUtils
+ PugsMesh)
diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt
index 703f19857cfb99c0053aecf1af0984dab371d72e..f8cf2f0b445882a672d54e60dbc75858bbb246c7 100644
--- a/src/language/modules/CMakeLists.txt
+++ b/src/language/modules/CMakeLists.txt
@@ -11,5 +11,6 @@ add_library(PugsLanguageModules
add_dependencies(PugsLanguageModules
+ PugsLanguageAlgorithms
PugsUtils
PugsMesh)
diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index b9a65729637128081a7206ad6fa5d98ad95e9978..89e678c3fb4095350842b19a9145502cd7f558a8 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -9,7 +9,7 @@
#include <language/utils/TypeDescriptor.hpp>
#include <mesh/CartesianMeshBuilder.hpp>
#include <mesh/Connectivity.hpp>
-#include <mesh/DiamondDualMeshBuilder.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
#include <mesh/GmshReader.hpp>
#include <mesh/Mesh.hpp>
#include <utils/Exceptions.hpp>
@@ -184,9 +184,30 @@ MeshModule::MeshModule()
std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IMesh>(
const std::shared_ptr<const IMesh>&)>>(
- [](const std::shared_ptr<const IMesh>& p_mesh) -> std::shared_ptr<const IMesh> {
- DiamondDualMeshBuilder builder{p_mesh};
- return builder.mesh();
+ [](const std::shared_ptr<const IMesh>& i_mesh) -> std::shared_ptr<const IMesh> {
+ switch (i_mesh->dimension()) {
+ case 1: {
+ using MeshType = Mesh<Connectivity<1>>;
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ return DiamondDualMeshManager::instance().getDiamondDualMesh(p_mesh);
+ }
+ case 2: {
+ using MeshType = Mesh<Connectivity<2>>;
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ return DiamondDualMeshManager::instance().getDiamondDualMesh(p_mesh);
+ }
+ case 3: {
+ using MeshType = Mesh<Connectivity<3>>;
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ return DiamondDualMeshManager::instance().getDiamondDualMesh(p_mesh);
+ }
+ default: {
+ throw UnexpectedError("invalid dimension");
+ }
+ }
}
));
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index c04f29ad41df122df252cdc9e7009b49241e8ed2..542690a4ffcfb73ba174f0014b358548208d830f 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -1,215 +1,21 @@
#include <language/modules/SchemeModule.hpp>
+#include <language/algorithms/AcousticSolverAlgorithm.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/TypeDescriptor.hpp>
#include <mesh/Mesh.hpp>
-#include <scheme/AcousticSolver.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/FourierBoundaryConditionDescriptor.hpp>
+#include <scheme/FreeBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryDescriptor.hpp>
#include <scheme/NamedBoundaryDescriptor.hpp>
+#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/NumberedBoundaryDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <memory>
-/////////// TEMPORARY
-
-#include <language/utils/PugsFunctionAdapter.hpp>
-#include <output/VTKWriter.hpp>
-
-template <typename T>
-class InterpolateItemValue;
-template <typename OutputType, typename InputType>
-class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
-{
- static constexpr size_t Dimension = OutputType::Dimension;
- using Adapter = PugsFunctionAdapter<OutputType(InputType)>;
-
- public:
- template <ItemType item_type>
- static inline ItemValue<OutputType, item_type>
- interpolate(const FunctionSymbolId& function_symbol_id, const ItemValue<const InputType, item_type>& position)
- {
- auto& expression = Adapter::getFunctionExpression(function_symbol_id);
- auto convert_result = Adapter::getResultConverter(expression.m_data_type);
-
- Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
-
- using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
- Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
- const IConnectivity& connectivity = *position.connectivity_ptr();
-
- ItemValue<OutputType, item_type> value(connectivity);
- using ItemId = ItemIdT<item_type>;
-
- parallel_for(connectivity.template numberOf<item_type>(), [=, &expression, &tokens](ItemId i) {
- const int32_t t = tokens.acquire();
-
- auto& execution_policy = context_list[t];
-
- Adapter::convertArgs(execution_policy.currentContext(), position[i]);
- auto result = expression.execute(execution_policy);
- value[i] = convert_result(std::move(result));
-
- tokens.release(t);
- });
-
- return value;
- }
-};
-
-template <size_t Dimension>
-struct GlaceScheme
-{
- using ConnectivityType = Connectivity<Dimension>;
- using MeshType = Mesh<ConnectivityType>;
- using MeshDataType = MeshData<Dimension>;
- using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>;
-
- std::shared_ptr<const MeshType> m_mesh;
-
- GlaceScheme(std::shared_ptr<const IMesh> i_mesh,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
- const FunctionSymbolId& rho_id,
- const FunctionSymbolId& u_id,
- const FunctionSymbolId& p_id)
- : m_mesh{std::dynamic_pointer_cast<const MeshType>(i_mesh)}
- {
- std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
-
- std::vector<BoundaryConditionHandler> bc_list;
- {
- constexpr ItemType FaceType = [] {
- if constexpr (Dimension > 1) {
- return ItemType::face;
- } else {
- return ItemType::node;
- }
- }();
-
- for (const auto& bc_descriptor : bc_descriptor_list) {
- switch (bc_descriptor->type()) {
- case IBoundaryConditionDescriptor::Type::symmetry: {
- const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
- dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
- for (size_t i_ref_face_list = 0;
- i_ref_face_list < m_mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
- const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
- const RefId& ref = ref_face_list.refId();
- if (ref == sym_bc_descriptor.boundaryDescriptor()) {
- SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
- new SymmetryBoundaryCondition<MeshType::Dimension>(
- MeshFlatNodeBoundary<MeshType::Dimension>(m_mesh, ref_face_list));
- std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
- bc_list.push_back(BoundaryConditionHandler(bc));
- }
- }
- break;
- }
- default: {
- throw UnexpectedError("Unknown BCDescription\n");
- }
- }
- }
- }
-
- UnknownsType unknowns(*m_mesh);
-
- {
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
- unknowns.rhoj() =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id,
- mesh_data.xj());
-
- unknowns.pj() =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj());
-
- unknowns.uj() = InterpolateItemValue<TinyVector<Dimension>(
- TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id, mesh_data.xj());
- }
- unknowns.gammaj().fill(1.4);
-
- AcousticSolver acoustic_solver(m_mesh, bc_list);
-
- const double tmax = 0.2;
- double t = 0;
-
- int itermax = std::numeric_limits<int>::max();
- int iteration = 0;
-
- CellValue<double>& rhoj = unknowns.rhoj();
- CellValue<double>& ej = unknowns.ej();
- CellValue<double>& pj = unknowns.pj();
- CellValue<double>& gammaj = unknowns.gammaj();
- CellValue<double>& cj = unknowns.cj();
- CellValue<TinyVector<Dimension>>& uj = unknowns.uj();
- CellValue<double>& Ej = unknowns.Ej();
- CellValue<double>& mj = unknowns.mj();
- CellValue<double>& inv_mj = unknowns.invMj();
-
- BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
- block_eos.updateEandCFromRhoP();
-
- {
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
- const CellValue<const double> Vj = mesh_data.Vj();
-
- parallel_for(
- m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
-
- parallel_for(
- m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; });
-
- parallel_for(
- m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
- }
-
- VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
-
- while ((t < tmax) and (iteration < itermax)) {
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
- vtk_writer.write(m_mesh,
- {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
- NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
- t);
- double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
- if (t + dt > tmax) {
- dt = tmax - t;
- }
-
- std::cout.setf(std::cout.scientific);
- std::cout << "iteration " << rang::fg::cyan << std::setw(4) << iteration << rang::style::reset
- << " time=" << rang::fg::green << t << rang::style::reset << " dt=" << rang::fgB::blue << dt
- << rang::style::reset << '\n';
-
- m_mesh = acoustic_solver.computeNextStep(dt, unknowns);
-
- block_eos.updatePandCFromRhoE();
-
- t += dt;
- ++iteration;
- }
- std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after "
- << rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations"
- << rang::style::reset << '\n';
- {
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
- vtk_writer.write(m_mesh,
- {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
- NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
- NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
- t, true); // forces last output
- }
- }
-};
-
SchemeModule::SchemeModule()
{
this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryDescriptor>>);
@@ -247,6 +53,34 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("pressure",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+ const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
+ const FunctionSymbolId&)>>(
+
+ [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+ const FunctionSymbolId& pressure_id)
+ -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+ return std::make_shared<DirichletBoundaryConditionDescriptor>("pressure", boundary,
+ pressure_id);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("velocity",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+ const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
+ const FunctionSymbolId&)>>(
+
+ [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+ const FunctionSymbolId& velocity_id)
+ -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+ return std::make_shared<DirichletBoundaryConditionDescriptor>("velocity", boundary,
+ velocity_id);
+ }
+
+ ));
+
this->_addBuiltinFunction("glace",
std::make_shared<BuiltinFunctionEmbedder<
void(std::shared_ptr<const IMesh>,
@@ -260,15 +94,77 @@ SchemeModule::SchemeModule()
const FunctionSymbolId& p_id) -> void {
switch (p_mesh->dimension()) {
case 1: {
- GlaceScheme<1>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+ AcousticSolverAlgorithm<1>{AcousticSolverType::Glace,
+ p_mesh,
+ bc_descriptor_list,
+ rho_id,
+ u_id,
+ p_id};
+ break;
+ }
+ case 2: {
+ AcousticSolverAlgorithm<2>{AcousticSolverType::Glace,
+ p_mesh,
+ bc_descriptor_list,
+ rho_id,
+ u_id,
+ p_id};
+ break;
+ }
+ case 3: {
+ AcousticSolverAlgorithm<3>{AcousticSolverType::Glace,
+ p_mesh,
+ bc_descriptor_list,
+ rho_id,
+ u_id,
+ p_id};
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+ }
+
+ ));
+
+ this->_addBuiltinFunction("eucclhyd",
+ std::make_shared<BuiltinFunctionEmbedder<
+ void(std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
+
+ [](std::shared_ptr<const IMesh> p_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id,
+ const FunctionSymbolId& p_id) -> void {
+ switch (p_mesh->dimension()) {
+ case 1: {
+ AcousticSolverAlgorithm<1>{AcousticSolverType::Eucclhyd,
+ p_mesh,
+ bc_descriptor_list,
+ rho_id,
+ u_id,
+ p_id};
break;
}
case 2: {
- GlaceScheme<2>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+ AcousticSolverAlgorithm<2>{AcousticSolverType::Eucclhyd,
+ p_mesh,
+ bc_descriptor_list,
+ rho_id,
+ u_id,
+ p_id};
break;
}
case 3: {
- GlaceScheme<3>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+ AcousticSolverAlgorithm<3>{AcousticSolverType::Eucclhyd,
+ p_mesh,
+ bc_descriptor_list,
+ rho_id,
+ u_id,
+ p_id};
break;
}
default: {
diff --git a/src/language/modules/VTKModule.cpp b/src/language/modules/VTKModule.cpp
index f815a78e763b387ee5707d9ac65ec1f038a36b14..ce84141ee866f7ac6fdbb8448db78d0d1b64de88 100644
--- a/src/language/modules/VTKModule.cpp
+++ b/src/language/modules/VTKModule.cpp
@@ -24,7 +24,10 @@ VTKModule::VTKModule()
const std::shared_ptr<const MeshType> mesh =
std::dynamic_pointer_cast<const MeshType>(p_mesh);
- writer.write(mesh, OutputNamedItemValueSet{}, time, true);
+ writer.write(mesh,
+ OutputNamedItemValueSet{
+ NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+ time, true);
break;
}
case 2: {
@@ -32,7 +35,10 @@ VTKModule::VTKModule()
const std::shared_ptr<const MeshType> mesh =
std::dynamic_pointer_cast<const MeshType>(p_mesh);
- writer.write(mesh, OutputNamedItemValueSet{}, time, true);
+ writer.write(mesh,
+ OutputNamedItemValueSet{
+ NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+ time, true);
break;
}
case 3: {
@@ -40,7 +46,10 @@ VTKModule::VTKModule()
const std::shared_ptr<const MeshType> mesh =
std::dynamic_pointer_cast<const MeshType>(p_mesh);
- writer.write(mesh, OutputNamedItemValueSet{}, time, true);
+ writer.write(mesh,
+ OutputNamedItemValueSet{
+ NamedItemValue{"cell_number", mesh->connectivity().cellNumber()}},
+ time, true);
break;
}
}
diff --git a/src/language/utils/InterpolateItemValue.hpp b/src/language/utils/InterpolateItemValue.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..16a87b4952c7cb0c2907833705e8142d9c0a7178
--- /dev/null
+++ b/src/language/utils/InterpolateItemValue.hpp
@@ -0,0 +1,82 @@
+#ifndef INTERPOLATE_ITEM_VALUE_HPP
+#define INTERPOLATE_ITEM_VALUE_HPP
+
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <mesh/ItemType.hpp>
+#include <mesh/ItemValue.hpp>
+
+template <typename T>
+class InterpolateItemValue;
+template <typename OutputType, typename InputType>
+class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
+{
+ static constexpr size_t Dimension = OutputType::Dimension;
+ using Adapter = PugsFunctionAdapter<OutputType(InputType)>;
+
+ public:
+ template <ItemType item_type>
+ static inline ItemValue<OutputType, item_type>
+ interpolate(const FunctionSymbolId& function_symbol_id, const ItemValue<const InputType, item_type>& position)
+ {
+ auto& expression = Adapter::getFunctionExpression(function_symbol_id);
+ auto convert_result = Adapter::getResultConverter(expression.m_data_type);
+
+ Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
+
+ using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+ Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+ const IConnectivity& connectivity = *position.connectivity_ptr();
+
+ ItemValue<OutputType, item_type> value(connectivity);
+ using ItemId = ItemIdT<item_type>;
+
+ parallel_for(connectivity.template numberOf<item_type>(), [=, &expression, &tokens](ItemId i) {
+ const int32_t t = tokens.acquire();
+
+ auto& execution_policy = context_list[t];
+
+ Adapter::convertArgs(execution_policy.currentContext(), position[i]);
+ auto result = expression.execute(execution_policy);
+ value[i] = convert_result(std::move(result));
+
+ tokens.release(t);
+ });
+
+ return value;
+ }
+
+ template <ItemType item_type>
+ static inline Array<OutputType>
+ interpolate(const FunctionSymbolId& function_symbol_id,
+ const ItemValue<const InputType, item_type>& position,
+ const Array<const ItemIdT<item_type>>& list_of_items)
+ {
+ auto& expression = Adapter::getFunctionExpression(function_symbol_id);
+ auto convert_result = Adapter::getResultConverter(expression.m_data_type);
+
+ Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
+
+ using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+ Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+
+ Array<OutputType> value{list_of_items.size()};
+ using ItemId = ItemIdT<item_type>;
+
+ parallel_for(list_of_items.size(), [=, &expression, &tokens](size_t i_item) {
+ ItemId item_id = list_of_items[i_item];
+ const int32_t t = tokens.acquire();
+
+ auto& execution_policy = context_list[t];
+
+ Adapter::convertArgs(execution_policy.currentContext(), position[item_id]);
+ auto result = expression.execute(execution_policy);
+ value[i_item] = convert_result(std::move(result));
+
+ tokens.release(t);
+ });
+
+ return value;
+ }
+};
+
+#endif // INTERPOLATE_ITEM_VALUE_HPP
diff --git a/src/main.cpp b/src/main.cpp
index 9b853028133fd9035f67f743c62ee2a0a298202a..aad67d58a5d228cdc502e9b024cfc497fef23291 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,6 +1,8 @@
#include <utils/PugsUtils.hpp>
#include <language/PugsParser.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/SynchronizerManager.hpp>
@@ -11,9 +13,13 @@ main(int argc, char* argv[])
SynchronizerManager::create();
MeshDataManager::create();
+ DiamondDualConnectivityManager::create();
+ DiamondDualMeshManager::create();
parser(filename);
+ DiamondDualMeshManager::destroy();
+ DiamondDualConnectivityManager::destroy();
MeshDataManager::destroy();
SynchronizerManager::destroy();
finalize();
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 92c865c599767db0d5e5fed84e7d16d0cbec8579..6c0088831c86023e31914b953c9ea137e5411d21 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -8,8 +8,11 @@ add_library(
ConnectivityComputer.cpp
ConnectivityDispatcher.cpp
DiamondDualConnectivityBuilder.cpp
+ DiamondDualConnectivityManager.cpp
DiamondDualMeshBuilder.cpp
+ DiamondDualMeshManager.cpp
GmshReader.cpp
+ IConnectivity.cpp
IMesh.cpp
LogicalConnectivityBuilder.cpp
MeshBuilderBase.cpp
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 5a888190841e7886d4548a5405dc3126223318ac..3bf5cd4201c574749e5b6a67a2b6c4d667bd88f1 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -12,7 +12,6 @@
#include <mesh/RefId.hpp>
#include <mesh/RefItemList.hpp>
#include <mesh/SubItemValuePerItem.hpp>
-#include <mesh/SynchronizerManager.hpp>
#include <utils/CSRGraph.hpp>
#include <utils/Exceptions.hpp>
#include <utils/PugsAssert.hpp>
@@ -651,11 +650,7 @@ class Connectivity final : public IConnectivity
void _buildFrom(const ConnectivityDescriptor& descriptor);
public:
- ~Connectivity()
- {
- auto& manager = SynchronizerManager::instance();
- manager.deleteConnectivitySynchronizer(this);
- }
+ ~Connectivity() = default;
};
template <size_t Dim>
diff --git a/src/mesh/ConnectivityBuilderBase.cpp b/src/mesh/ConnectivityBuilderBase.cpp
index f5d47e25003a710e1e2c18ba6b772ffd04bb532f..60c38bfb84fb975c0e2cdc3b07d1d97e39eb0021 100644
--- a/src/mesh/ConnectivityBuilderBase.cpp
+++ b/src/mesh/ConnectivityBuilderBase.cpp
@@ -2,12 +2,13 @@
#include <mesh/Connectivity.hpp>
#include <mesh/ConnectivityDescriptor.hpp>
-#include <mesh/ConnectivityDispatcher.hpp>
#include <mesh/ItemId.hpp>
#include <mesh/Mesh.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
+#include <map>
+#include <unordered_map>
#include <vector>
template <size_t Dimension>
diff --git a/src/mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp b/src/mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8d35d2ff9541fbfab7e8e6e30976bd5f88970f56
--- /dev/null
+++ b/src/mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp
@@ -0,0 +1,153 @@
+#ifndef CONNECTIVITY_TO_DIAMOND_DUAL_CONNECTIVITY_DATA_MAPPER_HPP
+#define CONNECTIVITY_TO_DIAMOND_DUAL_CONNECTIVITY_DATA_MAPPER_HPP
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemIdToItemIdMap.hpp>
+#include <mesh/ItemValue.hpp>
+#include <utils/Array.hpp>
+#include <utils/PugsAssert.hpp>
+
+class IConnectivityToDiamondDualConnectivityDataMapper
+{
+ public:
+ IConnectivityToDiamondDualConnectivityDataMapper(const IConnectivityToDiamondDualConnectivityDataMapper&) = delete;
+ IConnectivityToDiamondDualConnectivityDataMapper(IConnectivityToDiamondDualConnectivityDataMapper&&) = delete;
+
+ IConnectivityToDiamondDualConnectivityDataMapper() = default;
+ virtual ~IConnectivityToDiamondDualConnectivityDataMapper() = default;
+};
+
+template <size_t Dimension>
+class ConnectivityToDiamondDualConnectivityDataMapper : public IConnectivityToDiamondDualConnectivityDataMapper
+{
+ private:
+ const IConnectivity* m_primal_connectivity;
+ const IConnectivity* m_dual_connectivity;
+
+ NodeIdToNodeIdMap m_primal_node_to_dual_node_map;
+ CellIdToNodeIdMap m_primal_cell_to_dual_node_map;
+ FaceIdToCellIdMap m_primal_face_to_dual_cell_map;
+
+ public:
+ template <typename OriginDataType1, typename OriginDataType2, typename DestinationDataType>
+ void
+ toDualNode(const NodeValue<OriginDataType1>& primal_node_value,
+ const CellValue<OriginDataType2>& primal_cell_value,
+ const NodeValue<DestinationDataType>& dual_node_value) const
+ {
+ static_assert(not std::is_const_v<DestinationDataType>, "destination data type must not be constant");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType1>, DestinationDataType>, "incompatible types");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType2>, DestinationDataType>, "incompatible types");
+
+ Assert(m_primal_connectivity == primal_cell_value.connectivity_ptr().get());
+ Assert(m_primal_connectivity == primal_node_value.connectivity_ptr().get());
+ Assert(m_dual_connectivity == dual_node_value.connectivity_ptr().get());
+
+ parallel_for(
+ m_primal_node_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_node_id, dual_node_id] = m_primal_node_to_dual_node_map[i];
+
+ dual_node_value[dual_node_id] = primal_node_value[primal_node_id];
+ });
+
+ parallel_for(
+ m_primal_cell_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_cell_id, dual_node_id] = m_primal_cell_to_dual_node_map[i];
+ dual_node_value[dual_node_id] = primal_cell_value[primal_cell_id];
+ });
+ }
+
+ template <typename OriginDataType, typename DestinationDataType1, typename DestinationDataType2>
+ void
+ fromDualNode(const NodeValue<OriginDataType>& dual_node_value,
+ const NodeValue<DestinationDataType1>& primal_node_value,
+ const CellValue<DestinationDataType2>& primal_cell_value) const
+ {
+ static_assert(not std::is_const_v<DestinationDataType1>, "destination data type must not be constant");
+ static_assert(not std::is_const_v<DestinationDataType2>, "destination data type must not be constant");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType1>, "incompatible types");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType2>, "incompatible types");
+
+ Assert(m_primal_connectivity == primal_cell_value.connectivity_ptr().get());
+ Assert(m_primal_connectivity == primal_node_value.connectivity_ptr().get());
+ Assert(m_dual_connectivity == dual_node_value.connectivity_ptr().get());
+
+ parallel_for(
+ m_primal_node_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_node_id, dual_node_id] = m_primal_node_to_dual_node_map[i];
+
+ primal_node_value[primal_node_id] = dual_node_value[dual_node_id];
+ });
+
+ parallel_for(
+ m_primal_cell_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_cell_id, dual_node_id] = m_primal_cell_to_dual_node_map[i];
+ primal_cell_value[primal_cell_id] = dual_node_value[dual_node_id];
+ });
+ }
+
+ template <typename OriginDataType, typename DestinationDataType>
+ void
+ toDualCell(const FaceValue<OriginDataType>& primal_face_value,
+ const CellValue<DestinationDataType>& dual_cell_value) const
+ {
+ static_assert(not std::is_const_v<DestinationDataType>, "destination data type must not be constant");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType>, "incompatible types");
+
+ Assert(m_primal_connectivity == primal_face_value.connectivity_ptr().get());
+ Assert(m_dual_connectivity == dual_cell_value.connectivity_ptr().get());
+
+ parallel_for(
+ m_primal_face_to_dual_cell_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+
+ dual_cell_value[dual_cell_id] = primal_face_value[primal_face_id];
+ });
+
+ parallel_for(
+ m_primal_cell_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+
+ dual_cell_value[dual_cell_id] = primal_face_value[primal_face_id];
+ });
+ }
+
+ template <typename OriginDataType, typename DestinationDataType>
+ void
+ fromDualCell(const CellValue<DestinationDataType>& dual_cell_value,
+ const FaceValue<OriginDataType>& primal_face_value) const
+ {
+ static_assert(not std::is_const_v<DestinationDataType>, "destination data type must not be constant");
+ static_assert(std::is_same_v<std::remove_const_t<OriginDataType>, DestinationDataType>, "incompatible types");
+
+ Assert(m_primal_connectivity == primal_face_value.connectivity_ptr().get());
+ Assert(m_dual_connectivity == dual_cell_value.connectivity_ptr().get());
+
+ parallel_for(
+ m_primal_face_to_dual_cell_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+
+ primal_face_value[primal_face_id] = dual_cell_value[dual_cell_id];
+ });
+
+ parallel_for(
+ m_primal_cell_to_dual_node_map.size(), PUGS_LAMBDA(size_t i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+ primal_face_value[primal_face_id] = dual_cell_value[dual_cell_id];
+ });
+ }
+
+ ConnectivityToDiamondDualConnectivityDataMapper(const Connectivity<Dimension>& primal_connectivity,
+ const Connectivity<Dimension>& dual_connectivity,
+ const NodeIdToNodeIdMap& primal_node_to_dual_node_map,
+ const CellIdToNodeIdMap& primal_cell_to_dual_node_map,
+ const FaceIdToCellIdMap& primal_face_to_dual_cell_map)
+ : m_primal_connectivity{&primal_connectivity},
+ m_dual_connectivity{&dual_connectivity},
+ m_primal_node_to_dual_node_map{primal_node_to_dual_node_map},
+ m_primal_cell_to_dual_node_map{primal_cell_to_dual_node_map},
+ m_primal_face_to_dual_cell_map{primal_face_to_dual_cell_map}
+ {}
+};
+
+#endif // CONNECTIVITY_TO_DIAMOND_DUAL_CONNECTIVITY_DATA_MAPPER_HPP
diff --git a/src/mesh/DiamondDualConnectivityBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp
index 8bcbe9bad8e15cbef1dfa449c0d6defc115f126f..31d54be8997e9571958399f17c44dbcede1c35e4 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.cpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.cpp
@@ -3,12 +3,16 @@
#include <mesh/Connectivity.hpp>
#include <mesh/ConnectivityDescriptor.hpp>
#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/RefId.hpp>
#include <utils/Array.hpp>
+#include <utils/ArrayUtils.hpp>
#include <utils/Messenger.hpp>
+#include <vector>
+
template <size_t Dimension>
void
DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<Dimension>& primal_connectivity,
@@ -18,32 +22,56 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
+ const auto& primal_node_number = primal_connectivity.nodeNumber();
+ const auto& primal_cell_number = primal_connectivity.cellNumber();
+
const size_t diamond_number_of_nodes = primal_number_of_cells + primal_number_of_nodes;
+ const size_t diamond_number_of_cells = primal_number_of_faces;
- diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);
+ {
+ m_primal_node_to_dual_node_map = NodeIdToNodeIdMap{primal_number_of_nodes};
- const auto& primal_node_number = primal_connectivity.nodeNumber();
+ NodeId diamond_node_id = 0;
+ for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) {
+ m_primal_node_to_dual_node_map[primal_node_id] = std::make_pair(primal_node_id, diamond_node_id++);
+ }
- for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
- diamond_descriptor.node_number_vector[primal_node_id] = primal_node_number[primal_node_id];
+ m_primal_cell_to_dual_node_map = CellIdToNodeIdMap{primal_number_of_cells};
+ for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
+ m_primal_cell_to_dual_node_map[primal_cell_id] = std::make_pair(primal_cell_id, diamond_node_id++);
+ }
}
- const auto& primal_cell_number = primal_connectivity.cellNumber();
+ diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);
- const size_t max_node_number = max(primal_node_number);
+ parallel_for(m_primal_node_to_dual_node_map.size(), [&](size_t i) {
+ const auto [primal_node_id, diamond_dual_node_id] = m_primal_node_to_dual_node_map[i];
- for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
- diamond_descriptor.node_number_vector[primal_number_of_nodes + primal_cell_id] =
- primal_cell_number[primal_cell_id] + max_node_number;
+ diamond_descriptor.node_number_vector[diamond_dual_node_id] = primal_node_number[primal_node_id];
+ });
+
+ const size_t cell_number_shift = max(primal_node_number) + 1;
+ parallel_for(primal_number_of_cells, [&](size_t i) {
+ const auto [primal_cell_id, diamond_dual_node_id] = m_primal_cell_to_dual_node_map[i];
+
+ diamond_descriptor.node_number_vector[diamond_dual_node_id] =
+ primal_cell_number[primal_cell_id] + cell_number_shift;
+ });
+
+ {
+ m_primal_face_to_dual_cell_map = FaceIdToCellIdMap{primal_number_of_faces};
+ CellId diamond_cell_id = 0;
+ for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_faces; ++primal_face_id) {
+ m_primal_face_to_dual_cell_map[primal_face_id] = std::make_pair(primal_face_id, diamond_cell_id++);
+ }
}
- const size_t diamond_number_of_cells = primal_number_of_faces;
diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells);
-
const auto& primal_face_number = primal_connectivity.faceNumber();
- for (FaceId i_primal_face = 0; i_primal_face < primal_number_of_faces; ++i_primal_face) {
- diamond_descriptor.cell_number_vector[i_primal_face] = primal_face_number[i_primal_face];
- }
+ parallel_for(diamond_number_of_cells, [&](size_t i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+ diamond_descriptor.cell_number_vector[dual_cell_id] = primal_face_number[primal_face_id];
+ });
if constexpr (Dimension == 3) {
const size_t number_of_edges = diamond_descriptor.edge_to_node_vector.size();
@@ -60,9 +88,9 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();
- for (FaceId i_face = 0; i_face < primal_connectivity.numberOfFaces(); ++i_face) {
- const size_t i_cell = i_face;
- const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face];
+ parallel_for(primal_number_of_faces, [&](FaceId face_id) {
+ const size_t i_cell = face_id;
+ const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
if (primal_face_cell_list.size() == 1) {
diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle;
@@ -79,22 +107,22 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
diamond_descriptor.cell_type_vector[i_cell] = CellType::Diamond;
}
}
- }
+ });
diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
const auto& primal_face_local_number_in_their_cells = primal_connectivity.faceLocalNumbersInTheirCells();
const auto& cell_face_is_reversed = primal_connectivity.cellFaceIsReversed();
- for (FaceId i_face = 0; i_face < primal_connectivity.numberOfFaces(); ++i_face) {
- const size_t& i_diamond_cell = i_face;
- const auto& primal_face_cell_list = primal_face_to_cell_matrix[i_face];
- const auto& primal_face_node_list = primal_face_to_node_matrix[i_face];
+ parallel_for(primal_number_of_faces, [&](FaceId face_id) {
+ const size_t& i_diamond_cell = face_id;
+ const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
+ const auto& primal_face_node_list = primal_face_to_node_matrix[face_id];
if (primal_face_cell_list.size() == 1) {
diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(primal_face_node_list.size() + 1);
const CellId cell_id = primal_face_cell_list[0];
- const auto i_face_in_cell = primal_face_local_number_in_their_cells(i_face, 0);
+ const auto i_face_in_cell = primal_face_local_number_in_their_cells(face_id, 0);
for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node] = primal_face_node_list[i_node];
@@ -121,7 +149,7 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
const CellId cell0_id = primal_face_cell_list[0];
const CellId cell1_id = primal_face_cell_list[1];
- const auto i_face_in_cell0 = primal_face_local_number_in_their_cells(i_face, 0);
+ const auto i_face_in_cell0 = primal_face_local_number_in_their_cells(face_id, 0);
if constexpr (Dimension == 2) {
Assert(primal_face_node_list.size() == 2);
@@ -148,7 +176,7 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
}
}
}
- }
+ });
}
template <>
@@ -162,23 +190,34 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor<1>(const Con
const size_t diamond_number_of_nodes = 2 * primal_number_of_nodes - primal_number_of_cells;
- diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);
-
- const auto& primal_node_number = primal_connectivity.nodeNumber();
-
const size_t diamond_number_of_cells = primal_number_of_faces;
const size_t number_of_kept_nodes = 2 * (diamond_number_of_nodes - diamond_number_of_cells);
const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
size_t next_kept_node_id = 0;
- for (NodeId node_id = 0; node_id < primal_connectivity.numberOfNodes(); ++node_id) {
- const auto& primal_node_cell_list = primal_node_to_cell_matrix[node_id];
+ diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);
+
+ const auto& primal_node_number = primal_connectivity.nodeNumber();
+
+ for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
+ const auto& primal_node_cell_list = primal_node_to_cell_matrix[primal_node_id];
if (primal_node_cell_list.size() == 1) {
- diamond_descriptor.node_number_vector[next_kept_node_id++] = primal_node_number[node_id];
+ diamond_descriptor.node_number_vector[next_kept_node_id++] = primal_node_number[primal_node_id];
}
}
+ const size_t primal_number_of_kept_nodes = next_kept_node_id;
+
+ const auto& primal_cell_number = primal_connectivity.cellNumber();
+
+ const size_t cell_number_shift = max(primal_node_number) + 1;
+
+ for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
+ diamond_descriptor.node_number_vector[primal_number_of_kept_nodes + primal_cell_id] =
+ primal_cell_number[primal_cell_id] + cell_number_shift;
+ }
+
if (number_of_kept_nodes != next_kept_node_id) {
throw UnexpectedError("unexpected number of kept node" + std::to_string(next_kept_node_id) +
" != " + std::to_string(number_of_kept_nodes));
@@ -200,8 +239,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor<1>(const Con
diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
- ;
-
const auto& primal_node_local_number_in_their_cells = primal_connectivity.nodeLocalNumbersInTheirCells();
{
@@ -278,7 +315,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
node_array[i] = diamond_node_list[i];
}
diamond_descriptor.addRefItemList(RefNodeList{primal_ref_node_list.refId(), node_array});
- std::cout << "stored " << primal_ref_node_list.refId() << '\n';
}
}
}
@@ -332,7 +368,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
face_array[i] = diamond_face_list[i];
}
diamond_descriptor.addRefItemList(RefFaceList{primal_ref_face_list.refId(), face_array});
- std::cout << "stored " << primal_ref_face_list.refId() << '\n';
}
}
}
@@ -384,7 +419,6 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
edge_array[i] = diamond_edge_list[i];
}
diamond_descriptor.addRefItemList(RefEdgeList{primal_ref_edge_list.refId(), edge_array});
- std::cout << "stored " << primal_ref_edge_list.refId() << '\n';
}
}
}
@@ -466,15 +500,60 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
diamond_descriptor.edge_owner_vector.resize(edge_cell_owner.size());
for (size_t i_edge = 0; i_edge < edge_cell_owner.size(); ++i_edge) {
- diamond_descriptor.face_owner_vector[i_edge] = diamond_descriptor.cell_owner_vector[edge_cell_owner[i_edge]];
+ diamond_descriptor.edge_owner_vector[i_edge] = diamond_descriptor.cell_owner_vector[edge_cell_owner[i_edge]];
}
}
m_connectivity = ConnectivityType::build(diamond_descriptor);
+
+ if constexpr (Dimension == 1) {
+ const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
+
+ NodeId dual_node_id = 0;
+ m_primal_node_to_dual_node_map = [&]() {
+ std::vector<std::pair<NodeId, NodeId>> primal_node_to_dual_node_vector;
+ for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
+ if (node_to_cell_matrix[primal_node_id].size() == 1) {
+ primal_node_to_dual_node_vector.push_back(std::make_pair(primal_node_id, dual_node_id++));
+ }
+ }
+ return convert_to_array(primal_node_to_dual_node_vector);
+ }();
+
+ m_primal_cell_to_dual_node_map = [&]() {
+ CellIdToNodeIdMap primal_cell_to_dual_node_map{primal_number_of_cells};
+ for (CellId primal_cell_id = 0; primal_cell_id < primal_cell_to_dual_node_map.size(); ++primal_cell_id) {
+ primal_cell_to_dual_node_map[primal_cell_id] = std::make_pair(primal_cell_id, dual_node_id++);
+ }
+ return primal_cell_to_dual_node_map;
+ }();
+
+ m_primal_face_to_dual_cell_map = [&]() {
+ FaceIdToCellIdMap primal_face_to_dual_cell_map{primal_connectivity.numberOfFaces()};
+ for (size_t id = 0; id < primal_face_to_dual_cell_map.size(); ++id) {
+ const CellId dual_cell_id = id;
+ const FaceId primal_face_id = id;
+
+ primal_face_to_dual_cell_map[id] = std::make_pair(primal_face_id, dual_cell_id);
+ }
+ return primal_face_to_dual_cell_map;
+ }();
+ }
+
+ m_mapper =
+ std::make_shared<ConnectivityToDiamondDualConnectivityDataMapper<Dimension>>(primal_connectivity,
+ dynamic_cast<const ConnectivityType&>(
+ *m_connectivity),
+ m_primal_node_to_dual_node_map,
+ m_primal_cell_to_dual_node_map,
+ m_primal_face_to_dual_cell_map);
}
DiamondDualConnectivityBuilder::DiamondDualConnectivityBuilder(const IConnectivity& connectivity)
{
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Construction of diamond dual mesh is not implemented in parallel");
+ }
switch (connectivity.dimension()) {
case 1: {
this->_buildDiamondConnectivityFrom<1>(connectivity);
diff --git a/src/mesh/DiamondDualConnectivityBuilder.hpp b/src/mesh/DiamondDualConnectivityBuilder.hpp
index 670372df9e53d3bfa8dd3f1b03545d9264b41dc4..36b120811ad6478af90b89112740d272014b2e7a 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.hpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.hpp
@@ -2,6 +2,7 @@
#define DIAMOND_DUAL_CONNECTIVITY_BUILDER_HPP
#include <mesh/ConnectivityBuilderBase.hpp>
+#include <mesh/ItemIdToItemIdMap.hpp>
#include <memory>
@@ -9,17 +10,33 @@ template <size_t>
class Connectivity;
class ConnectivityDescriptor;
+class IConnectivityToDiamondDualConnectivityDataMapper;
+
class DiamondDualConnectivityBuilder : public ConnectivityBuilderBase
{
private:
+ NodeIdToNodeIdMap m_primal_node_to_dual_node_map;
+ CellIdToNodeIdMap m_primal_cell_to_dual_node_map;
+ FaceIdToCellIdMap m_primal_face_to_dual_cell_map;
+
+ std::shared_ptr<IConnectivityToDiamondDualConnectivityDataMapper> m_mapper;
+
template <size_t Dimension>
void _buildDiamondConnectivityDescriptor(const Connectivity<Dimension>&, ConnectivityDescriptor&);
template <size_t Dimension>
void _buildDiamondConnectivityFrom(const IConnectivity&);
- public:
+ friend class DiamondDualConnectivityManager;
DiamondDualConnectivityBuilder(const IConnectivity&);
+
+ public:
+ std::shared_ptr<IConnectivityToDiamondDualConnectivityDataMapper>
+ mapper() const
+ {
+ return m_mapper;
+ }
+
~DiamondDualConnectivityBuilder() = default;
};
diff --git a/src/mesh/DiamondDualConnectivityManager.cpp b/src/mesh/DiamondDualConnectivityManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9459010b69e89acaf3f247866d2f6e81ef00cd81
--- /dev/null
+++ b/src/mesh/DiamondDualConnectivityManager.cpp
@@ -0,0 +1,96 @@
+#include <utils/PugsAssert.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/DiamondDualConnectivityBuilder.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <utils/Exceptions.hpp>
+
+#include <sstream>
+
+DiamondDualConnectivityManager* DiamondDualConnectivityManager::m_instance{nullptr};
+
+void
+DiamondDualConnectivityManager::create()
+{
+ Assert(m_instance == nullptr, "DiamondDualConnectivityManager is already created");
+ m_instance = new DiamondDualConnectivityManager;
+}
+
+void
+DiamondDualConnectivityManager::destroy()
+{
+ Assert(m_instance != nullptr, "DiamondDualConnectivityManager was not created!");
+
+ if (m_instance->m_connectivity_to_diamond_dual_connectivity_info_map.size() > 0) {
+ std::stringstream error;
+ error << ": some connectivities are still registered\n";
+ for (const auto& [connectivity, diamond_dual_connectivity_info] :
+ m_instance->m_connectivity_to_diamond_dual_connectivity_info_map) {
+ error << " - connectivity " << rang::fgB::magenta << connectivity << rang::style::reset << '\n';
+ }
+ throw UnexpectedError(error.str());
+ }
+ delete m_instance;
+ m_instance = nullptr;
+}
+
+void
+DiamondDualConnectivityManager::deleteConnectivity(const IConnectivity* p_connectivity)
+{
+ m_connectivity_to_diamond_dual_connectivity_info_map.erase(p_connectivity);
+}
+
+DiamondDualConnectivityManager::DiamondDualConnectivityInfo
+DiamondDualConnectivityManager::_getDiamondDualConnectivityInfo(const IConnectivity& connectivity)
+{
+ const IConnectivity* p_connectivity = &connectivity;
+
+ if (auto i_connectivity = m_connectivity_to_diamond_dual_connectivity_info_map.find(p_connectivity);
+ i_connectivity != m_connectivity_to_diamond_dual_connectivity_info_map.end()) {
+ return i_connectivity->second;
+ } else {
+ DiamondDualConnectivityBuilder builder{connectivity};
+
+ DiamondDualConnectivityInfo connectivity_info{builder.connectivity(), builder.mapper()};
+
+ m_connectivity_to_diamond_dual_connectivity_info_map[p_connectivity] = connectivity_info;
+
+ return connectivity_info;
+ }
+}
+
+template <size_t Dimension>
+std::shared_ptr<const Connectivity<Dimension>>
+DiamondDualConnectivityManager::getDiamondDualConnectivity(const Connectivity<Dimension>& connectivity)
+{
+ return std::dynamic_pointer_cast<const Connectivity<Dimension>>(
+ this->_getDiamondDualConnectivityInfo(connectivity).diamondDualConnectivity());
+}
+
+template <size_t Dimension>
+std::shared_ptr<const ConnectivityToDiamondDualConnectivityDataMapper<Dimension>>
+DiamondDualConnectivityManager::getConnectivityToDiamondDualConnectivityDataMapper(
+ const Connectivity<Dimension>& connectivity)
+{
+ return std::dynamic_pointer_cast<const ConnectivityToDiamondDualConnectivityDataMapper<Dimension>>(
+ this->_getDiamondDualConnectivityInfo(connectivity).connectivityToDiamondDualConnectivityDataMapper());
+}
+
+template std::shared_ptr<const Connectivity<1>> DiamondDualConnectivityManager::getDiamondDualConnectivity(
+ const Connectivity<1>& connectivity);
+
+template std::shared_ptr<const Connectivity<2>> DiamondDualConnectivityManager::getDiamondDualConnectivity(
+ const Connectivity<2>& connectivity);
+
+template std::shared_ptr<const Connectivity<3>> DiamondDualConnectivityManager::getDiamondDualConnectivity(
+ const Connectivity<3>& connectivity);
+
+template std::shared_ptr<const ConnectivityToDiamondDualConnectivityDataMapper<1>>
+DiamondDualConnectivityManager::getConnectivityToDiamondDualConnectivityDataMapper(const Connectivity<1>&);
+
+template std::shared_ptr<const ConnectivityToDiamondDualConnectivityDataMapper<2>>
+DiamondDualConnectivityManager::getConnectivityToDiamondDualConnectivityDataMapper(const Connectivity<2>&);
+
+template std::shared_ptr<const ConnectivityToDiamondDualConnectivityDataMapper<3>>
+DiamondDualConnectivityManager::getConnectivityToDiamondDualConnectivityDataMapper(const Connectivity<3>&);
diff --git a/src/mesh/DiamondDualConnectivityManager.hpp b/src/mesh/DiamondDualConnectivityManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..33dfbaf4cf7cc8e59eebfc095624d43c5bf038b3
--- /dev/null
+++ b/src/mesh/DiamondDualConnectivityManager.hpp
@@ -0,0 +1,94 @@
+#ifndef DIAMOND_DUAL_CONNECTIVITY_MANAGER_HPP
+#define DIAMOND_DUAL_CONNECTIVITY_MANAGER_HPP
+
+#include <mesh/IConnectivity.hpp>
+
+#include <memory>
+#include <unordered_map>
+
+template <size_t Dimension>
+class Connectivity;
+
+class IConnectivityToDiamondDualConnectivityDataMapper;
+
+template <size_t Dimension>
+class ConnectivityToDiamondDualConnectivityDataMapper;
+
+class DiamondDualConnectivityManager
+{
+ private:
+ class DiamondDualConnectivityInfo
+ {
+ private:
+ std::shared_ptr<const IConnectivity> m_diamond_dual_connectivity;
+ std::shared_ptr<const IConnectivityToDiamondDualConnectivityDataMapper>
+ m_connectivity_to_diamond_dual_connectivity_data_mapper;
+
+ public:
+ PUGS_INLINE
+ std::shared_ptr<const IConnectivity>
+ diamondDualConnectivity() const
+ {
+ return m_diamond_dual_connectivity;
+ }
+
+ PUGS_INLINE
+ std::shared_ptr<const IConnectivityToDiamondDualConnectivityDataMapper>
+ connectivityToDiamondDualConnectivityDataMapper()
+ {
+ return m_connectivity_to_diamond_dual_connectivity_data_mapper;
+ }
+
+ DiamondDualConnectivityInfo& operator=(const DiamondDualConnectivityInfo&) = default;
+ DiamondDualConnectivityInfo& operator=(DiamondDualConnectivityInfo&&) = default;
+
+ DiamondDualConnectivityInfo() = default;
+ DiamondDualConnectivityInfo(const DiamondDualConnectivityInfo&) = default;
+ DiamondDualConnectivityInfo(DiamondDualConnectivityInfo&&) = default;
+
+ DiamondDualConnectivityInfo(const std::shared_ptr<const IConnectivity>& diamond_dual_connectivity,
+ const std::shared_ptr<const IConnectivityToDiamondDualConnectivityDataMapper>&
+ connectivity_to_diamond_dual_connectivity_data_mapper)
+ : m_diamond_dual_connectivity{diamond_dual_connectivity},
+ m_connectivity_to_diamond_dual_connectivity_data_mapper{connectivity_to_diamond_dual_connectivity_data_mapper}
+ {}
+
+ ~DiamondDualConnectivityInfo() = default;
+ };
+
+ DiamondDualConnectivityInfo _getDiamondDualConnectivityInfo(const IConnectivity& connectivity);
+
+ std::unordered_map<const IConnectivity*, DiamondDualConnectivityInfo>
+ m_connectivity_to_diamond_dual_connectivity_info_map;
+
+ static DiamondDualConnectivityManager* m_instance;
+
+ DiamondDualConnectivityManager(const DiamondDualConnectivityManager&) = delete;
+ DiamondDualConnectivityManager(DiamondDualConnectivityManager&&) = delete;
+
+ DiamondDualConnectivityManager() = default;
+ ~DiamondDualConnectivityManager() = default;
+
+ public:
+ static void create();
+ static void destroy();
+
+ PUGS_INLINE
+ static DiamondDualConnectivityManager&
+ instance()
+ {
+ Assert(m_instance != nullptr, "DiamondDualConnectivityManager was not created!");
+ return *m_instance;
+ }
+
+ void deleteConnectivity(const IConnectivity*);
+
+ template <size_t Dimension>
+ std::shared_ptr<const Connectivity<Dimension>> getDiamondDualConnectivity(const Connectivity<Dimension>&);
+
+ template <size_t Dimension>
+ std::shared_ptr<const ConnectivityToDiamondDualConnectivityDataMapper<Dimension>>
+ getConnectivityToDiamondDualConnectivityDataMapper(const Connectivity<Dimension>&);
+};
+
+#endif // DIAMOND_DUAL_CONNECTIVITY_MANAGER_HPP
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 2218541b338e196a913b3c0859b882ad668d252e..df4157c9c19b377f77b754d3f33ac54141539f2f 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -1,117 +1,60 @@
#include <mesh/DiamondDualMeshBuilder.hpp>
-#include <mesh/DiamondDualConnectivityBuilder.hpp>
-
#include <mesh/Connectivity.hpp>
-#include <mesh/ConnectivityDescriptor.hpp>
-#include <mesh/ConnectivityDispatcher.hpp>
+#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/DiamondDualConnectivityBuilder.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
-#include <mesh/RefId.hpp>
-#include <utils/Array.hpp>
-#include <utils/Messenger.hpp>
template <size_t Dimension>
void
-DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(
- const Mesh<Connectivity<Dimension>>& primal_mesh,
- const std::shared_ptr<const Connectivity<Dimension>>& p_diamond_connectivity)
+DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const std::shared_ptr<const IMesh>& p_i_mesh)
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<Connectivity<Dimension>>;
- const ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
-
- NodeValue<TinyVector<Dimension>> diamond_xr{diamond_connectivity};
-
- const NodeValue<const TinyVector<Dimension>> primal_xr = primal_mesh.xr();
- MeshData<Dimension>& primal_mesh_data = MeshDataManager::instance().getMeshData(primal_mesh);
- const CellValue<const TinyVector<Dimension>> primal_xj = primal_mesh_data.xj();
-
- NodeId i_node = 0;
- for (; i_node < primal_mesh.numberOfNodes(); ++i_node) {
- diamond_xr[i_node] = primal_xr[i_node];
- }
+ std::shared_ptr p_primal_mesh = std::dynamic_pointer_cast<const MeshType>(p_i_mesh);
+ const MeshType& primal_mesh = *p_primal_mesh;
- for (CellId i_cell = 0; i_cell < primal_mesh.numberOfCells(); ++i_cell) {
- diamond_xr[i_node++] = primal_xj[i_cell];
- }
+ DiamondDualConnectivityManager& manager = DiamondDualConnectivityManager::instance();
- m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
-}
-
-template <>
-void
-DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const Mesh<Connectivity<1>>& primal_mesh,
- const std::shared_ptr<const Connectivity<1>>& p_diamond_connectivity)
-{
- using ConnectivityType = Connectivity<1>;
- using MeshType = Mesh<Connectivity<1>>;
-
- const ConnectivityType& primal_connectivity = primal_mesh.connectivity();
- const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
+ std::shared_ptr<const ConnectivityType> p_diamond_connectivity =
+ manager.getDiamondDualConnectivity(primal_mesh.connectivity());
const ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
- NodeValue<TinyVector<1>> diamond_xr{diamond_connectivity};
+ const NodeValue<const TinyVector<Dimension>> primal_xr = primal_mesh.xr();
- const NodeValue<const TinyVector<1>> primal_xr = primal_mesh.xr();
- MeshData<1>& primal_mesh_data = MeshDataManager::instance().getMeshData(primal_mesh);
- const CellValue<const TinyVector<1>> primal_xj = primal_mesh_data.xj();
+ MeshData<Dimension>& primal_mesh_data = MeshDataManager::instance().getMeshData(primal_mesh);
+ const CellValue<const TinyVector<Dimension>> primal_xj = primal_mesh_data.xj();
- NodeId next_node_id = 0;
- for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
- if (node_to_cell_matrix[primal_node_id].size() == 1) {
- diamond_xr[next_node_id++] = primal_xr[primal_node_id];
- }
- }
+ std::shared_ptr connectivity_to_diamond_dual_connectivity_data_mapper =
+ manager.getConnectivityToDiamondDualConnectivityDataMapper(primal_mesh.connectivity());
- for (CellId i_cell = 0; i_cell < primal_mesh.numberOfCells(); ++i_cell) {
- diamond_xr[next_node_id++] = primal_xj[i_cell];
- }
+ NodeValue<TinyVector<Dimension>> diamond_xr{diamond_connectivity};
+ connectivity_to_diamond_dual_connectivity_data_mapper->toDualNode(primal_xr, primal_xj, diamond_xr);
m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
}
DiamondDualMeshBuilder::DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& p_mesh)
{
+ std::cout << "building DiamondDualMesh\n";
+
switch (p_mesh->dimension()) {
case 1: {
- using ConnectivityType = Connectivity<1>;
- using MeshType = Mesh<ConnectivityType>;
-
- std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(p_mesh);
- DiamondDualConnectivityBuilder builder(mesh->connectivity());
-
- std::shared_ptr p_diamond_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
-
- this->_buildDualDiamondMeshFrom(*mesh, p_diamond_connectivity);
+ this->_buildDualDiamondMeshFrom<1>(p_mesh);
break;
}
case 2: {
- using ConnectivityType = Connectivity<2>;
- using MeshType = Mesh<ConnectivityType>;
-
- std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(p_mesh);
- DiamondDualConnectivityBuilder builder(mesh->connectivity());
-
- std::shared_ptr p_diamond_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
-
- this->_buildDualDiamondMeshFrom(*mesh, p_diamond_connectivity);
+ this->_buildDualDiamondMeshFrom<2>(p_mesh);
break;
}
case 3: {
- using ConnectivityType = Connectivity<3>;
- using MeshType = Mesh<ConnectivityType>;
-
- std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(p_mesh);
- DiamondDualConnectivityBuilder builder(mesh->connectivity());
-
- std::shared_ptr p_diamond_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
-
- this->_buildDualDiamondMeshFrom(*mesh, p_diamond_connectivity);
+ this->_buildDualDiamondMeshFrom<3>(p_mesh);
break;
}
default: {
diff --git a/src/mesh/DiamondDualMeshBuilder.hpp b/src/mesh/DiamondDualMeshBuilder.hpp
index ded3cfb0f40ca4d1b1bd9dc06f7981391f9401c6..2e47cfe21b4790a63e0e4aeba2fcfb9c59699286 100644
--- a/src/mesh/DiamondDualMeshBuilder.hpp
+++ b/src/mesh/DiamondDualMeshBuilder.hpp
@@ -5,21 +5,16 @@
#include <memory>
-template <size_t>
-class Connectivity;
-
-template <typename ConnectivityType>
-class Mesh;
-
class DiamondDualMeshBuilder : public MeshBuilderBase
{
private:
template <size_t Dimension>
- void _buildDualDiamondMeshFrom(const Mesh<Connectivity<Dimension>>&,
- const std::shared_ptr<const Connectivity<Dimension>>&);
+ void _buildDualDiamondMeshFrom(const std::shared_ptr<const IMesh>&);
- public:
+ friend class DiamondDualMeshManager;
DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>&);
+
+ public:
~DiamondDualMeshBuilder() = default;
};
diff --git a/src/mesh/DiamondDualMeshManager.cpp b/src/mesh/DiamondDualMeshManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0619736753f2b49a4229c60492ea6fd4c0de7212
--- /dev/null
+++ b/src/mesh/DiamondDualMeshManager.cpp
@@ -0,0 +1,65 @@
+#include <mesh/DiamondDualMeshManager.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/DiamondDualMeshBuilder.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Exceptions.hpp>
+#include <utils/PugsAssert.hpp>
+
+#include <sstream>
+
+DiamondDualMeshManager* DiamondDualMeshManager::m_instance{nullptr};
+
+void
+DiamondDualMeshManager::create()
+{
+ Assert(m_instance == nullptr, "DiamondDualMeshManager is already created");
+ m_instance = new DiamondDualMeshManager;
+}
+
+void
+DiamondDualMeshManager::destroy()
+{
+ Assert(m_instance != nullptr, "DiamondDualMeshManager was not created!");
+
+ if (m_instance->m_mesh_to_diamond_dual_mesh_map.size() > 0) {
+ std::stringstream error;
+ error << ": some meshes are still registered\n";
+ for (const auto& i_mesh_data : m_instance->m_mesh_to_diamond_dual_mesh_map) {
+ error << " - mesh " << rang::fgB::magenta << i_mesh_data.first << rang::style::reset << '\n';
+ }
+ throw UnexpectedError(error.str());
+ }
+ delete m_instance;
+ m_instance = nullptr;
+}
+
+void
+DiamondDualMeshManager::deleteMesh(const IMesh* p_mesh)
+{
+ m_mesh_to_diamond_dual_mesh_map.erase(p_mesh);
+}
+
+template <size_t Dimension>
+std::shared_ptr<const Mesh<Connectivity<Dimension>>>
+DiamondDualMeshManager::getDiamondDualMesh(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh)
+{
+ const IMesh* p_mesh = mesh.get();
+
+ if (auto i_mesh_data = m_mesh_to_diamond_dual_mesh_map.find(p_mesh);
+ i_mesh_data != m_mesh_to_diamond_dual_mesh_map.end()) {
+ return std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(i_mesh_data->second);
+ } else {
+ DiamondDualMeshBuilder builder{mesh};
+
+ m_mesh_to_diamond_dual_mesh_map[p_mesh] = builder.mesh();
+ return std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(builder.mesh());
+ }
+}
+
+template std::shared_ptr<const Mesh<Connectivity<1>>> DiamondDualMeshManager::getDiamondDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<1>>>);
+template std::shared_ptr<const Mesh<Connectivity<2>>> DiamondDualMeshManager::getDiamondDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<2>>>);
+template std::shared_ptr<const Mesh<Connectivity<3>>> DiamondDualMeshManager::getDiamondDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<3>>>);
diff --git a/src/mesh/DiamondDualMeshManager.hpp b/src/mesh/DiamondDualMeshManager.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9f802eb96d13704235f50e82ad75e84ba9be8a64
--- /dev/null
+++ b/src/mesh/DiamondDualMeshManager.hpp
@@ -0,0 +1,49 @@
+#ifndef DIAMOND_DUAL_MESH_MANAGER_HPP
+#define DIAMOND_DUAL_MESH_MANAGER_HPP
+
+#include <mesh/IMesh.hpp>
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <memory>
+#include <unordered_map>
+
+template <size_t Dimension>
+class Connectivity;
+
+template <typename ConnectivityType>
+class Mesh;
+
+class DiamondDualMeshManager
+{
+ private:
+ std::unordered_map<const IMesh*, std::shared_ptr<const IMesh>> m_mesh_to_diamond_dual_mesh_map;
+
+ static DiamondDualMeshManager* m_instance;
+
+ DiamondDualMeshManager(const DiamondDualMeshManager&) = delete;
+ DiamondDualMeshManager(DiamondDualMeshManager&&) = delete;
+
+ DiamondDualMeshManager() = default;
+ ~DiamondDualMeshManager() = default;
+
+ public:
+ static void create();
+ static void destroy();
+
+ PUGS_INLINE
+ static DiamondDualMeshManager&
+ instance()
+ {
+ Assert(m_instance != nullptr, "DiamondDualMeshManager was not created!");
+ return *m_instance;
+ }
+
+ void deleteMesh(const IMesh*);
+
+ template <size_t Dimension>
+ std::shared_ptr<const Mesh<Connectivity<Dimension>>> getDiamondDualMesh(
+ std::shared_ptr<const Mesh<Connectivity<Dimension>>>);
+};
+
+#endif // DIAMOND_DUAL_MESH_MANAGER_HPP
diff --git a/src/mesh/IConnectivity.cpp b/src/mesh/IConnectivity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..48cf57d55ebceb7818f755fd9214807a247f6575
--- /dev/null
+++ b/src/mesh/IConnectivity.cpp
@@ -0,0 +1,10 @@
+#include <mesh/IConnectivity.hpp>
+
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <mesh/SynchronizerManager.hpp>
+
+IConnectivity::~IConnectivity()
+{
+ SynchronizerManager::instance().deleteConnectivitySynchronizer(this);
+ DiamondDualConnectivityManager::instance().deleteConnectivity(this);
+}
diff --git a/src/mesh/IConnectivity.hpp b/src/mesh/IConnectivity.hpp
index 27be707cbcc4976b17e9db1e7bf6f37195c9d63a..9bfb8e5043088a18880440235cd72b1677c1f519 100644
--- a/src/mesh/IConnectivity.hpp
+++ b/src/mesh/IConnectivity.hpp
@@ -35,7 +35,7 @@ class IConnectivity : public std::enable_shared_from_this<IConnectivity>
IConnectivity() = default;
IConnectivity(IConnectivity&&) = delete;
IConnectivity(const IConnectivity&) = delete;
- ~IConnectivity() = default;
+ virtual ~IConnectivity();
};
template <>
diff --git a/src/mesh/IMesh.cpp b/src/mesh/IMesh.cpp
index cb2accc141ed09c778237587373494797fee726d..f315f4b74b01eb371c42eb4bf029f81e93ae6016 100644
--- a/src/mesh/IMesh.cpp
+++ b/src/mesh/IMesh.cpp
@@ -1,8 +1,10 @@
#include <mesh/IMesh.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
#include <mesh/MeshDataManager.hpp>
IMesh::~IMesh()
{
- MeshDataManager::instance().deleteMeshData(*this);
+ MeshDataManager::instance().deleteMeshData(this);
+ DiamondDualMeshManager::instance().deleteMesh(this);
}
diff --git a/src/mesh/ItemIdToItemIdMap.hpp b/src/mesh/ItemIdToItemIdMap.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0e92eb13ce3943c4e2487081ee2e48fc3deac8ee
--- /dev/null
+++ b/src/mesh/ItemIdToItemIdMap.hpp
@@ -0,0 +1,30 @@
+#ifndef ITEM_ID_TO_ITEM_ID_MAP_HPP
+#define ITEM_ID_TO_ITEM_ID_MAP_HPP
+
+#include <mesh/ItemId.hpp>
+#include <utils/Array.hpp>
+
+template <ItemType type1, ItemType type2>
+using ItemIdToItemIdMap = Array<std::pair<ItemIdT<type1>, ItemIdT<type2>>>;
+
+using NodeIdToNodeIdMap = ItemIdToItemIdMap<ItemType::node, ItemType::node>;
+using NodeIdToEdgeIdMap = ItemIdToItemIdMap<ItemType::node, ItemType::edge>;
+using NodeIdToFaceIdMap = ItemIdToItemIdMap<ItemType::node, ItemType::face>;
+using NodeIdToCellIdMap = ItemIdToItemIdMap<ItemType::node, ItemType::cell>;
+
+using EdgeIdToNodeIdMap = ItemIdToItemIdMap<ItemType::edge, ItemType::node>;
+using EdgeIdToEdgeIdMap = ItemIdToItemIdMap<ItemType::edge, ItemType::edge>;
+using EdgeIdToFaceIdMap = ItemIdToItemIdMap<ItemType::edge, ItemType::face>;
+using EdgeIdToCellIdMap = ItemIdToItemIdMap<ItemType::edge, ItemType::cell>;
+
+using FaceIdToNodeIdMap = ItemIdToItemIdMap<ItemType::face, ItemType::node>;
+using FaceIdToEdgeIdMap = ItemIdToItemIdMap<ItemType::face, ItemType::edge>;
+using FaceIdToFaceIdMap = ItemIdToItemIdMap<ItemType::face, ItemType::face>;
+using FaceIdToCellIdMap = ItemIdToItemIdMap<ItemType::face, ItemType::cell>;
+
+using CellIdToNodeIdMap = ItemIdToItemIdMap<ItemType::cell, ItemType::node>;
+using CellIdToEdgeIdMap = ItemIdToItemIdMap<ItemType::cell, ItemType::edge>;
+using CellIdToFaceIdMap = ItemIdToItemIdMap<ItemType::cell, ItemType::face>;
+using CellIdToCellIdMap = ItemIdToItemIdMap<ItemType::cell, ItemType::cell>;
+
+#endif // ITEM_ID_TO_ITEM_ID_MAP_HPP
diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
index e8d7071af5ae2c09b60c10478560087f5b2bdc86..1cb20d3ee595cf0963b95b4cc67f201cf2cddb5c 100644
--- a/src/mesh/ItemToItemMatrix.hpp
+++ b/src/mesh/ItemToItemMatrix.hpp
@@ -29,8 +29,10 @@ class ItemToItemMatrix
}
PUGS_INLINE
- TargetItemId operator[](size_t j) const
+ TargetItemId
+ operator[](size_t j) const
{
+ Assert(j < m_row.length);
return m_row(j);
}
@@ -73,15 +75,19 @@ class ItemToItemMatrix
}
PUGS_INLINE
- auto operator[](const SourceItemId& source_id) const
+ auto
+ operator[](const SourceItemId& source_id) const
{
+ Assert(source_id < m_connectivity_matrix.numRows());
using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
}
template <typename IndexType>
- PUGS_INLINE const auto& operator[](const IndexType& source_id) const
+ PUGS_INLINE const auto&
+ operator[](const IndexType& source_id) const
{
+ Assert(source_id < m_connectivity_matrix.numRows());
static_assert(std::is_same_v<IndexType, SourceItemId>, "ItemToItemMatrix must be indexed using correct ItemId");
using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 86e78f5f7fde25656882535258f6fb3779855d01..c82e9ce06ec8b7d43ce1c6cc2fec282f89dc1fee 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -32,36 +32,113 @@ class MeshData : public IMeshData
private:
const MeshType& m_mesh;
- std::shared_ptr<NodeValuePerCell<const Rd>> m_Cjr;
- std::shared_ptr<NodeValuePerCell<const double>> m_ljr;
- 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;
+ NodeValuePerFace<const Rd> m_Nlr;
+ NodeValuePerFace<const Rd> m_nlr;
+ NodeValuePerCell<const Rd> m_Cjr;
+ NodeValuePerCell<const double> m_ljr;
+ NodeValuePerCell<const Rd> m_njr;
+ FaceValue<const Rd> m_xl;
+ CellValue<const Rd> m_cell_centroid;
+ CellValue<const Rd> m_cell_iso_barycenter;
+ CellValue<const double> m_Vj;
+ FaceValue<const double> m_ll;
PUGS_INLINE
void
- _computeIsobarycenter()
+ _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
+ _compute_nlr()
+ {
+ if constexpr (Dimension == 1) {
+ static_assert(Dimension != 1, "nlr does not make sense in 1d");
+ } else {
+ const auto& Nlr = this->Nlr();
+
+ NodeValuePerFace<Rd> nlr{m_mesh.connectivity()};
+
+ parallel_for(
+ Nlr.numberOfValues(), PUGS_LAMBDA(size_t lr) {
+ double length = l2Norm(Nlr[lr]);
+ nlr[lr] = 1. / length * Nlr[lr];
+ });
+
+ m_nlr = nlr;
+ }
+ }
+
+ PUGS_INLINE
+ void
+ _computeFaceIsobarycenter()
+ { // Computes vertices isobarycenter
+ if constexpr (Dimension == 1) {
+ static_assert(Dimension != 1, "xl does not make sense in 1d");
+ } else {
+ const NodeValue<const Rd>& xr = m_mesh.xr();
+
+ const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+ FaceValue<Rd> xl(m_mesh.connectivity());
+ parallel_for(
+ m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ Rd X = zero;
+ const auto& face_nodes = face_to_node_matrix[face_id];
+ for (size_t R = 0; R < face_nodes.size(); ++R) {
+ X += xr[face_nodes[R]];
+ }
+ xl[face_id] = 1. / face_nodes.size() * X;
+ });
+ m_xl = xl;
+ }
+ }
+
+ PUGS_INLINE
+ void
+ _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 = std::make_shared<CellValue<const Rd>>(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;
@@ -69,9 +146,98 @@ 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 = std::make_shared<CellValue<const Rd>>(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;
+ }
}
}
@@ -95,7 +261,7 @@ class MeshData : public IMeshData
}
Vj[j] = inv_Dimension * sum_cjr_xr;
});
- m_Vj = std::make_shared<CellValue<const double>>(Vj);
+ m_Vj = Vj;
}
PUGS_INLINE
@@ -123,7 +289,7 @@ class MeshData : public IMeshData
Nlr(l, 0) = Nr;
Nlr(l, 1) = Nr;
});
- m_Nlr = std::make_shared<NodeValuePerFace<const Rd>>(Nlr);
+ m_Nlr = Nlr;
} else {
const NodeValue<const Rd>& xr = m_mesh.xr();
@@ -150,7 +316,7 @@ class MeshData : public IMeshData
Nlr(l, r) = Nr;
}
});
- m_Nlr = std::make_shared<NodeValuePerFace<const Rd>>(Nlr);
+ m_Nlr = Nlr;
}
}
@@ -165,7 +331,7 @@ class MeshData : public IMeshData
Cjr(j, 0) = -1;
Cjr(j, 1) = 1;
});
- m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
+ m_Cjr = Cjr;
} else if constexpr (Dimension == 2) {
const NodeValue<const Rd>& xr = m_mesh.xr();
const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
@@ -181,7 +347,7 @@ class MeshData : public IMeshData
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) {
auto Nlr = this->Nlr();
@@ -228,7 +394,7 @@ class MeshData : public IMeshData
}
});
- m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
+ m_Cjr = Cjr;
}
}
@@ -241,13 +407,13 @@ class MeshData : public IMeshData
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);
+ m_ljr = ljr;
} 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);
+ m_ljr = ljr;
}
}
@@ -265,19 +431,18 @@ class MeshData : public IMeshData
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);
+ m_njr = njr;
}
}
void
_checkCellVolume()
{
- Assert(m_Vj);
- auto& Vj = *m_Vj;
+ Assert(m_Vj.isBuilt());
bool is_valid = [&] {
for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
- if (Vj[j] <= 0) {
+ if (m_Vj[j] <= 0) {
return false;
}
}
@@ -297,65 +462,105 @@ class MeshData : public IMeshData
return m_mesh;
}
+ PUGS_INLINE
+ FaceValue<const double>
+ ll()
+ {
+ if (not m_ll.isBuilt()) {
+ this->_compute_ll();
+ }
+ return m_ll;
+ }
+
PUGS_INLINE
NodeValuePerFace<const Rd>
Nlr()
{
- if (not m_Nlr) {
+ if (not m_Nlr.isBuilt()) {
this->_computeNlr();
}
- return *m_Nlr;
+ return m_Nlr;
+ }
+
+ PUGS_INLINE
+ NodeValuePerFace<const Rd>
+ nlr()
+ {
+ if (not m_nlr.isBuilt()) {
+ this->_compute_nlr();
+ }
+ return m_nlr;
}
PUGS_INLINE
NodeValuePerCell<const Rd>
Cjr()
{
- if (not m_Cjr) {
+ if (not m_Cjr.isBuilt()) {
this->_computeCjr();
}
- return *m_Cjr;
+ return m_Cjr;
}
PUGS_INLINE
NodeValuePerCell<const double>
ljr()
{
- if (not m_ljr) {
+ if (not m_ljr.isBuilt()) {
this->_compute_ljr();
}
- return *m_ljr;
+ return m_ljr;
}
PUGS_INLINE
NodeValuePerCell<const Rd>
njr()
{
- if (not m_njr) {
+ if (not m_njr.isBuilt()) {
this->_compute_njr();
}
- return *m_njr;
+ return m_njr;
+ }
+
+ PUGS_INLINE
+ FaceValue<const Rd>
+ xl()
+ {
+ if (not m_xl.isBuilt()) {
+ this->_computeFaceIsobarycenter();
+ }
+ 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) {
- this->_computeIsobarycenter();
+ if (not m_cell_centroid.isBuilt()) {
+ this->_computeCellCentroid();
}
- return *m_xj;
+ return m_cell_centroid;
}
PUGS_INLINE
CellValue<const double>
Vj()
{
- if (not m_Vj) {
+ if (not m_Vj.isBuilt()) {
this->_computeCellVolume();
this->_checkCellVolume();
}
- return *m_Vj;
+ return m_Vj;
}
private:
diff --git a/src/mesh/MeshDataManager.cpp b/src/mesh/MeshDataManager.cpp
index 281e7517aa704e6c4352769e88f854b86cbc7366..8ca184a6013c083f5e8df7412d19f80b8f50f95a 100644
--- a/src/mesh/MeshDataManager.cpp
+++ b/src/mesh/MeshDataManager.cpp
@@ -35,9 +35,9 @@ MeshDataManager::destroy()
}
void
-MeshDataManager::deleteMeshData(const IMesh& mesh)
+MeshDataManager::deleteMeshData(const IMesh* p_mesh)
{
- m_mesh_mesh_data_map.erase(&mesh);
+ m_mesh_mesh_data_map.erase(p_mesh);
}
template <size_t Dimension>
diff --git a/src/mesh/MeshDataManager.hpp b/src/mesh/MeshDataManager.hpp
index 7b127148b889bc585082d626b79c435293fa4cd7..fd04101dea29a2d15e0feb2193f22f785ff65dc3 100644
--- a/src/mesh/MeshDataManager.hpp
+++ b/src/mesh/MeshDataManager.hpp
@@ -5,8 +5,8 @@
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
-#include <map>
#include <memory>
+#include <unordered_map>
class IMesh;
@@ -22,7 +22,7 @@ class MeshData;
class MeshDataManager
{
private:
- std::map<const IMesh*, std::shared_ptr<IMeshData>> m_mesh_mesh_data_map;
+ std::unordered_map<const IMesh*, std::shared_ptr<IMeshData>> m_mesh_mesh_data_map;
static MeshDataManager* m_instance;
@@ -44,7 +44,7 @@ class MeshDataManager
return *m_instance;
}
- void deleteMeshData(const IMesh&);
+ void deleteMeshData(const IMesh*);
template <size_t Dimension>
MeshData<Dimension>& getMeshData(const Mesh<Connectivity<Dimension>>&);
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index 9ce6369e07b63537ac698149703547aa107e83cd..66af694c7255da41ff7570d294dd4777acb0164f 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -54,16 +54,18 @@ class SubItemValuePerItem
const size_t m_size;
public:
- PUGS_INLINE
- const DataType& operator[](size_t i) const noexcept(NO_ASSERT)
+ template <typename IndexType>
+ PUGS_INLINE const DataType& operator[](IndexType i) const noexcept(NO_ASSERT)
{
+ static_assert(std::is_integral_v<IndexType>, "SubView is indexed by integral values");
Assert(i < m_size);
return m_sub_values[i];
}
- PUGS_FORCEINLINE
- DataType& operator[](size_t i) noexcept(NO_ASSERT)
+ template <typename IndexType>
+ PUGS_FORCEINLINE DataType& operator[](IndexType i) noexcept(NO_ASSERT)
{
+ static_assert(std::is_integral_v<IndexType>, "SubView is indexed by integral values");
Assert(i < m_size);
return m_sub_values[i];
}
@@ -96,23 +98,15 @@ class SubItemValuePerItem
return m_connectivity_ptr.use_count() != 0;
}
- // Following Kokkos logic, these classes are view and const view does allow
- // changes in data
- PUGS_FORCEINLINE
- DataType&
- operator()(ItemId j, size_t r) const noexcept(NO_ASSERT)
- {
- Assert(this->isBuilt());
- return m_values[m_host_row_map(size_t{j}) + r];
- }
-
- template <typename IndexType>
+ template <typename IndexType, typename SubIndexType>
PUGS_FORCEINLINE DataType&
- operator()(IndexType j, size_t r) const noexcept(NO_ASSERT)
+ operator()(IndexType item_id, SubIndexType i) const noexcept(NO_ASSERT)
{
- static_assert(std::is_same_v<IndexType, ItemId>, "SubItemValuePerItem indexed by ItemId");
+ static_assert(std::is_same_v<IndexType, ItemId>, "first index must be of the correct ItemId type");
+ static_assert(std::is_integral_v<SubIndexType>, "second index must be an integral type");
Assert(this->isBuilt());
- return m_values[m_host_row_map(size_t{j}) + r];
+ Assert(i + m_host_row_map(size_t{item_id}) < m_host_row_map(size_t{item_id} + 1));
+ return m_values[m_host_row_map(size_t{item_id}) + i];
}
PUGS_INLINE
@@ -125,19 +119,12 @@ class SubItemValuePerItem
// Following Kokkos logic, these classes are view and const view does allow
// changes in data
- PUGS_FORCEINLINE
- DataType& operator[](size_t i) const noexcept(NO_ASSERT)
+ template <typename ArrayIndexType>
+ DataType& operator[](const ArrayIndexType& i) const noexcept(NO_ASSERT)
{
+ static_assert(std::is_integral_v<ArrayIndexType>, "index must be an integral type");
Assert(this->isBuilt());
- return m_values[i];
- }
-
- template <typename IndexType>
- DataType& operator[](const IndexType& i) const noexcept(NO_ASSERT)
- {
- static_assert(std::is_same_v<IndexType, size_t>, "Access to SubItemValuePerItem's array must be indexed by "
- "size_t");
- Assert(this->isBuilt());
+ Assert(i < m_values.size());
return m_values[i];
}
@@ -150,34 +137,40 @@ class SubItemValuePerItem
return m_host_row_map.extent(0);
}
- PUGS_INLINE
- size_t
- numberOfSubValues(size_t i_cell) const noexcept(NO_ASSERT)
+ template <typename IndexType>
+ PUGS_INLINE size_t
+ numberOfSubValues(IndexType item_id) const noexcept(NO_ASSERT)
{
+ static_assert(std::is_same_v<IndexType, ItemId>, "index must be an ItemId");
Assert(this->isBuilt());
- return m_host_row_map(i_cell + 1) - m_host_row_map(i_cell);
+ Assert(item_id < m_host_row_map.extent(0));
+ return m_host_row_map(size_t{item_id} + 1) - m_host_row_map(size_t{item_id});
}
- PUGS_INLINE
- SubView
- itemValues(size_t i_cell) noexcept(NO_ASSERT)
+ template <typename IndexType>
+ PUGS_INLINE SubView
+ itemValues(IndexType item_id) noexcept(NO_ASSERT)
{
+ static_assert(std::is_same_v<IndexType, ItemId>, "index must be an ItemId");
Assert(this->isBuilt());
- const auto& cell_begin = m_host_row_map(i_cell);
- const auto& cell_end = m_host_row_map(i_cell + 1);
- return SubView(m_values, cell_begin, cell_end);
+ Assert(item_id < m_host_row_map.extent(0));
+ const auto& item_begin = m_host_row_map(size_t{item_id});
+ const auto& item_end = m_host_row_map(size_t{item_id} + 1);
+ return SubView(m_values, item_begin, item_end);
}
// Following Kokkos logic, these classes are view and const view does allow
// changes in data
- PUGS_INLINE
- SubView
- itemValues(size_t i_cell) const noexcept(NO_ASSERT)
+ template <typename IndexType>
+ PUGS_INLINE SubView
+ itemValues(IndexType item_id) const noexcept(NO_ASSERT)
{
+ static_assert(std::is_same_v<IndexType, ItemId>, "index must be an ItemId");
Assert(this->isBuilt());
- const auto& cell_begin = m_host_row_map(i_cell);
- const auto& cell_end = m_host_row_map(i_cell + 1);
- return SubView(m_values, cell_begin, cell_end);
+ Assert(item_id < m_host_row_map.extent(0));
+ const auto& item_begin = m_host_row_map(size_t{item_id});
+ const auto& item_end = m_host_row_map(size_t{item_id} + 1);
+ return SubView(m_values, item_begin, item_end);
}
template <typename DataType2, typename ConnectivityPtr2>
diff --git a/src/mesh/SynchronizerManager.hpp b/src/mesh/SynchronizerManager.hpp
index 79cb3db0a3473bdbe9f88d4d20836a835aca9260..8bcb441f05c4e841a71c103f9a252bb52c59dcb2 100644
--- a/src/mesh/SynchronizerManager.hpp
+++ b/src/mesh/SynchronizerManager.hpp
@@ -4,8 +4,8 @@
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
-#include <map>
#include <memory>
+#include <unordered_map>
class IConnectivity;
class Synchronizer;
@@ -13,7 +13,7 @@ class Synchronizer;
class SynchronizerManager
{
private:
- std::map<const IConnectivity*, std::shared_ptr<Synchronizer>> m_connectivity_synchronizer_map;
+ std::unordered_map<const IConnectivity*, std::shared_ptr<Synchronizer>> m_connectivity_synchronizer_map;
static SynchronizerManager* m_instance;
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 5d85b3630e2afc77ef7f2a3422b2f7d207abe5b7..f3c202ffc51633d291d36150d662fe950023ef58 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -8,12 +8,14 @@
#include <output/OutputNamedItemValueSet.hpp>
#include <utils/Exceptions.hpp>
#include <utils/Messenger.hpp>
+#include <utils/PugsAssert.hpp>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
+#include <unordered_map>
class VTKWriter
{
@@ -234,6 +236,10 @@ class VTKWriter
std::visit([&, name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, name, item_value); },
item_value_variant);
}
+ if constexpr (MeshType::Dimension == 3) {
+ fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faces\"/>\n";
+ fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\"/>\n";
+ }
fout << "</PCells>\n";
fout << "<PPointData>\n";
@@ -315,7 +321,9 @@ class VTKWriter
{
Array<int8_t> types(mesh->numberOfCells());
- const auto& cell_type = mesh->connectivity().cellType();
+ const auto& cell_type = mesh->connectivity().cellType();
+ const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
switch (cell_type[j]) {
@@ -336,7 +344,11 @@ class VTKWriter
break;
}
case CellType::Pyramid: {
- types[j] = 14;
+ if (cell_to_node_matrix[j].size() == 5) {
+ types[j] = 14; // quadrangle basis
+ } else {
+ types[j] = 41; // polygonal basis
+ }
break;
}
case CellType::Prism: {
@@ -347,6 +359,10 @@ class VTKWriter
types[j] = 12;
break;
}
+ case CellType::Diamond: {
+ types[j] = 41;
+ break;
+ }
default: {
std::ostringstream os;
os << __FILE__ << ':' << __LINE__ << ": unknown cell type";
@@ -355,6 +371,70 @@ class VTKWriter
}
});
_write_array(fout, "types", types);
+ if constexpr (MeshType::Dimension == 3) {
+ const bool has_general_polyhedron = [&] {
+ for (size_t i = 0; i < types.size(); ++i) {
+ if (types[i] == 41)
+ return true;
+ }
+ return false;
+ }();
+ if (has_general_polyhedron) {
+ const auto& cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
+ const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+ const auto& cell_face_is_reversed = mesh->connectivity().cellFaceIsReversed();
+
+ Array<size_t> faces_offsets(mesh->numberOfCells());
+ size_t next_offset = 0;
+ fout << "<DataArray type=\"Int64\" IdType=\"1\" Name=\"faces\">\n";
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ std::unordered_map<size_t, size_t> node_id_to_local_number_in_cell;
+ for (size_t i_cell_node = 0; i_cell_node < cell_nodes.size(); ++i_cell_node) {
+ node_id_to_local_number_in_cell[cell_nodes[i_cell_node]] = i_cell_node;
+ }
+
+ const auto& cell_faces = cell_to_face_matrix[cell_id];
+ fout << cell_faces.size() << '\n';
+ next_offset++;
+ for (size_t i_cell_face = 0; i_cell_face < cell_faces.size(); ++i_cell_face) {
+ const FaceId& face_id = cell_faces[i_cell_face];
+ const auto& face_nodes = face_to_node_matrix[face_id];
+ fout << face_nodes.size();
+ next_offset++;
+ Array<size_t> face_node_in_cell(face_nodes.size());
+
+ for (size_t i_face_node = 0; i_face_node < face_nodes.size(); ++i_face_node) {
+ const NodeId& node_id = face_nodes[i_face_node];
+ auto i_node_id_to_local_number_in_cell = node_id_to_local_number_in_cell.find(node_id);
+ Assert(i_node_id_to_local_number_in_cell != node_id_to_local_number_in_cell.end());
+ face_node_in_cell[i_face_node] = i_node_id_to_local_number_in_cell->second;
+ }
+
+ if (cell_face_is_reversed(cell_id, i_cell_face)) {
+ for (size_t i = 0; i < face_node_in_cell.size(); ++i) {
+ fout << ' ' << face_node_in_cell[face_node_in_cell.size() - 1 - i];
+ }
+ } else {
+ for (size_t i = 0; i < face_node_in_cell.size(); ++i) {
+ fout << ' ' << face_node_in_cell[i];
+ }
+ }
+
+ next_offset += face_nodes.size();
+
+ fout << '\n';
+ }
+ faces_offsets[cell_id] = next_offset;
+ }
+ fout << "</DataArray>\n";
+ fout << "<DataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\">\n";
+ for (size_t i_face_offsets = 0; i_face_offsets < faces_offsets.size(); ++i_face_offsets) {
+ fout << faces_offsets[i_face_offsets] << '\n';
+ }
+ fout << "</DataArray>\n";
+ }
+ }
}
fout << "</Cells>\n";
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index d0eb4d0b699709e0cddff74035e6b2263b7fffd2..3e3151d25429593955cb281c2dd814690f67bfa5 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -1,47 +1,170 @@
#ifndef ACOUSTIC_SOLVER_HPP
#define ACOUSTIC_SOLVER_HPP
-#include <rang.hpp>
-
-#include <utils/ArrayUtils.hpp>
-
-#include <scheme/BlockPerfectGas.hpp>
-#include <utils/PugsAssert.hpp>
-
-#include <scheme/BoundaryCondition.hpp>
-#include <scheme/FiniteVolumesEulerUnknowns.hpp>
-
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
-#include <mesh/MeshDataManager.hpp>
-
#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
-
#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/SubItemValuePerItem.hpp>
-
+#include <scheme/AcousticSolverType.hpp>
+#include <scheme/BlockPerfectGas.hpp>
+#include <scheme/FiniteVolumesEulerUnknowns.hpp>
+#include <utils/ArrayUtils.hpp>
#include <utils/Exceptions.hpp>
#include <utils/Messenger.hpp>
+#include <utils/PugsAssert.hpp>
+
+#include <rang.hpp>
#include <iostream>
template <typename MeshType>
class AcousticSolver
{
+ public:
+ class PressureBoundaryCondition;
+ class SymmetryBoundaryCondition;
+ class VelocityBoundaryCondition;
+
+ using BoundaryCondition =
+ std::variant<PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
constexpr static size_t Dimension = MeshType::Dimension;
+ private:
using MeshDataType = MeshData<Dimension>;
using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>;
std::shared_ptr<const MeshType> m_mesh;
const typename MeshType::Connectivity& m_connectivity;
- const std::vector<BoundaryConditionHandler>& m_boundary_condition_list;
using Rd = TinyVector<Dimension>;
using Rdd = TinyMatrix<Dimension>;
- private:
+ const BoundaryConditionList m_boundary_condition_list;
+
+ const AcousticSolverType m_solver_type;
+
+ void
+ _applyPressureBC()
+ {
+ for (const auto& boundary_condition : m_boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
+ MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+ if constexpr (Dimension == 1) {
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+ const auto& node_to_cell_matrix = m_connectivity.nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
+
+ const auto& node_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+ const auto& node_cell_list = node_to_cell_matrix[node_id];
+ Assert(node_cell_list.size() == 1);
+
+ CellId node_cell_id = node_cell_list[0];
+ size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
+
+ m_br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
+ });
+ } else {
+ const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+
+ const auto& face_to_cell_matrix = m_connectivity.faceToCellMatrix();
+ const auto& face_to_node_matrix = m_connectivity.faceToNodeMatrix();
+ const auto& face_local_numbers_in_their_cells = m_connectivity.faceLocalNumbersInTheirCells();
+ const auto& face_cell_is_reversed = m_connectivity.cellFaceIsReversed();
+
+ const auto& face_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const auto& face_cell_list = face_to_cell_matrix[face_id];
+ Assert(face_cell_list.size() == 1);
+
+ CellId face_cell_id = face_cell_list[0];
+ size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+ const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
+
+ const auto& face_nodes = face_to_node_matrix[face_id];
+
+ for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ NodeId node_id = face_nodes[i_node];
+ m_br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _applySymmetryBC()
+ {
+ for (const auto& boundary_condition : m_boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+ const Rd& n = bc.outgoingNormal();
+
+ const Rdd I = identity;
+ const Rdd nxn = tensorProduct(n, n);
+ const Rdd P = I - nxn;
+
+ const Array<const NodeId>& node_list = bc.nodeList();
+ parallel_for(
+ bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
+ const NodeId r = node_list[r_number];
+
+ m_Ar[r] = P * m_Ar[r] * P + nxn;
+ m_br[r] = P * m_br[r];
+ });
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _applyVelocityBC()
+ {
+ for (const auto& boundary_condition : m_boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+ const auto& value_list = bc.valueList();
+
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ NodeId node_id = node_list[i_node];
+ const auto& value = value_list[i_node];
+
+ m_Ar[node_id] = identity;
+ m_br[node_id] = value;
+ });
+ }
+ },
+ boundary_condition);
+ }
+ }
+
PUGS_INLINE
const CellValue<const double>
computeRhoCj(const CellValue<const double>& rhoj, const CellValue<const double>& cj)
@@ -53,11 +176,13 @@ class AcousticSolver
PUGS_INLINE
void
- computeAjr(const CellValue<const double>& rhocj,
- const NodeValuePerCell<const Rd>& Cjr,
- const NodeValuePerCell<const double>& /* ljr */,
- const NodeValuePerCell<const Rd>& njr)
+ _computeGlaceAjr(const CellValue<const double>& rhocj)
{
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+ const NodeValuePerCell<const Rd> njr = mesh_data.njr();
+
parallel_for(
m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
const size_t& nb_nodes = m_Ajr.numberOfSubValues(j);
@@ -68,6 +193,71 @@ class AcousticSolver
});
}
+ PUGS_INLINE
+ void
+ _computeEucclhydAjr(const CellValue<const double>& rhocj)
+ {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+ const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
+
+ const auto& face_to_node_matrix = m_connectivity.faceToNodeMatrix();
+ const auto& cell_to_node_matrix = m_connectivity.cellToNodeMatrix();
+ const auto& cell_to_face_matrix = m_connectivity.cellToFaceMatrix();
+
+ parallel_for(
+ m_Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { m_Ajr[jr] = zero; });
+
+ 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 double& rho_c = rhocj[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];
+
+ 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();
+ };
+
+ for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
+ m_Ajr(j, R) += tensorProduct(rho_c * Nlr(l, rl), nlr(l, rl));
+ }
+ }
+ });
+ }
+
+ PUGS_INLINE
+ void
+ computeAjr(const CellValue<const double>& rhocj)
+ {
+ switch (m_solver_type) {
+ case AcousticSolverType::Glace: {
+ this->_computeGlaceAjr(rhocj);
+ break;
+ }
+ case AcousticSolverType::Eucclhyd: {
+ if constexpr (Dimension > 1) {
+ this->_computeEucclhydAjr(rhocj);
+ } else {
+ this->_computeGlaceAjr(rhocj);
+ }
+ break;
+ }
+ }
+ }
+
PUGS_INLINE
const NodeValue<const Rdd>
computeAr(const NodeValuePerCell<const Rdd>& Ajr)
@@ -94,11 +284,13 @@ class AcousticSolver
PUGS_INLINE
const NodeValue<const Rd>
- computeBr(const NodeValuePerCell<Rdd>& Ajr,
- const NodeValuePerCell<const Rd>& Cjr,
- const CellValue<const Rd>& uj,
- const CellValue<const double>& pj)
+ computeBr(const CellValue<const Rd>& uj, const CellValue<const double>& pj)
{
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
+ const NodeValuePerCell<const Rdd>& Ajr = m_Ajr;
+
const auto& node_to_cell_matrix = m_connectivity.nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
@@ -121,58 +313,33 @@ class AcousticSolver
void
applyBoundaryConditions()
{
- for (const auto& handler : m_boundary_condition_list) {
- switch (handler.boundaryCondition().type()) {
- case BoundaryCondition::normal_velocity: {
- throw NotImplementedError("normal_velocity BC");
- }
- case BoundaryCondition::velocity: {
- throw NotImplementedError("velocity BC");
- }
- case BoundaryCondition::pressure: {
- throw NotImplementedError("pressure BC");
- }
- case BoundaryCondition::symmetry: {
- const SymmetryBoundaryCondition<Dimension>& symmetry_bc =
- dynamic_cast<const SymmetryBoundaryCondition<Dimension>&>(handler.boundaryCondition());
- const Rd& n = symmetry_bc.outgoingNormal();
-
- const Rdd I = identity;
- const Rdd nxn = tensorProduct(n, n);
- const Rdd P = I - nxn;
-
- const Array<const NodeId>& node_list = symmetry_bc.nodeList();
- parallel_for(
- symmetry_bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
- const NodeId r = node_list[r_number];
-
- m_Ar[r] = P * m_Ar[r] * P + nxn;
- m_br[r] = P * m_br[r];
- });
- break;
- }
- }
- }
+ this->_applyPressureBC();
+ this->_applySymmetryBC();
+ this->_applyVelocityBC();
}
- NodeValue<Rd>
- computeUr(const NodeValue<const Rdd>& Ar, const NodeValue<const Rd>& br)
+ void
+ computeUr()
{
+ const NodeValue<const Rdd> Ar = m_Ar;
+ const NodeValue<const Rd> br = m_br;
+
inverse(Ar, m_inv_Ar);
const NodeValue<const Rdd> invAr = m_inv_Ar;
parallel_for(
m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { m_ur[r] = invAr[r] * br[r]; });
-
- return m_ur;
}
void
- computeFjr(const NodeValuePerCell<Rdd>& Ajr,
- const NodeValue<const Rd>& ur,
- const NodeValuePerCell<const Rd>& Cjr,
- const CellValue<const Rd>& uj,
- const CellValue<const double>& pj)
+ computeFjr(const CellValue<const Rd>& uj, const CellValue<const double>& pj)
{
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+ const NodeValuePerCell<const Rdd> Ajr = m_Ajr;
+
+ const NodeValue<const Rd> ur = m_ur;
+
const auto& cell_to_node_matrix = m_mesh->connectivity().cellToNodeMatrix();
parallel_for(
@@ -197,26 +364,22 @@ class AcousticSolver
computeExplicitFluxes(const CellValue<const double>& rhoj,
const CellValue<const Rd>& uj,
const CellValue<const double>& pj,
- const CellValue<const double>& cj,
- const NodeValuePerCell<const Rd>& Cjr,
- const NodeValuePerCell<const double>& ljr,
- const NodeValuePerCell<const Rd>& njr)
+ const CellValue<const double>& cj)
{
const CellValue<const double> rhocj = computeRhoCj(rhoj, cj);
- computeAjr(rhocj, Cjr, ljr, njr);
+ computeAjr(rhocj);
NodeValuePerCell<const Rdd> Ajr = m_Ajr;
this->computeAr(Ajr);
- this->computeBr(m_Ajr, Cjr, uj, pj);
+ this->computeBr(uj, pj);
this->applyBoundaryConditions();
synchronize(m_Ar);
synchronize(m_br);
- NodeValue<Rd>& ur = m_ur;
- ur = computeUr(m_Ar, m_br);
- computeFjr(m_Ajr, ur, Cjr, uj, pj);
+ computeUr();
+ computeFjr(uj, pj);
}
NodeValue<Rd> m_br;
@@ -229,10 +392,13 @@ class AcousticSolver
CellValue<double> m_Vj_over_cj;
public:
- AcousticSolver(std::shared_ptr<const MeshType> p_mesh, const std::vector<BoundaryConditionHandler>& bc_list)
+ AcousticSolver(std::shared_ptr<const MeshType> p_mesh,
+ const BoundaryConditionList& bc_list,
+ const AcousticSolverType solver_type)
: m_mesh(p_mesh),
m_connectivity(m_mesh->connectivity()),
m_boundary_condition_list(bc_list),
+ m_solver_type(solver_type),
m_br(m_connectivity),
m_Ajr(m_connectivity),
m_Ar(m_connectivity),
@@ -281,12 +447,9 @@ class AcousticSolver
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();
- const NodeValuePerCell<const double> ljr = mesh_data.ljr();
- const NodeValuePerCell<const Rd> njr = mesh_data.njr();
+ const CellValue<const double> Vj = mesh_data.Vj();
- computeExplicitFluxes(rhoj, uj, pj, cj, Cjr, ljr, njr);
+ computeExplicitFluxes(rhoj, uj, pj, cj);
const NodeValuePerCell<Rd>& Fjr = m_Fjr;
const NodeValue<const Rd> ur = m_ur;
@@ -326,4 +489,124 @@ class AcousticSolver
}
};
+template <typename MeshType>
+class AcousticSolver<MeshType>::PressureBoundaryCondition
+{
+ private:
+ const Array<const double> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ PressureBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}
+ {}
+
+ ~PressureBoundaryCondition() = default;
+};
+
+template <>
+class AcousticSolver<Mesh<Connectivity<1>>>::PressureBoundaryCondition
+{
+ private:
+ const Array<const double> m_value_list;
+ const Array<const NodeId> m_face_list;
+
+ public:
+ const Array<const NodeId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ PressureBoundaryCondition(const Array<const NodeId>& face_list, const Array<const double>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}
+ {}
+
+ ~PressureBoundaryCondition() = default;
+};
+
+template <typename MeshType>
+class AcousticSolver<MeshType>::VelocityBoundaryCondition
+{
+ private:
+ const Array<const TinyVector<Dimension>> m_value_list;
+ const Array<const NodeId> m_node_list;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_node_list;
+ }
+
+ const Array<const TinyVector<Dimension>>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ VelocityBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list)
+ : m_value_list{value_list}, m_node_list{node_list}
+ {}
+
+ ~VelocityBoundaryCondition() = default;
+};
+
+template <typename MeshType>
+class AcousticSolver<MeshType>::SymmetryBoundaryCondition
+{
+ public:
+ static constexpr size_t Dimension = MeshType::Dimension;
+
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+
+ public:
+ const Rd&
+ outgoingNormal() const
+ {
+ return m_mesh_flat_node_boundary.outgoingNormal();
+ }
+
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_flat_node_boundary.nodeList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_flat_node_boundary.nodeList();
+ }
+
+ SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
+ : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+ {
+ ;
+ }
+
+ ~SymmetryBoundaryCondition() = default;
+};
+
#endif // ACOUSTIC_SOLVER_HPP
diff --git a/src/scheme/AcousticSolverType.hpp b/src/scheme/AcousticSolverType.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3095dcc864d5598e985290302ad94106adbd6b2e
--- /dev/null
+++ b/src/scheme/AcousticSolverType.hpp
@@ -0,0 +1,10 @@
+#ifndef ACOUSTIC_SOLVER_TYPE_HPP
+#define ACOUSTIC_SOLVER_TYPE_HPP
+
+enum class AcousticSolverType
+{
+ Eucclhyd,
+ Glace
+};
+
+#endif // ACOUSTIC_SOLVER_TYPE_HPP
diff --git a/src/scheme/BoundaryCondition.hpp b/src/scheme/BoundaryCondition.hpp
deleted file mode 100644
index 3d3332fc92cfa88792ecf556cb97781de007a104..0000000000000000000000000000000000000000
--- a/src/scheme/BoundaryCondition.hpp
+++ /dev/null
@@ -1,146 +0,0 @@
-#ifndef BOUNDARY_CONDITION_HANDLER_HPP
-#define BOUNDARY_CONDITION_HANDLER_HPP
-
-#include <utils/Array.hpp>
-
-#include <mesh/MeshNodeBoundary.hpp>
-#include <mesh/RefItemList.hpp>
-
-#include <memory>
-#include <vector>
-
-class BoundaryCondition // clazy:exclude=copyable-polymorphic
-{
- public:
- enum Type
- {
- velocity,
- normal_velocity,
- pressure,
- symmetry
- };
-
- const Type m_type;
-
- protected:
- BoundaryCondition(Type type) : m_type(type)
- {
- ;
- }
-
- public:
- const Type&
- type() const
- {
- return m_type;
- }
-
- virtual ~BoundaryCondition() = default;
-};
-
-class PressureBoundaryCondition final : public BoundaryCondition // clazy:exclude=copyable-polymorphic
-{
- private:
- const double m_value;
-
- const size_t m_number_of_faces;
- Array<const unsigned int> m_face_list;
-
- public:
- double
- value() const
- {
- return m_value;
- }
-
- const size_t&
- numberOfFaces() const
- {
- return m_number_of_faces;
- }
-
- const Array<const unsigned int>&
- faceList() const
- {
- return m_face_list;
- }
-
- PressureBoundaryCondition(double value, const std::vector<unsigned int>& faces)
- : BoundaryCondition(BoundaryCondition::pressure), m_value(value), m_number_of_faces(faces.size())
- {
- Array<unsigned int> face_list(faces.size());
- parallel_for(
- m_number_of_faces, PUGS_LAMBDA(int f) { face_list[f] = faces[f]; });
- m_face_list = face_list;
- }
-
- ~PressureBoundaryCondition() = default;
-};
-
-template <size_t dimension>
-class SymmetryBoundaryCondition : public BoundaryCondition
-{
- public:
- using Rd = TinyVector<dimension, double>;
-
- private:
- const MeshFlatNodeBoundary<dimension> m_mesh_flat_node_boundary;
-
- public:
- const Rd&
- outgoingNormal() const
- {
- return m_mesh_flat_node_boundary.outgoingNormal();
- }
-
- size_t
- numberOfNodes() const
- {
- return m_mesh_flat_node_boundary.nodeList().size();
- }
-
- const Array<const NodeId>&
- nodeList() const
- {
- return m_mesh_flat_node_boundary.nodeList();
- }
-
- SymmetryBoundaryCondition(const MeshFlatNodeBoundary<dimension>& mesh_flat_node_boundary)
- : BoundaryCondition(BoundaryCondition::symmetry), m_mesh_flat_node_boundary(mesh_flat_node_boundary)
- {
- ;
- }
-
- ~SymmetryBoundaryCondition() = default;
-};
-
-class BoundaryConditionHandler
-{
- private:
- std::shared_ptr<BoundaryCondition> m_boundary_condition;
-
- public:
- const BoundaryCondition&
- boundaryCondition() const
- {
- return *m_boundary_condition;
- }
-
- PUGS_INLINE
- BoundaryConditionHandler& operator=(BoundaryConditionHandler&) = default;
-
- PUGS_INLINE
- BoundaryConditionHandler(const BoundaryConditionHandler&) = default;
-
- PUGS_INLINE
- BoundaryConditionHandler(BoundaryConditionHandler&&) = default;
-
- PUGS_INLINE
- BoundaryConditionHandler(std::shared_ptr<BoundaryCondition> boundary_condition)
- : m_boundary_condition(boundary_condition)
- {}
-
- ~BoundaryConditionHandler() = default;
-};
-
-#endif // BOUNDARY_CONDITION_HANDLER_HPP
diff --git a/src/scheme/DirichletBoundaryConditionDescriptor.hpp b/src/scheme/DirichletBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5fa6c42f9d4b5333ad1ca5ef2fdb46a4f06b4242
--- /dev/null
+++ b/src/scheme/DirichletBoundaryConditionDescriptor.hpp
@@ -0,0 +1,64 @@
+#ifndef DIRICHLET_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define DIRICHLET_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class DirichletBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+ std::ostream&
+ _write(std::ostream& os) const final
+ {
+ os << "dirichlet(" << m_name << ',' << *m_boundary_descriptor << ")";
+ return os;
+ }
+
+ const std::string_view m_name;
+
+ std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+ const FunctionSymbolId m_rhs_symbol_id;
+
+ public:
+ std::string_view
+ name() const
+ {
+ return m_name;
+ }
+
+ FunctionSymbolId
+ rhsSymbolId() const
+ {
+ return m_rhs_symbol_id;
+ }
+
+ const IBoundaryDescriptor&
+ boundaryDescriptor() const
+ {
+ return *m_boundary_descriptor;
+ }
+
+ Type
+ type() const final
+ {
+ return Type::dirichlet;
+ }
+
+ DirichletBoundaryConditionDescriptor(const std::string_view name,
+ std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+ const FunctionSymbolId& rhs_symbol_id)
+ : m_name{name}, m_boundary_descriptor(boundary_descriptor), m_rhs_symbol_id{rhs_symbol_id}
+ {
+ ;
+ }
+
+ DirichletBoundaryConditionDescriptor(const DirichletBoundaryConditionDescriptor&) = delete;
+ DirichletBoundaryConditionDescriptor(DirichletBoundaryConditionDescriptor&&) = delete;
+
+ ~DirichletBoundaryConditionDescriptor() = default;
+};
+
+#endif // DIRICHLET_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/FourierBoundaryConditionDescriptor.hpp b/src/scheme/FourierBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e7aeb830c3447bcf9783060210bc97e59bf2ebdf
--- /dev/null
+++ b/src/scheme/FourierBoundaryConditionDescriptor.hpp
@@ -0,0 +1,75 @@
+#ifndef FOURIER_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define FOURIER_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class FourierBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+ std::ostream&
+ _write(std::ostream& os) const final
+ {
+ os << "fourier(" << m_name << ',' << *m_boundary_descriptor << ")";
+ return os;
+ }
+
+ const std::string_view m_name;
+
+ std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+ const FunctionSymbolId m_mass_symbol_id;
+ const FunctionSymbolId m_rhs_symbol_id;
+
+ public:
+ std::string_view
+ name() const
+ {
+ return m_name;
+ }
+
+ FunctionSymbolId
+ massSymbolId() const
+ {
+ return m_mass_symbol_id;
+ }
+
+ FunctionSymbolId
+ rhsSymbolId() const
+ {
+ return m_rhs_symbol_id;
+ }
+
+ const IBoundaryDescriptor&
+ boundaryDescriptor() const
+ {
+ return *m_boundary_descriptor;
+ }
+
+ Type
+ type() const final
+ {
+ return Type::fourier;
+ }
+
+ FourierBoundaryConditionDescriptor(const std::string_view name,
+ std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+ const FunctionSymbolId& mass_symbol_id,
+ const FunctionSymbolId& rhs_symbol_id)
+ : m_name{name},
+ m_boundary_descriptor(boundary_descriptor),
+ m_mass_symbol_id{mass_symbol_id},
+ m_rhs_symbol_id{rhs_symbol_id}
+ {
+ ;
+ }
+
+ FourierBoundaryConditionDescriptor(const FourierBoundaryConditionDescriptor&) = delete;
+ FourierBoundaryConditionDescriptor(FourierBoundaryConditionDescriptor&&) = delete;
+
+ ~FourierBoundaryConditionDescriptor() = default;
+};
+
+#endif // FOURIER_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/FreeBoundaryConditionDescriptor.hpp b/src/scheme/FreeBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b0eb2ceb218d77cae347a44219f13d278a25aee5
--- /dev/null
+++ b/src/scheme/FreeBoundaryConditionDescriptor.hpp
@@ -0,0 +1,46 @@
+#ifndef FREE_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define FREE_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class FreeBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+ std::ostream&
+ _write(std::ostream& os) const final
+ {
+ os << "free(" << *m_boundary_descriptor << ")";
+ return os;
+ }
+
+ std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+
+ public:
+ const IBoundaryDescriptor&
+ boundaryDescriptor() const
+ {
+ return *m_boundary_descriptor;
+ }
+
+ Type
+ type() const final
+ {
+ return Type::free;
+ }
+
+ FreeBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
+ : m_boundary_descriptor(boundary_descriptor)
+ {
+ ;
+ }
+
+ FreeBoundaryConditionDescriptor(const FreeBoundaryConditionDescriptor&) = delete;
+ FreeBoundaryConditionDescriptor(FreeBoundaryConditionDescriptor&&) = delete;
+
+ ~FreeBoundaryConditionDescriptor() = default;
+};
+
+#endif // FREE_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index 0814c9e8580f6c36ad1152c3143d5877bff45dbc..30a25db66f1278e72eb2bea0ce7a017d8022da65 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -8,6 +8,10 @@ class IBoundaryConditionDescriptor
public:
enum class Type
{
+ dirichlet,
+ fourier,
+ free,
+ neumann,
symmetry
};
diff --git a/src/scheme/NeumannBoundaryConditionDescriptor.hpp b/src/scheme/NeumannBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..edc1b1fb35b3f530c47adcd2f0073c9727707965
--- /dev/null
+++ b/src/scheme/NeumannBoundaryConditionDescriptor.hpp
@@ -0,0 +1,64 @@
+#ifndef NEUMANN_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define NEUMANN_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class NeumannBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+ std::ostream&
+ _write(std::ostream& os) const final
+ {
+ os << "neumann(" << m_name << ',' << *m_boundary_descriptor << ")";
+ return os;
+ }
+
+ const std::string_view m_name;
+
+ std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+ const FunctionSymbolId m_rhs_symbol_id;
+
+ public:
+ std::string_view
+ name() const
+ {
+ return m_name;
+ }
+
+ FunctionSymbolId
+ rhsSymbolId() const
+ {
+ return m_rhs_symbol_id;
+ }
+
+ const IBoundaryDescriptor&
+ boundaryDescriptor() const
+ {
+ return *m_boundary_descriptor;
+ }
+
+ Type
+ type() const final
+ {
+ return Type::neumann;
+ }
+
+ NeumannBoundaryConditionDescriptor(const std::string_view name,
+ std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+ const FunctionSymbolId& rhs_symbol_id)
+ : m_name{name}, m_boundary_descriptor(boundary_descriptor), m_rhs_symbol_id{rhs_symbol_id}
+ {
+ ;
+ }
+
+ NeumannBoundaryConditionDescriptor(const NeumannBoundaryConditionDescriptor&) = delete;
+ NeumannBoundaryConditionDescriptor(NeumannBoundaryConditionDescriptor&&) = delete;
+
+ ~NeumannBoundaryConditionDescriptor() = default;
+};
+
+#endif // NEUMANN_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 2dfd3234486b352e5b9f27e13e10001846008bc1..7cd1ed89023e619af478625b58441002ecae0af2 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -90,6 +90,7 @@ target_link_libraries (unit_tests
PugsLanguage
PugsLanguageAST
PugsLanguageModules
+ PugsLanguageAlgorithms
PugsLanguageUtils
PugsMesh
PugsUtils
diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp
index b7cb2d3bd8afe280a7d8c416a662016851e0d820..99b94515ba7795903434db761d509b9052fcaf71 100644
--- a/tests/test_BiCGStab.cpp
+++ b/tests/test_BiCGStab.cpp
@@ -45,7 +45,7 @@ TEST_CASE("BiCGStab", "[algebra]")
Vector<double> x{5};
x = 0;
- BiCGStab{b, A, x, 10, 1e-12};
+ BiCGStab<false>{b, A, x, 10, 1e-12};
Vector error = x - x_exact;
REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
}
diff --git a/tests/test_PCG.cpp b/tests/test_PCG.cpp
index 4b2f2df9bed4bec0ea31502be3c5480e1eec6989..d4564761d8c564082e4668360e5a2e2a8ecce008 100644
--- a/tests/test_PCG.cpp
+++ b/tests/test_PCG.cpp
@@ -62,7 +62,7 @@ TEST_CASE("PCG", "[algebra]")
Vector<double> x{5};
x = 0;
- PCG{b, A, A, x, 10, 1e-12};
+ PCG<false>{b, A, A, x, 10, 1e-12};
REQUIRE(std::sqrt((x, x)) == 0);
}
@@ -104,7 +104,7 @@ TEST_CASE("PCG", "[algebra]")
Vector<double> x{5};
x = 0;
- PCG{b, A, A, x, 1, 1e-12};
+ PCG<false>{b, A, A, x, 1, 1e-12};
Vector error = x - x_exact;
REQUIRE(std::sqrt((error, error)) > 1E-10 * std::sqrt((x, x)));
}
diff --git a/tests/test_SparseMatrixDescriptor.cpp b/tests/test_SparseMatrixDescriptor.cpp
index 4e8b90134060c786cc7f97bc8fc8cd00bfba3257..cd66380b1d1beb6bcb9ba1a9c0d9c495947a3235 100644
--- a/tests/test_SparseMatrixDescriptor.cpp
+++ b/tests/test_SparseMatrixDescriptor.cpp
@@ -168,4 +168,30 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
REQUIRE(value_array[6] == 1);
REQUIRE(value_array[7] == -2);
}
+
+ SECTION("output")
+ {
+ SparseMatrixDescriptor<int, uint8_t> S{5};
+ S(0, 2) = 3;
+ S(1, 2) = 11;
+ S(1, 1) = 1;
+ S(0, 2) += 2;
+ S(3, 3) = 5;
+ S(3, 1) = -3;
+ S(4, 1) = 1;
+ S(2, 2) = 4;
+ S(4, 4) = -2;
+
+ std::ostringstream output;
+ output << '\n' << S;
+
+ std::string expected_output = R"(
+0 | 2:5
+1 | 1:1 2:11
+2 | 2:4
+3 | 1:-3 3:5
+4 | 1:1 4:-2
+)";
+ REQUIRE(output.str() == expected_output);
+ }
}