Replace scalar diffusion coefficient with tensorial one

src/scheme/CMakeLists.txt

@@ -9,4 +9,5 @@ add_library(
   DiscreteFunctionVectorIntegrator.cpp
   DiscreteFunctionVectorInterpoler.cpp
   ScalarDiamondScheme.cpp
-  VectorDiamondScheme.cpp)
+  VectorDiamondScheme.cpp
+  ScalarNodalScheme.cpp)

src/scheme/ScalarNodalScheme.cpp

@@ -3,137 +3,6 @@
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
 
-template <size_t Dimension>
-class ScalarNodalSchemeHandler::InterpolationWeightsManager
-{
- private:
-  std::shared_ptr<const Mesh<Connectivity<Dimension>>> m_mesh;
-  FaceValue<bool> m_primal_face_is_on_boundary;
-  NodeValue<bool> m_primal_node_is_on_boundary;
-  CellValuePerNode<double> m_w_rj;
-  FaceValuePerNode<double> m_w_rl;
-
- public:
-  CellValuePerNode<double>&
-  wrj()
-  {
-    return m_w_rj;
-  }
-
-  FaceValuePerNode<double>&
-  wrl()
-  {
-    return m_w_rl;
-  }
-
-  void
-  compute()
-  {
-    using MeshDataType      = MeshData<Dimension>;
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-    const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
-
-    const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
-    const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
-    const auto& node_to_cell_matrix                  = m_mesh->connectivity().nodeToCellMatrix();
-    const auto& node_to_face_matrix                  = m_mesh->connectivity().nodeToFaceMatrix();
-    CellValuePerNode<double> w_rj{m_mesh->connectivity()};
-    FaceValuePerNode<double> w_rl{m_mesh->connectivity()};
-
-    for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
-      w_rl[i] = std::numeric_limits<double>::signaling_NaN();
-    }
-
-    for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) {
-      SmallVector<double> b{Dimension + 1};
-      b[0] = 1;
-      for (size_t i = 1; i < Dimension + 1; i++) {
-        b[i] = xr[i_node][i - 1];
-      }
-      const auto& node_to_cell = node_to_cell_matrix[i_node];
-
-      if (not m_primal_node_is_on_boundary[i_node]) {
-        SmallMatrix<double> A{Dimension + 1, node_to_cell.size()};
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          A(0, j) = 1;
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          for (size_t j = 0; j < node_to_cell.size(); j++) {
-            const CellId J = node_to_cell[j];
-            A(i, j)        = xj[J][i - 1];
-          }
-        }
-        SmallVector<double> x{node_to_cell.size()};
-        x = zero;
-
-        LeastSquareSolver ls_solver;
-        ls_solver.solveLocalSystem(A, x, b);
-
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          w_rj(i_node, j) = x[j];
-        }
-      } else {
-        int nb_face_used = 0;
-        for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-          FaceId face_id = node_to_face_matrix[i_node][i_face];
-          if (m_primal_face_is_on_boundary[face_id]) {
-            nb_face_used++;
-          }
-        }
-        SmallMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used};
-        for (size_t j = 0; j < node_to_cell.size() + nb_face_used; j++) {
-          A(0, j) = 1;
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          for (size_t j = 0; j < node_to_cell.size(); j++) {
-            const CellId J = node_to_cell[j];
-            A(i, j)        = xj[J][i - 1];
-          }
-        }
-        for (size_t i = 1; i < Dimension + 1; i++) {
-          int cpt_face = 0;
-          for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-            FaceId face_id = node_to_face_matrix[i_node][i_face];
-            if (m_primal_face_is_on_boundary[face_id]) {
-              A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
-              cpt_face++;
-            }
-          }
-        }
-
-        SmallVector<double> x{node_to_cell.size() + nb_face_used};
-        x = zero;
-
-        LeastSquareSolver ls_solver;
-        ls_solver.solveLocalSystem(A, x, b);
-
-        for (size_t j = 0; j < node_to_cell.size(); j++) {
-          w_rj(i_node, j) = x[j];
-        }
-        int cpt_face = node_to_cell.size();
-        for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
-          FaceId face_id = node_to_face_matrix[i_node][i_face];
-          if (m_primal_face_is_on_boundary[face_id]) {
-            w_rl(i_node, i_face) = x[cpt_face++];
-          }
-        }
-      }
-    }
-    m_w_rj = w_rj;
-    m_w_rl = w_rl;
-  }
-
-  InterpolationWeightsManager(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh,
-                              FaceValue<bool> primal_face_is_on_boundary,
-                              NodeValue<bool> primal_node_is_on_boundary)
-    : m_mesh(mesh),
-      m_primal_face_is_on_boundary(primal_face_is_on_boundary),
-      m_primal_node_is_on_boundary(primal_node_is_on_boundary)
-  {}
-
-  ~InterpolationWeightsManager() = default;
-};
 
 class ScalarNodalSchemeHandler::IScalarNodalScheme
 {
@@ -249,23 +118,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     ~FourierBoundaryCondition() = default;
   };
 
