diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
index 68994e2525c9069ba08d178b7b57a7fb8208319e..06db3fc99abf7154d2285c0e78171f2bd9fdef14 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -16,11 +16,13 @@
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/SubItemValuePerItem.hpp>
#include <output/VTKWriter.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/FourierBoundaryConditionDescriptor.hpp>
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/Timer.hpp>
template <size_t Dimension>
HeatDiamondScheme<Dimension>::HeatDiamondScheme(
@@ -390,18 +392,97 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
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()};
+ 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;
+ }();
+
+ 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()};
+
+ 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;
+ }();
+
+ 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;
+ }();
+ Timer my_timer;
+ 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;
+ }();
+ std::cout << my_timer.seconds() << "\n";
+
+ 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 = [&] {
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_kappaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
+ alpha_j[diamond_cell_id] = dual_kappaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
});
FaceValue<double> computed_alpha_l{mesh->connectivity()};
mapper->fromDualCell(alpha_j, computed_alpha_l);
-
VTKWriter vtk_writer("DD_" + std::to_string(Dimension), 0.01);
vtk_writer.write(diamond_mesh,
@@ -416,7 +497,8 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
SparseMatrixDescriptor S{number_of_dof};
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
- const double beta_l = 1. / Dimension * alpha_l[face_id] * mes_l[face_id];
+ // const double beta_l = 1. / Dimension * alpha_l[face_id] * mes_l[face_id];
+ const double beta_l = l2Norm(dualClj[face_id]) * alpha_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) {
@@ -433,57 +515,32 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
}
- FaceValue<const CellId> face_dual_cell_id = [=]() {
- FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
- CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
- 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;
- }();
-
- 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()};
-
- 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;
- }();
-
- const auto& dual_Cjr = diamond_mesh_data.Cjr();
-
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();
+ // 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 = alpha_l[face_id];
+ const double alpha_face_id = mes_l[face_id] * alpha_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);
+ // 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);
+ // return -primal_nlr(face_id, i_node);
+ // } else {
+ // 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/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index c82e9ce06ec8b7d43ce1c6cc2fec282f89dc1fee..6d767e6bee0a96dba8b2e383fc182e7c47d9ba49 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -59,13 +59,25 @@ class MeshData : public IMeshData
m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
const auto& face_nodes = face_to_node_matrix[face_id];
- double lenght = 0;
+ double lenght = 0;
+ TinyVector<Dimension, double> normal = zero;
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
- lenght += l2Norm(Nlr(face_id, i_node));
+ normal += Nlr(face_id, i_node);
}
-
+ lenght = l2Norm(normal);
ll[face_id] = lenght;
});
+ // parallel_for(
+ // m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ // const auto& face_nodes = face_to_node_matrix[face_id];
+
+ // double lenght = 0;
+ // for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ // lenght += l2Norm(Nlr(face_id, i_node));
+ // }
+
+ // ll[face_id] = lenght;
+ // });
m_ll = ll;
}