From cce9806f1c293bac7dfcd0e1ec0454b7163aed19 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Mon, 14 Oct 2013 19:19:08 +0000 Subject: [PATCH] Replace the first Threads by WorkStream in matrix_tools.cc git-svn-id: https://svn.dealii.org/trunk@31230 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/numerics/matrix_tools.cc | 581 ++++++++++++------------ 1 file changed, 297 insertions(+), 284 deletions(-) diff --git a/deal.II/source/numerics/matrix_tools.cc b/deal.II/source/numerics/matrix_tools.cc index af0e940cc8..6f8bff1a54 100644 --- a/deal.II/source/numerics/matrix_tools.cc +++ b/deal.II/source/numerics/matrix_tools.cc @@ -16,7 +16,6 @@ #include #include -#include #include #include #include @@ -603,6 +602,41 @@ namespace MatrixCreator data.dof_indices, *matrix); } + + + + namespace AssemblerBoundary + { + struct Scratch + { + Scratch() {} + }; + + template + struct CopyData + { + CopyData() {}; + + CopyData(CopyData const &data); + + unsigned int dofs_per_cell; + std::vector dofs; + std::vector > dof_is_on_face; + typename DoFHandler::active_cell_iterator cell; + std::vector > cell_matrix; + std::vector > cell_vector; + }; + + template + CopyData::CopyData(CopyData const &data) : + dofs_per_cell(data.dofs_per_cell), + dofs(data.dofs), + dof_is_on_face(data.dof_is_on_face), + cell(data.cell), + cell_matrix(data.cell_matrix), + cell_vector(data.cell_vector) + {} + } } } @@ -834,36 +868,30 @@ namespace MatrixCreator { template void - create_boundary_mass_matrix_1 (std_cxx1x::tuple &, - const DoFHandler &, - const Quadrature & > commons, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *const coefficient, - const std::vector &component_mapping, - const MatrixCreator::internal::IteratorRange > range, - Threads::Mutex &mutex) + create_boundary_mass_matrix_1 (typename DoFHandler::active_cell_iterator const &cell, + MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch, + MatrixCreator::internal::AssemblerBoundary::CopyData + ©_data, + Mapping const &mapping, + FiniteElement const &fe, + Quadrature const &q, + typename FunctionMap::type const &boundary_functions, + Function const *const coefficient, + std::vector const &component_mapping) + { // All assertions for this function // are in the calling function // before creating threads. - const Mapping &mapping = std_cxx1x::get<0>(commons); - const DoFHandler &dof = std_cxx1x::get<1>(commons); - const Quadrature& q = std_cxx1x::get<2>(commons); - - const FiniteElement &fe = dof.get_fe(); - const unsigned int n_components = fe.n_components(); + const unsigned int n_components = fe.n_components(); const unsigned int n_function_components = boundary_functions.begin()->second->n_components; - const bool fe_is_system = (n_components != 1); + const bool fe_is_system = (n_components != 1); const bool fe_is_primitive = fe.is_primitive(); - const unsigned int dofs_per_cell = fe.dofs_per_cell, - dofs_per_face = fe.dofs_per_face; - - FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); - Vector cell_vector(dofs_per_cell); + const unsigned int dofs_per_face = fe.dofs_per_face; + + copy_data.cell = cell; + copy_data.dofs_per_cell = fe.dofs_per_cell; UpdateFlags update_flags = UpdateFlags (update_values | @@ -883,236 +911,244 @@ namespace MatrixCreator std::vector rhs_values_scalar (fe_values.n_quadrature_points); std::vector > rhs_values_system (fe_values.n_quadrature_points, - Vector(n_function_components)); + Vector(n_function_components)); + + copy_data.dofs.resize(copy_data.dofs_per_cell); + cell->get_dof_indices (copy_data.dofs); - std::vector dofs (dofs_per_cell); std::vector dofs_on_face_vector (dofs_per_face); - // for each dof on the cell, have a - // flag whether it is on the face - std::vector dof_is_on_face(dofs_per_cell); + // Because CopyData objects are reused and that push_back is used, + // dof_is_on_face, cell_matrix, and cell_vector must be cleared before + // they are reused + copy_data.dof_is_on_face.clear(); + copy_data.cell_matrix.clear(); + copy_data.cell_vector.clear(); + + for (unsigned int face=0; face::faces_per_cell; ++face) + // check if this face is on that part of + // the boundary we are interested in + if (boundary_functions.find(cell->face(face)->boundary_indicator()) != + boundary_functions.end()) + { + copy_data.cell_matrix.push_back(FullMatrix (copy_data.dofs_per_cell, + copy_data.dofs_per_cell)); + copy_data.cell_vector.push_back(Vector (copy_data.dofs_per_cell)); + fe_values.reinit (cell, face); - typename DoFHandler::active_cell_iterator cell = range.first; - for (; cell!=range.second; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - // check if this face is on that part of - // the boundary we are interested in - if (boundary_functions.find(cell->face(face)->boundary_indicator()) != - boundary_functions.end()) - { - cell_matrix = 0; - cell_vector = 0; + if (fe_is_system) + // FE has several components + { + boundary_functions.find(cell->face(face)->boundary_indicator()) + ->second->vector_value_list (fe_values.get_quadrature_points(), + rhs_values_system); + + if (coefficient_is_vector) + // If coefficient is + // vector valued, fill + // all components + coefficient->vector_value_list (fe_values.get_quadrature_points(), + coefficient_vector_values); + else + { + // If a scalar + // function is + // given, update + // the values, if + // not, use the + // default one set + // in the + // constructor above + if (coefficient != 0) + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + // Copy scalar + // values into vector + for (unsigned int point=0; point > normal_adjustment(fe_values.n_quadrature_points, + std::vector(n_components, 1.)); - if (fe_is_system) - // FE has several components - { - boundary_functions.find(cell->face(face)->boundary_indicator()) - ->second->vector_value_list (fe_values.get_quadrature_points(), - rhs_values_system); + for (unsigned int comp = 0; comp &base = fe.base_element(fe.component_to_base_index(comp).first); + const unsigned int bcomp = fe.component_to_base_index(comp).second; - if (coefficient_is_vector) - // If coefficient is - // vector valued, fill - // all components - coefficient->vector_value_list (fe_values.get_quadrature_points(), - coefficient_vector_values); - else - { - // If a scalar - // function is - // geiven, update - // the values, if - // not, use the - // default one set - // in the - // constructor above - if (coefficient != 0) - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - // Copy scalar - // values into vector + if (!base.conforms(FiniteElementData::H1) && + base.conforms(FiniteElementData::Hdiv)) for (unsigned int point=0; point > normal_adjustment(fe_values.n_quadrature_points, - std::vector(n_components, 1.)); - - for (unsigned int comp = 0; comp &base = fe.base_element(fe.component_to_base_index(comp).first); - const unsigned int bcomp = fe.component_to_base_index(comp).second; - - if (!base.conforms(FiniteElementData::H1) && - base.conforms(FiniteElementData::Hdiv)) - for (unsigned int point=0; pointface(face)->boundary_indicator()) - ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); + normal_adjustment[point][comp] = fe_values.normal_vector(point)(bcomp) + * fe_values.normal_vector(point)(bcomp); + } - if (coefficient != 0) - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; pointget_dof_indices (dofs); - cell->face(face)->get_dof_indices (dofs_on_face_vector); - for (unsigned int i=0; iface(face)->boundary_indicator()) + ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); - // lock the matrix - Threads::Mutex::ScopedLock lock (mutex); - for (unsigned int i=0; ivalue_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; pointface(face)->get_dof_indices (dofs_on_face_vector); + // for each dof on the cell, have a + // flag whether it is on the face + copy_data.dof_is_on_face.push_back(std::vector(copy_data.dofs_per_cell)); + for (unsigned int i=0; i + void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData + const ©_data, + typename FunctionMap::type const &boundary_functions, + std::vector const &dof_to_boundary_mapping, + SparseMatrix &matrix, + Vector &rhs_vector) + { + // now transfer cell matrix and vector to the whole boundary matrix + // + // in the following: dof[i] holds the global index of the i-th degree of + // freedom on the present cell. If it is also a dof on the boundary, it + // must be a nonzero entry in the dof_to_boundary_mapping and then + // the boundary index of this dof is dof_to_boundary_mapping[dof[i]]. + // + // if dof[i] is not on the boundary, it should be zero on the boundary + // therefore on all quadrature points and finally all of its + // entries in the cell matrix and vector should be zero. If not, we + // throw an error (note: because of the evaluation of the shape + // functions only up to machine precision, the term "must be zero" + // really should mean: "should be very small". since this is only an + // assertion and not part of the code, we may choose "very small" + // quite arbitrarily) + // + // the main problem here is that the matrix or vector entry should also + // be zero if the degree of freedom dof[i] is on the boundary, but not + // on the present face, i.e. on another face of the same cell also + // on the boundary. We can therefore not rely on the + // dof_to_boundary_mapping[dof[i]] being !=-1, we really have to + // determine whether dof[i] is a dof on the present face. We do so + // by getting the dofs on the face into @p{dofs_on_face_vector}, + // a vector as always. Usually, searching in a vector is + // inefficient, so we copy the dofs into a set, which enables binary + // searches. + unsigned int pos(0); + for (unsigned int face=0; face::faces_per_cell; ++face) + { + // check if this face is on that part of + // the boundary we are interested in + if (boundary_functions.find(copy_data.cell->face(face)->boundary_indicator()) != + boundary_functions.end()) + { + for (unsigned int i=0; i void - create_boundary_mass_matrix_1<2,3> (std_cxx1x::tuple &, - const DoFHandler<2,3> &, - const Quadrature<1> & > , - SparseMatrix &, - const FunctionMap<3>::type &, - Vector &, - std::vector &, - const Function<3> *const , - const std::vector &, - const MatrixCreator::internal::IteratorRange > , - Threads::Mutex &) + create_boundary_mass_matrix_1<2,3> (typename DoFHandler<2,3>::active_cell_iterator const &cell, + MatrixCreator::internal::AssemblerBoundary::Scratch const + &scratch, + MatrixCreator::internal::AssemblerBoundary::CopyData<2,3> + ©_data, + Mapping<2,3> const &mapping, + FiniteElement<2,3> const &fe, + Quadrature<1> const &q, + typename FunctionMap<3>::type const &boundary_functions, + Function<3> const *const coefficient, + std::vector const &component_mapping) { Assert(false,ExcNotImplemented()); } @@ -1121,17 +1157,17 @@ namespace MatrixCreator template <> void - create_boundary_mass_matrix_1<1,3> (std_cxx1x::tuple &, - const DoFHandler<1,3> &, - const Quadrature<0> & > , - SparseMatrix &, - const FunctionMap<3>::type &, - Vector &, - std::vector &, - const Function<3> *const , - const std::vector &, - const MatrixCreator::internal::IteratorRange > , - Threads::Mutex &) + create_boundary_mass_matrix_1<1,3> (typename DoFHandler<1,3>::active_cell_iterator const &cell, + MatrixCreator::internal::AssemblerBoundary::Scratch const + &scratch, + MatrixCreator::internal::AssemblerBoundary::CopyData<1,3> + ©_data, + Mapping<1,3> const &mapping, + FiniteElement<1,3> const &fe, + Quadrature<0> const &q, + typename FunctionMap<3>::type const &boundary_functions, + Function<3> const *const coefficient, + std::vector const &component_mapping) { Assert(false,ExcNotImplemented()); } @@ -1184,47 +1220,24 @@ namespace MatrixCreator else AssertDimension (n_components, component_mapping.size()); - const unsigned int n_threads = multithread_info.n_threads(); - Threads::ThreadGroup<> threads; - - // define starting and end point - // for each thread - typedef typename DoFHandler::active_cell_iterator active_cell_iterator; - std::vector > thread_ranges - = Threads::split_range (dof.begin_active(), - dof.end(), n_threads); - - // mutex to synchronise access to - // the matrix - Threads::Mutex mutex; - - typedef std_cxx1x::tuple &, - const DoFHandler &, - const Quadrature&> Commons; - - // then assemble in parallel - typedef void (*create_boundary_mass_matrix_1_t) - (Commons, - SparseMatrix &matrix, - const typename FunctionMap::type &boundary_functions, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *const coefficient, - const std::vector &component_mapping, - const MatrixCreator::internal::IteratorRange > range, - Threads::Mutex &mutex); - create_boundary_mass_matrix_1_t p - = &create_boundary_mass_matrix_1; - -//TODO: Use WorkStream here - for (unsigned int thread=0; thread copy_data; + + WorkStream::run(dof.begin_active(),dof.end(), + static_cast::active_cell_iterator + const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, + MatrixCreator::internal::AssemblerBoundary::CopyData &)> > + (std_cxx1x::bind(create_boundary_mass_matrix_1,std_cxx1x::_1,std_cxx1x::_2, + std_cxx1x::_3, + std_cxx1x::cref(mapping),std_cxx1x::cref(fe),std_cxx1x::cref(q), + std_cxx1x::cref(boundary_functions),coefficient, + std_cxx1x::cref(component_mapping))), + static_cast const &)> > (std_cxx1x::bind( + copy_boundary_mass_matrix_1,std_cxx1x::_1, + std_cxx1x::cref(boundary_functions),std_cxx1x::cref(dof_to_boundary_mapping), + std_cxx1x::ref(matrix),std_cxx1x::ref(rhs_vector))), + scratch,copy_data); } -- 2.39.5