-  class SymmetryBoundaryCondition
-  {
-   private:
-    const Array<const FaceId> m_face_list;
-
-   public:
-    const Array<const FaceId>&
-    faceList() const
-    {
-      return m_face_list;
-    }
-
-   public:
-    SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
-
-    ~SymmetryBoundaryCondition() = default;
-  };
 
  public:
   std::shared_ptr<const IDiscreteFunction>
@@ -276,18 +128,14 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
 
   ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
                       const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& alpha,
-                      const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mub,
-                      const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mu,
+                      const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& nod_k_bound,
+                      const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& nod_k,
                       const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
                       const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
   {
-    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mu->mesh(),
-           "diffusion coefficient is not defined on the dual mesh!");
-    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mub->mesh(),
-           "boundary diffusion coefficient is not defined on the dual mesh!");
 
     using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition,
-                                           NeumannBoundaryCondition, SymmetryBoundaryCondition>;
+                                           NeumannBoundaryCondition>;
 
     using BoundaryConditionList = std::vector<BoundaryCondition>;
 
@@ -297,10 +145,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       bool is_valid_boundary_condition = true;
 
       switch (bc_descriptor->type()) {
-      case IBoundaryConditionDescriptor::Type::symmetry: {
-        throw NotImplementedError("NIY");
-        break;
-      }
       case IBoundaryConditionDescriptor::Type::dirichlet: {
         const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
           dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
@@ -381,7 +225,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
               using T = std::decay_t<decltype(bc)>;
               if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
                             (std::is_same_v<T, FourierBoundaryCondition>) or
-                            (std::is_same_v<T, SymmetryBoundaryCondition>) or
                             (std::is_same_v<T, DirichletBoundaryCondition>)) {
                 const auto& face_list = bc.faceList();
 
@@ -411,8 +254,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             [&](auto&& bc) {
               using T = std::decay_t<decltype(bc)>;
               if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
-                            (std::is_same_v<T, FourierBoundaryCondition>) or
-                            (std::is_same_v<T, SymmetryBoundaryCondition>)) {
+                            (std::is_same_v<T, FourierBoundaryCondition>)) {
                 const auto& face_list = bc.faceList();
 
                 for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
@@ -466,10 +308,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]));
       }
 
-      InterpolationWeightsManager iwm(mesh, primal_face_is_on_boundary, primal_node_is_on_boundary);
-      iwm.compute();
-      CellValuePerNode<double> w_rj = iwm.wrj();
-      FaceValuePerNode<double> w_rl = iwm.wrl();
 
       MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
 
@@ -485,8 +323,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         std::shared_ptr mapper =
           DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity());
 
-        CellValue<const double> dual_kappaj  = dual_mu->cellValues();
-        CellValue<const double> dual_kappajb = dual_mub->cellValues();
+        CellValue<const TinyMatrix<Dimension>> nod_kappar  = nod_k->cellValues();
+        CellValue<const TinyMatrix<Dimension>> nod_kapparb = nod_k_bound->cellValues();
 
         const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
 
@@ -588,29 +426,25 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         }();
 
         FaceValue<const double> alpha_l = [&] {
-          CellValue<double> alpha_j{diamond_mesh->connectivity()};
+          FaceValue<double> alpha_j{mesh->connectivity()};
 
           parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_j[diamond_cell_id] = dual_kappaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
+            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+              alpha_j[face_id] = 1;
             });
 
-          FaceValue<double> computed_alpha_l{mesh->connectivity()};
-          mapper->fromDualCell(alpha_j, computed_alpha_l);
-          return computed_alpha_l;
+          return alpha_j;
         }();
 
         FaceValue<const double> alphab_l = [&] {
-          CellValue<double> alpha_jb{diamond_mesh->connectivity()};
+          FaceValue<double> alpha_lb{mesh->connectivity()};
 
           parallel_for(
-            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
-              alpha_jb[diamond_cell_id] = dual_kappajb[diamond_cell_id] / dual_Vj[diamond_cell_id];
+            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+              alpha_lb[face_id] = 1; //Refaire
             });
 
-          FaceValue<double> computed_alpha_lb{mesh->connectivity()};
-          mapper->fromDualCell(alpha_jb, computed_alpha_lb);
-          return computed_alpha_lb;
+          return alpha_lb;
         }();
 
         const CellValue<const double> primal_Vj = mesh_data.Vj();
