Add transfer function

src/scheme/ImplicitAcousticSolver.hpp

@@ -85,6 +85,64 @@ class ImplicitAcousticSolver
       });
   }
 
+  PUGS_INLINE
+  const int
+  Transferfunction(const int& k, const CellId& j)
+  {
+    return (k * m_mesh.numberOfCells()) + j;
+  }
+
+  // PUGS_INLINE
+  // const SparseMatrixDescriptor<int, uint8_t>
+  // ComputeMatrix(const CellValue<const double>& rhocj,
+  //               const NodeValue<const double>& rhocr,
+  //               const double& dt,
+  //               const CellValue<const double>& inv_mj)
+  // {
+  //   SparseMatrixDescriptor<int, uint8_t> S{2 * m_mesh.numberOfCells()};
+  //   for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
+  //     int i                  = 0;
+  //     const auto& cell_nodes = cell_to_node_matrix[j];
+  //     const size_t nb_nodes  = cell_nodes.size();
+  //     S(i, i)                = 1;
+  //     for (size_t R = 0; R < nb_nodes; ++R) {
+  //       NodeId r = cell_nodes[R];
+  //       S(i, i) +=S( T(0,j),T(0,j))
+  //     }
+  //     i += 1;
+  //   }
+  //   for (size_t i = 0; i < m_mesh.numberOfCells() - 1; ++i) {
+  //     const CellId un = 1;
+  //     S(i, i)         = 1 + (rhocj[i] * rhocj[i] * dt * inv_mj[i] * 0.5) * (rhocr[i] + rhocr[i + un]);
+  //     S(i, i + un)    = -rhocj[i] * rhocj[i] * dt * inv_mj[i] * 0.5 * rhocr[i + 1];
+  //     S(i + 1, i)     = -rhocj[i + 1] * rhocj[i + 1] * dt * inv_mj[i + 1] * 0.5 * rhocr[i];
+  //   }
+  //   S(m_mesh.numberOfCells() - 1, m_mesh.numberOfCells() - 1) =
+  //     1 + (rhocj[m_mesh.numberOfCells() - 1] * rhocj[m_mesh.numberOfCells() - 1] * dt *
+  //          inv_mj[m_mesh.numberOfCells() - 1] * 0.5) *
+  //           (rhocr[m_mesh.numberOfCells() - 1] + rhocr[m_mesh.numberOfCells()]);
+
+  // for (size_t i = m_mesh.numberOfCells(); i < 2 * m_mesh.numberOfCells() - 1; ++i) {
+  //   S(i, i)     = 1 + (dt * inv_mj[i] * 0.5) * (rhocr[i + 1] + rhocr[i]);
+  //   S(i, i + 1) = -dt * inv_mj[i] * 0.5 * rhocr[i];
+  //   S(i + 1, i) = -dt * inv_mj[i + 1] * 0.5 * rhocr[i];
+  // }
+  // S(2 * m_mesh.numberOfCells() - 1, 2 * m_mesh.numberOfCells() - 1) =
+  //   1 + (dt * inv_mj[m_mesh.numberOfCells() - 1] * 0.5) *
+  //         (rhocr[m_mesh.numberOfCells() - 1] + rhocr[m_mesh.numberOfCells()]);
+
+  // for (size_t i = 0; i < m_mesh.numberOfCells() - 1; ++i) {
+  //   S(i, i + m_mesh.numberOfCells() + 1) = dt * inv_mj[i] * 0.5 * rhocj[i];
+  //   S(i + 1, m_mesh.numberOfCells() + i) = -dt * inv_mj[i + 1] * 0.5 * rhocj[i + 1];
+
+  //   for (size_t i = 0; i < m_mesh.numberOfCells() - 1; ++i) {
+  //     S(i + m_mesh.numberOfCells(), i + 1) = dt * 0.5 * inv_mj[m_mesh.numberOfCells() + i];
+  //     S(i + m_mesh.numberOfCells() + 1, i) = -dt * 0.5 * inv_mj[m_mesh.numberOfCells() + i + 1];
+  //   }
+
+  //   return S;
+  // }
+
   PUGS_INLINE
   const NodeValue<const Rdd>
   computeAr(const NodeValuePerCell<const Rdd>& Ajr)
@@ -320,6 +378,9 @@ class ImplicitAcousticSolver
     for (NodeId r{0}; r < m_mesh.numberOfNodes(); ++r) {
       std::cout << "rhocr[" << r << "]=" << rhocr[r] << '\n';
     }
+    CellId x = 2;
+
+    std::cout << "transferfunction" << Transferfunction(0, x) << '\n';
 
     std::cout << "travailler\n";
     std::exit(0);