From f2c3a45b78a067ffb004aeaf0314ebe67f75d069 Mon Sep 17 00:00:00 2001
From: Drouard <axelle.drouard@orange.fr>
Date: Tue, 7 Jun 2022 14:06:41 +0200
Subject: [PATCH] Change arguments of "nodalheat" for implicit scheme.

---
 src/language/modules/SchemeModule.cpp |  34 ++++----
 src/scheme/ScalarNodalScheme.cpp      | 119 +++++++++-----------------
 src/scheme/ScalarNodalScheme.hpp      |   7 +-
 3 files changed, 62 insertions(+), 98 deletions(-)

diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 8800b37b4..abb60970a 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -603,23 +603,25 @@ SchemeModule::SchemeModule()
       }
 
       ));
-  this->_addBuiltinFunction("nodalheat",
-                            std::make_shared<BuiltinFunctionEmbedder<
-                              std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
-                                                                       const std::shared_ptr<const IDiscreteFunction>&,
-                                                                       const std::shared_ptr<const IDiscreteFunction>&,
-                                                                       const std::vector<std::shared_ptr<
-                                                                         const IBoundaryConditionDescriptor>>&)>>(
-
-                              [](const std::shared_ptr<const IDiscreteFunction>& mub_dual,
-                                 const std::shared_ptr<const IDiscreteFunction>& mu_dual,
-                                 const std::shared_ptr<const IDiscreteFunction>& f,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list) -> const std::shared_ptr<const IDiscreteFunction> {
-                                return ScalarNodalSchemeHandler{mub_dual, mu_dual, f, bc_descriptor_list}.solution();
-                              }
+  this
+    ->_addBuiltinFunction("nodalheat",
+                          std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                            const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+                                                     const std::shared_ptr<const IDiscreteFunction>&,
+                                                     const std::vector<
+                                                       std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+                                                     const std::shared_ptr<const IDiscreteFunction>&, const double&)>>(
+
+                            [](const std::shared_ptr<const IDiscreteFunction>& kappa,
+                               const std::shared_ptr<const IDiscreteFunction>& f,
+                               const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                 bc_descriptor_list,
+                               const std::shared_ptr<const IDiscreteFunction>& P,
+                               const double& dt) -> const std::shared_ptr<const IDiscreteFunction> {
+                              return ScalarNodalSchemeHandler{kappa, f, bc_descriptor_list, P, dt}.solution();
+                            }
 
-                              ));
+                            ));
 
   this->_addBuiltinFunction("unsteadyelasticity",
                             std::make_shared<BuiltinFunctionEmbedder<
diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index 2837089cd..72936fe92 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -143,10 +143,11 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
   }
 
   ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
-                    const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_b,
                     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)
+                    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+                    const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& P,
+                    const double& dt)
   {
     using BoundaryCondition =
       std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, NeumannBoundaryCondition>;
@@ -408,8 +409,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     }();
 
     {
-      CellValue<const TinyMatrix<Dimension>> cell_kappaj  = cell_k->cellValues();
-      CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_b->cellValues();
+      CellValue<const TinyMatrix<Dimension>> cell_kappaj = cell_k->cellValues();
 
       const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
         NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
@@ -429,24 +429,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         return kappa;
       }();
 
-      const NodeValue<const TinyMatrix<Dimension>> node_kapparb = [&] {
-        NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
-        parallel_for(
-          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 i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
-              const CellId cell_id = node_to_cell[i_cell];
-              double local_weight  = 1. / l2Norm(xr[node_id] - xj[cell_id]);
-              kappa[node_id] += local_weight * cell_kappajb[cell_id];
-              weight += local_weight;
-            }
-            kappa[node_id] = 1. / weight * kappa[node_id];
-          });
-        return kappa;
-      }();
-
       const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] {
         NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
         parallel_for(
@@ -474,17 +456,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         return kappa_invBeta;
       }();
 
-      const NodeValue<const TinyMatrix<Dimension>> kapparb_invBetar = [&] {
-        NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
-        parallel_for(
-          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
-            if (not node_is_corner[node_id]) {
-              kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(node_betar[node_id]);
-            }
-          });
-        return kappa_invBeta;
-      }();
-
       const NodeValue<const TinyMatrix<Dimension>> corner_betar = [&] {
         NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
         parallel_for(
@@ -538,17 +509,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         return kappa_invBeta;
       }();
 
-      const NodeValue<const TinyMatrix<Dimension>> corner_kapparb_invBetar = [&] {
-        NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
-        parallel_for(
-          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
-            if (node_is_corner[node_id]) {
-              kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(corner_betar[node_id]);
-            }
-          });
-        return kappa_invBeta;
-      }();
-
       const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] {
         NodeValue<TinyVector<Dimension>> compute_sum_Cjr{mesh->connectivity()};
         parallel_for(
@@ -570,8 +530,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         parallel_for(
           mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
             if (is_boundary_node[node_id]) {
-              compute_v[node_id] = 1. / l2Norm(node_kapparb[node_id] * exterior_normal[node_id]) *
-                                   node_kapparb[node_id] * exterior_normal[node_id];
+              compute_v[node_id] = 1. / l2Norm(node_kappar[node_id] * exterior_normal[node_id]) * node_kappar[node_id] *
+                                   exterior_normal[node_id];
             }
           });
         return compute_v;
@@ -634,7 +594,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
           const TinyMatrix<Dimension> I   = identity;
           const TinyVector<Dimension> n   = exterior_normal[node_id];
           const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
-          const TinyMatrix<Dimension> P   = I - nxn;
+          const TinyMatrix<Dimension> Q   = I - nxn;
 
           for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
             const CellId cell_id_j = node_to_cell[i_cell_j];
@@ -647,7 +607,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
               S(cell_id_j, cell_id_k) +=
                 dot(kappar_invBetar[node_id] *
                       (Cjr(cell_id_k, i_node_k) - theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
-                    P * Cjr(cell_id_k, i_node_j));
+                    Q * Cjr(cell_id_k, i_node_j));
             }
           }
         } else if ((node_is_dirichlet[node_id]) && (not node_is_corner[node_id])) {
@@ -679,18 +639,25 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         }
       }
 
+      for (CellId cell_id_j = 0; cell_id_j < mesh->numberOfCells(); ++cell_id_j) {
+        for (CellId cell_id_k = 0; cell_id_k < mesh->numberOfCells(); ++cell_id_k) {
+          S(cell_id_j, cell_id_k) *= dt / Vj[cell_id_j];
+        }
+      };
+
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-        S(cell_id, cell_id) += Vj[cell_id];
+        S(cell_id, cell_id) += 1;
       }
 
       CellValue<const double> fj = f->cellValues();
+      CellValue<const double> Ph = P->cellValues();
 
       CRSMatrix A{S.getCRSMatrix()};
 
       Vector<double> b{mesh->numberOfCells()};
       b = zero;
       for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-        b[cell_id] = fj[cell_id] * Vj[cell_id];
+        b[cell_id] = fj[cell_id];
       };
 
       for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
@@ -699,16 +666,16 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         const TinyMatrix<Dimension> I              = identity;
         const TinyVector<Dimension> n              = exterior_normal[node_id];
         const TinyMatrix<Dimension> nxn            = tensorProduct(n, n);
-        const TinyMatrix<Dimension> P              = I - nxn;
+        const TinyMatrix<Dimension> Q              = I - nxn;
         if ((node_is_neumann[node_id]) && (not node_is_dirichlet[node_id])) {
           if ((is_boundary_node[node_id]) and (not node_is_corner[node_id])) {
             for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
               const CellId cell_id = node_to_cell[i_cell];
               const size_t i_node  = node_local_number_in_its_cells[i_cell];
 
-              b[cell_id] += node_boundary_values[node_id] *
-                            (1. / (sum_theta[node_id] * l2Norm(node_kapparb[node_id] * n)) *
-                               dot(kapparb_invBetar[node_id] * sum_Cjr[node_id], P * Cjr(cell_id, i_node)) +
+              b[cell_id] += 1. / Vj[cell_id] * node_boundary_values[node_id] *
+                            (1. / (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
+                               dot(kappar_invBetar[node_id] * sum_Cjr[node_id], Q * Cjr(cell_id, i_node)) +
                              dot(Cjr(cell_id, i_node), n));
             }
           } else if (node_is_corner[node_id]) {
@@ -719,7 +686,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
               FaceId face_id = node_to_face[i_face];
               sum_mes_l += mes_l[face_id];
             }
-            b[cell_id] += 0.5 * node_boundary_values[node_id] * sum_mes_l;
+            b[cell_id] += 1. / Vj[cell_id] * 0.5 * node_boundary_values[node_id] * sum_mes_l;
           }
         } else if (node_is_dirichlet[node_id]) {
           if (not node_is_corner[node_id]) {
@@ -731,9 +698,9 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                 const CellId cell_id_k = node_to_cell[i_cell_k];
                 const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
 
-                b[cell_id_j] +=
-                  dot(node_boundary_values[node_id] * kapparb_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
-                      Cjr(cell_id_j, i_node_j));
+                b[cell_id_j] += 1. / Vj[cell_id_j] *
+                                dot(node_boundary_values[node_id] * kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+                                    Cjr(cell_id_j, i_node_j));
               }
             }
 
@@ -747,7 +714,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                 const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
 
                 b[cell_id_j] +=
-                  dot(node_boundary_values[node_id] * corner_kapparb_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+                  1. / Vj[cell_id_j] *
+                  dot(node_boundary_values[node_id] * corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
                       Cjr(cell_id_j, i_node_j));
               }
             }
@@ -757,8 +725,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         Vector<double> T{mesh->numberOfCells()};
         T = zero;
 
+        Vector<double> B{mesh->numberOfCells()};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { B[cell_id] = dt * b[cell_id] + Ph[cell_id]; });
+
         LinearSolver solver;
-        solver.solveLocalSystem(A, T, b);
+        solver.solveLocalSystem(A, T, B);
 
         m_solution     = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
         auto& solution = *m_solution;
@@ -776,17 +748,18 @@ ScalarNodalSchemeHandler::solution() const
 }
 
 ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
-  const std::shared_ptr<const IDiscreteFunction>& cell_k_b,
   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)
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const std::shared_ptr<const IDiscreteFunction>& P,
+  const double& dt)
 {
   const std::shared_ptr i_mesh = getCommonMesh({cell_k, f});
   if (not i_mesh) {
     throw NormalError("primal discrete functions are not defined on the same mesh");
   }
 
-  checkDiscretizationType({cell_k_b, cell_k, f}, DiscreteFunctionType::P0);
+  checkDiscretizationType({cell_k, f}, DiscreteFunctionType::P0);
 
   switch (i_mesh->dimension()) {
   case 1: {
@@ -801,26 +774,14 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
     m_scheme =
-      std::make_unique<ScalarNodalScheme<2>>(mesh,
-                                             std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_b),
-                                             std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
+      std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
                                              std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
-                                             bc_descriptor_list);
+                                             bc_descriptor_list,
+                                             std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(P), dt);
     break;
   }
   case 3: {
-    using MeshType                   = Mesh<Connectivity<3>>;
-    using DiscreteTensorFunctionType = DiscreteFunctionP0<3, TinyMatrix<3>>;
-    using DiscreteScalarFunctionType = DiscreteFunctionP0<3, double>;
-
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    m_scheme =
-      std::make_unique<ScalarNodalScheme<3>>(mesh,
-                                             std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_b),
-                                             std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
-                                             std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
-                                             bc_descriptor_list);
+    throw NormalError("The scheme is not implemented in dimension 3");
     break;
   }
   default: {
diff --git a/src/scheme/ScalarNodalScheme.hpp b/src/scheme/ScalarNodalScheme.hpp
index 27ec1818e..a43ff3b4c 100644
--- a/src/scheme/ScalarNodalScheme.hpp
+++ b/src/scheme/ScalarNodalScheme.hpp
@@ -40,10 +40,11 @@ class ScalarNodalSchemeHandler
 
   std::shared_ptr<const IDiscreteFunction> solution() const;
 
-  ScalarNodalSchemeHandler(const std::shared_ptr<const IDiscreteFunction>& cell_k_bound,
-                           const std::shared_ptr<const IDiscreteFunction>& cell_k,
+  ScalarNodalSchemeHandler(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);
+                           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+                           const std::shared_ptr<const IDiscreteFunction>& P,
+                           const double& dt);
 
   ~ScalarNodalSchemeHandler();
 };
-- 
GitLab