solver works for totally implicit or totally explicit mesh

src/scheme/ImplicitAcousticSolver.cpp

@@ -166,6 +166,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
             Ajr(j, r) = tensorProduct(rho[j] * c[j] * Cjr(j, r), njr(j, r));
           }
         });
+
       break;
     }
     }
@@ -292,7 +293,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
                   const Rd coef = m_Ajr(id_cell_i, node_nb_in_i) * m_inv_Ar[node_id] * Cjr(id_cell_j, node_nb_in_j);
                   A(j_index_p, i_index_u) += coef[i_dimension];
                   A(i_index_u, j_index_p) -= coef[i_dimension];
-                  // il y a un pb je passe deux fois au noeud du bord
                 }
               }
             }
@@ -480,8 +480,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
           }
         }
       }
-    }
-  }
       for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
         if (m_is_implicit_cell[cell_id]) {
           for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
@@ -493,6 +491,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
           }
         }
       }
+    }
 
     for (const auto& boundary_condition : m_boundary_condition_list) {
       std::visit(
@@ -906,7 +905,9 @@ ImplicitAcousticSolver(const SolverType solver_type,
                          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
                          const std::vector<std::shared_ptr<const IZoneDescriptor>>& explicit_zone_list,
                          const double& dt)
-  : m_solver_type{solver_type}, m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)}, m_mesh{*p_mesh}
+    : m_solver_type{solver_type},
+      m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)},
+      m_mesh{*p_mesh}
   {
     const auto node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
     const auto cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
@@ -1054,7 +1055,7 @@ ImplicitAcousticSolver(const SolverType solver_type,
 
       std::cout << "ratio" << norm_inf_sol / norm_inf_Un << "\n";
 
-  } while ((norm_inf_sol > 1e-14 * norm_inf_Un) and (nb_iter < 10000));
+    } while ((norm_inf_sol > 1e-14 * norm_inf_Un) and (nb_iter < 10000) and (m_number_of_implicit_cells > 0));
 
     m_predicted_u = [&] {
       CellValue<Rd> predicted_u = copy(u.cellValues());
@@ -1195,6 +1196,7 @@ apply(const double& dt,
                   NodeId node_id    = node_list[i_node];
                   const auto& value = value_list[i_node];
                   b[node_id]        = value;
+                  // std::cout << "node" << i_node << "value" << b[node_id] << '\n';
                 });
             } else {
               throw UnexpectedError("boundary condition not handled");
@@ -1204,7 +1206,8 @@ apply(const double& dt,
       }
 
       // NodeValue<const Rdxd> Ar = _computeAr(m_Ajr);
-
+      // std::cout << "inverse Ar" << m_inv_Ar << '\n';
+      // std::cout << "br" << b << '\n';
       const NodeValue<Rd> computed_ur(m_mesh.connectivity());
       parallel_for(
         m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { computed_ur[r] = m_inv_Ar[r] * b[r]; });
@@ -1212,6 +1215,7 @@ apply(const double& dt,
       return computed_ur;
     }();
 
+    // std::cout << "ur" << ur << '\n';
     // std::cout << "flux vitesse" << ur << '\n';
 
     //  p_j+1/2
@@ -1234,10 +1238,12 @@ apply(const double& dt,
     // for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
     //   const size_t& nb_nodes = m_Ajr.numberOfSubValues(cell_id);
     //   for (size_t r = 0; r < nb_nodes; ++r) {
+    //     std::cout << "Ajr(" << cell_id << "," << r << ")=" << m_Ajr(cell_id, r) << '\n';
     //     std::cout << "Fjr(" << cell_id << "," << r << ")=" << Fjr(cell_id, r) << '\n';
     //   }
     // }
 
+    // std::exit(0);
     // time n+1
     const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
 
@@ -1334,8 +1340,7 @@ ImplicitAcousticSolver(const SolverType solver_type,
   ImplicitAcousticSolver()                         = default;
   ImplicitAcousticSolver(ImplicitAcousticSolver&&) = default;
   ~ImplicitAcousticSolver()                        = default;
-}
-;
+};
 
 template <size_t Dimension>
 class ImplicitAcousticSolverHandler::ImplicitAcousticSolver<Dimension>::PressureBoundaryCondition