From fab23cde9fcd605ba123d57b1f994254b8fc57a2 Mon Sep 17 00:00:00 2001 From: Scott Miller Date: Tue, 13 Aug 2013 02:50:00 +0000 Subject: [PATCH] step-51 use workstream for postprocessing the local solution git-svn-id: https://svn.dealii.org/trunk@30297 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-51/step-51.cc | 238 +++++++++++++++++++--------- 1 file changed, 165 insertions(+), 73 deletions(-) diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc index 662f3f25ad..bc4778c23a 100644 --- a/deal.II/examples/step-51/step-51.cc +++ b/deal.II/examples/step-51/step-51.cc @@ -301,6 +301,10 @@ private: struct PerTaskData; struct ScratchData; + + struct EmptyData; + struct PostProcessScratchData; + void NullFunction(const EmptyData &data); void setup_system (); void assemble_system (const bool reconstruct_trace = false); @@ -310,6 +314,9 @@ private: void copy_local_to_global(const PerTaskData &data); void solve (); void postprocess (); + void postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, + PostProcessScratchData &scratch, + EmptyData &task_data); void refine_grid (const unsigned int cylce); void output_results (const unsigned int cycle); @@ -508,6 +515,69 @@ struct Step51::ScratchData }; +template +struct Step51::EmptyData +{ + EmptyData(){} + void reset(){} +}; + +template +struct Step51::PostProcessScratchData +{ + FEValues fe_values_local; + FEValues fe_values; + + std::vector u_values; + std::vector > u_gradients; + FullMatrix cell_matrix; + + Vector cell_rhs; + Vector cell_sol; + + // Full constructor + PostProcessScratchData(const FiniteElement &fe, + const FiniteElement &fe_local, + const QGauss &quadrature_formula, + const UpdateFlags local_flags, + const UpdateFlags flags) + : + fe_values_local (fe_local, quadrature_formula, local_flags), + fe_values (fe, quadrature_formula, flags), + u_values (quadrature_formula.size()), + u_gradients (quadrature_formula.size()), + cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell), + cell_rhs (fe.dofs_per_cell), + cell_sol (fe.dofs_per_cell) + {} + + // Copy constructor + PostProcessScratchData(const PostProcessScratchData &sd) + : + fe_values_local (sd.fe_values_local.get_fe(), + sd.fe_values_local.get_quadrature(), + sd.fe_values_local.get_update_flags()), + fe_values (sd.fe_values.get_fe(), + sd.fe_values.get_quadrature(), + sd.fe_values.get_update_flags()), + u_values (sd.u_values), + u_gradients (sd.u_gradients), + cell_matrix (sd.cell_matrix), + cell_rhs (sd.cell_rhs), + cell_sol (sd.cell_sol) + {} + + void reset() { + cell_matrix = 0.; + cell_rhs = 0.; + cell_sol = 0.; + } + +}; + +template +void Step51::NullFunction(const EmptyData &data){} + template void Step51::copy_local_to_global(const PerTaskData &data) { @@ -796,7 +866,31 @@ template void Step51::postprocess() { - const unsigned int n_active_cells=triangulation.n_active_cells(); + // construct post-processed solution with (hopefully) higher order of + // accuracy + { + const QGauss quadrature_formula(fe_u_post.degree+1); + const UpdateFlags local_flags (update_values); + const UpdateFlags flags ( update_values | update_gradients | + update_JxW_values); + + EmptyData task_data; + + PostProcessScratchData scratch (fe_u_post, fe_local, + quadrature_formula, + local_flags, + flags); + + WorkStream::run(dof_handler_u_post.begin_active(), + dof_handler_u_post.end(), + *this, + &Step51::postprocess_one_cell, + &Step51::NullFunction, + scratch, + task_data); + } + +// Compute some convergence rates, etc., and add to a table Vector difference_per_cell (triangulation.n_active_cells()); ComponentSelectFunction value_select (dim, dim+1); @@ -820,86 +914,84 @@ Step51::postprocess() &gradient_select); const double grad_error = difference_per_cell.l2_norm(); - convergence_table.add_value("cells", n_active_cells); - convergence_table.add_value("dofs", dof_handler.n_dofs()); - convergence_table.add_value("val L2", L2_error); - convergence_table.add_value("grad L2", grad_error); - - // construct post-processed solution with (hopefully) higher order of - // accuracy - QGauss quadrature(fe_u_post.degree+1); - FEValues fe_values(fe_u_post, quadrature, - update_values | update_JxW_values | - update_gradients); - - const unsigned int n_q_points = quadrature.size(); - std::vector u_values(n_q_points); - std::vector > u_gradients(n_q_points); - FEValuesExtractors::Vector fluxes(0); - FEValuesExtractors::Scalar scalar(dim); - FEValues fe_values_local(fe_local, quadrature, update_values); - FullMatrix cell_matrix(fe_u_post.dofs_per_cell, - fe_u_post.dofs_per_cell); - Vector cell_rhs(fe_u_post.dofs_per_cell); - Vector cell_sol(fe_u_post.dofs_per_cell); - - typename DoFHandler::active_cell_iterator - cell_loc = dof_handler_local.begin_active(), - cell = dof_handler_u_post.begin_active(), - endc = dof_handler_u_post.end(); - for ( ; cell != endc; ++cell, ++cell_loc) - { - fe_values.reinit(cell); - fe_values_local.reinit(cell_loc); - - fe_values_local[scalar].get_function_values(solution_local, u_values); - fe_values_local[fluxes].get_function_values(solution_local, u_gradients); - for (unsigned int i=1; idistribute_local_to_global(cell_sol, solution_u_post); - } - VectorTools::integrate_difference (dof_handler_u_post, solution_u_post, Solution(), difference_per_cell, QGauss(fe.degree+3), VectorTools::L2_norm); - double post_error = difference_per_cell.l2_norm(); - convergence_table.add_value("val L2-post", post_error); + const double post_error = difference_per_cell.l2_norm(); + + convergence_table.add_value("cells", triangulation.n_active_cells()); + convergence_table.add_value("dofs", dof_handler.n_dofs()); + convergence_table.add_value("val L2", L2_error); + convergence_table.add_value("grad L2", grad_error); + convergence_table.add_value("val L2-post", post_error); } +template +void +Step51::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, + PostProcessScratchData &scratch, + EmptyData &task_data) +{ + scratch.reset(); + + typename DoFHandler::active_cell_iterator + loc_cell (&triangulation, + cell->level(), + cell->index(), + &dof_handler_local); + + scratch.fe_values_local.reinit (loc_cell); + scratch.fe_values.reinit(cell); + + FEValuesExtractors::Vector fluxes(0); + FEValuesExtractors::Scalar scalar(dim); + + const unsigned int n_q_points = scratch.fe_values.get_quadrature().size(); + const unsigned int dofs_per_cell = scratch.fe_values.dofs_per_cell; + + scratch.fe_values_local[scalar].get_function_values(solution_local, scratch.u_values); + scratch.fe_values_local[fluxes].get_function_values(solution_local, scratch.u_gradients); + + double sum = 0; + for (unsigned int i=1; idistribute_local_to_global(scratch.cell_sol, solution_u_post); +} template -- 2.39.5