From: turcksin Date: Thu, 17 Oct 2013 17:22:16 +0000 (+0000) Subject: Use WorkStream in Step-14. Need to update the documentation and to use TBB instead... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b9130d85fa4688bd907ed4664747ab8dcb9ec2be;p=dealii-svn.git Use WorkStream in Step-14. Need to update the documentation and to use TBB instead of new_thread. git-svn-id: https://svn.dealii.org/trunk@31286 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-13/step-13.cc b/deal.II/examples/step-13/step-13.cc index 32ee1a99e5..eea40b8811 100644 --- a/deal.II/examples/step-13/step-13.cc +++ b/deal.II/examples/step-13/step-13.cc @@ -666,12 +666,12 @@ namespace Step13 void assemble_matrix (const typename DoFHandler::active_cell_iterator &cell, Assembler::Scratch &scratch, - Assembler::CopyData ©_data); + Assembler::CopyData ©_data) const; void copy_local_to_global(Assembler::CopyData const ©_data, - LinearSystem &linear_system); + LinearSystem &linear_system) const; }; @@ -818,7 +818,7 @@ namespace Step13 void Solver::assemble_matrix (const typename DoFHandler::active_cell_iterator &cell, Assembler::Scratch &scratch, - Assembler::CopyData ©_data) + Assembler::CopyData ©_data) const { FEValues fe_values (*fe, *quadrature, update_gradients | update_JxW_values); @@ -847,7 +847,7 @@ namespace Step13 template void Solver::copy_local_to_global(Assembler::CopyData const ©_data, - LinearSystem &linear_system) + LinearSystem &linear_system) const { // In the step-9 program, we have shown that you have to use the // mutex to lock the matrix when copying the elements from the local diff --git a/deal.II/examples/step-14/step-14.cc b/deal.II/examples/step-14/step-14.cc index 58359c9e97..934e0a9dd3 100644 --- a/deal.II/examples/step-14/step-14.cc +++ b/deal.II/examples/step-14/step-14.cc @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -59,6 +60,23 @@ namespace Step14 { using namespace dealii; + namespace Assembler + { + struct Scratch + { + Scratch() {} + }; + + struct CopyData + { + CopyData() {} + + unsigned int dofs_per_cell; + FullMatrix cell_matrix; + std::vector local_dof_indices; + }; + } + // @sect3{Evaluating the solution} // As mentioned in the introduction, significant parts of the program have @@ -474,10 +492,14 @@ namespace Step14 assemble_linear_system (LinearSystem &linear_system); void - assemble_matrix (LinearSystem &linear_system, - const typename DoFHandler::active_cell_iterator &begin_cell, - const typename DoFHandler::active_cell_iterator &end_cell, - Threads::Mutex &mutex) const; + assemble_matrix (const typename DoFHandler::active_cell_iterator &cell, + Assembler::Scratch &scratch, + Assembler::CopyData ©_data) const; + + + void + copy_local_to_global(Assembler::CopyData const ©_data, + LinearSystem &linear_system) const; }; @@ -550,15 +572,13 @@ namespace Step14 dof_handler.end (), n_threads); - Threads::Mutex mutex; - Threads::ThreadGroup<> threads; - for (unsigned int thread=0; thread::assemble_matrix, - *this, - linear_system, - thread_ranges[thread].first, - thread_ranges[thread].second, - mutex); + Assembler::Scratch scratch; + Assembler::CopyData copy_data; + WorkStream::run(dof_handler.begin_active(),dof_handler.end(), + std::bind(&Solver::assemble_matrix,this,std_cxx1x::_1,std_cxx1x::_2,std_cxx1x::_3), + std::bind(&Solver::copy_local_to_global,this,std_cxx1x::_1,std_cxx1x::ref(linear_system)), + scratch,copy_data); + assemble_rhs (linear_system.rhs); linear_system.hanging_node_constraints.condense (linear_system.rhs); @@ -569,7 +589,6 @@ namespace Step14 *boundary_values, boundary_value_map); - threads.join_all (); linear_system.hanging_node_constraints.condense (linear_system.matrix); MatrixTools::apply_boundary_values (boundary_value_map, @@ -581,47 +600,77 @@ namespace Step14 template void - Solver::assemble_matrix (LinearSystem &linear_system, - const typename DoFHandler::active_cell_iterator &begin_cell, - const typename DoFHandler::active_cell_iterator &end_cell, - Threads::Mutex &mutex) const + Solver::assemble_matrix (const typename DoFHandler::active_cell_iterator &cell, + Assembler::Scratch &scratch, + Assembler::CopyData ©_data) const { FEValues fe_values (*fe, *quadrature, update_gradients | update_JxW_values); - const unsigned int dofs_per_cell = fe->dofs_per_cell; + copy_data.dofs_per_cell = fe->dofs_per_cell; const unsigned int n_q_points = quadrature->size(); - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + copy_data.cell_matrix = FullMatrix (copy_data.dofs_per_cell, copy_data.dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + copy_data.local_dof_indices.resize(copy_data.dofs_per_cell); - for (typename DoFHandler::active_cell_iterator cell=begin_cell; - cell!=end_cell; ++cell) - { - cell_matrix = 0; + fe_values.reinit (cell); - fe_values.reinit (cell); + for (unsigned int q_point=0; q_pointget_dof_indices (copy_data.local_dof_indices); + } - cell->get_dof_indices (local_dof_indices); - Threads::Mutex::ScopedLock lock (mutex); - for (unsigned int i=0; i + void + Solver::copy_local_to_global(Assembler::CopyData const ©_data, + LinearSystem &linear_system) const + { + for (unsigned int i=0; iThreads::ThreadGroup class here, but rather use + // the one created thread object directly to wait for this particular + // thread's exit. + // + // Note that taking up the address of the + // DoFTools::make_hanging_node_constraints function is a + // little tricky, since there are actually three of them, one for each + // supported space dimension. Taking addresses of overloaded functions is + // somewhat complicated in C++, since the address-of operator + // & in that case returns more like a set of values (the + // addresses of all functions with that name), and selecting the right one + // is then the next step. If the context dictates which one to take (for + // example by assigning to a function pointer of known type), then the + // compiler can do that by itself, but if this set of pointers shall be + // given as the argument to a function that takes a template, the compiler + // could choose all without having a preference for one. We therefore have + // to make it clear to the compiler which one we would like to have; for + // this, we could use a cast, but for more clarity, we assign it to a + // temporary mhnc_p (short for pointer to + // make_hanging_node_constraints) with the right type, and using + // this pointer instead. template Solver::LinearSystem:: LinearSystem (const DoFHandler &dof_handler)