From 4267419ad3bdbd07365b807c06c8c9ff4705b89e Mon Sep 17 00:00:00 2001
From: Julie PATELA <julie.patela.ocre@cea.fr>
Date: Thu, 30 Jul 2020 13:23:52 +0200
Subject: [PATCH] Add internal node's matrix to HeatDiamondAlgorithm; Pass
 second member in parameter

---
 .../algorithms/Heat5PointsAlgorithm.cpp       |  27 ++-
 .../algorithms/Heat5PointsAlgorithm.hpp       |   3 +-
 .../algorithms/HeatDiamondAlgorithm.cpp       | 156 ++++++++++++++----
 .../algorithms/HeatDiamondAlgorithm.hpp       |   3 +-
 src/language/modules/SchemeModule.cpp         |  22 +--
 5 files changed, 151 insertions(+), 60 deletions(-)

diff --git a/src/language/algorithms/Heat5PointsAlgorithm.cpp b/src/language/algorithms/Heat5PointsAlgorithm.cpp
index b46f0739f..f3ce0a3a8 100644
--- a/src/language/algorithms/Heat5PointsAlgorithm.cpp
+++ b/src/language/algorithms/Heat5PointsAlgorithm.cpp
@@ -21,7 +21,8 @@ Heat5PointsAlgorithm<Dimension>::Heat5PointsAlgorithm(
   std::shared_ptr<const IMesh> i_mesh,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
   const FunctionSymbolId& T_id,
-  const FunctionSymbolId& kappa_id)
+  const FunctionSymbolId& kappa_id,
+  const FunctionSymbolId& f_id)
 {
   using ConnectivityType = Connectivity<Dimension>;
   using MeshType         = Mesh<ConnectivityType>;
@@ -176,24 +177,18 @@ Heat5PointsAlgorithm<Dimension>::Heat5PointsAlgorithm(
           }
         }
       }
+      CellValue<double> fj =
+        InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
+
       const CellValue<const double> primal_Vj = mesh_data.Vj();
       CRSMatrix A{S};
       Vector<double> b{mesh->numberOfCells()};
-      double pi = 3.14159265;
-      for (CellId i_cell = 0; i_cell < mesh->numberOfCells(); ++i_cell) {
-        if constexpr (Dimension == 2) {
-          b[i_cell] =
-            -2 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) * primal_Vj[i_cell];
-
-        } else if constexpr (Dimension == 3) {
-          b[i_cell] = -3 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) *
-                      sin(pi * xj[i_cell][2]) * primal_Vj[i_cell];
-        }
-      }
+      parallel_for(
+        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { b[cell_id] = fj[cell_id] * primal_Vj[cell_id]; });
 
       Vector<double> T{mesh->numberOfCells()};
       T = 0;
-      PCG{b, A, A, T, 100, 1e-12};
+      PCG{b, A, A, T, 1000, 1e-12};
 
       CellValue<double> Temperature{mesh->connectivity()};
 
@@ -202,7 +197,8 @@ Heat5PointsAlgorithm<Dimension>::Heat5PointsAlgorithm(
 
       Vector<double> error{mesh->numberOfCells()};
       parallel_for(
-        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { error[cell_id] = Temperature[cell_id] - Tj[cell_id]; });
+        mesh->numberOfCells(),
+        PUGS_LAMBDA(CellId cell_id) { error[cell_id] = (Temperature[cell_id] - Tj[cell_id]) * primal_Vj[cell_id]; });
 
       std::cout << "||Error||_2 = " << std::sqrt((error, error)) << "\n";
 
