diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index 38ef1f5e86bceff24452643f35e39e1d2f5743e6..8f0a00407ffc5b5c3a93034a4ac8185fb91dde35 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -8,8 +8,6 @@ class ScalarHybridSchemeHandler::IScalarHybridScheme
public:
virtual std::shared_ptr<const IDiscreteFunction> getSolution() const = 0;
- // virtual std::shared_ptr<const IDiscreteFunction> getSolution() const = 0;
-
IScalarHybridScheme() = default;
virtual ~IScalarHybridScheme() = default;
};
@@ -23,7 +21,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
using MeshDataType = MeshData<Dimension>;
std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_cell_temperature;
- // std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_solution;
class DirichletBoundaryCondition
{
@@ -156,11 +153,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
{
return {m_cell_temperature};
}
- // std::shared_ptr<const IDiscreteFunction>
- // getSolution() const final
- // {
- // return m_solution;
- // }
ScalarHybridScheme(const std::shared_ptr<const MeshType>& mesh,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k,
@@ -332,7 +324,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
throw NormalError(error_msg.str());
}
}
- // std::cout << "Test 1.2 \n";
+
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
@@ -346,7 +338,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
const auto is_boundary_node = mesh->connectivity().isBoundaryNode();
const auto is_boundary_face = mesh->connectivity().isBoundaryFace();
- // std::cout << "Test 1.3 \n";
const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
@@ -359,7 +350,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
const auto& nl = mesh_data.nl();
- // std::cout << "Test 1.4 \n";
const NodeValue<const bool> node_is_neumann = [&] {
NodeValue<bool> compute_node_is_neumann{mesh->connectivity()};
compute_node_is_neumann.fill(false);
@@ -383,7 +373,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_node_is_neumann;
}();
- // std::cout << "Test 1.5 \n";
const FaceValue<const bool> face_is_neumann = [&] {
FaceValue<bool> compute_face_is_neumann{mesh->connectivity()};
@@ -404,7 +393,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_face_is_neumann;
}();
- // std::cout << "Test 1.6 \n";
const NodeValue<const bool> node_is_dirichlet = [&] {
NodeValue<bool> compute_node_is_dirichlet{mesh->connectivity()};
@@ -427,7 +415,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_node_is_dirichlet;
}();
- // std::cout << "Test 1.7 \n";
const FaceValue<const bool> face_is_dirichlet = [&] {
FaceValue<bool> compute_face_is_dirichlet{mesh->connectivity()};
@@ -448,7 +435,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_face_is_dirichlet;
}();
- // std::cout << "Test 1.8 \n";
const NodeValue<const bool> node_is_corner = [&] {
NodeValue<bool> compute_node_is_corner{mesh->connectivity()};
@@ -481,7 +467,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
});
return compute_node_is_edge;
}();
- // std::cout << "Test 1.9 \n";
const NodeValue<const bool> node_is_angle = [&] {
NodeValue<bool> compute_node_is_angle{mesh->connectivity()};
@@ -507,12 +492,10 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
if (l2Norm(n1 - n2) > (1.E-15) and l2Norm(n1 + n2) > (1.E-15) and not node_is_corner[node_id]) {
compute_node_is_angle[node_id] = true;
}
- // std ::cout << node_id << " " << n1 << " " << n2 << " " << compute_node_is_angle[node_id] << "\n";
}
});
return compute_node_is_angle;
}();
- // std::cout << "Test 1.10 \n";
const NodeValue<const TinyVector<Dimension>> exterior_normal = [&] {
NodeValue<TinyVector<Dimension>> compute_exterior_normal{mesh->connectivity()};
@@ -528,16 +511,10 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
const double norm_exterior_normal = l2Norm(compute_exterior_normal[node_id]);
compute_exterior_normal[node_id] *= 1. / norm_exterior_normal;
- //// std::cout << node_id << " " << compute_exterior_normal[node_id] << "\n";
}
}
return compute_exterior_normal;
}();
- // std::cout << "Test 1.11 \n";
-
- // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
- // // std::cout << node_id << " " << exterior_normal[node_id] << " " << node_is_angle[node_id] << "\n";
- //};
const FaceValue<const double> mes_l = [&] {
if constexpr (Dimension == 1) {
@@ -548,7 +525,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
return mesh_data.ll();
}
}();
- // std::cout << "Test 1.12 \n";
const NodeValue<const TinyVector<Dimension>> normal_angle_base = [&] {
NodeValue<TinyVector<Dimension>> compute_normal_base{mesh->connectivity()};
@@ -588,7 +564,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_normal_base;
}();
- // std::cout << "Test 1.13 \n";
const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
@@ -608,7 +583,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
});
return kappa;
}();
- // std::cout << "Test 1.14 \n";
const FaceValue<const TinyMatrix<Dimension>> face_kappal = [&] {
FaceValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
@@ -624,7 +598,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
});
return kappa;
}();
- // std::cout << "Test 1.15 \n";
const FaceValuePerCell<const TinyVector<Dimension>> xj1_xj2 = [&] {
FaceValuePerCell<TinyVector<Dimension>> compute_xj1_xj2{mesh->connectivity()};
@@ -660,14 +633,12 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_xj1_xj2;
}();
- // std::cout << "Test 1.16 \n";
const FaceValuePerCell<const TinyVector<Dimension>> xj1_xj2_orth = [&] {
FaceValuePerCell<TinyVector<Dimension>> compute_xj1_xj2_orth{mesh->connectivity()};
compute_xj1_xj2_orth.fill(zero);
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
const auto& cell_to_face = cell_to_face_matrix[cell_id];
- //// std::cout << "Test 1.16.1 \n";
for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
const FaceId face_id = cell_to_face[i_face];
@@ -677,7 +648,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_xj1_xj2_orth;
}();
- // std::cout << "Test 1.17 \n";
const FaceValuePerCell<const TinyVector<Dimension>> w1 = [&] {
FaceValuePerCell<TinyVector<Dimension>> compute_w1{mesh->connectivity()};
@@ -694,8 +664,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
if (pow(dot(v, xj1_xj2(cell_id, i_face)), 2) == pow(l2Norm(v) * l2Norm(xj1_xj2(cell_id, i_face)), 2)) {
v[2] = 1;
}
-
- // compute_w1(cell_id, i_face) = crossProduct(v, xj1_xj2(cell_id, i_face));
compute_w1(cell_id, i_face)[0] = v[1] * xj1_xj2(cell_id, i_face)[2] - v[2] * xj1_xj2(cell_id, i_face)[1];
compute_w1(cell_id, i_face)[1] = v[2] * xj1_xj2(cell_id, i_face)[0] - v[0] * xj1_xj2(cell_id, i_face)[2];
compute_w1(cell_id, i_face)[2] = v[0] * xj1_xj2(cell_id, i_face)[1] - v[1] * xj1_xj2(cell_id, i_face)[0];
@@ -706,7 +674,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_w1;
}();
- // std::cout << "Test 1.17.1 \n";
const FaceValuePerCell<const TinyVector<Dimension>> w2 = [&] {
FaceValuePerCell<TinyVector<Dimension>> compute_w2{mesh->connectivity()};
@@ -715,8 +682,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
const auto& cell_to_face = cell_to_face_matrix[cell_id];
for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
const FaceId face_id = cell_to_face[i_face];
-
- // compute_w2(cell_id, i_face) = crossProduct(xj1_xj2(cell_id, i_face), w1(cell_id, i_face));
if (Dimension == 3) {
compute_w2(cell_id, i_face)[0] = w1(cell_id, i_face)[1] * xj1_xj2(cell_id, i_face)[2] -
w1(cell_id, i_face)[2] * xj1_xj2(cell_id, i_face)[1];
@@ -731,7 +696,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_w2;
}();
- // std::cout << "Test 1.17.2 \n";
const FaceValuePerCell<const TinyVector<Dimension>> normal_face_base = [&] {
FaceValuePerCell<TinyVector<Dimension>> compute_normal_base{mesh->connectivity()};
@@ -741,26 +705,14 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
const FaceId face_id = cell_to_face[i_face];
- // TinyMatrix<Dimension> A;
TinyVector<Dimension> face_normal;
- // A(0, 0) = xj1_xj2_orth(cell_id, i_face)[0];
- // A(0, 1) = xj1_xj2(cell_id, i_face)[0];
- // A(1, 0) = xj1_xj2_orth(cell_id, i_face)[1];
- // A(1, 1) = xj1_xj2(cell_id, i_face)[1];
-
if (dot(nl[face_id], xj1_xj2(cell_id, i_face)) > 0) {
face_normal = nl[face_id];
} else {
face_normal = -nl[face_id];
}
- // compute_normal_base(cell_id, i_face) = inverse(A) * (face_kappal[face_id] * face_normal);
-
- // std::cout << face_kappal[face_id] << "\n";
- // std::cout << xj1_xj2_orth(cell_id, i_face) << "\n";
- // std::cout << face_normal << "\n";
-
if (Dimension == 2) {
compute_normal_base(cell_id, i_face)[1] =
dot(face_kappal[face_id] * face_normal, xj1_xj2_orth(cell_id, i_face)) /
@@ -777,15 +729,11 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
dot(face_kappal[face_id] * face_normal, xj1_xj2(cell_id, i_face)) /
pow(l2Norm(xj1_xj2(cell_id, i_face)), 2);
}
-
- // std::cout << cell_id << " " << face_id << "; alpha_l = " << compute_normal_base(cell_id, i_face)[0]
- // << "; delta_l = " << compute_normal_base(cell_id, i_face)[1] << "\n";
}
}
return compute_normal_base;
}();
- // std::cout << "Test 1.18 \n";
const FaceValuePerCell<const TinyVector<Dimension>> base_xj1_xj2_orth = [&] {
FaceValuePerCell<TinyVector<Dimension>> compute_w{mesh->connectivity()};
@@ -818,16 +766,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
bool l1_OK = false;
if (node_is_edge[node_id]) {
- // const auto& node_to_cell = node_to_cell_matrix[node_id];
- // CellId cell_id = node_to_cell[0];
- // const auto& cell_to_node = cell_to_node_matrix[cell_id];
- // for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
- // const NodeId node_id_in_cell = cell_to_node[i_node];
- // if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) {
- // cell_id = node_to_cell[1];
- // }
- // }
-
const auto& cell_to_face = cell_to_face_matrix[cell_id];
for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
const FaceId face_id = cell_to_face[i_face];
@@ -862,16 +800,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
if (node_is_edge[node_id]) {
const FaceId l1 = normal_edge_index(cell_id, i_node)[0];
- // const auto& node_to_cell = node_to_cell_matrix[node_id];
- // CellId cell_id = node_to_cell[0];
- // const auto& cell_to_node = cell_to_node_matrix[cell_id];
- // for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
- // const NodeId node_id_in_cell = cell_to_node[i_node];
- // if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) {
- // cell_id = node_to_cell[1];
- // }
- // }
-
if (dot(nl[l1], xl[l1] - xj[cell_id]) <= 0) {
compute_normal1_edge(cell_id, i_node) = -nl[l1];
} else {
@@ -893,16 +821,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
const NodeId node_id = cell_to_node[i_node];
if (node_is_edge[node_id]) {
const FaceId l2 = normal_edge_index(cell_id, i_node)[1];
-
- // const auto& node_to_cell = node_to_cell_matrix[node_id];
- // CellId cell_id = node_to_cell[0];
- // const auto& cell_to_node = cell_to_node_matrix[cell_id];
- // for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
- // const NodeId node_id_in_cell = cell_to_node[i_node];
- // if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) {
- // cell_id = node_to_cell[1];
- // }
- // }
if (dot(nl[l2], xl[l2] - xj[cell_id]) <= 0) {
compute_normal2_edge(cell_id, i_node) = -nl[l2];
} else {
@@ -972,7 +890,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_face_boundary_values;
}();
- // std::cout << "Test 1.19 \n";
const FaceValue<const double> face_boundary_prev_values = [&] {
FaceValue<double> compute_face_boundary_values{mesh->connectivity()};
@@ -1007,7 +924,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_face_boundary_values;
}();
- // std::cout << "Test 1.20 \n";
const NodeValue<const double> node_boundary_values = [&] {
NodeValue<double> compute_node_boundary_values{mesh->connectivity()};
@@ -1062,11 +978,9 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id]) and not node_is_angle[node_id]) {
compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
}
- //// std::cout << node_id << " " << compute_node_boundary_values[node_id] << "\n";
}
return compute_node_boundary_values;
}();
- // std::cout << "Test 1.21 \n";
const NodeValue<const double> node_boundary_prev_values = [&] {
NodeValue<double> compute_node_boundary_values{mesh->connectivity()};
@@ -1123,7 +1037,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
return compute_node_boundary_values;
}();
- // std::cout << "Test 1.22 \n";
const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] {
NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
@@ -1141,21 +1054,17 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
});
return beta;
}();
- // std::cout << "Test 1.23 \n";
const NodeValue<const TinyMatrix<Dimension>> kappar_invBetar = [&] {
NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
- // std::cout << node_id << " " << node_is_corner[node_id] << "\n";
if (not node_is_corner[node_id]) {
- // std::cout << node_id << " " << node_betar[node_id];
kappa_invBeta[node_id] = node_kappar[node_id] * inverse(node_betar[node_id]);
}
});
return kappa_invBeta;
}();
- // std::cout << "Test 1.24 \n";
const NodeValue<const TinyMatrix<Dimension>> corner_betar = [&] {
NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
@@ -1168,7 +1077,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
size_t i_cell = 0;
const CellId cell_id = node_to_cell[i_cell];
if (Dimension == 2) {
- // std::cout << "Test 1.24.1 \n";
const unsigned int i_node = node_local_number_in_its_cell[i_cell];
const auto& cell_nodes = cell_to_node_matrix[cell_id];
@@ -1192,15 +1100,10 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
Cjr2[1] = 0.5 * (xrp2[0] - xrp1[0]);
beta[node_id] = tensorProduct(Cjr1, (xr[node_id] - xj1)) + tensorProduct(Cjr2, (xr[node_id] - xj2));
} else if (Dimension == 3 && not node_is_edge[node_id]) {
- // std::cout << "Test 1.24.2 \n";
- // std::cout << "coin : " << node_id << " " << Vj[cell_id] << "\n";
const TinyMatrix<Dimension> I = identity;
beta[node_id] = 0.125 * Vj[cell_id] * I;
- // throw NotImplementedError("The scheme is not implemented in this dimension.");
} else if (Dimension == 3 && node_is_edge[node_id]) {
- // std::cout << "Test 1.24.3 \n";
- // std::cout << " bord droit : " << node_id << " " << Vj[cell_id] << "\n";
const TinyMatrix<Dimension> I = identity;
beta[node_id] += 0.125 * Vj[cell_id] * I;
}
@@ -1208,7 +1111,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
});
return beta;
}();
- // std::cout << "Test 1.25 \n";
const NodeValue<const TinyMatrix<Dimension>> corner_kappar_invBetar = [&] {
NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
@@ -1220,7 +1122,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
});
return kappa_invBeta;
}();
- // std::cout << "Test 1.26 \n";
const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] {
NodeValue<TinyVector<Dimension>> compute_sum_Cjr{mesh->connectivity()};
@@ -1238,7 +1139,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
});
return compute_sum_Cjr;
}();
- // std::cout << "Test 1.27 \n";
const NodeValue<const TinyVector<Dimension>> v = [&] {
NodeValue<TinyVector<Dimension>> compute_v{mesh->connectivity()};
@@ -1251,7 +1151,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
});
return compute_v;
}();
- // std::cout << "Test 1.28 \n";
const NodeValuePerCell<const double> theta = [&] {
NodeValuePerCell<double> compute_theta{mesh->connectivity()};
@@ -1267,7 +1166,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
});
return compute_theta;
}();
- // std::cout << "Test 1.29 \n";
const NodeValue<const double> sum_theta = [&] {
NodeValue<double> compute_sum_theta{mesh->connectivity()};
@@ -1287,11 +1185,9 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
});
return compute_sum_theta;
}();
- // std::cout << "Test 1.30 \n";
- // std::cout << "Test 1.1 \n";
const Array<int> non_zeros{mesh->numberOfCells()};
- non_zeros.fill(4); // Modif pour optimiser
+ non_zeros.fill(Dimension);
CRSMatrixDescriptor<double> S1(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
CRSMatrixDescriptor<double> S2(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
@@ -1468,7 +1364,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
pow(l2Norm(node_kappar[node_id] * t), 2)) *
node_kappar[node_id] * t,
base_xj1_xj2_orth(cell_id_j, i_face));
- // std::cout << "Test 1.3.3 \n";
}
}
}
@@ -1593,8 +1488,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
}
};
- // std::cout << "Test 1.31" << "\n";
- //// std::cout << "S1 = " << S1 << "\n";
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
S1(cell_id, cell_id) += Vj[cell_id];
@@ -1611,8 +1504,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
Ph[cell_id] = Pj[cell_id];
};
- // std::cout << "P = " << Ph << "\n";
-
CRSMatrix A1{S1.getCRSMatrix()};
CRSMatrix A2{S2.getCRSMatrix()};
@@ -1663,7 +1554,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
lambda * mes_l[face_id] / numberOfNodesPerFace * node_boundary_values[node_id] /
(sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], base_xj1_xj2_orth(cell_id, i_face));
- // std::cout << "xl - xj2 = " << xl[face_id] -
+
b_prev[cell_id] +=
lambda * mes_l[face_id] / numberOfNodesPerFace * node_boundary_prev_values[node_id] /
(sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
@@ -1676,8 +1567,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
} else if (node_is_corner[node_id] && not node_is_edge[node_id]) {
const auto& node_to_face = node_to_face_matrix[node_id];
- // const CellId cell_id = node_to_cell[0];
- const CellId cell_id = node_to_cell[0];
+ const CellId cell_id = node_to_cell[0];
for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
FaceId face_id = node_to_face[i_face];
const auto& face_to_node = face_to_node_matrix[face_id];
@@ -1689,7 +1579,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
} else if (node_is_edge[node_id]) {
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
- // const auto& node_to_cell = node_to_cell_matrix[node_id];
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
@@ -1795,7 +1684,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
b[cell_id_j] +=
lambda * mes_l[face_id] / numberOfNodesPerFace * node_boundary_values[node_id] *
dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], base_xj1_xj2_orth(cell_id_j, i_face));
- // std::cout << "xl - xj2 = " << xl[face_id] -
b_prev[cell_id_j] +=
lambda * mes_l[face_id] / numberOfNodesPerFace * node_boundary_prev_values[node_id] *
dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], base_xj1_xj2_orth(cell_id_j, i_face));
@@ -1871,7 +1759,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
Vector<double> T{mesh->numberOfCells()};
T = zero;
- // Array<double> Ph2 = Ph;
Vector<double> A2P = A2 * Ph;
Vector<double> B{mesh->numberOfCells()};
@@ -1892,7 +1779,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
ur[node_id] += Pj[cell_id] * Cjr(cell_id, i_node);
}
ur[node_id] = -inverse(node_betar[node_id]) * ur[node_id];
- // std::cout << node_id << "; ur = " << ur[node_id] << "\n";
} else if (is_boundary_node[node_id] and node_is_dirichlet[node_id]) {
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
@@ -1904,7 +1790,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
} else {
ur[node_id] = inverse(corner_betar[node_id]) * ur[node_id];
}
- // std::cout << "bord : " << node_id << "; ur = " << ur[node_id] << "\n";
}
}
@@ -1916,7 +1801,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
compute_primal_node_temperature.fill(0);
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
- // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); node_id++) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
double sum_Vj = 0;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
@@ -1929,9 +1813,8 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
return compute_primal_node_temperature;
}();
- // m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
m_cell_temperature = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
- // auto& solution = *m_solution;
+
auto& cell_temperature = *m_cell_temperature;
parallel_for(
@@ -1939,11 +1822,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
};
-// std::shared_ptr<const IDiscreteFunction>
-// ScalarHybridSchemeHandler::solution() const
-// {
-// return m_scheme->getSolution();
-// }
std::shared_ptr<const IDiscreteFunction>
ScalarHybridSchemeHandler::solution() const
{
diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index 0897f6d8323c2cccd368f5cedb23fcdf33b3c7a3..4953c8f31ef5f00a0b5b12990b843c2e0123a874 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -9,8 +9,6 @@ class ScalarNodalSchemeHandler::IScalarNodalScheme
virtual std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> getSolution()
const = 0;
- // virtual std::shared_ptr<const IDiscreteFunction> getSolution() const = 0;
-
IScalarNodalScheme() = default;
virtual ~IScalarNodalScheme() = default;
};
@@ -25,7 +23,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_cell_temperature;
std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_node_temperature;
- // std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_solution;
class DirichletBoundaryCondition
{
@@ -158,11 +155,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
{
return {m_cell_temperature, m_node_temperature};
}
- // std::shared_ptr<const IDiscreteFunction>
- // getSolution() const final
- // {
- // return m_solution;
- // }
ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k,
@@ -351,7 +343,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const auto is_boundary_node = mesh->connectivity().isBoundaryNode();
const auto is_boundary_face = mesh->connectivity().isBoundaryFace();
- // const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
@@ -444,7 +435,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const NodeValue<const bool> node_is_angle = [&] {
NodeValue<bool> compute_node_is_angle{mesh->connectivity()};
compute_node_is_angle.fill(false);
- // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
if (is_boundary_node[node_id]) {
@@ -466,7 +456,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
if (l2Norm(n1 - n2) > (1.E-15) and l2Norm(n1 + n2) > (1.E-15) and not node_is_corner[node_id]) {
compute_node_is_angle[node_id] = true;
}
- // std::cout << node_id << " " << n1 << " " << n2 << " " << compute_node_is_angle[node_id] << "\n";
}
});
return compute_node_is_angle;
@@ -486,16 +475,11 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
const double norm_exterior_normal = l2Norm(compute_exterior_normal[node_id]);
compute_exterior_normal[node_id] *= 1. / norm_exterior_normal;
- // std::cout << node_id << " " << compute_exterior_normal[node_id] << "\n";
}
}
return compute_exterior_normal;
}();
- // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
- // std::cout << node_id << " " << exterior_normal[node_id] << " " << node_is_angle[node_id] << "\n";
- //};
-
const FaceValue<const double> mes_l = [&] {
if constexpr (Dimension == 1) {
FaceValue<double> compute_mes_l{mesh->connectivity()};
@@ -598,7 +582,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id]) and not node_is_angle[node_id]) {
compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
}
- // std::cout << node_id << " " << compute_node_boundary_values[node_id] << "\n";
}
return compute_node_boundary_values;
}();
@@ -825,10 +808,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return compute_sum_theta;
}();
- // const double lambda = 0;
-
const Array<int> non_zeros{mesh->numberOfCells()};
- non_zeros.fill(4); // Modif pour optimiser
+ non_zeros.fill(Dimension);
CRSMatrixDescriptor<double> S1(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
CRSMatrixDescriptor<double> S2(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
@@ -945,8 +926,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
CRSMatrix A1{S1.getCRSMatrix()};
CRSMatrix A2{S2.getCRSMatrix()};
- // std::cout << A1 << "\n";
-
Vector<double> b{mesh->numberOfCells()};
b = zero;
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
@@ -1034,7 +1013,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
Vector<double> T{mesh->numberOfCells()};
T = zero;
- // Array<double> Ph2 = Ph;
Vector<double> A2P = A2 * Ph;
Vector<double> B{mesh->numberOfCells()};
@@ -1043,8 +1021,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
B[cell_id] = dt * ((1. - cn_coeff / 2.) * b[cell_id] + cn_coeff / 2. * b_prev[cell_id]) + A2P[cell_id];
});
- // std::cout << "g = " << node_boundary_values << "\n";
-
NodeValue<TinyVector<Dimension>> ur{mesh->connectivity()};
ur.fill(zero);
for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
@@ -1057,7 +1033,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
ur[node_id] += Pj[cell_id] * Cjr(cell_id, i_node);
}
ur[node_id] = -inverse(node_betar[node_id]) * ur[node_id];
- // std::cout << node_id << "; ur = " << ur[node_id] << "\n";
} else if (is_boundary_node[node_id] and node_is_dirichlet[node_id]) {
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
@@ -1069,7 +1044,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
} else {
ur[node_id] = inverse(corner_betar[node_id]) * ur[node_id];
}
- // std::cout << "bord : " << node_id << "; ur = " << ur[node_id] << "\n";
}
}
@@ -1081,7 +1055,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
compute_primal_node_temperature.fill(0);
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
- // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); node_id++) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
double sum_Vj = 0;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
@@ -1100,10 +1073,9 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return my_dual_temperature;
}();
- // m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
m_cell_temperature = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
m_node_temperature = std::make_shared<DiscreteFunctionP0<Dimension, double>>(dual_mesh);
- // auto& solution = *m_solution;
+
auto& cell_temperature = *m_cell_temperature;
auto& node_temperature = *m_node_temperature;
parallel_for(
@@ -1115,11 +1087,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
};
-// std::shared_ptr<const IDiscreteFunction>
-// ScalarNodalSchemeHandler::solution() const
-// {
-// return m_scheme->getSolution();
-// }
std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>>
ScalarNodalSchemeHandler::solution() const
{