From 43648826e310fe73eb46124063a41efac9c29817 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 17 Jul 2015 14:53:34 -0500 Subject: [PATCH] Do the work of compute_fill() in parallel over the base elements. --- include/deal.II/fe/fe_system.h | 24 ++- source/fe/fe_system.cc | 295 ++++++++++++++++++--------------- 2 files changed, 186 insertions(+), 133 deletions(-) diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h index 7919bcf2f4..1514b1e0ed 100644 --- a/include/deal.II/fe/fe_system.h +++ b/include/deal.II/fe/fe_system.h @@ -840,7 +840,7 @@ protected: const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fe_data, - FEValuesData &data) const ; + FEValuesData &data) const; /** * Implementation of the same function in FiniteElement. @@ -856,7 +856,7 @@ protected: const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fe_data, - FEValuesData &data) const ; + FEValuesData &data) const; /** @@ -878,7 +878,7 @@ protected: CellSimilarity::Similarity cell_similarity, typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fe_data, - FEValuesData &data) const ; + FEValuesData &data) const; private: @@ -1037,6 +1037,24 @@ private: std::vector > hp_object_dof_identities (const FiniteElement &fe_other) const; + /** + * Compute the equivalent of compute_fill(), but only for the base element + * specified by the second-to-last argument. + */ + template + void + compute_fill_one_base ( + const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + const unsigned int base_element, + FEValuesData &data) const; + /** * Usually: Fields of cell-independent data. * diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index e460f80e0b..b9e8f74a32 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -1084,6 +1084,156 @@ FESystem::fill_fe_subface_values ( +template +template +void +FESystem::compute_fill_one_base ( + const Mapping &mapping, + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + CellSimilarity::Similarity cell_similarity, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + const unsigned int base_no, + FEValuesData &data) const +{ + InternalData &fe_data = static_cast (fedata); + const unsigned int n_q_points = quadrature.size(); + + const FiniteElement & + base_fe = base_element(base_no); + typename FiniteElement::InternalDataBase & + base_fe_data = fe_data.get_fe_data(base_no); + FEValuesData & + base_data = fe_data.get_fe_values_data(base_no); + + // fill_fe_face_values needs argument Quadrature for both cases + // dim_1==dim-1 and dim_1=dim. Hence the following workaround + const Quadrature *cell_quadrature = 0; + const Quadrature *face_quadrature = 0; + + // static cast to the common base class of quadrature being either + // Quadrature or Quadrature: + const Subscriptor *quadrature_base_pointer = &quadrature; + + if (face_no==invalid_face_number) + { + Assert(dim_1==dim, ExcDimensionMismatch(dim_1,dim)); + Assert (dynamic_cast *>(quadrature_base_pointer) != 0, + ExcInternalError()); + + cell_quadrature + = static_cast *>(quadrature_base_pointer); + } + else + { + Assert(dim_1==dim-1, ExcDimensionMismatch(dim_1,dim-1)); + Assert (dynamic_cast *>(quadrature_base_pointer) != 0, + ExcInternalError()); + + face_quadrature + = static_cast *>(quadrature_base_pointer); + } + + + // Copy all of the things that the mapping set in the FEValuesData + // that we store here into the corresponding objects we pass down + // to the various base elements. some of these arrays may be empty, + // in which case copying is cheap + base_data.JxW_values = data.JxW_values; + base_data.jacobians = data.jacobians; + base_data.jacobian_grads = data.jacobian_grads; + base_data.inverse_jacobians = data.inverse_jacobians; + base_data.quadrature_points = data.quadrature_points; + base_data.normal_vectors = data.normal_vectors; + base_data.boundary_forms = data.boundary_forms; + + // Make sure that in the case of fill_fe_values the data is only + // copied from base_data to data if base_data is changed. therefore + // use fe_fe_data.current_update_flags() + // + // for the case of fill_fe_(sub)face_values the data needs to be + // copied from base_data to data on each face, therefore use + // base_fe_data.update_flags. + if (face_no==invalid_face_number) + base_fe.fill_fe_values(mapping, cell, *cell_quadrature, mapping_data, + base_fe_data, base_data, cell_similarity); + else if (sub_no==invalid_face_number) + base_fe.fill_fe_face_values(mapping, cell, face_no, + *face_quadrature, mapping_data, base_fe_data, base_data); + else + base_fe.fill_fe_subface_values(mapping, cell, face_no, sub_no, + *face_quadrature, mapping_data, base_fe_data, base_data); + + // now data has been generated, so copy it. we used to work by + // looping over all base elements (i.e. this outer loop), then over + // multiplicity, then over the shape functions from that base + // element, but that requires that we can infer the global number of + // a shape function from its number in the base element. for that we + // had the component_to_system_table. + // + // however, this does of course no longer work since we have + // non-primitive elements. so we go the other way round: loop over + // all shape functions of the composed element, and here only treat + // those shape functions that belong to a given base element + //TODO: Introduce the needed table and loop only over base element shape functions. This here is not efficient at all AND very bad style + const UpdateFlags base_flags(dim_1==dim ? + base_fe_data.current_update_flags() : + base_fe_data.update_flags); + + // if the current cell is just a translation of the previous one, + // the underlying data has not changed, and we don't even need to + // enter this section + if (cell_similarity != CellSimilarity::translation) + for (unsigned int system_index=0; system_indexdofs_per_cell; + ++system_index) + if (this->system_to_base_table[system_index].first.first == base_no) + { + const unsigned int + base_index = this->system_to_base_table[system_index].second; + Assert (base_indexn_nonzero_components(i); + unsigned int in_index = 0; + for (unsigned int i=0; in_nonzero_components(system_index) == + base_fe.n_nonzero_components(base_index), + ExcInternalError()); + for (unsigned int s=0; sn_nonzero_components(system_index); ++s) + { + if (base_flags & update_values) + for (unsigned int q=0; q template void @@ -1142,138 +1292,23 @@ FESystem::compute_fill ( } } - // fill_fe_face_values needs argument Quadrature for both cases - // dim_1==dim-1 and dim_1=dim. Hence the following workaround - const Quadrature *cell_quadrature = 0; - const Quadrature *face_quadrature = 0; - - // static cast to the common base class of quadrature being either - // Quadrature or Quadrature: - const Subscriptor *quadrature_base_pointer = &quadrature; - - if (face_no==invalid_face_number) - { - Assert(dim_1==dim, ExcDimensionMismatch(dim_1,dim)); - Assert (dynamic_cast *>(quadrature_base_pointer) != 0, - ExcInternalError()); - - cell_quadrature - = static_cast *>(quadrature_base_pointer); - } - else - { - Assert(dim_1==dim-1, ExcDimensionMismatch(dim_1,dim-1)); - Assert (dynamic_cast *>(quadrature_base_pointer) != 0, - ExcInternalError()); - - face_quadrature - = static_cast *>(quadrature_base_pointer); - } - // let base elements update the necessary data + Threads::TaskGroup<> task_group; for (unsigned int base_no=0; base_non_base_elements(); ++base_no) - { - const FiniteElement & - base_fe = base_element(base_no); - typename FiniteElement::InternalDataBase & - base_fe_data = fe_data.get_fe_data(base_no); - FEValuesData & - base_data = fe_data.get_fe_values_data(base_no); - - // Copy all of the things that the mapping set in the FEValuesData - // that we store here into the corresponding objects we pass down - // to the various base elements. some of these arrays may be empty, - // in which case copying is cheap - base_data.JxW_values = data.JxW_values; - base_data.jacobians = data.jacobians; - base_data.jacobian_grads = data.jacobian_grads; - base_data.inverse_jacobians = data.inverse_jacobians; - base_data.quadrature_points = data.quadrature_points; - base_data.normal_vectors = data.normal_vectors; - base_data.boundary_forms = data.boundary_forms; - - - // Make sure that in the case of fill_fe_values the data is only - // copied from base_data to data if base_data is changed. therefore - // use fe_fe_data.current_update_flags() - // - // for the case of fill_fe_(sub)face_values the data needs to be - // copied from base_data to data on each face, therefore use - // base_fe_data.update_flags. - if (face_no==invalid_face_number) - base_fe.fill_fe_values(mapping, cell, *cell_quadrature, mapping_data, - base_fe_data, base_data, cell_similarity); - else if (sub_no==invalid_face_number) - base_fe.fill_fe_face_values(mapping, cell, face_no, - *face_quadrature, mapping_data, base_fe_data, base_data); - else - base_fe.fill_fe_subface_values(mapping, cell, face_no, sub_no, - *face_quadrature, mapping_data, base_fe_data, base_data); - - // now data has been generated, so copy it. we used to work by - // looping over all base elements (i.e. this outer loop), then over - // multiplicity, then over the shape functions from that base - // element, but that requires that we can infer the global number of - // a shape function from its number in the base element. for that we - // had the component_to_system_table. - // - // however, this does of course no longer work since we have - // non-primitive elements. so we go the other way round: loop over - // all shape functions of the composed element, and here only treat - // those shape functions that belong to a given base element -//TODO: Introduce the needed table and loop only over base element shape functions. This here is not efficient at all AND very bad style - const UpdateFlags base_flags(dim_1==dim ? - base_fe_data.current_update_flags() : - base_fe_data.update_flags); - - // if the current cell is just a translation of the previous one, - // the underlying data has not changed, and we don't even need to - // enter this section - if (cell_similarity != CellSimilarity::translation) - for (unsigned int system_index=0; system_indexdofs_per_cell; - ++system_index) - if (this->system_to_base_table[system_index].first.first == base_no) - { - const unsigned int - base_index = this->system_to_base_table[system_index].second; - Assert (base_indexn_nonzero_components(i); - unsigned int in_index = 0; - for (unsigned int i=0; in_nonzero_components(system_index) == - base_fe.n_nonzero_components(base_index), - ExcInternalError()); - for (unsigned int s=0; sn_nonzero_components(system_index); ++s) - { - if (base_flags & update_values) - for (unsigned int q=0; q(std_cxx11::bind (&FESystem::template compute_fill_one_base, + this, + std_cxx11::cref(mapping), + std_cxx11::cref(cell), + face_no, + sub_no, + std_cxx11::cref(quadrature), + CellSimilarity::none, + std_cxx11::ref(mapping_data), + std_cxx11::ref(fedata), + base_no, + std_cxx11::ref(data)))); + task_group.join_all(); } if (fe_data.compute_hessians && cell_similarity != CellSimilarity::translation) -- 2.39.5