diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 92d354fc52ddfd3e35ab69a8d892edb32f6432de..5846781f5b0cdb8d471e9cd50a66c3ef0fe911f4 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -1,23 +1,16 @@
#ifndef MESH_NODE_BOUNDARY_HPP
#define MESH_NODE_BOUNDARY_HPP
-#include <utils/Array.hpp>
-
+#include <Kokkos_Vector.hpp>
#include <algebra/TinyVector.hpp>
-
-#include <mesh/ItemValue.hpp>
-#include <mesh/RefItemList.hpp>
-
#include <mesh/ConnectivityMatrix.hpp>
#include <mesh/IConnectivity.hpp>
-
+#include <mesh/ItemValue.hpp>
+#include <mesh/RefItemList.hpp>
+#include <utils/Array.hpp>
#include <utils/Exceptions.hpp>
#include <utils/Messenger.hpp>
-#include <Kokkos_Vector.hpp>
-
-#include <iostream>
-
template <size_t Dimension>
class MeshNodeBoundary // clazy:exclude=copyable-polymorphic
{
@@ -35,10 +28,10 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic
}
template <typename MeshType>
- MeshNodeBoundary(const std::shared_ptr<const MeshType>& mesh, const RefFaceList& ref_face_list)
+ MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
{
static_assert(Dimension == MeshType::Dimension);
- const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
+ const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
const Array<const FaceId>& face_list = ref_face_list.list();
parallel_for(
@@ -52,7 +45,7 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic
Kokkos::vector<unsigned int> node_ids;
// not enough but should reduce significantly the number of resizing
node_ids.reserve(Dimension * face_list.size());
- const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+ const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
for (size_t l = 0; l < face_list.size(); ++l) {
const FaceId face_number = face_list[l];
@@ -73,8 +66,7 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic
}
template <typename MeshType>
- MeshNodeBoundary(const std::shared_ptr<const MeshType>&, const RefNodeList& ref_node_list)
- : m_node_list(ref_node_list.list())
+ MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list) : m_node_list(ref_node_list.list())
{
static_assert(Dimension == MeshType::Dimension);
}
@@ -97,16 +89,16 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension> // clazy:exclu
const Rd m_outgoing_normal;
template <typename MeshType>
- PUGS_INLINE Rd _getNormal(const std::shared_ptr<const MeshType>& mesh);
+ PUGS_INLINE Rd _getNormal(const MeshType& mesh);
template <typename MeshType>
PUGS_INLINE void _checkBoundaryIsFlat(TinyVector<2, double> normal,
TinyVector<2, double> xmin,
TinyVector<2, double> xmax,
- const std::shared_ptr<const MeshType>& mesh) const;
+ const MeshType& mesh) const;
template <typename MeshType>
- PUGS_INLINE Rd _getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh);
+ PUGS_INLINE Rd _getOutgoingNormal(const MeshType& mesh);
public:
const Rd&
@@ -119,14 +111,14 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension> // clazy:exclu
MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
template <typename MeshType>
- MeshFlatNodeBoundary(const std::shared_ptr<const MeshType>& mesh, const RefFaceList& ref_face_list)
+ MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
: MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh))
{
;
}
template <typename MeshType>
- MeshFlatNodeBoundary(const std::shared_ptr<const MeshType>& mesh, const RefNodeList& ref_node_list)
+ MeshFlatNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
: MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh))
{
;
@@ -144,7 +136,7 @@ void
MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(TinyVector<2, double> normal,
TinyVector<2, double> xmin,
TinyVector<2, double> xmax,
- const std::shared_ptr<const MeshType>& mesh) const
+ const MeshType& mesh) const
{
static_assert(MeshType::Dimension == 2);
using R2 = TinyVector<2, double>;
@@ -152,7 +144,7 @@ MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(TinyVector<2, double> normal,
const R2 origin = 0.5 * (xmin + xmax);
const double length = l2Norm(xmax - xmin);
- const NodeValue<const R2>& xr = mesh->xr();
+ const NodeValue<const R2>& xr = mesh.xr();
parallel_for(
m_node_list.size(), PUGS_LAMBDA(size_t r) {
@@ -166,14 +158,14 @@ MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(TinyVector<2, double> normal,
template <>
template <typename MeshType>
PUGS_INLINE TinyVector<1, double>
-MeshFlatNodeBoundary<1>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
{
static_assert(MeshType::Dimension == 1);
using R = TinyVector<1, double>;
const size_t number_of_bc_nodes = [&]() {
size_t number_of_bc_nodes = 0;
- auto node_is_owned = mesh->connectivity().nodeIsOwned();
+ auto node_is_owned = mesh.connectivity().nodeIsOwned();
for (size_t i_node = 0; i_node < m_node_list.size(); ++i_node) {
number_of_bc_nodes += (node_is_owned[m_node_list[i_node]]);
}
@@ -190,12 +182,12 @@ MeshFlatNodeBoundary<1>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
template <>
template <typename MeshType>
PUGS_INLINE TinyVector<2, double>
-MeshFlatNodeBoundary<2>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
{
static_assert(MeshType::Dimension == 2);
using R2 = TinyVector<2, double>;
- const NodeValue<const R2>& xr = mesh->xr();
+ const NodeValue<const R2>& xr = mesh.xr();
R2 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
@@ -245,7 +237,7 @@ MeshFlatNodeBoundary<2>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
template <>
template <typename MeshType>
PUGS_INLINE TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
{
static_assert(MeshType::Dimension == 3);
using R3 = TinyVector<3, double>;
@@ -258,7 +250,7 @@ MeshFlatNodeBoundary<3>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
R3 ymax = xmax;
R3 zmax = xmax;
- const NodeValue<const R3>& xr = mesh->xr();
+ const NodeValue<const R3>& xr = mesh.xr();
for (size_t r = 0; r < m_node_list.size(); ++r) {
const R3& x = xr[m_node_list[r]];
@@ -365,7 +357,7 @@ MeshFlatNodeBoundary<3>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
template <>
template <typename MeshType>
PUGS_INLINE TinyVector<1, double>
-MeshFlatNodeBoundary<1>::_getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<1>::_getOutgoingNormal(const MeshType& mesh)
{
static_assert(MeshType::Dimension == 1);
using R = TinyVector<1, double>;
@@ -375,10 +367,10 @@ MeshFlatNodeBoundary<1>::_getOutgoingNormal(const std::shared_ptr<const MeshType
double max_height = 0;
if (m_node_list.size() > 0) {
- const NodeValue<const R>& xr = mesh->xr();
- const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ const NodeValue<const R>& xr = mesh.xr();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const NodeId r0 = m_node_list[0];
const CellId j0 = node_to_cell_matrix[r0][0];
@@ -410,7 +402,7 @@ MeshFlatNodeBoundary<1>::_getOutgoingNormal(const std::shared_ptr<const MeshType
template <>
template <typename MeshType>
PUGS_INLINE TinyVector<2, double>
-MeshFlatNodeBoundary<2>::_getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<2>::_getOutgoingNormal(const MeshType& mesh)
{
static_assert(MeshType::Dimension == 2);
using R2 = TinyVector<2, double>;
@@ -420,10 +412,10 @@ MeshFlatNodeBoundary<2>::_getOutgoingNormal(const std::shared_ptr<const MeshType
double max_height = 0;
if (m_node_list.size() > 0) {
- const NodeValue<const R2>& xr = mesh->xr();
- const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ const NodeValue<const R2>& xr = mesh.xr();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const NodeId r0 = m_node_list[0];
const CellId j0 = node_to_cell_matrix[r0][0];
@@ -454,7 +446,7 @@ MeshFlatNodeBoundary<2>::_getOutgoingNormal(const std::shared_ptr<const MeshType
template <>
template <typename MeshType>
PUGS_INLINE TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<3>::_getOutgoingNormal(const MeshType& mesh)
{
static_assert(MeshType::Dimension == 3);
using R3 = TinyVector<3, double>;
@@ -464,10 +456,10 @@ MeshFlatNodeBoundary<3>::_getOutgoingNormal(const std::shared_ptr<const MeshType
double max_height = 0;
if (m_node_list.size() > 0) {
- const NodeValue<const R3>& xr = mesh->xr();
- const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ const NodeValue<const R3>& xr = mesh.xr();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const NodeId r0 = m_node_list[0];
const CellId j0 = node_to_cell_matrix[r0][0];
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 6973037571b23f861cec5c70ed28fbe67c036fcf..0ad9b7901ceb6dd39cab12d817201e307d035308 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -18,15 +18,14 @@ template <size_t Dimension>
double
acoustic_dt(const DiscreteFunctionP0<Dimension, double>& c)
{
- const std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh =
- std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(c.mesh());
+ const Mesh<Connectivity<Dimension>>& mesh = dynamic_cast<const Mesh<Connectivity<Dimension>>&>(*c.mesh());
- const auto Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
- const auto Sj = MeshDataManager::instance().getMeshData(*mesh).sumOverRLjr();
+ const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+ const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
- CellValue<double> local_dt{mesh->connectivity()};
+ CellValue<double> local_dt{mesh.connectivity()};
parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); });
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); });
return min(local_dt);
}
@@ -89,11 +88,11 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
{
Assert(rho.mesh() == c.mesh());
- std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*rho.mesh());
- CellValue<double> rhoc{mesh->connectivity()};
+ CellValue<double> rhoc{mesh.connectivity()};
parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rhoc[cell_id] = rho[cell_id] * c[cell_id]; });
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rhoc[cell_id] = rho[cell_id] * c[cell_id]; });
return rhoc;
}
@@ -247,7 +246,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
}
BoundaryConditionList
- _getBCList(const std::shared_ptr<const MeshType>& mesh,
+ _getBCList(const MeshType& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
BoundaryConditionList bc_list;
@@ -267,9 +266,9 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
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 < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
- const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ 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()) {
SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)};
@@ -285,8 +284,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
if (dirichlet_bc_descriptor.name() == "velocity") {
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);
+ 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};
@@ -296,15 +295,15 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
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);
+ 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") {
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);
+ 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();
@@ -314,9 +313,9 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
Array<const double> face_values = [&] {
if constexpr (Dimension == 1) {
return InterpolateItemValue<double(
- TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh->xr(), face_list);
+ TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh.xr(), face_list);
} else {
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
return InterpolateItemValue<double(
TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(), face_list);
@@ -394,16 +393,14 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
}
AcousticSolver(SolverType solver_type,
- const std::shared_ptr<const MeshType>& p_mesh,
+ const MeshType& mesh,
const DiscreteScalarFunction& rho,
const DiscreteScalarFunction& c,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
- : m_solver_type{solver_type}, m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)}
+ : m_solver_type{solver_type}, m_boundary_condition_list{this->_getBCList(mesh, bc_descriptor_list)}
{
- const MeshType& mesh = *p_mesh;
-
CellValue<const double> rhoc = this->_getRhoC(rho, c);
NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(mesh, rhoc);
@@ -425,27 +422,27 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
std::shared_ptr<const DiscreteFunctionP0<Dimension, Rd>>,
std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>>
apply(const double& dt,
- const std::shared_ptr<const MeshType>& mesh,
+ const MeshType& mesh,
const DiscreteFunctionP0<Dimension, double>& rho,
const DiscreteFunctionP0<Dimension, Rd>& u,
const DiscreteFunctionP0<Dimension, double>& E) const
{
- const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- NodeValue<Rd> new_xr = copy(mesh->xr());
+ NodeValue<Rd> new_xr = copy(mesh.xr());
parallel_for(
- mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * m_ur[r]; });
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * m_ur[r]; });
- std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr);
+ std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
- CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+ CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
CellValue<double> new_rho = copy(rho.cellValues());
CellValue<Rd> new_u = copy(u.cellValues());
CellValue<double> new_E = copy(E.cellValues());
parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
@@ -463,7 +460,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
return {new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_rho),
std::make_shared<DiscreteVectorFunction>(new_mesh, new_u),
@@ -488,7 +485,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
throw NormalError("acoustic solver expects P0 functions");
}
- return this->apply(dt, std::dynamic_pointer_cast<const MeshType>(i_mesh),
+ return this->apply(dt, dynamic_cast<const MeshType&>(*i_mesh),
*std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho),
*std::dynamic_pointer_cast<const DiscreteVectorFunction>(u),
*std::dynamic_pointer_cast<const DiscreteScalarFunction>(E));
@@ -502,7 +499,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
const std::shared_ptr<const IDiscreteFunction>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
: AcousticSolver{solver_type,
- std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(mesh),
+ dynamic_cast<const MeshType&>(*mesh),
*std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho),
*std::dynamic_pointer_cast<const DiscreteScalarFunction>(c),
*std::dynamic_pointer_cast<const DiscreteVectorFunction>(u),