@@ -222,16 +218,19 @@ template Heat5PointsAlgorithm<1>::Heat5PointsAlgorithm(
   std::shared_ptr<const IMesh>,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
   const FunctionSymbolId&,
+  const FunctionSymbolId&,
   const FunctionSymbolId&);
 
 template Heat5PointsAlgorithm<2>::Heat5PointsAlgorithm(
   std::shared_ptr<const IMesh>,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
   const FunctionSymbolId&,
+  const FunctionSymbolId&,
   const FunctionSymbolId&);
 
 template Heat5PointsAlgorithm<3>::Heat5PointsAlgorithm(
   std::shared_ptr<const IMesh>,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
   const FunctionSymbolId&,
+  const FunctionSymbolId&,
   const FunctionSymbolId&);
diff --git a/src/language/algorithms/Heat5PointsAlgorithm.hpp b/src/language/algorithms/Heat5PointsAlgorithm.hpp
index 3271ee37f..a061972d6 100644
--- a/src/language/algorithms/Heat5PointsAlgorithm.hpp
+++ b/src/language/algorithms/Heat5PointsAlgorithm.hpp
@@ -14,7 +14,8 @@ struct Heat5PointsAlgorithm
   Heat5PointsAlgorithm(std::shared_ptr<const IMesh> i_mesh,
                        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
                        const FunctionSymbolId& T_id,
-                       const FunctionSymbolId& kappa_id);
+                       const FunctionSymbolId& kappa_id,
+                       const FunctionSymbolId& f_id);
 };
 
 #endif   // HEAT_5POINTS_ALGORITHM_HPP
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
index 4c2e5d692..6a6ff4076 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -21,7 +21,8 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
   std::shared_ptr<const IMesh> i_mesh,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
   const FunctionSymbolId& T_id,
-  const FunctionSymbolId& kappa_id)
+  const FunctionSymbolId& kappa_id,
+  const FunctionSymbolId& f_id)
 {
   using ConnectivityType = Connectivity<Dimension>;
   using MeshType         = Mesh<ConnectivityType>;
@@ -115,26 +116,15 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
                        0, true);   // forces last output
 
       const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
-      const auto& face_to_node_matrix       = m_mesh->connectivity().faceToNodeMatrix();
 
-      const FaceValue<const double> mes_l = [=] {
-        FaceValue<double> compute_mes_l{m_mesh->connectivity()};
-        const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
+      const FaceValue<const double> mes_l = [&] {
         if constexpr (Dimension == 1) {
+          FaceValue<double> compute_mes_l{m_mesh->connectivity()};
           compute_mes_l.fill(1);
-        } else if constexpr (Dimension == 2) {
-          parallel_for(
-            m_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
-              const auto& face_to_node              = face_to_node_matrix[face_id];
-              const NodeId node_id1                 = face_to_node[0];
-              const NodeId node_id2                 = face_to_node[1];
-              const TinyVector<Dimension, double> r = xr[node_id1] - xr[node_id2];
-              compute_mes_l[face_id]                = l2Norm(r);
-            });
+          return compute_mes_l;
         } else {
-          throw NotImplementedError("Not implemented in 3D");
+          return mesh_data.ll();
         }
-        return compute_mes_l;
       }();
 
       const CellValue<const double> dual_mes_l_j = [=] {
@@ -167,16 +157,30 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
         return computed_alpha_l;
       }();
 
+      const auto& primal_face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix();
+      const auto& face_to_cell_matrix        = m_mesh->connectivity().faceToCellMatrix();
+
+      NodeValue<bool> primal_node_is_on_boundary(m_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 < m_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;
+          }
+        }
+      }
+
       SparseMatrixDescriptor S{m_mesh->numberOfCells()};
-      const auto& face_to_cell_matrix = m_mesh->connectivity().faceToCellMatrix();
       for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) {
         const auto& primal_face_to_cell                          = face_to_cell_matrix[face_id];
         const double beta_l                                      = 0.5 * alpha_l[face_id] * mes_l[face_id];
         const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
 
-        //          CellValue<double> dual_cell_id{diamond_mesh->connectivity()};
-        //          mapper->toDualCell(face_id, dual_cell_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) {
@@ -186,24 +190,104 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
             } else {
               S(cell_id1, cell_id2) += beta_l;
             }
