Add the calculation of the error

src/language/algorithms/ElasticityDiamondAlgorithm.cpp

@@ -628,16 +628,10 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
       Vector<double> U{number_of_dof * Dimension};
       U = 0;
 
-      BiCGStab{b, A, U, 1000, 1e-15};
+      BiCGStab{b, A, U, 10000, 1e-15};
 
       Vector r = A * U - b;
 
-      // for (CellId j = 0; j < number_of_dof; ++j) {
-      //   for (size_t l = 0; l < Dimension; ++l) {
-      //     std::cout << "Ucalc : " << U[(cell_dof_number[j] * Dimension) + l] << "\n";
-      //   }
-      // }
-
       std::cout << "real residu = " << std::sqrt((r, r)) << '\n';
 
       CellValue<TinyVector<Dimension>> Speed{mesh->connectivity()};
@@ -653,32 +647,16 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
           Uexacte[(cell_dof_number[j] * Dimension) + l] = Uj[j][l];
         }
       }
-      // std::cout << "A*Uex - b =\n";
-      // for (CellId i_id = 0; i_id < number_of_dof; ++i_id) {
-      //   for (size_t j = 0; j < Dimension; ++j) {
-      //     std::cout << (A * Uexacte - b)[(cell_dof_number[i_id] * Dimension) + j] << "\n";
-      //   }
-      // }
-
-      // std::cout << "Vecteur AU - b: \n";
-      // for (CellId i = 0; i < number_of_dof; ++i) {
-      //   // double temp = 0;
-      //   for (size_t k = 0; k < Dimension; ++k) {
-      //     double temp = 0;
-      //     for (CellId j = 0; j < number_of_dof; ++j) {
-      //       for (size_t l = 0; l < Dimension; ++l) {
-      //         temp += S((cell_dof_number[i] * Dimension) + k, (cell_dof_number[j] * Dimension) + l) * Uj[j][l];
-      //       }
-      //     }
-      //     std::cout << temp << "; " << b[(cell_dof_number[i] * Dimension) + k] << "=> "
-      //               << temp - b[(cell_dof_number[i] * Dimension) + k] << "\n";
-      //   }
-      // }
-
-      // std::cout << "Vecteur Uexacte : \n";
-      // for (CellId i = 0; i < number_of_dof; ++i) {
-      //   std::cout << Uj[i] << "\n";
-      // }
+
+      Vector<double> error{number_of_dof * Dimension};
+      for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+        for (size_t i = 0; i < Dimension; ++i) {
+          error[(cell_id * Dimension) + i] = (Speed[cell_id][i] - Uj[cell_id][i]) * sqrt(primal_Vj[cell_id]);
+        }
+      }
+
+      std::cout << "||Error||_2 = " << std::sqrt((error, error)) << "\n";
+
       VTKWriter vtk_writer("Speed_" + std::to_string(Dimension), 0.01);
 
       vtk_writer.write(mesh, {NamedItemValue{"U", Speed}, NamedItemValue{"Uexacte", Uj}}, 0,