From: scott.miller Date: Fri, 9 Aug 2013 19:17:54 +0000 (+0000) Subject: use meshworker for step-51 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c6304575c6f9938484c11c7ddde4fbc0aaf33662;p=dealii-svn.git use meshworker for step-51 git-svn-id: https://svn.dealii.org/trunk@30268 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc index 44599da697..0a097f89ab 100644 --- a/deal.II/examples/step-51/step-51.cc +++ b/deal.II/examples/step-51/step-51.cc @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -295,16 +296,22 @@ public: void run (); private: + + struct PerTaskData; + struct ScratchData; + void setup_system (); void assemble_system (const bool reconstruct_trace = false); + void assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch, + PerTaskData &task_data); + void copy_local_to_global(const PerTaskData &data); void solve (); void postprocess (); void refine_grid (const unsigned int cylce); void output_results (const unsigned int cycle); Triangulation triangulation; - - const MappingQ mapping; // The 'local' solutions are interior to each element. These // represent the primal solution field $u$ as well as the auxiliary @@ -363,7 +370,6 @@ private: template Step51::Step51 (const unsigned int degree, const RefinementMode refinement_mode) : - mapping (3), fe_local (FE_DGQ(degree), dim, FE_DGQ(degree), 1), dof_handler_local (triangulation), @@ -374,7 +380,178 @@ Step51::Step51 (const unsigned int degree, refinement_mode (refinement_mode) {} +template +struct Step51::PerTaskData +{ + FullMatrix cell_matrix; + Vector cell_vector; + std::vector dof_indices; + + bool trace_reconstruct; + + PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct) + : cell_matrix(n_dofs, n_dofs), + cell_vector(n_dofs), + dof_indices(n_dofs), + trace_reconstruct(trace_reconstruct) + {} + + void reset(){ + cell_matrix = 0.0; + cell_vector = 0.0; + } +}; + +template +struct Step51::ScratchData +{ + const Vector &solution; + Vector &solution_local; + + FEValues fe_values_local; + FEFaceValues fe_face_values_local; + FEFaceValues fe_face_values; + + FullMatrix ll_matrix; + FullMatrix lf_matrix; + FullMatrix fl_matrix; + FullMatrix tmp_matrix; + FullMatrix ff_matrix; + Vector l_rhs; + Vector f_rhs; + Vector tmp_rhs; + + std::vector > q_phi; + std::vector q_phi_div; + std::vector u_phi; + std::vector > u_phi_grad; + std::vector tr_phi; + std::vector trace_values; + + std::vector > fe_local_support_on_face; + std::vector > fe_support_on_face; + + bool trace_reconstruct; + + ConvectionVelocity convection_velocity; + RightHandSide right_hand_side; + const Solution exact_solution; + + // Full constructor + ScratchData(const Vector &solution, + Vector &solution_local, + const FiniteElement &fe, + const FiniteElement &fe_local, + const QGauss &quadrature_formula, + const QGauss &face_quadrature_formula, + const UpdateFlags local_flags, + const UpdateFlags local_face_flags, + const UpdateFlags flags, + const bool trace_reconstruct) + : + solution(solution), + solution_local(solution_local), + fe_values_local (fe_local, quadrature_formula, local_flags), + fe_face_values_local (fe_local, face_quadrature_formula, local_face_flags), + fe_face_values (fe, face_quadrature_formula, flags), + ll_matrix (fe_local.dofs_per_cell, fe_local.dofs_per_cell), + lf_matrix (fe_local.dofs_per_cell, fe.dofs_per_cell), + fl_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell), + tmp_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell), + ff_matrix (fe.dofs_per_cell, fe.dofs_per_cell), + l_rhs (fe_local.dofs_per_cell), + f_rhs (fe.dofs_per_cell), + tmp_rhs (fe_local.dofs_per_cell), + q_phi (fe_local.dofs_per_cell), + q_phi_div (fe_local.dofs_per_cell), + u_phi (fe_local.dofs_per_cell), + u_phi_grad (fe_local.dofs_per_cell), + tr_phi (fe.dofs_per_cell), + trace_values(face_quadrature_formula.size()), + fe_local_support_on_face(GeometryInfo::faces_per_cell), + fe_support_on_face(GeometryInfo::faces_per_cell), + trace_reconstruct (trace_reconstruct) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int i=0; i::faces_per_cell; ++face) + for (unsigned int i=0; i +void Step51::copy_local_to_global(const PerTaskData &data) +{ + if(data.trace_reconstruct == false) + constraints.distribute_local_to_global (data.cell_matrix, + data.cell_vector, + data.dof_indices, + system_matrix, system_rhs); +} template void @@ -397,9 +574,9 @@ Step51::setup_system () constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, constraints); typename FunctionMap::type boundary_functions; - Solution solution; - boundary_functions[0] = &solution; - VectorTools::project_boundary_values (mapping, dof_handler, + Solution solution_function; + boundary_functions[0] = &solution_function; + VectorTools::project_boundary_values (dof_handler, boundary_functions, QGauss(fe.degree+1), constraints); @@ -414,177 +591,166 @@ Step51::setup_system () system_matrix.reinit (sparsity_pattern); } +template +void +Step51::assemble_system (const bool trace_reconstruct) +{ + const QGauss quadrature_formula(fe.degree+1); + const QGauss face_quadrature_formula(fe.degree+1); + + const UpdateFlags local_flags (update_values | update_gradients | + update_JxW_values | update_quadrature_points); + + const UpdateFlags local_face_flags (update_values); + + const UpdateFlags flags ( update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + PerTaskData task_data (fe.dofs_per_cell, + trace_reconstruct); + ScratchData scratch (solution, + solution_local, + fe, fe_local, + quadrature_formula, + face_quadrature_formula, + local_flags, + local_face_flags, + flags, + trace_reconstruct); + + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &Step51::assemble_system_one_cell, + &Step51::copy_local_to_global, + scratch, + task_data); +} template void -Step51::assemble_system (const bool trace_reconstruct) +Step51::assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch, + PerTaskData &task_data) { - QGauss quadrature_formula(fe.degree+1); - QGauss face_quadrature_formula(fe.degree+1); - - FEValues fe_values_local (mapping, fe_local, quadrature_formula, - update_values | update_gradients | - update_JxW_values | update_quadrature_points); - FEFaceValues fe_face_values (mapping, fe, face_quadrature_formula, - update_values | update_normal_vectors | - update_quadrature_points | - update_JxW_values); - FEFaceValues fe_face_values_local (mapping, fe_local, - face_quadrature_formula, - update_values); - - const unsigned int n_q_points = quadrature_formula.size(); - const unsigned int n_face_q_points = face_quadrature_formula.size(); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int loc_dofs_per_cell = fe_local.dofs_per_cell; - - FullMatrix ll_matrix (loc_dofs_per_cell, loc_dofs_per_cell); - FullMatrix lf_matrix (loc_dofs_per_cell, dofs_per_cell); - FullMatrix fl_matrix (dofs_per_cell, loc_dofs_per_cell); - FullMatrix tmp_matrix (dofs_per_cell, loc_dofs_per_cell); - FullMatrix ff_matrix (dofs_per_cell, dofs_per_cell); - Vector l_rhs (loc_dofs_per_cell); - Vector f_rhs (dofs_per_cell); - Vector tmp_rhs (loc_dofs_per_cell); - - std::vector dof_indices (dofs_per_cell); - std::vector loc_dof_indices (loc_dofs_per_cell); - - std::vector > q_phi (loc_dofs_per_cell); - std::vector q_phi_div (loc_dofs_per_cell); - std::vector u_phi (loc_dofs_per_cell); - std::vector > u_phi_grad (loc_dofs_per_cell); - std::vector tr_phi (dofs_per_cell); - - std::vector trace_values(n_face_q_points); + // Construct iterator for dof_handler_local + typename DoFHandler::active_cell_iterator + loc_cell (&triangulation, + cell->level(), + cell->index(), + &dof_handler_local); + + const unsigned int n_q_points = scratch.fe_values_local.get_quadrature().size(); + const unsigned int n_face_q_points = scratch.fe_face_values_local.get_quadrature().size(); + + // const unsigned int dofs_per_cell = scratch.fe_face_values.get_fe().dofs_per_cell; + const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell; // Choose stabilization parameter to be 5 * diffusion = 5 const double tau_stab_diffusion = 5.; - ConvectionVelocity convection_velocity; - RightHandSide right_hand_side; - const Solution exact_solution; - const FEValuesExtractors::Vector fluxes (0); const FEValuesExtractors::Scalar scalar (dim); - std::vector > - fe_local_support_on_face(GeometryInfo::faces_per_cell); - for (unsigned int face=0; face::faces_per_cell; ++face) - for (unsigned int i=0; i > - fe_support_on_face(GeometryInfo::faces_per_cell); - for (unsigned int face=0; face::faces_per_cell; ++face) - for (unsigned int i=0; i::active_cell_iterator - cell = dof_handler.begin_active(), - loc_cell = dof_handler_local.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell, ++loc_cell) - { - ll_matrix = 0; - l_rhs = 0; - if (!trace_reconstruct) + scratch.ll_matrix = 0; + scratch.l_rhs = 0; + if (!scratch.trace_reconstruct) { - lf_matrix = 0; - fl_matrix = 0; - ff_matrix = 0; - f_rhs = 0; + scratch.lf_matrix = 0; + scratch.fl_matrix = 0; + scratch.ff_matrix = 0; + scratch.f_rhs = 0; } - fe_values_local.reinit (loc_cell); + scratch.fe_values_local.reinit (loc_cell); for (unsigned int q=0; q convection - = convection_velocity.value(fe_values_local.quadrature_point(q)); - const double JxW = fe_values_local.JxW(q); + = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q)); + const double JxW = scratch.fe_values_local.JxW(q); for (unsigned int k=0; k::faces_per_cell; ++face) { - fe_face_values_local.reinit(loc_cell, face); - fe_face_values.reinit(cell, face); - if (trace_reconstruct) - fe_face_values.get_function_values (solution, trace_values); + scratch.fe_face_values_local.reinit(loc_cell, face); + scratch.fe_face_values.reinit(cell, face); + if (scratch.trace_reconstruct) + scratch.fe_face_values.get_function_values (scratch.solution, scratch.trace_values); for (unsigned int q=0; q normal = fe_face_values.normal_vector(q); + const double JxW = scratch.fe_face_values.JxW(q); + const Point normal = scratch.fe_face_values.normal_vector(q); const Tensor<1,dim> convection - = convection_velocity.value(fe_face_values.quadrature_point(q)); + = scratch.convection_velocity.value(scratch.fe_face_values.quadrature_point(q)); const double tau_stab = (tau_stab_diffusion + std::abs(convection * normal)); - for (unsigned int k=0; k::assemble_system (const bool trace_reconstruct) (cell->face(face)->boundary_indicator() == 1)) { const double neumann_value = - exact_solution.value(fe_face_values.quadrature_point(q)); - for (unsigned int i=0; iget_dof_indices(dof_indices); - constraints.distribute_local_to_global (ff_matrix, f_rhs, - dof_indices, - system_matrix, system_rhs); + scratch.fl_matrix.mmult(scratch.tmp_matrix, scratch.ll_matrix); + scratch.tmp_matrix.vmult_add(scratch.f_rhs, scratch.l_rhs); + scratch.tmp_matrix.mmult(scratch.ff_matrix, scratch.lf_matrix, true); + cell->get_dof_indices(task_data.dof_indices); + task_data.cell_matrix = scratch.ff_matrix; + task_data.cell_vector = scratch.f_rhs; } else { - ll_matrix.vmult(tmp_rhs, l_rhs); - loc_cell->set_dof_values(tmp_rhs, solution_local); + scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs); + loc_cell->set_dof_values(scratch.tmp_rhs, scratch.solution_local); } - } } @@ -658,6 +822,7 @@ void Step51::solve () system_matrix.clear(); sparsity_pattern.reinit(0,0,0,1); + constraints.distribute(solution); // update local values @@ -674,7 +839,7 @@ Step51::postprocess() Vector difference_per_cell (triangulation.n_active_cells()); ComponentSelectFunction value_select (dim, dim+1); - VectorTools::integrate_difference (mapping, dof_handler_local, + VectorTools::integrate_difference (dof_handler_local, solution_local, SolutionAndGradient(), difference_per_cell, @@ -685,7 +850,7 @@ Step51::postprocess() ComponentSelectFunction gradient_select (std::pair(0, dim), dim+1); - VectorTools::integrate_difference (mapping, dof_handler_local, + VectorTools::integrate_difference (dof_handler_local, solution_local, SolutionAndGradient(), difference_per_cell, @@ -702,7 +867,7 @@ Step51::postprocess() // construct post-processed solution with (hopefully) higher order of // accuracy QGauss quadrature(fe_u_post.degree+1); - FEValues fe_values(mapping, fe_u_post, quadrature, + FEValues fe_values(fe_u_post, quadrature, update_values | update_JxW_values | update_gradients); @@ -711,7 +876,7 @@ Step51::postprocess() std::vector > u_gradients(n_q_points); FEValuesExtractors::Vector fluxes(0); FEValuesExtractors::Scalar scalar(dim); - FEValues fe_values_local(mapping, fe_local, quadrature, update_values); + 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); @@ -764,7 +929,7 @@ Step51::postprocess() cell->distribute_local_to_global(cell_sol, solution_u_post); } - VectorTools::integrate_difference (mapping, dof_handler_u_post, + VectorTools::integrate_difference (dof_handler_u_post, solution_u_post, Solution(), difference_per_cell, @@ -846,32 +1011,18 @@ void Step51::output_results (const unsigned int cycle) template void Step51::refine_grid (const unsigned int cycle) { - const bool do_cube = true; if (cycle == 0) { - if (!do_cube) - { - GridGenerator::hyper_ball (triangulation); - static const HyperBallBoundary boundary; - triangulation.set_boundary(0, boundary); - triangulation.refine_global(6-2*dim); - } - else - GridGenerator::subdivided_hyper_cube (triangulation, 2, -1, 1); + GridGenerator::subdivided_hyper_cube (triangulation, 2, -1, 1); } else switch (refinement_mode) { case global_refinement: { - if (do_cube) - { - triangulation.clear(); - GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1); - triangulation.refine_global(3-dim+cycle/2); - } - else - triangulation.refine_global (1); + triangulation.clear(); + GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1); + triangulation.refine_global(3-dim+cycle/2); break; }