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

Replace std::shared<const MeshType> -> const MeshType 

It was actually unnecessary to use a shared pointer at this point and
was quite disturbing at use since it enforced use of shared pointers
in caller methods (which was a bad design)
parent 557a1192
No related branches found
No related tags found
1 merge request!96Add random number engine encapsulation.
Show whitespace changes
Inline Side-by-side
src/mesh/MeshNodeBoundary.hpp
+ 33
− 41
View file @ 0db09f3c
#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];
......
src/scheme/AcousticSolver.cpp
+ 31
− 34
View file @ 0db09f3c
......@@ -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),
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment