From: Timo Heister Date: Tue, 22 Feb 2011 10:32:28 +0000 (+0000) Subject: step-32: calculate and output l2 divergence error for testing X-Git-Tag: v8.0.0~4291 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1a11757276f53df9570d77ecd1acae4a770a9dc5;p=dealii.git step-32: calculate and output l2 divergence error for testing git-svn-id: https://svn.dealii.org/trunk@23411 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index 228afda708..f05dbd6b69 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -3207,6 +3207,57 @@ void BoussinesqFlowProblem::output_results () { computing_timer.enter_section ("Postprocessing"); + //calculate l2 divergence error + { + double my_cells_error = 0.0; + QGauss<1> q_base(stokes_degree+1); + QIterated err_quadrature(q_base, 2); + + const unsigned int n_q_points = err_quadrature.size(); + FEValues fe_values (stokes_fe, err_quadrature, update_JxW_values | update_gradients); + const unsigned int dofs_per_cell = fe_values.get_fe().dofs_per_cell; + const FEValuesExtractors::Vector velocities (0); + + std::vector local_dof_indices (fe_values.dofs_per_cell); + + + typename DoFHandler::active_cell_iterator + cell = stokes_dof_handler.begin_active(), + endc = stokes_dof_handler.end(); + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == + Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)) + { + fe_values.reinit (cell); + cell->get_dof_indices(local_dof_indices); + + double cell_error = 0.0; + for (unsigned int q = 0; q < n_q_points; ++q) + { + double tmp = 0; + for (unsigned int i=0; i < dofs_per_cell; ++i) + { + tmp += + fe_values[velocities].divergence(i, q) + * stokes_solution(local_dof_indices[i]); + } + cell_error += tmp * tmp * fe_values.JxW(q); + } + my_cells_error += cell_error; // = sqrt(cell_error)^2 + } + + double div_error = 0.0; +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + MPI_Allreduce (&my_cells_error, &div_error, 1, MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD); +#else + div_error = my_cells_error; +#endif + + div_error = std::sqrt(div_error); + pcout << "> divergence error:" << div_error << std::endl; + } + const FESystem joint_fe (stokes_fe, 1, temperature_fe, 1);