diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
index 64e7837fa6d3a7bcf8a77d2853aa5a01e2e111ef..2d4c3b74a82f39f3e1db0a29d8a68d3162f34d1c 100644
--- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
@@ -172,35 +172,426 @@ class RoeViscousFormEulerianCompositeSolver_v2
}
public:
- void _applyNeumannBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
-
- void _applyNeumannflatBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
-
- void _applyInflowBoundaryCondition(const BoundaryConditionList& bc_list,
+ void
+ _applyOutflowBoundaryCondition(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<OutflowBoundaryCondition, T>) {
+ std::cout << " Traitement Outflow \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 xj = mesh.xj();
+ const auto xr = mesh.xr();
+ const auto xf = mesh.xl();
+ const auto xe = mesh.xe();
+
+ 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);
+
+ 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];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateNode[node_cell_id][node_local_number_in_cell][dim] += vectorSym[dim];
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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);
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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();
+
+ 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];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ 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];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+ }
+
+ // throw NormalError("Not implemented");
+ }
+ */
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _applySymmetricBoundaryCondition(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<SymmetryBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement SYMMETRY \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();
+
+ 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);
+
+ 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];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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);
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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();
+
+ 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];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ 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];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _applyNeumannflatBoundaryCondition(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
+ 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<NeumannflatBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement WALL \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();
+
+ 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);
+
+ 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);
+
+ 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();
+
+ 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];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ 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,
+ 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<InflowListBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement INFLOW \n";
- void _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
- void _applySymmetricBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
+ 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& face_array_list = bc.faceArrayList();
+ const auto& node_array_list = bc.nodeArrayList();
+
+ 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);
+
+ 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];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateNode[node_cell_id][node_local_number_in_cell][dim] = 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);
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][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& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+ const auto& edge_list = bc.edgeList();
+
+ const auto& edge_array_list = bc.edgeArrayList();
+
+ 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);
+
+ // Assert(face_cell_list.size() == 1);
+
+ 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];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
public:
double
@@ -280,6 +671,24 @@ class RoeViscousFormEulerianCompositeSolver_v2
double maxabsVpNormal;
// bool changingSign;
double MaxErrorProd;
+
+ public:
+ JacobianInformations(const Rpxp& Jacobiand,
+ const Rpxp& LeftEigenVectorsd,
+ const Rpxp& RightEigenVectorsd,
+ const Rp& EigenValuesd,
+ const Rpxp& AbsJacobiand,
+ const double maxabsVpNormald,
+ const double MaxErrorProdd)
+ : Jacobian(Jacobiand),
+ LeftEigenVectors(LeftEigenVectorsd),
+ RightEigenVectors(RightEigenVectorsd),
+ EigenValues(EigenValuesd),
+ AbsJacobian(AbsJacobiand),
+ maxabsVpNormal(maxabsVpNormald),
+ MaxErrorProd(MaxErrorProdd)
+ {}
+ ~JacobianInformations() {}
};
Rp
@@ -478,8 +887,8 @@ class RoeViscousFormEulerianCompositeSolver_v2
const double c = normal[2];
if ((a == b) and (b == c)) {
- static constexpr double invsqrt2 = 1. / sqrt(2.);
- static constexpr double invsqrt6 = 1. / sqrt(6.);
+ static double invsqrt2 = 1. / sqrt(2.);
+ static double invsqrt6 = 1. / sqrt(6.);
ortho[0] = Rd{invsqrt2, -invsqrt2, 0};
ortho[1] = Rd{invsqrt6, invsqrt6, -2 * invsqrt6}; // au signe pres
@@ -564,8 +973,8 @@ class RoeViscousFormEulerianCompositeSolver_v2
Rpxp ProdLeft = Right * EigenAbs;
Rpxp JacobianR = ProdLeft * Left;
- Rpxp id = identity;
- maxErr = 0;
+ // Rpxp id = identity;
+ maxErr = 0;
// double maxErrLeftRightId = 0;
for (size_t i = 0; i < Dimension + 2; ++i)
for (size_t j = 0; j < Dimension + 2; ++j) {
@@ -1949,431 +2358,3 @@ roeViscousFormEulerianCompositeSolver_v2(
},
mesh_v->variant());
}
-
-template <MeshConcept MeshType>
-void
-RoeViscousFormEulerianCompositeSolver_v2<MeshType>::_applyOutflowBoundaryCondition(
- 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<OutflowBoundaryCondition, T>) {
- std::cout << " Traitement Outflow \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 xj = mesh.xj();
- const auto xr = mesh.xr();
- const auto xf = mesh.xl();
- const auto xe = mesh.xe();
-
- 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);
-
- 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];
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateNode[node_cell_id][node_local_number_in_cell][dim] += vectorSym[dim];
-
- Rd vectorSym(zero);
- for (size_t dim = 0; dim < Dimension; ++dim)
- vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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);
-
- Rd vectorSym(zero);
- for (size_t dim = 0; dim < Dimension; ++dim)
- vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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();
-
- 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];
- // Assert(face_cell_list.size() == 1);
- const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
-
- 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];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- for (size_t dim = 0; dim < Dimension; ++dim)
- stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
- }
- }
-
- // throw NormalError("Not implemented");
- }
- */
- }
- },
- boundary_condition);
- }
-}
-
-template <MeshConcept MeshType>
-void
-RoeViscousFormEulerianCompositeSolver_v2<MeshType>::_applySymmetricBoundaryCondition(
- 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<SymmetryBoundaryCondition, T>) {
- // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- std::cout << " Traitement SYMMETRY \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();
-
- 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);
-
- 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];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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);
-
- Rd vectorSym(zero);
- for (size_t dim = 0; dim < Dimension; ++dim)
- vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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();
-
- 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];
- // Assert(face_cell_list.size() == 1);
- const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
-
- 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];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- for (size_t dim = 0; dim < Dimension; ++dim)
- stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
- }
- }
- }
- }
- },
- boundary_condition);
- }
-}
-
-template <MeshConcept MeshType>
-void
-RoeViscousFormEulerianCompositeSolver_v2<MeshType>::_applyNeumannflatBoundaryCondition(
- 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<NeumannflatBoundaryCondition, T>) {
- // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- std::cout << " Traitement WALL \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();
-
- 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);
-
- 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);
-
- 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();
-
- 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];
- // Assert(face_cell_list.size() == 1);
- const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
-
- 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);
- }
-}
-
-template <MeshConcept MeshType>
-void
-RoeViscousFormEulerianCompositeSolver_v2<MeshType>::_applyInflowBoundaryCondition(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<InflowListBoundaryCondition, T>) {
- // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- std::cout << " Traitement INFLOW \n";
-
- 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& face_array_list = bc.faceArrayList();
- const auto& node_array_list = bc.nodeArrayList();
-
- 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);
-
- 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];
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateNode[node_cell_id][node_local_number_in_cell][dim] = 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);
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][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& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
-
- const auto& edge_list = bc.edgeList();
-
- const auto& edge_array_list = bc.edgeArrayList();
-
- 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);
-
- // Assert(face_cell_list.size() == 1);
-
- 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];
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim];
- }
- }
- }
- }
- },
- boundary_condition);
- }
-}
diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp
index 8f189d32ff348063c69cac24ff89f56ea444bb25..4192bacc4a677921c7151e801f17d39d02b65032 100644
--- a/src/scheme/RusanovEulerianCompositeSolver.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver.cpp
@@ -170,35 +170,426 @@ class RusanovEulerianCompositeSolver
}
public:
- void _applyNeumannBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
-
- void _applyNeumannflatBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
-
- void _applyInflowBoundaryCondition(const BoundaryConditionList& bc_list,
+ void
+ _applyOutflowBoundaryCondition(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<OutflowBoundaryCondition, T>) {
+ std::cout << " Traitement Outflow \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 xj = mesh.xj();
+ const auto xr = mesh.xr();
+ const auto xf = mesh.xl();
+ const auto xe = mesh.xe();
+
+ 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);
+
+ 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];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateNode[node_cell_id][node_local_number_in_cell][dim] += vectorSym[dim];
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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);
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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();
+
+ 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];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ 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];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+ }
+
+ // throw NormalError("Not implemented");
+ }
+ */// throw NormalError("Not implemented");
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _applySymmetricBoundaryCondition(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<SymmetryBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement SYMMETRY \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();
+
+ 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);
+
+ 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];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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);
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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();
+
+ 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];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ 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];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _applyNeumannflatBoundaryCondition(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
+ 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<NeumannflatBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement WALL \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();
+
+ 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);
+
+ 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);
+
+ 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();
+
+ 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];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ 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,
+ 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<InflowListBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement INFLOW \n";
- void _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
- void _applySymmetricBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
+ 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& face_array_list = bc.faceArrayList();
+ const auto& node_array_list = bc.nodeArrayList();
+
+ 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);
+
+ 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];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateNode[node_cell_id][node_local_number_in_cell][dim] = 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);
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][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& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+ const auto& edge_list = bc.edgeList();
+
+ const auto& edge_array_list = bc.edgeArrayList();
+
+ 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);
+
+ // Assert(face_cell_list.size() == 1);
+
+ 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];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
public:
double
@@ -1539,428 +1930,3 @@ rusanovEulerianCompositeSolver(
},
mesh_v->variant());
}
-
-template <MeshConcept MeshType>
-void
-RusanovEulerianCompositeSolver<MeshType>::_applyOutflowBoundaryCondition(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<OutflowBoundaryCondition, T>) {
- std::cout << " Traitement Outflow \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 xj = mesh.xj();
- const auto xr = mesh.xr();
- const auto xf = mesh.xl();
- const auto xe = mesh.xe();
-
- 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);
-
- 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];
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateNode[node_cell_id][node_local_number_in_cell][dim] += vectorSym[dim];
-
- Rd vectorSym(zero);
- for (size_t dim = 0; dim < Dimension; ++dim)
- vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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);
-
- Rd vectorSym(zero);
- for (size_t dim = 0; dim < Dimension; ++dim)
- vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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();
-
- 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];
- // Assert(face_cell_list.size() == 1);
- const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
-
- 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];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- for (size_t dim = 0; dim < Dimension; ++dim)
- stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
- }
- }
-
- // throw NormalError("Not implemented");
- }
- */// throw NormalError("Not implemented");
- }
- },
- boundary_condition);
- }
-}
-
-template <MeshConcept MeshType>
-void
-RusanovEulerianCompositeSolver<MeshType>::_applySymmetricBoundaryCondition(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<SymmetryBoundaryCondition, T>) {
- // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- std::cout << " Traitement SYMMETRY \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();
-
- 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);
-
- 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];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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);
-
- Rd vectorSym(zero);
- for (size_t dim = 0; dim < Dimension; ++dim)
- vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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();
-
- 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];
- // Assert(face_cell_list.size() == 1);
- const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
-
- 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];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- for (size_t dim = 0; dim < Dimension; ++dim)
- stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
- }
- }
- }
- }
- },
- boundary_condition);
- }
-}
-
-template <MeshConcept MeshType>
-void
-RusanovEulerianCompositeSolver<MeshType>::_applyNeumannflatBoundaryCondition(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<NeumannflatBoundaryCondition, T>) {
- // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- std::cout << " Traitement WALL \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();
-
- 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);
-
- 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);
-
- 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();
-
- 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];
- // Assert(face_cell_list.size() == 1);
- const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
-
- 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);
- }
-}
-
-template <MeshConcept MeshType>
-void
-RusanovEulerianCompositeSolver<MeshType>::_applyInflowBoundaryCondition(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<InflowListBoundaryCondition, T>) {
- // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- std::cout << " Traitement INFLOW \n";
-
- 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& face_array_list = bc.faceArrayList();
- const auto& node_array_list = bc.nodeArrayList();
-
- 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);
-
- 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];
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateNode[node_cell_id][node_local_number_in_cell][dim] = 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);
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][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& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
-
- const auto& edge_list = bc.edgeList();
-
- const auto& edge_array_list = bc.edgeArrayList();
-
- 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);
-
- // Assert(face_cell_list.size() == 1);
-
- 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];
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim];
- }
- }
- }
- }
- },
- boundary_condition);
- }
-}
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
index cd3d407cb8840d594a8341ca196759a40e57b9bb..1a003ade61c35f99fc8bbb6b251b478999f14dc1 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
@@ -170,35 +170,426 @@ class RusanovEulerianCompositeSolver_v2
}
public:
- void _applyNeumannBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
-
- void _applyNeumannflatBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
-
- void _applyInflowBoundaryCondition(const BoundaryConditionList& bc_list,
+ void
+ _applyOutflowBoundaryCondition(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<OutflowBoundaryCondition, T>) {
+ std::cout << " Traitement Outflow \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 xj = mesh.xj();
+ const auto xr = mesh.xr();
+ const auto xf = mesh.xl();
+ const auto xe = mesh.xe();
+
+ 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);
+
+ 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];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateNode[node_cell_id][node_local_number_in_cell][dim] += vectorSym[dim];
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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);
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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();
+
+ 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];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ 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];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+ }
+
+ // throw NormalError("Not implemented");
+ }
+ */
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _applySymmetricBoundaryCondition(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<SymmetryBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement SYMMETRY \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();
+
+ 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);
+
+ 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];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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);
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ 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();
+
+ 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];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ 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];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _applyNeumannflatBoundaryCondition(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
+ 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<NeumannflatBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement WALL \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();
+
+ 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);
+
+ 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];
- void _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
+ const auto& face_cell_list = face_to_cell_matrix[face_id];
+ Assert(face_cell_list.size() == 1);
- void _applySymmetricBoundaryCondition(const BoundaryConditionList& bc_list,
- const MeshType& mesh,
- NodeValuePerCell<Rp>& stateNode,
- EdgeValuePerCell<Rp>& stateEdge,
- FaceValuePerCell<Rp>& stateFace) const;
+ CellId face_cell_id = face_cell_list[0];
+ size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+ 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();
+
+ 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];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ 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,
+ 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<InflowListBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement INFLOW \n";
+
+ 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& face_array_list = bc.faceArrayList();
+ const auto& node_array_list = bc.nodeArrayList();
+
+ 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);
+
+ 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];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateNode[node_cell_id][node_local_number_in_cell][dim] = 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);
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][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& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+ const auto& edge_list = bc.edgeList();
+
+ const auto& edge_array_list = bc.edgeArrayList();
+
+ 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);
+
+ // Assert(face_cell_list.size() == 1);
+
+ 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];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
public:
double
@@ -1442,428 +1833,3 @@ rusanovEulerianCompositeSolver_v2(
},
mesh_v->variant());
}
-
-template <MeshConcept MeshType>
-void
-RusanovEulerianCompositeSolver_v2<MeshType>::_applyOutflowBoundaryCondition(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<OutflowBoundaryCondition, T>) {
- std::cout << " Traitement Outflow \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 xj = mesh.xj();
- const auto xr = mesh.xr();
- const auto xf = mesh.xl();
- const auto xe = mesh.xe();
-
- 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);
-
- 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];
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateNode[node_cell_id][node_local_number_in_cell][dim] += vectorSym[dim];
-
- Rd vectorSym(zero);
- for (size_t dim = 0; dim < Dimension; ++dim)
- vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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);
-
- Rd vectorSym(zero);
- for (size_t dim = 0; dim < Dimension; ++dim)
- vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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();
-
- 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];
- // Assert(face_cell_list.size() == 1);
- const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
-
- 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];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- for (size_t dim = 0; dim < Dimension; ++dim)
- stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
- }
- }
-
- // throw NormalError("Not implemented");
- }
- */
- }
- },
- boundary_condition);
- }
-}
-
-template <MeshConcept MeshType>
-void
-RusanovEulerianCompositeSolver_v2<MeshType>::_applySymmetricBoundaryCondition(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<SymmetryBoundaryCondition, T>) {
- // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- std::cout << " Traitement SYMMETRY \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();
-
- 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);
-
- 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];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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);
-
- Rd vectorSym(zero);
- for (size_t dim = 0; dim < Dimension; ++dim)
- vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- 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();
-
- 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];
- // Assert(face_cell_list.size() == 1);
- const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
-
- 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];
-
- Rdxd MatriceProj(identity);
- MatriceProj -= tensorProduct(normal, normal);
- vectorSym = MatriceProj * vectorSym;
-
- for (size_t dim = 0; dim < Dimension; ++dim)
- stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
- }
- }
- }
- }
- },
- boundary_condition);
- }
-}
-
-template <MeshConcept MeshType>
-void
-RusanovEulerianCompositeSolver_v2<MeshType>::_applyNeumannflatBoundaryCondition(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<NeumannflatBoundaryCondition, T>) {
- // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- std::cout << " Traitement WALL \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();
-
- 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);
-
- 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);
-
- 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();
-
- 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];
- // Assert(face_cell_list.size() == 1);
- const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
-
- 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);
- }
-}
-
-template <MeshConcept MeshType>
-void
-RusanovEulerianCompositeSolver_v2<MeshType>::_applyInflowBoundaryCondition(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<InflowListBoundaryCondition, T>) {
- // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- std::cout << " Traitement INFLOW \n";
-
- 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& face_array_list = bc.faceArrayList();
- const auto& node_array_list = bc.nodeArrayList();
-
- 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);
-
- 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];
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateNode[node_cell_id][node_local_number_in_cell][dim] = 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);
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][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& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
-
- const auto& edge_list = bc.edgeList();
-
- const auto& edge_array_list = bc.edgeArrayList();
-
- 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);
-
- // Assert(face_cell_list.size() == 1);
-
- 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];
-
- for (size_t dim = 0; dim < Dimension + 2; ++dim)
- stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim];
- }
- }
- }
- }
- },
- boundary_condition);
- }
-}