diff --git a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp index e7e8de72edcae22d73ba20f7aadd280801b4ca18..e5c7d3c2786d874e00fa3f4a4247cce725fc693e 100644 --- a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp +++ b/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,