diff --git a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp
index 949397d46fb7502a16b0c8f88e7abb1ae5e55585..0ddfa777494f0ac15991a1db9d47a89668b9a5f4 100644
--- a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp
+++ b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp
@@ -422,36 +422,7 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
return compute_mes_j;
}();
- const CellValue<const double> primal_Vj = mesh_data.Vj();
-
- FaceValue<const double> alpha_lambda_l = [&] {
- CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
- parallel_for(
- diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
- alpha_j[diamond_cell_id] =
- dual_mes_l_j[diamond_cell_id] * dual_lambdaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
- });
-
- FaceValue<double> computed_alpha_l{mesh->connectivity()};
- mapper->fromDualCell(alpha_j, computed_alpha_l);
- return computed_alpha_l;
- }();
-
- FaceValue<const double> alpha_mu_l = [&] {
- CellValue<double> alpha_j{diamond_mesh->connectivity()};
-
- parallel_for(
- diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
- alpha_j[diamond_cell_id] =
- dual_mes_l_j[diamond_cell_id] * dual_muj[diamond_cell_id] / dual_Vj[diamond_cell_id];
- });
-
- FaceValue<double> computed_alpha_l{mesh->connectivity()};
- mapper->fromDualCell(alpha_j, computed_alpha_l);
- return computed_alpha_l;
- }();
-
+ const CellValue<const double> primal_Vj = mesh_data.Vj();
FaceValue<const CellId> face_dual_cell_id = [=]() {
FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
@@ -479,24 +450,100 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
return computed_dual_node_primal_node_id;
}();
+ 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();
+ 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);
+ }
+ }
+ }
+ }
+ }
+ });
+ return computedClj;
+ }();
+
+ 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_lambda_l = [&] {
+ CellValue<double> alpha_j{diamond_mesh->connectivity()};
+
+ parallel_for(
+ diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
+ alpha_j[diamond_cell_id] = dual_lambdaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
+ });
+
+ FaceValue<double> computed_alpha_l{mesh->connectivity()};
+ mapper->fromDualCell(alpha_j, computed_alpha_l);
+ return computed_alpha_l;
+ }();
+
+ FaceValue<const double> alpha_mu_l = [&] {
+ CellValue<double> alpha_j{diamond_mesh->connectivity()};
+
+ parallel_for(
+ diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
+ alpha_j[diamond_cell_id] = dual_muj[diamond_cell_id] / dual_Vj[diamond_cell_id];
+ });
+
+ FaceValue<double> computed_alpha_l{mesh->connectivity()};
+ mapper->fromDualCell(alpha_j, computed_alpha_l);
+ return computed_alpha_l;
+ }();
+
TinyMatrix<Dimension, double> eye = identity;
- const auto& primal_face_cell_is_reversed = mesh->connectivity().cellFaceIsReversed();
- const auto& face_local_numbers_in_their_cells = mesh->connectivity().faceLocalNumbersInTheirCells();
SparseMatrixDescriptor S{number_of_dof * Dimension};
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- const double beta_mu_l = (1. / Dimension) * alpha_mu_l[face_id] * mes_l[face_id];
- const double beta_lambda_l = (1. / Dimension) * alpha_lambda_l[face_id] * mes_l[face_id];
+ const double beta_mu_l = l2Norm(dualClj[face_id]) * alpha_mu_l[face_id] * mes_l[face_id];
+ const double beta_lambda_l = l2Norm(dualClj[face_id]) * alpha_lambda_l[face_id] * mes_l[face_id];
const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
- const CellId i_id = primal_face_to_cell[i_cell];
- const bool is_face_reversed_for_cell_i =
- primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_cell));
+ const CellId i_id = primal_face_to_cell[i_cell];
+ const bool is_face_reversed_for_cell_i = ((dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
+
const TinyVector<Dimension> nil = [&] {
if (is_face_reversed_for_cell_i) {
- return -primal_nlr(face_id, 0);
+ return -nlj[face_id];
} else {
- return primal_nlr(face_id, 0);
+ return nlj[face_id];
}
}();
for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
@@ -531,27 +578,24 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
}
}
- const auto& dual_Cjr = diamond_mesh_data.Cjr();
const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- const double alpha_mu_face_id = alpha_mu_l[face_id];
- const double alpha_lambda_face_id = alpha_lambda_l[face_id];
+ const double alpha_mu_face_id = mes_l[face_id] * alpha_mu_l[face_id];
+ const double alpha_lambda_face_id = mes_l[face_id] * alpha_lambda_l[face_id];
for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
- CellId i_id = face_to_cell_matrix[face_id][i_face_cell];
- const bool is_face_reversed_for_cell_i =
- primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_face_cell));
+ CellId i_id = face_to_cell_matrix[face_id][i_face_cell];
+ const bool is_face_reversed_for_cell_i = ((dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
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];
- // if (not primal_node_is_on_boundary[node_id]) {
const TinyVector<Dimension> nil = [&] {
if (is_face_reversed_for_cell_i) {
- return -primal_nlr(face_id, i_node);
+ return -nlj[face_id];
} else {
- return primal_nlr(face_id, i_node);
+ return nlj[face_id];
}
}();
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
index 3a30d04856e21eb5099456d39729d96ca4f9bcdb..5f66eb48e79378b25ebab8191f39b88625a8973d 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -514,9 +514,6 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
- const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
- // const auto& primal_face_cell_is_reversed = mesh->connectivity().cellFaceIsReversed();
- // const auto& face_local_numbers_in_their_cells = mesh->connectivity().faceLocalNumbersInTheirCells();
const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
const double alpha_face_id = mes_l[face_id] * alpha_l[face_id];
@@ -524,17 +521,12 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
CellId i_id = face_to_cell_matrix[face_id][i_face_cell];
const bool is_face_reversed_for_cell_i = ((dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
- // primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_face_cell));
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];
const TinyVector<Dimension> nil = [&] {
if (is_face_reversed_for_cell_i) {
- // return -primal_nlr(face_id, i_node);
- // } else {
- // return primal_nlr(face_id, i_node);
- // }
return -nlj[face_id];
} else {
return nlj[face_id];