Apply boundary conditions on matrix Ar, put inv_Ar as a global variable

src/scheme/ImplicitAcousticSolver.cpp

@@ -43,6 +43,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
   const MeshType& m_mesh;
 
   NodeValuePerCell<const Rdxd> m_Ajr;
+  NodeValue<const Rdxd> m_inv_Ar;
   NodeValue<const bool> m_is_boundary_node;
 
   CellValue<const Rd> m_predicted_u;
@@ -138,9 +139,55 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
         Ar[r] = sum;
       });
 
+    for (const auto& boundary_condition : m_boundary_condition_list) {
+      std::visit(
+        [&](auto&& bc) {
+          using T = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
+            const auto& node_list = bc.nodeList();
+
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(size_t i_node) {
+                NodeId node_id = node_list[i_node];
+
+                Ar[node_id] = identity;
+              });
+          } else if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+            const Rd& n = bc.outgoingNormal();
+
+            const Rdxd I   = identity;
+            const Rdxd nxn = tensorProduct(n, n);
+            const Rdxd P   = I - nxn;
+
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
+                const NodeId r = node_list[r_number];
+
+                Ar[r] = P * Ar[r] * P + nxn;
+              });
+          } else {
+            throw UnexpectedError("boundary condition not handled");
+          }
+        },
+        boundary_condition);
+    }
+
     return Ar;
   }
 
+  NodeValue<const Rdxd>
+  _computeInvAr(const NodeValue<const Rdxd>& Ar) const
+  {
+    NodeValue<Rdxd> inv_Ar{m_mesh.connectivity()};
+
+    for (NodeId node_id = 0; node_id < m_mesh.numberOfNodes(); ++node_id) {
+      inv_Ar[node_id] = inverse(Ar[node_id]);
+    }
+
+    return inv_Ar;
+  }
+
   int
   mapP(CellId j) const
   {
@@ -156,7 +203,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
   }
 
   CRSMatrixDescriptor<double>
-  _getA(const NodeValuePerCell<const Rdxd>& Ajr, const NodeValue<const Rdxd>& Ar) const
+  _getA(const NodeValuePerCell<const Rdxd>& Ajr) const
   {
     Array<int> non_zeros{(Dimension + 1) * m_number_of_implicit_cells};
     non_zeros.fill(3);
@@ -169,7 +216,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
 
     for (NodeId node_id = 0; node_id < m_mesh.numberOfNodes(); ++node_id) {
       if (not m_is_boundary_node[node_id]) {
-        const Rdxd& inv_Ar                         = inverse(Ar[node_id]);
         const auto& node_cell                      = node_to_cell_matrix[node_id];
         const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemValues(node_id);
 
@@ -182,7 +228,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
             const size_t node_nb_in_j = node_local_number_in_its_cells[cell_j];
             const int j_index_p       = mapP(id_cell_j);
 
-            const double coef = (Ajr(id_cell_i, node_nb_in_i) * inv_Ar * Cjr(id_cell_j, node_nb_in_j))[0];
+            const double coef = (Ajr(id_cell_i, node_nb_in_i) * m_inv_Ar[node_id] * Cjr(id_cell_j, node_nb_in_j))[0];
             A(j_index_p, i_index_u) += coef;
             A(i_index_u, j_index_p) -= coef;
           }
@@ -763,10 +809,23 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
 
     m_Ajr                                = this->_computeAjr(rho, c);
     NodeValue<const Rdxd> Ar             = this->_computeAr(m_Ajr);
-    const CRSMatrixDescriptor<double>& A = this->_getA(m_Ajr, Ar);
+    m_inv_Ar                             = this->_computeInvAr(Ar);
+    const CRSMatrixDescriptor<double>& A = this->_getA(m_Ajr);
     Vector<double> Un                    = this->_getU(p, u);
     Vector<double> Uk                    = copy(Un);
 
+    std::cout << "Ajr" << '\n';
+    for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
+      const size_t& nb_nodes = m_Ajr.numberOfSubValues(j);
+      for (size_t r = 0; r < nb_nodes; ++r) {
+        std::cout << m_Ajr(j, r) << '\n';
+      }
+    }
+    // std::cout << "Ar" << '\n';
+    // std::cout << Ar << '\n';
+    // std::cout << m_inv_Ar << '\n';
+    // std::exit(0);
+
     int nb_iter = 0;
     double norm_inf_sol;