diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
index ce30fc3520ff3d9239a8cbbec3f539fcdfaf6466..8752bbecf4d17c54183babdde52addb24bc7e362 100644
--- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
@@ -352,7 +352,7 @@ class RoeViscousFormEulerianCompositeSolver_v2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
@@ -459,7 +459,7 @@ class RoeViscousFormEulerianCompositeSolver_v2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
@@ -519,34 +519,62 @@ 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
+ // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
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);
+ int nbnormal = 0;
+
+ for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
+ FaceId node_face_id = node_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 (node_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. / node_cell_list.size();
+
+ if (nbnormal == 0)
+ continue;
+ normal *= 1. / nbnormal;
+
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];
@@ -579,7 +607,7 @@ class RoeViscousFormEulerianCompositeSolver_v2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
@@ -592,27 +620,52 @@ 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 Cje = mesh_data.Cje();
+ const auto& edge_local_numbers_in_their_faces = mesh.connectivity().edgeLocalNumbersInTheirFaces();
+
+ 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();
+
+ 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];
diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
index d22e940a9c3d53575747cf3d09bb1160377f209a..b164663491b3753fb4a1b49dac507f01c554ce27 100644
--- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
@@ -360,7 +360,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
@@ -467,7 +467,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
@@ -527,35 +527,60 @@ 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
+ // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
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);
+ int nbnormal = 0;
+
+ for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
+ FaceId node_face_id = node_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 (node_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. / node_cell_list.size();
- normal *= 1. / l2Norm(normal);
+ if (nbnormal == 0)
+ continue;
+ normal *= 1. / nbnormal;
+ 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];
@@ -587,7 +612,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
@@ -600,27 +625,52 @@ 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();
+
+ 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];
@@ -944,7 +994,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2
const double gammaK,
const double cK,
const double pK,
-*/
+ */
const Rd& normal,
const bool check = false) const
{
diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp
index a13b34fb94582505b86683af9eed5efdca6f8af3..8374ae71679e9383c97ddb29546c954b5b352556 100644
--- a/src/scheme/RusanovEulerianCompositeSolver.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver.cpp
@@ -185,34 +185,61 @@ 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
+ // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
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);
+ int nbnormal = 0;
+
+ for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
+ FaceId node_face_id = node_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 (node_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. / node_cell_list.size();
+
+ if (nbnormal == 0)
+ continue;
+ normal *= 1. / nbnormal;
+
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];
@@ -245,7 +272,7 @@ class RusanovEulerianCompositeSolver
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
@@ -258,27 +285,51 @@ class RusanovEulerianCompositeSolver
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 Cje = mesh_data.Cje();
+ const auto& edge_local_numbers_in_their_faces = mesh.connectivity().edgeLocalNumbersInTheirFaces();
- 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_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];
@@ -354,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];
}
}
@@ -466,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];
}
}
@@ -481,7 +534,7 @@ class RusanovEulerianCompositeSolver
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
@@ -573,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];
}
}
@@ -588,7 +642,7 @@ class RusanovEulerianCompositeSolver
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
@@ -1215,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 867e6c4e5f7c2f9865820e954c26c219800f2cb9..e17e9dc7145e8db3808b972558d3376c572f8f8e 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp
@@ -193,34 +193,61 @@ 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
+ // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
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);
+ int nbnormal = 0;
+
+ for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
+ FaceId node_face_id = node_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 (node_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. / node_cell_list.size();
+
+ if (nbnormal == 0)
+ continue;
+ normal *= 1. / nbnormal;
+
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];
@@ -253,7 +280,7 @@ class RusanovEulerianCompositeSolver_o2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
@@ -266,27 +293,52 @@ 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();
+
+ 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];
@@ -489,7 +541,7 @@ class RusanovEulerianCompositeSolver_o2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
@@ -596,7 +648,7 @@ class RusanovEulerianCompositeSolver_o2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
@@ -795,7 +847,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 cde234af1fcfe35fa66938cbc30b5a37942f3392..a68b66b90edaf331b23a3d16218d013856f9c1f2 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
@@ -185,35 +185,59 @@ 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
+ // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
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);
+ int nbnormal = 0;
+ for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
+ FaceId node_face_id = node_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 (node_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;
- normal *= 1. / node_cell_list.size();
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];
@@ -246,7 +270,7 @@ class RusanovEulerianCompositeSolver_v2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
@@ -260,27 +284,51 @@ class RusanovEulerianCompositeSolver_v2
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];
@@ -302,6 +350,7 @@ class RusanovEulerianCompositeSolver_v2
boundary_condition);
}
}
+
void
_applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list,
const MeshType& mesh,
@@ -483,7 +532,7 @@ class RusanovEulerianCompositeSolver_v2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
@@ -590,7 +639,7 @@ class RusanovEulerianCompositeSolver_v2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
index 7386f4e84451e2d1b99d284df97eb5f9f44693d1..e1c4f46ab1d354587f7e23d54ca430cd0b937b43 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
@@ -193,34 +193,60 @@ 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
+ // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
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);
+ int nbnormal = 0;
+
+ for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
+ FaceId node_face_id = node_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 (node_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. / node_cell_list.size();
+
+ if (nbnormal == 0)
+ continue;
+ normal *= 1. / nbnormal;
+
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];
@@ -253,7 +279,7 @@ class RusanovEulerianCompositeSolver_v2_o2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
@@ -266,27 +292,51 @@ class RusanovEulerianCompositeSolver_v2_o2
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];
@@ -488,7 +538,7 @@ class RusanovEulerianCompositeSolver_v2_o2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
@@ -595,7 +645,7 @@ class RusanovEulerianCompositeSolver_v2_o2
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[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;