diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp index 397eb85e7d67a020138516acc5f602652bc063a1..fd856ffbc8fd9ce547e4666fb0eb471dfe413686 100644 --- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp +++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp @@ -620,26 +620,49 @@ class RoeViscousFormEulerianCompositeSolver_v2 const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); - const auto& edge_list = bc.edgeList(); + const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix(); + + const auto& edge_local_numbers_in_their_faces = mesh.connectivity().edgeLocalNumbersInTheirFaces(); - const auto Cje = mesh_data.Cje(); + 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]; - - const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + const auto& edge_face_list = edge_to_face_matrix[edge_id]; - // Normal locale approchée + // on va chercher les normale d'une face issue du edge de CL et contenue dans le faceList Rd normal(zero); - for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { - CellId edge_cell_id = edge_cell_list[i_cell]; - size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; - normal += Cje(edge_cell_id, edge_local_number_in_cell); + int nbnormal = 0; + for (size_t i_face = 0; i_face < edge_face_list.size(); ++i_face) { + FaceId edge_face_id = edge_face_list[i_face]; + + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (edge_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_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(facebc_id, 0); + + // Normal locale approchée + Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell); + normalloc *= 1. / l2Norm(normalloc); + normal += normalloc; + ++nbnormal; + break; + } + } } - normal *= 1. / edge_cell_list.size(); - normal *= 1. / l2Norm(normal); + + if (nbnormal == 0) + continue; + normal *= 1. / nbnormal; + + 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); for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { CellId edge_cell_id = edge_cell_list[i_cell]; diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp index 2c9e1711c951a3fb18f43aa2a1e73ec9094666aa..934c39b9674a6b258bf1a564eccee21146b30a1a 100644 --- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp +++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp @@ -625,26 +625,48 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); - const auto& edge_list = bc.edgeList(); + const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix(); + + const auto& edge_local_numbers_in_their_faces = mesh.connectivity().edgeLocalNumbersInTheirFaces(); - const auto Cje = mesh_data.Cje(); + 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]; + const auto& edge_face_list = edge_to_face_matrix[edge_id]; - const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); - - // Normal locale approchée + // on va chercher les normale d'une face issue du edge de CL et contenue dans le faceList Rd normal(zero); - for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { - CellId edge_cell_id = edge_cell_list[i_cell]; - size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; - normal += Cje(edge_cell_id, edge_local_number_in_cell); + int nbnormal = 0; + for (size_t i_face = 0; i_face < edge_face_list.size(); ++i_face) { + FaceId edge_face_id = edge_face_list[i_face]; + + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (edge_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_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(facebc_id, 0); + + // Normal locale approchée + Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell); + normalloc *= 1. / l2Norm(normalloc); + normal += normalloc; + ++nbnormal; + break; + } + } } - normal *= 1. / edge_cell_list.size(); - normal *= 1. / l2Norm(normal); + + if (nbnormal == 0) + continue; + normal *= 1. / nbnormal; + 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); for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { CellId edge_cell_id = edge_cell_list[i_cell]; diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp index 2bedc88795d33628eb3c410e22543e1408c8230a..8374ae71679e9383c97ddb29546c954b5b352556 100644 --- a/src/scheme/RusanovEulerianCompositeSolver.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver.cpp @@ -285,27 +285,51 @@ class RusanovEulerianCompositeSolver const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); - const auto& edge_list = bc.edgeList(); - - const auto Cje = mesh_data.Cje(); + const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix(); - for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { - const EdgeId edge_id = edge_list[i_edge]; + const auto& edge_local_numbers_in_their_faces = mesh.connectivity().edgeLocalNumbersInTheirFaces(); - const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; + const auto& edge_list = bc.edgeList(); - const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { + const EdgeId edge_id = edge_list[i_edge]; + const auto& edge_face_list = edge_to_face_matrix[edge_id]; - // Normal locale approchée + // on va chercher les normale d'une face issue du edge de CL et contenue dans le faceList Rd normal(zero); - for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { - CellId edge_cell_id = edge_cell_list[i_cell]; - size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; - normal += Cje(edge_cell_id, edge_local_number_in_cell); + int nbnormal = 0; + for (size_t i_face = 0; i_face < edge_face_list.size(); ++i_face) { + FaceId edge_face_id = edge_face_list[i_face]; + + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (edge_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_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(facebc_id, 0); + + // Normal locale approchée + Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell); + normalloc *= 1. / l2Norm(normalloc); + normal += normalloc; + ++nbnormal; + break; + } + } } - normal *= 1. / edge_cell_list.size(); + + if (nbnormal == 0) + continue; + normal *= 1. / nbnormal; + normal *= 1. / l2Norm(normal); + 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); + 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]; @@ -381,7 +405,8 @@ class RusanovEulerianCompositeSolver 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]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // + node_array_list[i_node][dim]; } } @@ -493,7 +518,8 @@ class RusanovEulerianCompositeSolver 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]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // + // node_array_list[i_node][dim]; } } @@ -600,7 +626,8 @@ class RusanovEulerianCompositeSolver 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]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // + // node_array_list[i_node][dim]; } } @@ -1242,7 +1269,8 @@ class RusanovEulerianCompositeSolver } /* std::cout << " maxErr face " << xl[l] << " err " << l2Norm(SumGjf) << " size face cell" - << face_to_cell.size() << " Cell " << face_to_cell[0] << " / " << face_to_cell[1] << " Cjf0 " + << face_to_cell.size() << " Cell " << face_to_cell[0] << " / " << face_to_cell[1] << " Cjf0 + " << Cjf(face_to_cell[0], 0) << " / Cjf1 " << Cjf(face_to_cell[1], 1) << "\n"; if (l2Norm(SumGjf) > 1e-3) diff --git a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp index 6fbaafb7a2813fc4f0a0e153d4001b9734c3f397..9e59e875db07c820e988da528302ff680cb9c1be 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp @@ -293,26 +293,49 @@ class RusanovEulerianCompositeSolver_o2 const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); - const auto& edge_list = bc.edgeList(); + const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix(); + + const auto& edge_local_numbers_in_their_faces = mesh.connectivity().edgeLocalNumbersInTheirFaces(); - const auto Cje = mesh_data.Cje(); + 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]; - - const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + const auto& edge_face_list = edge_to_face_matrix[edge_id]; - // Normal locale approchée + // on va chercher les normale d'une face issue du edge de CL et contenue dans le faceList Rd normal(zero); - for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { - CellId edge_cell_id = edge_cell_list[i_cell]; - size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; - normal += Cje(edge_cell_id, edge_local_number_in_cell); + int nbnormal = 0; + for (size_t i_face = 0; i_face < edge_face_list.size(); ++i_face) { + FaceId edge_face_id = edge_face_list[i_face]; + + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (edge_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_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(facebc_id, 0); + + // Normal locale approchée + Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell); + normalloc *= 1. / l2Norm(normalloc); + normal += normalloc; + ++nbnormal; + break; + } + } } - normal *= 1. / edge_cell_list.size(); - normal *= 1. / l2Norm(normal); + + if (nbnormal == 0) + continue; + normal *= 1. / nbnormal; + + 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); for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { CellId edge_cell_id = edge_cell_list[i_cell]; diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp index af7c70e379caec3ff033ab3d91edb6a4b384b3e1..52269034acfe3179ba76851ad5a627b74d013874 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp @@ -284,27 +284,49 @@ class RusanovEulerianCompositeSolver_v2 const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); - const auto& edge_list = bc.edgeList(); + const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix(); + + const auto& edge_local_numbers_in_their_faces = mesh.connectivity().edgeLocalNumbersInTheirFaces(); - const auto Cje = mesh_data.Cje(); + 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 EdgeId edge_id = edge_list[i_edge]; + const auto& edge_face_list = edge_to_face_matrix[edge_id]; + + // on va chercher les normale d'une face issue du edge de CL et contenue dans le faceList + Rd normal(zero); + int nbnormal = 0; + for (size_t i_face = 0; i_face < edge_face_list.size(); ++i_face) { + FaceId edge_face_id = edge_face_list[i_face]; + + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (edge_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_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(facebc_id, 0); + + // Normal locale approchée + Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell); + normalloc *= 1. / l2Norm(normalloc); + normal += normalloc; + ++nbnormal; + break; + } + } + } + + if (nbnormal == 0) + continue; + normal *= 1. / nbnormal; const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); - // Normal locale approchée .. 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]; - size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; - normal += Cje(edge_cell_id, edge_local_number_in_cell); - } - normal *= 1. / edge_cell_list.size(); - normal *= 1. / l2Norm(normal); - for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { CellId edge_cell_id = edge_cell_list[i_cell]; size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp index be4c0495efd775e73e684b2effe2c23d1437a130..399c0c230d7c5884353ac085675392468481286e 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp @@ -292,27 +292,49 @@ class RusanovEulerianCompositeSolver_v2_o2 const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); - const auto& edge_list = bc.edgeList(); + const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix(); + + const auto& edge_local_numbers_in_their_faces = mesh.connectivity().edgeLocalNumbersInTheirFaces(); - const auto Cje = mesh_data.Cje(); + 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 EdgeId edge_id = edge_list[i_edge]; + const auto& edge_face_list = edge_to_face_matrix[edge_id]; + + // on va chercher les normale d'une face issue du edge de CL et contenue dans le faceList + Rd normal(zero); + int nbnormal = 0; + for (size_t i_face = 0; i_face < edge_face_list.size(); ++i_face) { + FaceId edge_face_id = edge_face_list[i_face]; + + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (edge_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_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(facebc_id, 0); + + // Normal locale approchée + Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell); + normalloc *= 1. / l2Norm(normalloc); + normal += normalloc; + ++nbnormal; + break; + } + } + } + + if (nbnormal == 0) + continue; + normal *= 1. / nbnormal; const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); - // Normal locale approchée - Rd normal(zero); - for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { - CellId edge_cell_id = edge_cell_list[i_cell]; - size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; - normal += Cje(edge_cell_id, edge_local_number_in_cell); - } - normal *= 1. / edge_cell_list.size(); - normal *= 1. / l2Norm(normal); - for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { CellId edge_cell_id = edge_cell_list[i_cell]; size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];