From 54b65f6fd8368f385b099c957a2174151a598230 Mon Sep 17 00:00:00 2001 From: HOCH PHILIPPE <philippe.hoch@gmail.com> Date: Thu, 14 Nov 2024 21:53:32 +0100 Subject: [PATCH] merge of current dev --- src/language/modules/SchemeModule.cpp | 10 + src/scheme/IBoundaryConditionDescriptor.hpp | 3 +- ...eViscousFormEulerianCompositeSolver_v2.cpp | 173 +++++++++++++++-- ...scousFormEulerianCompositeSolver_v2_o2.cpp | 174 +++++++++++++++--- src/scheme/RusanovEulerianCompositeSolver.cpp | 168 +++++++++++++++-- .../RusanovEulerianCompositeSolver_o2.cpp | 164 +++++++++++++++-- .../RusanovEulerianCompositeSolver_v2.cpp | 163 ++++++++++++++-- .../RusanovEulerianCompositeSolver_v2_o2.cpp | 163 ++++++++++++++-- .../WallBoundaryConditionDescriptor.hpp | 40 ++++ 9 files changed, 946 insertions(+), 112 deletions(-) create mode 100644 src/scheme/WallBoundaryConditionDescriptor.hpp diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index 19beaad93..047ce78cd 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -51,6 +51,8 @@ #include <scheme/RusanovEulerianCompositeSolver_v2_o2.hpp> #include <scheme/SymmetryBoundaryConditionDescriptor.hpp> #include <scheme/VariableBCDescriptor.hpp> +#include <scheme/WallBoundaryConditionDescriptor.hpp> +#include <scheme/test_reconstruction.hpp> #include <utils/Socket.hpp> #include <memory> @@ -341,6 +343,14 @@ SchemeModule::SchemeModule() } )); + this->_addBuiltinFunction("wall", std::function( + + [](std::shared_ptr<const IBoundaryDescriptor> boundary) + -> std::shared_ptr<const IBoundaryConditionDescriptor> { + return std::make_shared<WallBoundaryConditionDescriptor>(boundary); + } + + )); this->_addBuiltinFunction("pressure", std::function( diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp index 14967ae3f..fc9ff381b 100644 --- a/src/scheme/IBoundaryConditionDescriptor.hpp +++ b/src/scheme/IBoundaryConditionDescriptor.hpp @@ -21,7 +21,8 @@ class IBoundaryConditionDescriptor inflow_list, neumann, outflow, - symmetry + symmetry, + wall }; protected: diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp index 3490fa4f1..dab114826 100644 --- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp +++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp @@ -36,14 +36,14 @@ class RoeViscousFormEulerianCompositeSolver_v2 class SymmetryBoundaryCondition; class InflowListBoundaryCondition; class OutflowBoundaryCondition; - class NeumannBoundaryCondition; + class WallBoundaryCondition; class NeumannflatBoundaryCondition; using BoundaryCondition = std::variant<SymmetryBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition, NeumannflatBoundaryCondition, - NeumannBoundaryCondition>; + WallBoundaryCondition>; using BoundaryConditionList = std::vector<BoundaryCondition>; @@ -57,17 +57,15 @@ class RoeViscousFormEulerianCompositeSolver_v2 bool is_valid_boundary_condition = true; switch (bc_descriptor->type()) { - case IBoundaryConditionDescriptor::Type::neumann: { + case IBoundaryConditionDescriptor::Type::wall: { if constexpr (Dimension == 2) { - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } else { static_assert(Dimension == 3); - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } break; } @@ -163,7 +161,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 } if (not is_valid_boundary_condition) { std::ostringstream error_msg; - error_msg << *bc_descriptor << " is an invalid boundary condition for Rusanov v2 Eulerian Composite solver"; + error_msg << *bc_descriptor << " is an invalid boundary condition for Roe v2 Eulerian Composite solver"; throw NormalError(error_msg.str()); } } @@ -503,6 +501,138 @@ class RoeViscousFormEulerianCompositeSolver_v2 } } + void + _applyWallBoundaryCondition(const BoundaryConditionList& bc_list, + const MeshType& mesh, + NodeValuePerCell<Rp>& stateNode, + EdgeValuePerCell<Rp>& stateEdge, + FaceValuePerCell<Rp>& stateFace) const + { + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<WallBoundaryCondition, T>) { + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + std::cout << " Traitement WALL local (non flat) \n"; + // const Rd& normal = bc.outgoingNormal(); + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); + // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& face_list = bc.faceList(); + const auto& node_list = bc.nodeList(); + + const auto Cjr = mesh_data.Cjr(); + const auto Cjf = mesh_data.Cjf(); + + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + const NodeId node_id = node_list[i_node]; + + const auto& node_cell_list = node_to_cell_matrix[node_id]; + // Assert(face_cell_list.size() == 1); + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + normal += Cjr(node_cell_id, node_local_number_in_cell); + } + normal *= 1. / node_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim]; + } + } + + 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); + + // Normal locale approchée + Rd normal(Cjf(face_cell_id, face_local_number_in_cell)); + normal *= 1. / l2Norm(normal); + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); + + const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); + + const auto& edge_list = bc.edgeList(); + + const auto Cje = mesh_data.Cje(); + + for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { + const EdgeId edge_id = edge_list[i_edge]; + + const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; + + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + normal += Cje(edge_cell_id, edge_local_number_in_cell); + } + normal *= 1. / edge_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + } + } + } + }, + boundary_condition); + } + } + void _applyInflowBoundaryCondition(const BoundaryConditionList& bc_list, const MeshType& mesh, @@ -1031,7 +1161,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 _applyInflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); //_applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); _applySymmetricBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); - _applyNeumannflatBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + _applyWallBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); // // Construction du flux .. ok pour l'ordre 1 @@ -1630,8 +1760,11 @@ class RoeViscousFormEulerianCompositeSolver_v2 double theta = 2. / 3.; //.5; double eta = 1. / 6.; //.2; if constexpr (Dimension == 2) { - eta = 0; + eta = 0; + theta = 1; // pour schema aux aretes + // theta=0; //pour schema aux noeuds } + // else { // theta = 1. / 3.; // eta = 1. / 3.; @@ -1685,12 +1818,12 @@ class RoeViscousFormEulerianCompositeSolver_v2 }; template <MeshConcept MeshType> -class RoeViscousFormEulerianCompositeSolver_v2<MeshType>::NeumannBoundaryCondition +class RoeViscousFormEulerianCompositeSolver_v2<MeshType>::WallBoundaryCondition { }; template <> -class RoeViscousFormEulerianCompositeSolver_v2<Mesh<2>>::NeumannBoundaryCondition +class RoeViscousFormEulerianCompositeSolver_v2<Mesh<2>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1725,7 +1858,7 @@ class RoeViscousFormEulerianCompositeSolver_v2<Mesh<2>>::NeumannBoundaryConditio return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary) { ; @@ -1733,7 +1866,7 @@ class RoeViscousFormEulerianCompositeSolver_v2<Mesh<2>>::NeumannBoundaryConditio }; template <> -class RoeViscousFormEulerianCompositeSolver_v2<Mesh<3>>::NeumannBoundaryCondition +class RoeViscousFormEulerianCompositeSolver_v2<Mesh<3>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1778,9 +1911,9 @@ class RoeViscousFormEulerianCompositeSolver_v2<Mesh<3>>::NeumannBoundaryConditio return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, - const MeshEdgeBoundary& mesh_edge_boundary, - const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, + const MeshEdgeBoundary& mesh_edge_boundary, + const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_edge_boundary(mesh_edge_boundary), diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp index 84ec1dd2b..371129c55 100644 --- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp +++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp @@ -40,14 +40,15 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 class SymmetryBoundaryCondition; class InflowListBoundaryCondition; class OutflowBoundaryCondition; - class NeumannBoundaryCondition; + // class NeumannBoundaryCondition; + class WallBoundaryCondition; class NeumannflatBoundaryCondition; using BoundaryCondition = std::variant<SymmetryBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition, NeumannflatBoundaryCondition, - NeumannBoundaryCondition>; + WallBoundaryCondition>; using BoundaryConditionList = std::vector<BoundaryCondition>; @@ -61,17 +62,15 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 bool is_valid_boundary_condition = true; switch (bc_descriptor->type()) { - case IBoundaryConditionDescriptor::Type::neumann: { + case IBoundaryConditionDescriptor::Type::wall: { if constexpr (Dimension == 2) { - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } else { static_assert(Dimension == 3); - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } break; } @@ -167,7 +166,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 } if (not is_valid_boundary_condition) { std::ostringstream error_msg; - error_msg << *bc_descriptor << " is an invalid boundary condition for Rusanov v2 Eulerian Composite solver"; + error_msg << *bc_descriptor << " is an invalid boundary condition for Roe v2 Eulerian Composite solver"; throw NormalError(error_msg.str()); } } @@ -191,7 +190,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 [&](auto&& bc) { using T = std::decay_t<decltype(bc)>; if constexpr (std::is_same_v<OutflowBoundaryCondition, T>) { - std::cout << " Traitement Outflow \n"; + std::cout << " Traitement OUTFLOW \n"; // const Rd& normal = bc.outgoingNormal(); /* const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); @@ -510,6 +509,138 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 } } + void + _applyWallBoundaryCondition(const BoundaryConditionList& bc_list, + const MeshType& mesh, + NodeValuePerCell<Rp>& stateNode, + EdgeValuePerCell<Rp>& stateEdge, + FaceValuePerCell<Rp>& stateFace) const + { + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<WallBoundaryCondition, T>) { + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + std::cout << " Traitement WALL local (non flat) \n"; + // const Rd& normal = bc.outgoingNormal(); + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); + // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& face_list = bc.faceList(); + const auto& node_list = bc.nodeList(); + + const auto Cjr = mesh_data.Cjr(); + const auto Cjf = mesh_data.Cjf(); + + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + const NodeId node_id = node_list[i_node]; + + const auto& node_cell_list = node_to_cell_matrix[node_id]; + // Assert(face_cell_list.size() == 1); + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + normal += Cjr(node_cell_id, node_local_number_in_cell); + } + normal *= 1. / node_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim]; + } + } + + 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); + + // Normal locale approchée + Rd normal(Cjf(face_cell_id, face_local_number_in_cell)); + normal *= 1. / l2Norm(normal); + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); + + const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); + + const auto& edge_list = bc.edgeList(); + + const auto Cje = mesh_data.Cje(); + + for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { + const EdgeId edge_id = edge_list[i_edge]; + + const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; + + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + normal += Cje(edge_cell_id, edge_local_number_in_cell); + } + normal *= 1. / edge_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + } + } + } + }, + boundary_condition); + } + } + void _applyInflowBoundaryCondition(const BoundaryConditionList& bc_list, const MeshType& mesh, @@ -1039,7 +1170,8 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 DiscreteFunctionDPk eps_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>(); // DiscreteFunctionDPk rho_u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>(); - // DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>(); + // DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const + // double>>(); CellValue<double> Limitor_rho(p_mesh->connectivity()); CellValue<double> Limitor_u(p_mesh->connectivity()); @@ -1203,7 +1335,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 _applyInflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); //_applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); _applySymmetricBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); - _applyNeumannflatBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + _applyWallBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); // // Construction du flux .. ok pour l'ordre 1 @@ -1857,12 +1989,12 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 }; template <MeshConcept MeshType> -class RoeViscousFormEulerianCompositeSolver_v2_o2<MeshType>::NeumannBoundaryCondition +class RoeViscousFormEulerianCompositeSolver_v2_o2<MeshType>::WallBoundaryCondition { }; template <> -class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<2>>::NeumannBoundaryCondition +class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<2>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1897,7 +2029,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<2>>::NeumannBoundaryCondi return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary) { ; @@ -1905,7 +2037,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<2>>::NeumannBoundaryCondi }; template <> -class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<3>>::NeumannBoundaryCondition +class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<3>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1950,9 +2082,9 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<3>>::NeumannBoundaryCondi return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, - const MeshEdgeBoundary& mesh_edge_boundary, - const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, + const MeshEdgeBoundary& mesh_edge_boundary, + const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_edge_boundary(mesh_edge_boundary), diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp index 05ce4a9aa..ce0d8355c 100644 --- a/src/scheme/RusanovEulerianCompositeSolver.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver.cpp @@ -34,14 +34,14 @@ class RusanovEulerianCompositeSolver class SymmetryBoundaryCondition; class InflowListBoundaryCondition; class OutflowBoundaryCondition; - class NeumannBoundaryCondition; + class WallBoundaryCondition; class NeumannflatBoundaryCondition; using BoundaryCondition = std::variant<SymmetryBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition, NeumannflatBoundaryCondition, - NeumannBoundaryCondition>; + WallBoundaryCondition>; using BoundaryConditionList = std::vector<BoundaryCondition>; @@ -55,17 +55,15 @@ class RusanovEulerianCompositeSolver bool is_valid_boundary_condition = true; switch (bc_descriptor->type()) { - case IBoundaryConditionDescriptor::Type::neumann: { + case IBoundaryConditionDescriptor::Type::wall: { if constexpr (Dimension == 2) { - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } else { static_assert(Dimension == 3); - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } break; } @@ -170,6 +168,138 @@ class RusanovEulerianCompositeSolver } public: + void + _applyWallBoundaryCondition(const BoundaryConditionList& bc_list, + const MeshType& mesh, + NodeValuePerCell<Rp>& stateNode, + EdgeValuePerCell<Rp>& stateEdge, + FaceValuePerCell<Rp>& stateFace) const + { + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<WallBoundaryCondition, T>) { + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + std::cout << " Traitement WALL local (non flat) \n"; + // const Rd& normal = bc.outgoingNormal(); + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); + // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& face_list = bc.faceList(); + const auto& node_list = bc.nodeList(); + + const auto Cjr = mesh_data.Cjr(); + const auto Cjf = mesh_data.Cjf(); + + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + const NodeId node_id = node_list[i_node]; + + const auto& node_cell_list = node_to_cell_matrix[node_id]; + // Assert(face_cell_list.size() == 1); + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + normal += Cjr(node_cell_id, node_local_number_in_cell); + } + normal *= 1. / node_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim]; + } + } + + 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); + + // Normal locale approchée + Rd normal(Cjf(face_cell_id, face_local_number_in_cell)); + normal *= 1. / l2Norm(normal); + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); + + const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); + + const auto& edge_list = bc.edgeList(); + + const auto Cje = mesh_data.Cje(); + + for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { + const EdgeId edge_id = edge_list[i_edge]; + + const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; + + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + normal += Cje(edge_cell_id, edge_local_number_in_cell); + } + normal *= 1. / edge_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + } + } + } + }, + boundary_condition); + } + } + void _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list, const MeshType& mesh, @@ -654,8 +784,9 @@ class RusanovEulerianCompositeSolver // Conditions aux limites des etats conservatifs _applyInflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); - //_applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + _applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); _applySymmetricBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + _applyWallBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); // // Construction du flux .. ok pour l'ordre 1 @@ -1260,14 +1391,13 @@ class RusanovEulerianCompositeSolver RusanovEulerianCompositeSolver() = default; ~RusanovEulerianCompositeSolver() = default; }; - template <MeshConcept MeshType> -class RusanovEulerianCompositeSolver<MeshType>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver<MeshType>::WallBoundaryCondition { }; template <> -class RusanovEulerianCompositeSolver<Mesh<2>>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver<Mesh<2>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1300,7 +1430,7 @@ class RusanovEulerianCompositeSolver<Mesh<2>>::NeumannBoundaryCondition return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary) { ; @@ -1308,7 +1438,7 @@ class RusanovEulerianCompositeSolver<Mesh<2>>::NeumannBoundaryCondition }; template <> -class RusanovEulerianCompositeSolver<Mesh<3>>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver<Mesh<3>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1353,9 +1483,9 @@ class RusanovEulerianCompositeSolver<Mesh<3>>::NeumannBoundaryCondition return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, - const MeshEdgeBoundary& mesh_edge_boundary, - const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, + const MeshEdgeBoundary& mesh_edge_boundary, + const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_edge_boundary(mesh_edge_boundary), diff --git a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp index 954dd8784..35c5dec36 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp @@ -38,14 +38,14 @@ class RusanovEulerianCompositeSolver_o2 class SymmetryBoundaryCondition; class InflowListBoundaryCondition; class OutflowBoundaryCondition; - class NeumannBoundaryCondition; + class WallBoundaryCondition; class NeumannflatBoundaryCondition; using BoundaryCondition = std::variant<SymmetryBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition, NeumannflatBoundaryCondition, - NeumannBoundaryCondition>; + WallBoundaryCondition>; using BoundaryConditionList = std::vector<BoundaryCondition>; @@ -59,17 +59,15 @@ class RusanovEulerianCompositeSolver_o2 bool is_valid_boundary_condition = true; switch (bc_descriptor->type()) { - case IBoundaryConditionDescriptor::Type::neumann: { + case IBoundaryConditionDescriptor::Type::wall: { if constexpr (Dimension == 2) { - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } else { static_assert(Dimension == 3); - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } break; } @@ -177,6 +175,138 @@ class RusanovEulerianCompositeSolver_o2 CellByCellLimitation<MeshType> Limitor; public: + void + _applyWallBoundaryCondition(const BoundaryConditionList& bc_list, + const MeshType& mesh, + NodeValuePerCell<Rp>& stateNode, + EdgeValuePerCell<Rp>& stateEdge, + FaceValuePerCell<Rp>& stateFace) const + { + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<WallBoundaryCondition, T>) { + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + std::cout << " Traitement WALL local (non flat) \n"; + // const Rd& normal = bc.outgoingNormal(); + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); + // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& face_list = bc.faceList(); + const auto& node_list = bc.nodeList(); + + const auto Cjr = mesh_data.Cjr(); + const auto Cjf = mesh_data.Cjf(); + + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + const NodeId node_id = node_list[i_node]; + + const auto& node_cell_list = node_to_cell_matrix[node_id]; + // Assert(face_cell_list.size() == 1); + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + normal += Cjr(node_cell_id, node_local_number_in_cell); + } + normal *= 1. / node_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim]; + } + } + + 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); + + // Normal locale approchée + Rd normal(Cjf(face_cell_id, face_local_number_in_cell)); + normal *= 1. / l2Norm(normal); + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); + + const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); + + const auto& edge_list = bc.edgeList(); + + const auto Cje = mesh_data.Cje(); + + for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { + const EdgeId edge_id = edge_list[i_edge]; + + const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; + + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + normal += Cje(edge_cell_id, edge_local_number_in_cell); + } + normal *= 1. / edge_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + } + } + } + }, + boundary_condition); + } + } + void _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list, const MeshType& mesh, @@ -1431,12 +1561,12 @@ class RusanovEulerianCompositeSolver_o2 }; template <MeshConcept MeshType> -class RusanovEulerianCompositeSolver_o2<MeshType>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver_o2<MeshType>::WallBoundaryCondition { }; template <> -class RusanovEulerianCompositeSolver_o2<Mesh<2>>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver_o2<Mesh<2>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1469,7 +1599,7 @@ class RusanovEulerianCompositeSolver_o2<Mesh<2>>::NeumannBoundaryCondition return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary) { ; @@ -1477,7 +1607,7 @@ class RusanovEulerianCompositeSolver_o2<Mesh<2>>::NeumannBoundaryCondition }; template <> -class RusanovEulerianCompositeSolver_o2<Mesh<3>>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver_o2<Mesh<3>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1522,9 +1652,9 @@ class RusanovEulerianCompositeSolver_o2<Mesh<3>>::NeumannBoundaryCondition return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, - const MeshEdgeBoundary& mesh_edge_boundary, - const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, + const MeshEdgeBoundary& mesh_edge_boundary, + const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_edge_boundary(mesh_edge_boundary), diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp index c00f2311b..ea8ca8d01 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp @@ -33,14 +33,14 @@ class RusanovEulerianCompositeSolver_v2 class SymmetryBoundaryCondition; class InflowListBoundaryCondition; class OutflowBoundaryCondition; - class NeumannBoundaryCondition; + class WallBoundaryCondition; class NeumannflatBoundaryCondition; using BoundaryCondition = std::variant<SymmetryBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition, NeumannflatBoundaryCondition, - NeumannBoundaryCondition>; + WallBoundaryCondition>; using BoundaryConditionList = std::vector<BoundaryCondition>; @@ -54,17 +54,15 @@ class RusanovEulerianCompositeSolver_v2 bool is_valid_boundary_condition = true; switch (bc_descriptor->type()) { - case IBoundaryConditionDescriptor::Type::neumann: { + case IBoundaryConditionDescriptor::Type::wall: { if constexpr (Dimension == 2) { - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } else { static_assert(Dimension == 3); - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } break; } @@ -169,6 +167,137 @@ class RusanovEulerianCompositeSolver_v2 } public: + void + _applyWallBoundaryCondition(const BoundaryConditionList& bc_list, + const MeshType& mesh, + NodeValuePerCell<Rp>& stateNode, + EdgeValuePerCell<Rp>& stateEdge, + FaceValuePerCell<Rp>& stateFace) const + { + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<WallBoundaryCondition, T>) { + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + std::cout << " Traitement WALL local (non flat) \n"; + // const Rd& normal = bc.outgoingNormal(); + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); + // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& face_list = bc.faceList(); + const auto& node_list = bc.nodeList(); + + const auto Cjr = mesh_data.Cjr(); + const auto Cjf = mesh_data.Cjf(); + + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + const NodeId node_id = node_list[i_node]; + + const auto& node_cell_list = node_to_cell_matrix[node_id]; + // Assert(face_cell_list.size() == 1); + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + normal += Cjr(node_cell_id, node_local_number_in_cell); + } + normal *= 1. / node_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim]; + } + } + + 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); + + // Normal locale approchée + Rd normal(Cjf(face_cell_id, face_local_number_in_cell)); + normal *= 1. / l2Norm(normal); + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); + + const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); + + const auto& edge_list = bc.edgeList(); + + const auto Cje = mesh_data.Cje(); + + for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { + const EdgeId edge_id = edge_list[i_edge]; + + const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; + + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + normal += Cje(edge_cell_id, edge_local_number_in_cell); + } + normal *= 1. / edge_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + } + } + } + }, + boundary_condition); + } + } void _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list, const MeshType& mesh, @@ -1167,12 +1296,12 @@ class RusanovEulerianCompositeSolver_v2 }; template <MeshConcept MeshType> -class RusanovEulerianCompositeSolver_v2<MeshType>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver_v2<MeshType>::WallBoundaryCondition { }; template <> -class RusanovEulerianCompositeSolver_v2<Mesh<2>>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver_v2<Mesh<2>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1207,7 +1336,7 @@ class RusanovEulerianCompositeSolver_v2<Mesh<2>>::NeumannBoundaryCondition return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary) { ; @@ -1215,7 +1344,7 @@ class RusanovEulerianCompositeSolver_v2<Mesh<2>>::NeumannBoundaryCondition }; template <> -class RusanovEulerianCompositeSolver_v2<Mesh<3>>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver_v2<Mesh<3>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1260,9 +1389,9 @@ class RusanovEulerianCompositeSolver_v2<Mesh<3>>::NeumannBoundaryCondition return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, - const MeshEdgeBoundary& mesh_edge_boundary, - const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, + const MeshEdgeBoundary& mesh_edge_boundary, + const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_edge_boundary(mesh_edge_boundary), diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp index 2826f8de7..9e38262d5 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp @@ -38,14 +38,14 @@ class RusanovEulerianCompositeSolver_v2_o2 class SymmetryBoundaryCondition; class InflowListBoundaryCondition; class OutflowBoundaryCondition; - class NeumannBoundaryCondition; + class WallBoundaryCondition; class NeumannflatBoundaryCondition; using BoundaryCondition = std::variant<SymmetryBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition, NeumannflatBoundaryCondition, - NeumannBoundaryCondition>; + WallBoundaryCondition>; using BoundaryConditionList = std::vector<BoundaryCondition>; @@ -59,17 +59,15 @@ class RusanovEulerianCompositeSolver_v2_o2 bool is_valid_boundary_condition = true; switch (bc_descriptor->type()) { - case IBoundaryConditionDescriptor::Type::neumann: { + case IBoundaryConditionDescriptor::Type::wall: { if constexpr (Dimension == 2) { - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } else { static_assert(Dimension == 3); - bc_list.emplace_back( - NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), - getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); } break; } @@ -177,6 +175,137 @@ class RusanovEulerianCompositeSolver_v2_o2 CellByCellLimitation<MeshType> Limitor; public: + void + _applyWallBoundaryCondition(const BoundaryConditionList& bc_list, + const MeshType& mesh, + NodeValuePerCell<Rp>& stateNode, + EdgeValuePerCell<Rp>& stateEdge, + FaceValuePerCell<Rp>& stateFace) const + { + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<WallBoundaryCondition, T>) { + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + std::cout << " Traitement WALL local (non flat) \n"; + // const Rd& normal = bc.outgoingNormal(); + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); + // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& face_list = bc.faceList(); + const auto& node_list = bc.nodeList(); + + const auto Cjr = mesh_data.Cjr(); + const auto Cjf = mesh_data.Cjf(); + + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + const NodeId node_id = node_list[i_node]; + + const auto& node_cell_list = node_to_cell_matrix[node_id]; + // Assert(face_cell_list.size() == 1); + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + normal += Cjr(node_cell_id, node_local_number_in_cell); + } + normal *= 1. / node_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim]; + } + } + + 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); + + // Normal locale approchée + Rd normal(Cjf(face_cell_id, face_local_number_in_cell)); + normal *= 1. / l2Norm(normal); + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); + + const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); + + const auto& edge_list = bc.edgeList(); + + const auto Cje = mesh_data.Cje(); + + for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { + const EdgeId edge_id = edge_list[i_edge]; + + const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; + + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + + // Normal locale approchée + Rd normal(zero); + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + normal += Cje(edge_cell_id, edge_local_number_in_cell); + } + normal *= 1. / edge_cell_list.size(); + normal *= 1. / l2Norm(normal); + + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + } + } + } + }, + boundary_condition); + } + } void _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list, const MeshType& mesh, @@ -1333,12 +1462,12 @@ class RusanovEulerianCompositeSolver_v2_o2 }; template <MeshConcept MeshType> -class RusanovEulerianCompositeSolver_v2_o2<MeshType>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver_v2_o2<MeshType>::WallBoundaryCondition { }; template <> -class RusanovEulerianCompositeSolver_v2_o2<Mesh<2>>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver_v2_o2<Mesh<2>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1373,7 +1502,7 @@ class RusanovEulerianCompositeSolver_v2_o2<Mesh<2>>::NeumannBoundaryCondition return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary) { ; @@ -1381,7 +1510,7 @@ class RusanovEulerianCompositeSolver_v2_o2<Mesh<2>>::NeumannBoundaryCondition }; template <> -class RusanovEulerianCompositeSolver_v2_o2<Mesh<3>>::NeumannBoundaryCondition +class RusanovEulerianCompositeSolver_v2_o2<Mesh<3>>::WallBoundaryCondition { using Rd = TinyVector<Dimension, double>; @@ -1426,9 +1555,9 @@ class RusanovEulerianCompositeSolver_v2_o2<Mesh<3>>::NeumannBoundaryCondition return m_mesh_face_boundary.faceList(); } - NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, - const MeshEdgeBoundary& mesh_edge_boundary, - const MeshFaceBoundary& mesh_face_boundary) + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, + const MeshEdgeBoundary& mesh_edge_boundary, + const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_edge_boundary(mesh_edge_boundary), diff --git a/src/scheme/WallBoundaryConditionDescriptor.hpp b/src/scheme/WallBoundaryConditionDescriptor.hpp new file mode 100644 index 000000000..cf16fba7d --- /dev/null +++ b/src/scheme/WallBoundaryConditionDescriptor.hpp @@ -0,0 +1,40 @@ +#ifndef WALL_BOUNDARY_CONDITION_DESCRIPTOR_HPP +#define WALL_BOUNDARY_CONDITION_DESCRIPTOR_HPP + +#include <language/utils/FunctionSymbolId.hpp> +#include <mesh/IBoundaryDescriptor.hpp> +#include <scheme/BoundaryConditionDescriptorBase.hpp> + +#include <memory> + +class WallBoundaryConditionDescriptor : public BoundaryConditionDescriptorBase +{ + private: + std::ostream& + _write(std::ostream& os) const final + { + os << "wall(" << *m_boundary_descriptor << ")"; + return os; + } + + public: + Type + type() const final + { + return Type::wall; + } + + WallBoundaryConditionDescriptor(const std::shared_ptr<const IBoundaryDescriptor>& boundary_descriptor) + + : BoundaryConditionDescriptorBase{boundary_descriptor} + { + ; + } + + WallBoundaryConditionDescriptor(const WallBoundaryConditionDescriptor&) = delete; + WallBoundaryConditionDescriptor(WallBoundaryConditionDescriptor&&) = delete; + + ~WallBoundaryConditionDescriptor() = default; +}; + +#endif // WALL_BOUNDARY_CONDITION_DESCRIPTOR_HPP -- GitLab