-            // S(cell_id1, cell_id2)+= alpha_l[face_id]*(,Cjr[dual_cell_id,dual_node_j]);
           }
         }
       }
+
+      FaceValue<const CellId> face_dual_cell_id = [=]() {
+        FaceValue<CellId> computed_face_dual_cell_id{m_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{m_mesh->connectivity()};
+        cell_ignored_id.fill(NodeId{std::numeric_limits<NodeId>::max()});
+        NodeValue<NodeId> node_primal_id{m_mesh->connectivity()};
+
+        parallel_for(
+          m_mesh->numberOfCells(), 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                       = m_mesh->connectivity().cellFaceIsReversed();
+      const auto& face_local_numbers_in_their_cells = m_mesh->connectivity().faceLocalNumbersInTheirCells();
+      const auto& primal_node_to_cell_matrix        = m_mesh->connectivity().nodeToCellMatrix();
+
+      for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) {
+        if (face_to_cell_matrix[face_id].size() > 1) {
+          const double alpha_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));
+
+            for (size_t j_face_cell = 0; j_face_cell < face_to_cell_matrix[face_id].size(); ++j_face_cell) {
+              CellId j_id         = face_to_cell_matrix[face_id][j_face_cell];
+              CellId dual_cell_id = face_dual_cell_id[face_id];
+
+              for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++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);
+                  }
+                }();
+
+                NodeId face_node_id             = primal_face_to_node_matrix[face_id][i_node];
+                const auto& primal_node_to_cell = primal_node_to_cell_matrix[face_node_id];
+                if (not primal_node_is_on_boundary[face_node_id]) {
+                  const double weight_rj = [&] {
+                    for (size_t i_cell = 0; i_cell < primal_node_to_cell.size(); ++i_cell) {
+                      if (j_id == primal_node_to_cell[i_cell]) {
+                        return w_rj(face_node_id, i_cell);
+                      }
+                    }
+                    Assert(false, "could not determine the weight");
+                    return std::numeric_limits<double>::signaling_NaN();
+                  }();
+                  for (size_t i_dual_node = 0; i_dual_node < dual_Cjr.numberOfSubValues(dual_cell_id); ++i_dual_node) {
+                    const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
+                    if (dual_node_id == face_node_id) {
+                      const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
+                      S(i_id, j_id) += weight_rj * alpha_face_id * (nil, Clr);
+                    }
+                  }
+                } else {
+                  std::cout << "*** ignore node " << face_node_id << '\n';
+                }
+              }
+            }
+          }
+        }
+      }
+
+      CellValue<double> fj =
+        InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
+
       const CellValue<const double> primal_Vj = mesh_data.Vj();
       CRSMatrix A{S};
       Vector<double> b{m_mesh->numberOfCells()};
-      double pi = 3.14159265;
-      for (CellId i_cell = 0; i_cell < m_mesh->numberOfCells(); ++i_cell) {
-        if constexpr (Dimension == 2) {
-          b[i_cell] =
-            -2 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) * primal_Vj[i_cell];
-
-        } else if constexpr (Dimension == 3) {
-          b[i_cell] = -3 * kappaj[i_cell] * pi * pi * sin(pi * xj[i_cell][0]) * sin(pi * xj[i_cell][1]) *
-                      sin(pi * xj[i_cell][2]) * primal_Vj[i_cell];
-        }
-      }
+      parallel_for(
+        m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { b[cell_id] = fj[cell_id] * primal_Vj[cell_id]; });
 
       Vector<double> T{m_mesh->numberOfCells()};
       T = 0;
@@ -216,7 +300,8 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
 
       Vector<double> error{m_mesh->numberOfCells()};
       parallel_for(
-        m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { error[cell_id] = Temperature[cell_id] - Tj[cell_id]; });
+        m_mesh->numberOfCells(),
+        PUGS_LAMBDA(CellId cell_id) { error[cell_id] = (Temperature[cell_id] - Tj[cell_id]) * primal_Vj[cell_id]; });
 
       std::cout << "||Error||_2 = " << std::sqrt((error, error)) << "\n";
 
