modification of the mesure of the face in 3D and some geometrical improvments...

modification of the mesure of the face in 3D and some geometrical improvments in HeatDiamond (ongoing)

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_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));
+          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];
               }
             }();
 

src/mesh/MeshData.hpp

@@ -60,12 +60,24 @@ class MeshData : public IMeshData
           const auto& face_nodes = face_to_node_matrix[face_id];
 
           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;
     }