From 66a727cc9e2bfedbe5c85c781a40c26bd9dfb011 Mon Sep 17 00:00:00 2001 From: HOCH PHILIPPE <philippe.hoch@gmail.com> Date: Mon, 31 Mar 2025 23:15:54 +0200 Subject: [PATCH] New implementation for Wall boundary at node in 2D .. still need coding in 3D --- ...eViscousFormEulerianCompositeSolver_v2.cpp | 65 +++++++++++------ ...scousFormEulerianCompositeSolver_v2_o2.cpp | 67 +++++++++++------- src/scheme/RusanovEulerianCompositeSolver.cpp | 65 +++++++++++------ .../RusanovEulerianCompositeSolver_o2.cpp | 69 ++++++++++++------- .../RusanovEulerianCompositeSolver_v2.cpp | 67 +++++++++++------- .../RusanovEulerianCompositeSolver_v2_o2.cpp | 68 +++++++++++------- 6 files changed, 262 insertions(+), 139 deletions(-) diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp index d56dc2a5..447bab9c 100644 --- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp +++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp @@ -519,48 +519,69 @@ class RoeViscousFormEulerianCompositeSolver_v2 // const Rd& normal = bc.outgoingNormal(); const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& node_to_face_matrix = mesh.connectivity().nodeToFaceMatrix(); const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& node_local_numbers_in_their_faces = mesh.connectivity().nodeLocalNumbersInTheirFaces(); const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); const auto& face_list = bc.faceList(); const auto& node_list = bc.nodeList(); - const auto Cjr = mesh_data.Cjr(); + // const auto Cjr = mesh_data.Cjr(); const auto Cjf = mesh_data.Cjf(); for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { const NodeId node_id = node_list[i_node]; - const auto& node_cell_list = node_to_cell_matrix[node_id]; + const auto& node_face_list = node_to_face_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); + // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id); - // Normal locale approchée - Rd normal(zero); - for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { - CellId node_cell_id = node_cell_list[i_cell]; - size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; - normal += Cjr(node_cell_id, node_local_number_in_cell); - } - normal *= 1. / node_cell_list.size(); - normal *= 1. / l2Norm(normal); + // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList - 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 i_face = 0; i_face < node_face_list.size(); ++i_face) { + FaceId node_face_id = node_face_list[i_face]; + Rd normal(zero); - Rd vectorSym(zero); - for (size_t dim = 0; dim < Dimension; ++dim) - vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (node_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_id]; + Assert(face_cell_list.size() == 1); - vectorSym -= dot(vectorSym, normal) * normal; + CellId face_cell_id = face_cell_list[0]; + size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0); - 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]; + // Normal locale approchée + normal = Cjf(face_cell_id, face_local_number_in_cell); + break; + } + } + if (l2Norm(normal) == 0.) + continue; + + normal *= 1. / l2Norm(normal); + 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]; + } } } diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp index e707e3f8..a2d09206 100644 --- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp +++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp @@ -527,48 +527,67 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 // const Rd& normal = bc.outgoingNormal(); const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& node_to_face_matrix = mesh.connectivity().nodeToFaceMatrix(); const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& node_local_numbers_in_their_faces = mesh.connectivity().nodeLocalNumbersInTheirFaces(); const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); const auto& face_list = bc.faceList(); const auto& node_list = bc.nodeList(); - const auto Cjr = mesh_data.Cjr(); + // const auto Cjr = mesh_data.Cjr(); const auto Cjf = mesh_data.Cjf(); for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - const NodeId node_id = node_list[i_node]; - - const auto& node_cell_list = node_to_cell_matrix[node_id]; + const NodeId node_id = node_list[i_node]; + const auto& node_face_list = node_to_face_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); + // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id); - // Normal locale approchée - Rd normal(zero); - for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { - CellId node_cell_id = node_cell_list[i_cell]; - size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; - normal += Cjr(node_cell_id, node_local_number_in_cell); - } - normal *= 1. / node_cell_list.size(); - normal *= 1. / l2Norm(normal); + // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList - 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 i_face = 0; i_face < node_face_list.size(); ++i_face) { + FaceId node_face_id = node_face_list[i_face]; + Rd normal(zero); - Rd vectorSym(zero); - for (size_t dim = 0; dim < Dimension; ++dim) - vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (node_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_id]; + Assert(face_cell_list.size() == 1); - vectorSym -= dot(vectorSym, normal) * normal; + CellId face_cell_id = face_cell_list[0]; + size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0); - 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]; + // Normal locale approchée + normal = Cjf(face_cell_id, face_local_number_in_cell); + break; + } + } + if (l2Norm(normal) == 0.) + continue; + + normal *= 1. / l2Norm(normal); + 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]; + } } } diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp index cda7e625..f8d7311e 100644 --- a/src/scheme/RusanovEulerianCompositeSolver.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver.cpp @@ -185,48 +185,69 @@ class RusanovEulerianCompositeSolver // const Rd& normal = bc.outgoingNormal(); const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& node_to_face_matrix = mesh.connectivity().nodeToFaceMatrix(); const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& node_local_numbers_in_their_faces = mesh.connectivity().nodeLocalNumbersInTheirFaces(); const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); const auto& face_list = bc.faceList(); const auto& node_list = bc.nodeList(); - const auto Cjr = mesh_data.Cjr(); + // const auto Cjr = mesh_data.Cjr(); const auto Cjf = mesh_data.Cjf(); for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { const NodeId node_id = node_list[i_node]; - const auto& node_cell_list = node_to_cell_matrix[node_id]; + const auto& node_face_list = node_to_face_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); + // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id); - // Normal locale approchée - Rd normal(zero); - for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { - CellId node_cell_id = node_cell_list[i_cell]; - size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; - normal += Cjr(node_cell_id, node_local_number_in_cell); - } - normal *= 1. / node_cell_list.size(); - normal *= 1. / l2Norm(normal); + // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList - 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 i_face = 0; i_face < node_face_list.size(); ++i_face) { + FaceId node_face_id = node_face_list[i_face]; + Rd normal(zero); - Rd vectorSym(zero); - for (size_t dim = 0; dim < Dimension; ++dim) - vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (node_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_id]; + Assert(face_cell_list.size() == 1); - vectorSym -= dot(vectorSym, normal) * normal; + CellId face_cell_id = face_cell_list[0]; + size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0); - 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]; + // Normal locale approchée + normal = Cjf(face_cell_id, face_local_number_in_cell); + break; + } + } + if (l2Norm(normal) == 0.) + continue; + + normal *= 1. / l2Norm(normal); + 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]; + } } } diff --git a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp index 7d53c2f8..4383a001 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp @@ -193,51 +193,71 @@ class RusanovEulerianCompositeSolver_o2 // const Rd& normal = bc.outgoingNormal(); const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& node_to_face_matrix = mesh.connectivity().nodeToFaceMatrix(); const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& node_local_numbers_in_their_faces = mesh.connectivity().nodeLocalNumbersInTheirFaces(); const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); const auto& face_list = bc.faceList(); const auto& node_list = bc.nodeList(); - const auto Cjr = mesh_data.Cjr(); + // const auto Cjr = mesh_data.Cjr(); const auto Cjf = mesh_data.Cjf(); for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { const NodeId node_id = node_list[i_node]; - const auto& node_cell_list = node_to_cell_matrix[node_id]; + const auto& node_face_list = node_to_face_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); + // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id); - // Normal locale approchée - Rd normal(zero); - for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { - CellId node_cell_id = node_cell_list[i_cell]; - size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; - normal += Cjr(node_cell_id, node_local_number_in_cell); - } - normal *= 1. / node_cell_list.size(); - normal *= 1. / l2Norm(normal); + // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList - 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 i_face = 0; i_face < node_face_list.size(); ++i_face) { + FaceId node_face_id = node_face_list[i_face]; + Rd normal(zero); - Rd vectorSym(zero); - for (size_t dim = 0; dim < Dimension; ++dim) - vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (node_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_id]; + Assert(face_cell_list.size() == 1); - vectorSym -= dot(vectorSym, normal) * normal; + CellId face_cell_id = face_cell_list[0]; + size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0); - 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]; + // Normal locale approchée + normal = Cjf(face_cell_id, face_local_number_in_cell); + break; + } + } + if (l2Norm(normal) == 0.) + continue; + + normal *= 1. / l2Norm(normal); + 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]; @@ -795,7 +815,8 @@ class RusanovEulerianCompositeSolver_o2 DiscreteFunctionDPk eps_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>(); // DiscreteFunctionDPk rho_u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>(); - // DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>(); + // DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const + // double>>(); CellValue<double> Limitor_rho(p_mesh->connectivity()); CellValue<double> Limitor_u(p_mesh->connectivity()); diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp index f367cb16..0fa8469d 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp @@ -185,49 +185,69 @@ class RusanovEulerianCompositeSolver_v2 // const Rd& normal = bc.outgoingNormal(); const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& node_to_face_matrix = mesh.connectivity().nodeToFaceMatrix(); const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& node_local_numbers_in_their_faces = mesh.connectivity().nodeLocalNumbersInTheirFaces(); const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); const auto& face_list = bc.faceList(); const auto& node_list = bc.nodeList(); - const auto Cjr = mesh_data.Cjr(); + // const auto Cjr = mesh_data.Cjr(); const auto Cjf = mesh_data.Cjf(); for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { const NodeId node_id = node_list[i_node]; - const auto& node_cell_list = node_to_cell_matrix[node_id]; + const auto& node_face_list = node_to_face_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); + // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id); - // Normal locale approchée - Rd normal(zero); - for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { - CellId node_cell_id = node_cell_list[i_cell]; - size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; - normal += Cjr(node_cell_id, node_local_number_in_cell); - } + // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList - normal *= 1. / node_cell_list.size(); - normal *= 1. / l2Norm(normal); + for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) { + FaceId node_face_id = node_face_list[i_face]; + Rd normal(zero); - for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { - CellId node_cell_id = node_cell_list[i_cell]; - size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (node_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_id]; + Assert(face_cell_list.size() == 1); - Rd vectorSym(zero); - for (size_t dim = 0; dim < Dimension; ++dim) - vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + CellId face_cell_id = face_cell_list[0]; + size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0); - vectorSym -= dot(vectorSym, normal) * normal; + // Normal locale approchée + normal = Cjf(face_cell_id, face_local_number_in_cell); + break; + } + } + if (l2Norm(normal) == 0.) + continue; - 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 + 1] = 0; + normal *= 1. / l2Norm(normal); + 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 + 1] = 0; + } } } @@ -271,7 +291,7 @@ class RusanovEulerianCompositeSolver_v2 const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); - // Normal locale approchée + // Normal locale approchée .. need better def Rd normal(zero); for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { CellId edge_cell_id = edge_cell_list[i_cell]; @@ -302,6 +322,7 @@ class RusanovEulerianCompositeSolver_v2 boundary_condition); } } + void _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list, const MeshType& mesh, diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp index 7a22ea9c..42cb12e5 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp @@ -193,48 +193,68 @@ class RusanovEulerianCompositeSolver_v2_o2 // const Rd& normal = bc.outgoingNormal(); const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& node_to_face_matrix = mesh.connectivity().nodeToFaceMatrix(); const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& node_local_numbers_in_their_faces = mesh.connectivity().nodeLocalNumbersInTheirFaces(); const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); const auto& face_list = bc.faceList(); const auto& node_list = bc.nodeList(); - const auto Cjr = mesh_data.Cjr(); + // const auto Cjr = mesh_data.Cjr(); const auto Cjf = mesh_data.Cjf(); for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - const NodeId node_id = node_list[i_node]; - - const auto& node_cell_list = node_to_cell_matrix[node_id]; + const NodeId node_id = node_list[i_node]; + const auto& node_face_list = node_to_face_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); + // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id); - // Normal locale approchée - Rd normal(zero); - for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { - CellId node_cell_id = node_cell_list[i_cell]; - size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; - normal += Cjr(node_cell_id, node_local_number_in_cell); - } - normal *= 1. / node_cell_list.size(); - normal *= 1. / l2Norm(normal); + // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList - 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 i_face = 0; i_face < node_face_list.size(); ++i_face) { + FaceId node_face_id = node_face_list[i_face]; + Rd normal(zero); - Rd vectorSym(zero); - for (size_t dim = 0; dim < Dimension; ++dim) - vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (node_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_id]; + Assert(face_cell_list.size() == 1); - vectorSym -= dot(vectorSym, normal) * normal; + CellId face_cell_id = face_cell_list[0]; + size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0); - 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]; + // Normal locale approchée + normal = Cjf(face_cell_id, face_local_number_in_cell); + break; + } + } + if (l2Norm(normal) == 0.) + continue; + + normal *= 1. / l2Norm(normal); + 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]; + } } } -- GitLab