@@ -228,7 +313,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
       }
     }
   } else {
-    throw NotImplementedError("not done in 3d");
+    throw NotImplementedError("not done in 1d");
   }
 }
 
@@ -236,16 +321,19 @@ template HeatDiamondScheme<1>::HeatDiamondScheme(
   std::shared_ptr<const IMesh>,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
   const FunctionSymbolId&,
+  const FunctionSymbolId&,
   const FunctionSymbolId&);
 
 template HeatDiamondScheme<2>::HeatDiamondScheme(
   std::shared_ptr<const IMesh>,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
   const FunctionSymbolId&,
+  const FunctionSymbolId&,
   const FunctionSymbolId&);
 
 template HeatDiamondScheme<3>::HeatDiamondScheme(
   std::shared_ptr<const IMesh>,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
   const FunctionSymbolId&,
+  const FunctionSymbolId&,
   const FunctionSymbolId&);
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.hpp b/src/language/algorithms/HeatDiamondAlgorithm.hpp
index 9516b8277..b17538951 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.hpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.hpp
@@ -14,7 +14,8 @@ struct HeatDiamondScheme
   HeatDiamondScheme(std::shared_ptr<const IMesh> i_mesh,
                     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
                     const FunctionSymbolId& T_id,
-                    const FunctionSymbolId& kappa_id);
+                    const FunctionSymbolId& kappa_id,
+                    const FunctionSymbolId& f_id);
 };
 
 #endif   // HEAT_DIAMOND_ALGORITHM_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index af11321e2..64c8fa03b 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -176,23 +176,24 @@ SchemeModule::SchemeModule()
   this->_addBuiltinFunction("heat", std::make_shared<BuiltinFunctionEmbedder<
                                       void(std::shared_ptr<const IMesh>,
                                            const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-                                           const FunctionSymbolId&, const FunctionSymbolId&)>>(
+                                           const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
 
                                       [](std::shared_ptr<const IMesh> p_mesh,
                                          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                            bc_descriptor_list,
-                                         const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id) -> void {
+                                         const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
+                                         const FunctionSymbolId& f_id) -> void {
                                         switch (p_mesh->dimension()) {
                                         case 1: {
-                                          HeatDiamondScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id};
+                                          HeatDiamondScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                           break;
                                         }
                                         case 2: {
-                                          HeatDiamondScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id};
+                                          HeatDiamondScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                           break;
                                         }
                                         case 3: {
-                                          HeatDiamondScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id};
+                                          HeatDiamondScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                           break;
                                         }
                                         default: {
@@ -207,23 +208,24 @@ SchemeModule::SchemeModule()
                             std::make_shared<BuiltinFunctionEmbedder<
                               void(std::shared_ptr<const IMesh>,
                                    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-                                   const FunctionSymbolId&, const FunctionSymbolId&)>>(
+                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId)>>(
 
                               [](std::shared_ptr<const IMesh> p_mesh,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list,
-                                 const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id) -> void {
+                                 const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
+                                 const FunctionSymbolId& f_id) -> void {
                                 switch (p_mesh->dimension()) {
                                 case 1: {
-                                  Heat5PointsAlgorithm<1>{p_mesh, bc_descriptor_list, T_id, kappa_id};
+                                  Heat5PointsAlgorithm<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                   break;
                                 }
                                 case 2: {
-                                  Heat5PointsAlgorithm<2>{p_mesh, bc_descriptor_list, T_id, kappa_id};
+                                  Heat5PointsAlgorithm<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                   break;
                                 }
                                 case 3: {
-                                  Heat5PointsAlgorithm<3>{p_mesh, bc_descriptor_list, T_id, kappa_id};
+                                  Heat5PointsAlgorithm<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id};
                                   break;
                                 }
                                 default: {
-- 
GitLab