@@ -644,61 +478,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
           S(j, j) += (*alpha)[cell_id] * primal_Vj[cell_id];
         }
 
-        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_face_id  = mes_l[face_id] * alpha_l[face_id];
-          const double alphab_face_id = mes_l[face_id] * alphab_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 = (dot(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];
-
-              const TinyVector<Dimension> nil = [&] {
-                if (is_face_reversed_for_cell_i) {
-                  return -nlj[face_id];
-                } else {
-                  return nlj[face_id];
-                }
-              }();
-
-              CellId dual_cell_id = face_dual_cell_id[face_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 dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                if (dual_node_primal_node_id[dual_node_id] == node_id) {
-                  const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
-
-                  const double a_ir  = alpha_face_id * dot(nil, Clr);
-                  const double ab_ir = alphab_face_id * dot(nil, Clr);
-
-                  for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
-                    CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
-                    S(cell_dof_number[i_id], cell_dof_number[j_id]) -= w_rj(node_id, j_cell) * a_ir;
-                    if (primal_face_is_neumann[face_id]) {
-                      S(face_dof_number[face_id], cell_dof_number[j_id]) += w_rj(node_id, j_cell) * ab_ir;
-                    }
-                  }
-                  if (primal_node_is_on_boundary[node_id]) {
-                    for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
-                      FaceId l_id = node_to_face_matrix[node_id][l_face];
-                      if (primal_face_is_on_boundary[l_id]) {
-                        S(cell_dof_number[i_id], face_dof_number[l_id]) -= w_rl(node_id, l_face) * a_ir;
-                        if (primal_face_is_neumann[face_id]) {
-                          S(face_dof_number[face_id], face_dof_number[l_id]) += w_rl(node_id, l_face) * ab_ir;
-                        }
-                      }
-                    }
-                  }
-                }
-              }
-            }
-          }
-        }
         for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
           if (primal_face_is_dirichlet[face_id]) {
             S(face_dof_number[face_id], face_dof_number[face_id]) += 1;
@@ -772,8 +552,8 @@ ScalarNodalSchemeHandler::solution() const
 
 ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
   const std::shared_ptr<const IDiscreteFunction>& alpha,
-  const std::shared_ptr<const IDiscreteFunction>& dual_mub,
-  const std::shared_ptr<const IDiscreteFunction>& dual_mu,
+  const std::shared_ptr<const IDiscreteFunction>& nod_k_bound,
+  const std::shared_ptr<const IDiscreteFunction>& nod_k,
   const std::shared_ptr<const IDiscreteFunction>& f,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
 {
@@ -781,15 +561,16 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
   if (not i_mesh) {
     throw NormalError("primal discrete functions are not defined on the same mesh");
   }
-  const std::shared_ptr i_dual_mesh = getCommonMesh({dual_mub, dual_mu});
+  const std::shared_ptr i_dual_mesh = getCommonMesh({nod_k_bound, nod_k});
   if (not i_dual_mesh) {
     throw NormalError("dual discrete functions are not defined on the same mesh");
   }
-  checkDiscretizationType({alpha, dual_mub, dual_mu, f}, DiscreteFunctionType::P0);
+  checkDiscretizationType({alpha, nod_k_bound, nod_k, f}, DiscreteFunctionType::P0);
 
   switch (i_mesh->dimension()) {
   case 1: {
     using MeshType                   = Mesh<Connectivity<1>>;
+    using DiscreteTensorFunctionType = DiscreteFunctionP0<1, TinyMatrix<1>>;
     using DiscreteScalarFunctionType = DiscreteFunctionP0<1, double>;
 
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
@@ -800,14 +581,15 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
 
     m_scheme =
       std::make_unique<ScalarNodalScheme<1>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
-                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
-                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
+                                               std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k_bound),
+                                               std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
                                                bc_descriptor_list);
     break;
   }
   case 2: {
     using MeshType                   = Mesh<Connectivity<2>>;
+    using DiscreteTensorFunctionType = DiscreteFunctionP0<2, TinyMatrix<2>>;
     using DiscreteScalarFunctionType = DiscreteFunctionP0<2, double>;
 
     std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
@@ -818,14 +600,15 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
 
     m_scheme =
       std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
-                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
-                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
+                                               std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k_bound),
+                                               std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
                                                bc_descriptor_list);
     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);
@@ -836,8 +619,8 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
 
     m_scheme =
       std::make_unique<ScalarNodalScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
-                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
-                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
+                                               std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k_bound),
+                                               std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k),
                                                std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
                                                bc_descriptor_list);
     break;