From c63fce4aaf0123f70e6a8cda3cacf9d4e0514979 Mon Sep 17 00:00:00 2001 From: turcksin Date: Mon, 14 Oct 2013 21:37:39 +0000 Subject: [PATCH] Replace the second Threads by WorkStream in matrix_tools.cc git-svn-id: https://svn.dealii.org/trunk@31232 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/numerics/matrix_tools.cc | 657 +++++++++++------------- 1 file changed, 300 insertions(+), 357 deletions(-) diff --git a/deal.II/source/numerics/matrix_tools.cc b/deal.II/source/numerics/matrix_tools.cc index 6f8bff1a54..d4f762140d 100644 --- a/deal.II/source/numerics/matrix_tools.cc +++ b/deal.II/source/numerics/matrix_tools.cc @@ -612,7 +612,7 @@ namespace MatrixCreator Scratch() {} }; - template + template struct CopyData { CopyData() {}; @@ -622,13 +622,13 @@ namespace MatrixCreator unsigned int dofs_per_cell; std::vector dofs; std::vector > dof_is_on_face; - typename DoFHandler::active_cell_iterator cell; + typename DH::active_cell_iterator cell; std::vector > cell_matrix; std::vector > cell_vector; }; - template - CopyData::CopyData(CopyData const &data) : + 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), @@ -870,8 +870,8 @@ namespace MatrixCreator void create_boundary_mass_matrix_1 (typename DoFHandler::active_cell_iterator const &cell, MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch, - MatrixCreator::internal::AssemblerBoundary::CopyData - ©_data, + MatrixCreator::internal::AssemblerBoundary::CopyData > ©_data, Mapping const &mapping, FiniteElement const &fe, Quadrature const &q, @@ -893,7 +893,6 @@ namespace MatrixCreator copy_data.cell = cell; copy_data.dofs_per_cell = fe.dofs_per_cell; - UpdateFlags update_flags = UpdateFlags (update_values | update_JxW_values | update_normal_vectors | @@ -1058,7 +1057,9 @@ namespace MatrixCreator cell->face(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)); + copy_data.dof_is_on_face.push_back(std::vector (copy_data.dofs_per_cell)); + // check for each of the dofs on this cell + // whether it is on the face for (unsigned int i=0; i - void copy_boundary_mass_matrix_1(MatrixCreator::internal::AssemblerBoundary::CopyData - const ©_data, + 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, @@ -1141,8 +1142,8 @@ namespace MatrixCreator 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, + MatrixCreator::internal::AssemblerBoundary::CopyData > ©_data, Mapping<2,3> const &mapping, FiniteElement<2,3> const &fe, Quadrature<1> const &q, @@ -1160,8 +1161,8 @@ namespace MatrixCreator 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, + MatrixCreator::internal::AssemblerBoundary::CopyData > ©_data, Mapping<1,3> const &mapping, FiniteElement<1,3> const &fe, Quadrature<0> const &q, @@ -1221,19 +1222,19 @@ namespace MatrixCreator AssertDimension (n_components, component_mapping.size()); MatrixCreator::internal::AssemblerBoundary::Scratch scratch; - MatrixCreator::internal::AssemblerBoundary::CopyData copy_data; + MatrixCreator::internal::AssemblerBoundary::CopyData > copy_data; WorkStream::run(dof.begin_active(),dof.end(), static_cast::active_cell_iterator const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &, - MatrixCreator::internal::AssemblerBoundary::CopyData &)> > + 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( + ::CopyData > 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))), @@ -1247,37 +1248,28 @@ namespace MatrixCreator template void - create_boundary_mass_matrix_1 (std_cxx1x::tuple &, - const hp::DoFHandler &, - const hp::QCollection &> 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_hp_boundary_mass_matrix_1 (typename hp::DoFHandler::active_cell_iterator const + &cell, + MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch, + MatrixCreator::internal::AssemblerBoundary + ::CopyData > ©_data, + hp::MappingCollection const &mapping, + hp::FECollection const &fe_collection, + hp::QCollection const &q, + const typename FunctionMap::type &boundary_functions, + Function const *const coefficient, + std::vector const &component_mapping) { - const hp::MappingCollection &mapping = std_cxx1x::get<0>(commons); - const hp::DoFHandler &dof = std_cxx1x::get<1>(commons); - const hp::QCollection& q = std_cxx1x::get<2>(commons); - const hp::FECollection &fe_collection = dof.get_fe(); const unsigned int n_components = fe_collection.n_components(); const unsigned int n_function_components = boundary_functions.begin()->second->n_components; const bool fe_is_system = (n_components != 1); -#ifdef DEBUG - if (true) - { - types::global_dof_index max_element = static_cast(0); - for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); - i!=dof_to_boundary_mapping.end(); ++i) - if ((*i != hp::DoFHandler::invalid_dof_index) && - (*i > max_element)) - max_element = *i; - Assert (max_element == matrix.n()-1, ExcInternalError()); - }; -#endif + const FiniteElement &fe = cell->get_fe(); + const unsigned int dofs_per_face = fe.dofs_per_face; + + copy_data.cell = cell; + copy_data.dofs_per_cell = fe.dofs_per_cell; + copy_data.dofs.resize(copy_data.dofs_per_cell); + cell->get_dof_indices (copy_data.dofs); const unsigned int max_dofs_per_cell = fe_collection.max_dofs_per_cell(), max_dofs_per_face = fe_collection.max_dofs_per_face(); @@ -1300,148 +1292,131 @@ namespace MatrixCreator std::vector rhs_values_scalar; std::vector > rhs_values_system; - std::vector dofs (max_dofs_per_cell); - std::vector dofs_on_face_vector (max_dofs_per_face); + 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(max_dofs_per_cell); + copy_data.dofs.resize(copy_data.dofs_per_cell); + cell->get_dof_indices (copy_data.dofs); + // 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(); - typename hp::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()) - { - x_fe_values.reinit (cell, face); - const FEFaceValues &fe_values = x_fe_values.get_present_fe_values (); + 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()) + { + x_fe_values.reinit (cell, face); - const FiniteElement &fe = cell->get_fe(); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int dofs_per_face = fe.dofs_per_face; + const FEFaceValues &fe_values = x_fe_values.get_present_fe_values (); - cell_matrix.reinit (dofs_per_cell, dofs_per_cell); - cell_vector.reinit (dofs_per_cell); - cell_matrix = 0; - cell_vector = 0; + 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)); - if (fe_is_system) - // FE has several components - { - rhs_values_system.resize (fe_values.n_quadrature_points, - Vector(n_function_components)); - boundary_functions.find(cell->face(face)->boundary_indicator()) - ->second->vector_value_list (fe_values.get_quadrature_points(), - rhs_values_system); + if (fe_is_system) + // FE has several components + { + rhs_values_system.resize (fe_values.n_quadrature_points, + Vector(n_function_components)); + boundary_functions.find(cell->face(face)->boundary_indicator()) + ->second->vector_value_list (fe_values.get_quadrature_points(), + rhs_values_system); - if (coefficient != 0) - { - if (coefficient->n_components==1) - { - coefficient_values.resize (fe_values.n_quadrature_points); - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; point(n_components)); - coefficient->vector_value_list (fe_values.get_quadrature_points(), - coefficient_vector_values); - for (unsigned int point=0; pointn_components==1) { - const double weight = fe_values.JxW(point); - for (unsigned int i=0; ivalue_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; pointface(face)->boundary_indicator()) - ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); - - if (coefficient != 0) + else + { + coefficient_vector_values.resize (fe_values.n_quadrature_points, + Vector(n_components)); + coefficient->vector_value_list (fe_values.get_quadrature_points(), + coefficient_vector_values); + for (unsigned int point=0; pointvalue_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; pointface(face)->boundary_indicator()) + ->second->value_list (fe_values.get_quadrature_points(), rhs_values_scalar); + + if (coefficient != 0) + { + coefficient_values.resize (fe_values.n_quadrature_points); + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); for (unsigned int point=0; pointface(face)->get_dof_indices (dofs_on_face_vector, + cell->active_fe_index()); + // 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)); + // check for each of the dofs on this cell + // whether it is on the face + for (unsigned int i=0; iget_dof_indices (dofs); - cell->face(face)->get_dof_indices (dofs_on_face_vector, - cell->active_fe_index()); - - // check for each of the - // dofs on this cell - // whether it is on the - // face - for (unsigned int i=0; i + void copy_hp_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()) + { #ifdef DEBUG - double max_diag_entry = 0; - for (unsigned int i=0; i max_diag_entry) - max_diag_entry = std::fabs(cell_matrix(i,i)); + // in debug mode: compute an element in the matrix which is + // guaranteed to belong to a boundary dof. We do this to check that the + // entries in the cell matrix are guaranteed to be zero if the + // respective dof is not on the boundary. Since because of + // round-off, the actual value of the matrix entry may be + // only close to zero, we assert that it is small relative to an element + // which is guaranteed to be nonzero. (absolute smallness does not + // suffice since the size of the domain scales in here) + // + // for this purpose we seek the diagonal of the matrix, where there + // must be an element belonging to the boundary. we take the maximum + // diagonal entry. + types::global_dof_index max_element = static_cast(0); + for (std::vector::const_iterator i=dof_to_boundary_mapping.begin(); + i!=dof_to_boundary_mapping.end(); ++i) + if ((*i != hp::DoFHandler::invalid_dof_index) && + (*i > max_element)) + max_element = *i; + Assert (max_element == matrix.n()-1, ExcInternalError()); + + double max_diag_entry = 0; + for (unsigned int i=0; i max_diag_entry) + max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i)); #endif - // lock the matrix - Threads::Mutex::ScopedLock lock (mutex); - for (unsigned int i=0; i threads; - - // define starting and end point - // for each thread - typedef typename hp::DoFHandler::active_cell_iterator active_cell_iterator; - std::vector > thread_ranges - = Threads::split_range (dof.begin_active(), - dof.end(), n_threads); - - typedef std_cxx1x::tuple &, - const hp::DoFHandler &, - const hp::QCollection&> Commons; - - // mutex to synchronise access to - // the matrix - Threads::Mutex mutex; - - // 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_hp_boundary_mass_matrix_1,std_cxx1x::_1,std_cxx1x::_2, + std_cxx1x::_3, + std_cxx1x::cref(mapping),std_cxx1x::cref(fe_collection),std_cxx1x::cref(q), + std_cxx1x::cref(boundary_functions),coefficient, + std_cxx1x::cref(component_mapping))), + static_cast > const &)> > (std_cxx1x::bind( + copy_hp_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