diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index fd4e2d5bd210e0aa3e61759c982770ce47685124..791e1c8ff7b94c39663bea036faff5ecc5a08b51 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -128,8 +128,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& alpha,
- const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& nod_k_bound,
- const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& nod_k,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_bound,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
@@ -246,7 +246,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return compute_face_dof_number;
}();
- const FaceValue<const bool> primal_face_is_neumann = [&] {
+ const FaceValue<const bool> face_is_neumann = [&] {
FaceValue<bool> face_is_neumann{mesh->connectivity()};
face_is_neumann.fill(false);
for (const auto& boundary_condition : boundary_condition_list) {
@@ -269,64 +269,27 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return face_is_neumann;
}();
- const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
- const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
- NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
- if (parallel::size() > 1) {
- throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
- }
-
- primal_node_is_on_boundary.fill(false);
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- if (face_to_cell_matrix[face_id].size() == 1) {
- for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
- NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
- primal_node_is_on_boundary[node_id] = true;
- }
- }
- }
-
- FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
- if (parallel::size() > 1) {
- throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
- }
-
- primal_face_is_on_boundary.fill(false);
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- if (face_to_cell_matrix[face_id].size() == 1) {
- primal_face_is_on_boundary[face_id] = true;
- }
- }
-
- FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
- if (parallel::size() > 1) {
- throw NotImplementedError("Calculation of face_is_neumann is incorrect");
- }
-
- primal_face_is_dirichlet.fill(false);
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]));
- }
-
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
- const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
+ const NodeValue<const TinyVector<Dimension>>& xr = mesh_data.xr();
- {
- std::shared_ptr diamond_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
+ const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
- MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
+ const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
- std::shared_ptr mapper =
- DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity());
+ {
- CellValue<const TinyMatrix<Dimension>> nod_kappar = nod_k->cellValues();
- CellValue<const TinyMatrix<Dimension>> nod_kapparb = nod_k_bound->cellValues();
- const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
+ CellValue<const TinyMatrix<Dimension>> cell_kappaj = cell_k->cellValues();
+ CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_bound->cellValues();
+
+ const CellValue<const double> Vj = mesh_data.Vj();
+ const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
const FaceValue<const double> mes_l = [&] {
if constexpr (Dimension == 1) {
@@ -338,92 +301,59 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
}();
- const CellValue<const double> dual_mes_l_j = [=] {
- CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
- mapper->toDualCell(mes_l, compute_mes_j);
-
- return compute_mes_j;
- }();
-
- FaceValue<const CellId> face_dual_cell_id = [=]() {
- FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
- CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
+ const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
+ NodeValue<const TinyMatrix<Dimension>> kappa;
parallel_for(
- diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
-
- mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id);
-
- return computed_face_dual_cell_id;
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ kappa[node_id] = zero;
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ double weight = 0;
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ double local_weight = 1 / l2Norm(xr[node_id]-xj[J]);
+ kappa[node_id] += cell_kappaj[J] * local_weight;
+ weight += local_weight;
+ }
+ kappa[node_id] /= weight;
+ }
+ );
+ return kappa;
}();
- NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
- CellValue<NodeId> cell_ignored_id{mesh->connectivity()};
- cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
- NodeValue<NodeId> node_primal_id{mesh->connectivity()};
-
+ const NodeValue<const TinyMatrix<Dimension>> node_kapparb = [&] {
+ NodeValue<const TinyMatrix<Dimension>> kappa;
parallel_for(
- mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
-
- NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
-
- mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
-
- return computed_dual_node_primal_node_id;
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ kappa[node_id] = zero;
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ double weight = 0;
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ double local_weight = 1 / l2Norm(xr[node_id]-xj[J]);
+ kappa[node_id] += cell_kappajb[J] * local_weight;
+ weight += local_weight;
+ }
+ kappa[node_id] /= weight;
+ });
+ return kappa;
}();
- CellValue<NodeId> primal_cell_dual_node_id = [=]() {
- CellValue<NodeId> cell_id{mesh->connectivity()};
- NodeValue<NodeId> node_ignored_id{mesh->connectivity()};
- node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
-
- NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()};
-
- parallel_for(
- diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; });
-
- CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()};
-
- mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id);
-
- return cell_id;
- }();
- const auto& dual_Cjr = diamond_mesh_data.Cjr();
- FaceValue<TinyVector<Dimension>> dualClj = [&] {
- FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()};
- const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix();
- const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
+ const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] {
+ NodeValue<const TinyMatrix<Dimension>> beta;
parallel_for(
- mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
- const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
- for (size_t i = 0; i < primal_face_to_cell.size(); i++) {
- CellId cell_id = primal_face_to_cell[i];
- const NodeId dual_node_id = primal_cell_dual_node_id[cell_id];
- for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size();
- i_dual_cell++) {
- const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell];
- if (face_dual_cell_id[face_id] == dual_cell_id) {
- for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
- i_dual_node++) {
- const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
- if (final_dual_node_id == dual_node_id) {
- computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node);
- }
- }
- }
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId 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];
+ beta[node_id] = zero;
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cell[j];
+ beta[node_id] += tensorProduct(Cjr(J,R),xr[node_id]-xj[J]);
}
- }
});
- return computedClj;
+ return beta;
}();
- FaceValue<TinyVector<Dimension>> nlj = [&] {
- FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()};
- parallel_for(
- mesh->numberOfFaces(),
- PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; });
- return computedNlj;
- }();
FaceValue<const double> alpha_l = [&] {
FaceValue<double> alpha_j{mesh->connectivity()};
@@ -447,23 +377,22 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return alpha_lb;
}();
- const CellValue<const double> primal_Vj = mesh_data.Vj();
const Array<int> non_zeros{number_of_dof};
- non_zeros.fill(Dimension);
+ non_zeros.fill(Dimension); //Modif pour optimiser
CRSMatrixDescriptor<double> S(number_of_dof, number_of_dof, non_zeros);
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
+ const auto& face_to_cell = face_to_cell_matrix[face_id];
const double beta_l = l2Norm(dualClj[face_id]) * alpha_l[face_id] * mes_l[face_id];
const double betab_l = l2Norm(dualClj[face_id]) * alphab_l[face_id] * mes_l[face_id];
- for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
- const CellId cell_id1 = primal_face_to_cell[i_cell];
- for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
- const CellId cell_id2 = primal_face_to_cell[j_cell];
+ for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+ const CellId cell_id1 = face_to_cell[i_cell];
+ for (size_t j_cell = 0; j_cell < face_to_cell.size(); ++j_cell) {
+ const CellId cell_id2 = face_to_cell[j_cell];
if (i_cell == j_cell) {
S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
- if (primal_face_is_neumann[face_id]) {
+ if (face_is_neumann[face_id]) {
S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= betab_l;
}
} else {
@@ -475,12 +404,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
const size_t j = cell_dof_number[cell_id];
- S(j, j) += (*alpha)[cell_id] * primal_Vj[cell_id];
+ S(j, j) += (*alpha)[cell_id] * Vj[cell_id];
}
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- if (primal_face_is_dirichlet[face_id]) {
+ if (face_is_dirichlet[face_id]) {
S(face_dof_number[face_id], face_dof_number[face_id]) += 1;
}
}
@@ -491,7 +420,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
Vector<double> b{number_of_dof};
b = zero;
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- b[cell_dof_number[cell_id]] = fj[cell_id] * primal_Vj[cell_id];
+ b[cell_dof_number[cell_id]] = fj[cell_id] * Vj[cell_id];
}
// Dirichlet on b^L_D
{
@@ -552,8 +481,8 @@ ScalarNodalSchemeHandler::solution() const
ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
const std::shared_ptr<const IDiscreteFunction>& alpha,
- const std::shared_ptr<const IDiscreteFunction>& nod_k_bound,
- const std::shared_ptr<const IDiscreteFunction>& nod_k,
+ const std::shared_ptr<const IDiscreteFunction>& k_bound,
+ const std::shared_ptr<const IDiscreteFunction>& cell_k,
const std::shared_ptr<const IDiscreteFunction>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
@@ -561,11 +490,11 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
if (not i_mesh) {
throw NormalError("primal discrete functions are not defined on the same mesh");
}
- const std::shared_ptr i_dual_mesh = getCommonMesh({nod_k_bound, nod_k});
+ const std::shared_ptr i_dual_mesh = getCommonMesh({cell_k_bound, cell_k});
if (not i_dual_mesh) {
throw NormalError("dual discrete functions are not defined on the same mesh");
}
- checkDiscretizationType({alpha, nod_k_bound, nod_k, f}, DiscreteFunctionType::P0);
+ checkDiscretizationType({alpha, cell_k_bound, cell_k, f}, DiscreteFunctionType::P0);
switch (i_mesh->dimension()) {
case 1: {
@@ -581,8 +510,8 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
m_scheme =
std::make_unique<ScalarNodalScheme<1>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k_bound),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k),
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
break;
@@ -600,8 +529,8 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
m_scheme =
std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k_bound),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k),
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
break;
@@ -619,8 +548,8 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
m_scheme =
std::make_unique<ScalarNodalScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k_bound),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k),
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
break;