From 9f2f4df14f4283d381d2a979943a4b88416d5ca4 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 12 Nov 2020 14:01:23 +0100 Subject: [PATCH] Make MatrixFree compile for arbitrary dim quadrature rules --- include/deal.II/hp/q_collection.h | 54 ++++++++++++++----- include/deal.II/matrix_free/fe_evaluation.h | 2 +- include/deal.II/matrix_free/mapping_info.h | 17 +++--- .../matrix_free/mapping_info.templates.h | 51 ++++++++++++++---- include/deal.II/matrix_free/matrix_free.h | 16 +++--- .../matrix_free/matrix_free.templates.h | 6 +-- include/deal.II/matrix_free/shape_info.h | 12 ++--- .../matrix_free/shape_info.templates.h | 7 ++- source/matrix_free/matrix_free.inst.in | 4 +- source/matrix_free/shape_info.inst.in | 16 ++++++ 10 files changed, 133 insertions(+), 52 deletions(-) diff --git a/include/deal.II/hp/q_collection.h b/include/deal.II/hp/q_collection.h index 0ead13dc43..073bac7968 100644 --- a/include/deal.II/hp/q_collection.h +++ b/include/deal.II/hp/q_collection.h @@ -52,13 +52,20 @@ namespace hp */ QCollection() = default; + /** + * Copy constructor. + */ + template + QCollection(const QCollection &other); + /** * Conversion constructor. This constructor creates a QCollection from a * single quadrature rule. More quadrature formulas can be added with * push_back(), if desired, though it would probably be clearer to add all * mappings the same way. */ - explicit QCollection(const Quadrature &quadrature); + template + explicit QCollection(const Quadrature &quadrature); /** * Constructor. This constructor creates a QCollection from one or @@ -93,8 +100,9 @@ namespace hp * is later destroyed by this object upon destruction of the entire * collection. */ + template void - push_back(const Quadrature &new_quadrature); + push_back(const Quadrature &new_quadrature); /** * Return a reference to the quadrature rule specified by the argument. @@ -150,21 +158,41 @@ namespace hp /* --------------- inline functions ------------------- */ + template + template + QCollection::QCollection(const QCollection &other) + { + for (unsigned int i = 0; i < other.size(); ++i) + push_back(other[i]); + } + + + template template QCollection::QCollection(const QTypes &... quadrature_objects) { - static_assert(is_base_of_all, QTypes...>::value, - "Not all of the input arguments of this function " - "are derived from Quadrature!"); - // loop over all of the given arguments and add the quadrature objects to // this collection. Inlining the definition of q_pointers causes internal // compiler errors on GCC 7.1.1 so we define it separately: - const auto q_pointers = { - (static_cast *>(&quadrature_objects))...}; - for (const auto p : q_pointers) - push_back(*p); + if (is_base_of_all, QTypes...>::value) + { + const auto q_pointers = { + (reinterpret_cast *>(&quadrature_objects))...}; + for (const auto p : q_pointers) + push_back(*p); + } + else if (is_base_of_all, QTypes...>::value) + { + const auto q_pointers = { + (reinterpret_cast *>(&quadrature_objects))...}; + for (const auto p : q_pointers) + push_back(*p); + } + else + { + Assert(false, ExcNotImplemented()); + } } @@ -223,7 +251,8 @@ namespace hp template - inline QCollection::QCollection(const Quadrature &quadrature) + template + inline QCollection::QCollection(const Quadrature &quadrature) { quadratures.push_back(std::make_shared>(quadrature)); } @@ -239,8 +268,9 @@ namespace hp template + template inline void - QCollection::push_back(const Quadrature &new_quadrature) + QCollection::push_back(const Quadrature &new_quadrature) { quadratures.push_back( std::make_shared>(new_quadrature)); diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 6446d407dd..f37f136841 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -3331,7 +3331,7 @@ inline FEEvaluationBaseData:: , // select the correct base element from the given FE component data(new internal::MatrixFreeFunctions::ShapeInfo( - quadrature, + Quadrature(quadrature), fe, fe.component_to_base_index(first_selected_component).first)) , jacobian(nullptr) diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h index b2ba970de7..e4852ea9eb 100644 --- a/include/deal.II/matrix_free/mapping_info.h +++ b/include/deal.II/matrix_free/mapping_info.h @@ -320,17 +320,17 @@ namespace internal * CellIterator::level() and CellIterator::index(), in order to allow * for different kinds of iterators, e.g. standard DoFHandler, * multigrid, etc.) on a fixed Triangulation. In addition, a mapping - * and several quadrature formulas are given. + * and several 1D quadrature formulas are given. */ void initialize( const dealii::Triangulation & tria, const std::vector> &cells, const FaceInfo & faces, - const std::vector & active_fe_index, - const Mapping & mapping, - const std::vector> &quad, - const UpdateFlags update_flags_cells, + const std::vector & active_fe_index, + const Mapping & mapping, + const std::vector> &quad, + const UpdateFlags update_flags_cells, const UpdateFlags update_flags_boundary_faces, const UpdateFlags update_flags_inner_faces, const UpdateFlags update_flags_faces_by_cells); @@ -506,9 +506,10 @@ namespace internal * internal functions to initialize all data as requested by the user. */ static UpdateFlags - compute_update_flags(const UpdateFlags update_flags, - const std::vector> &quad = - std::vector>()); + compute_update_flags( + const UpdateFlags update_flags, + const std::vector> &quad = + std::vector>()); }; diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 85d747d18c..629b1d6166 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -252,8 +252,8 @@ namespace internal template UpdateFlags MappingInfo::compute_update_flags( - const UpdateFlags update_flags, - const std::vector> &quad) + const UpdateFlags update_flags, + const std::vector> &quad) { // this class is build around the evaluation of jacobians, so compute // them in any case. The Jacobians will be inverted manually. Since we @@ -304,7 +304,7 @@ namespace internal const FaceInfo & face_info, const std::vector & active_fe_index, const Mapping & mapping, - const std::vector> & quad, + const std::vector> & quad, const UpdateFlags update_flags_cells, const UpdateFlags update_flags_boundary_faces, const UpdateFlags update_flags_inner_faces, @@ -337,18 +337,49 @@ namespace internal AssertIndexRange(0, n_hp_quads); cell_data[my_q].descriptor.resize(n_hp_quads); for (unsigned int q = 0; q < n_hp_quads; ++q) - cell_data[my_q].descriptor[q].initialize(quad[my_q][q], - update_default); + { + Assert(quad[my_q][q].is_tensor_product(), ExcNotImplemented()); + for (unsigned int i = 1; i < dim; ++i) + { + Assert(quad[my_q][q].get_tensor_basis()[0] == + quad[my_q][q].get_tensor_basis()[i], + ExcNotImplemented()); + } + + cell_data[my_q].descriptor[q].initialize( + quad[my_q][q].get_tensor_basis()[0], update_default); + } face_data[my_q].descriptor.resize(n_hp_quads); for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq) - face_data[my_q].descriptor[hpq].initialize( - quad[my_q][hpq], update_flags_boundary_faces); + { + Assert(quad[my_q][hpq].is_tensor_product(), ExcNotImplemented()); + for (unsigned int i = 1; i < dim; ++i) + { + Assert(quad[my_q][hpq].get_tensor_basis()[0] == + quad[my_q][hpq].get_tensor_basis()[i], + ExcNotImplemented()); + } + + face_data[my_q].descriptor[hpq].initialize( + quad[my_q][hpq].get_tensor_basis()[0], + update_flags_boundary_faces); + } face_data_by_cells[my_q].descriptor.resize(n_hp_quads); for (unsigned int hpq = 0; hpq < n_hp_quads; ++hpq) - face_data_by_cells[my_q].descriptor[hpq].initialize(quad[my_q][hpq], - update_default); + { + Assert(quad[my_q][hpq].is_tensor_product(), ExcNotImplemented()); + for (unsigned int i = 1; i < dim; ++i) + { + Assert(quad[my_q][hpq].get_tensor_basis()[0] == + quad[my_q][hpq].get_tensor_basis()[i], + ExcNotImplemented()); + } + + face_data_by_cells[my_q].descriptor[hpq].initialize( + quad[my_q][hpq].get_tensor_basis()[0], update_default); + } } // In case we have no hp adaptivity (active_fe_index is empty), we have @@ -2588,7 +2619,7 @@ namespace internal { FE_DGQ fe_geometry(mapping_degree); for (unsigned int my_q = 0; my_q < cell_data.size(); ++my_q) - shape_infos[my_q].reinit(cell_data[my_q].descriptor[0].quadrature_1d, + shape_infos[my_q].reinit(cell_data[my_q].descriptor[0].quadrature, fe_geometry); } diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 7dadca5886..ea4576eb49 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -2025,14 +2025,14 @@ private: * This is the actual reinit function that sets up the indices for the * DoFHandler case. */ - template + template void internal_reinit( const Mapping & mapping, const std::vector *> & dof_handlers, const std::vector *> &constraint, const std::vector & locally_owned_set, - const std::vector> & quad, + const std::vector> & quad, const AdditionalData & additional_data); /** @@ -2843,7 +2843,7 @@ MatrixFree::reinit( internal::MatrixFreeImplementation::extract_locally_owned_index_sets( dof_handlers, additional_data.mg_level); - std::vector> quad_hp; + std::vector> quad_hp; quad_hp.emplace_back(quad); internal_reinit(StaticMappingQ1::mapping, @@ -2877,7 +2877,7 @@ MatrixFree::reinit( internal::MatrixFreeImplementation::extract_locally_owned_index_sets( dof_handlers, additional_data.mg_level); - std::vector> quad_hp; + std::vector> quad_hp; quad_hp.emplace_back(quad); internal_reinit(mapping, @@ -2903,7 +2903,7 @@ MatrixFree::reinit( std::vector locally_owned_set = internal::MatrixFreeImplementation::extract_locally_owned_index_sets( dof_handler, additional_data.mg_level); - std::vector> quad_hp; + std::vector> quad_hp; for (unsigned int q = 0; q < quad.size(); ++q) quad_hp.emplace_back(quad[q]); internal_reinit(StaticMappingQ1::mapping, @@ -2949,7 +2949,7 @@ MatrixFree::reinit( std::vector locally_owned_set = internal::MatrixFreeImplementation::extract_locally_owned_index_sets( dof_handler, additional_data.mg_level); - std::vector> quad_hp; + std::vector> quad_hp; quad_hp.emplace_back(quad); internal_reinit(StaticMappingQ1::mapping, dof_handler, @@ -2994,7 +2994,7 @@ MatrixFree::reinit( std::vector locally_owned_set = internal::MatrixFreeImplementation::extract_locally_owned_index_sets( dof_handler, additional_data.mg_level); - std::vector> quad_hp; + std::vector> quad_hp; quad_hp.emplace_back(quad); internal_reinit(mapping, dof_handler, @@ -3020,7 +3020,7 @@ MatrixFree::reinit( std::vector locally_owned_set = internal::MatrixFreeImplementation::extract_locally_owned_index_sets( dof_handler, additional_data.mg_level); - std::vector> quad_hp; + std::vector> quad_hp; for (unsigned int q = 0; q < quad.size(); ++q) quad_hp.emplace_back(quad[q]); internal_reinit(mapping, diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 03f224613c..fdd4febe56 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -275,14 +275,14 @@ MatrixFree::copy_from( template -template +template void MatrixFree::internal_reinit( const Mapping & mapping, const std::vector *> & dof_handler, const std::vector *> &constraint, const std::vector & locally_owned_dofs, - const std::vector> & quad, + const std::vector> & quad, const typename MatrixFree::AdditionalData &additional_data) { @@ -1293,7 +1293,7 @@ MatrixFree::initialize_indices( Table<2, internal::MatrixFreeFunctions::ShapeInfo> shape_info_dummy( shape_info.size(0), shape_info.size(2)); { - QGauss<1> quad(1); + QGauss quad(1); for (unsigned int no = 0, c = 0; no < dof_handlers.size(); no++) for (unsigned int b = 0; b < dof_handlers[no]->get_fe(0).n_base_elements(); diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index 64daf31106..bb896a9c54 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -312,8 +312,8 @@ namespace internal /** * Constructor that initializes the data fields using the reinit method. */ - template - ShapeInfo(const Quadrature<1> & quad, + template + ShapeInfo(const Quadrature & quad, const FiniteElement &fe, const unsigned int base_element = 0); @@ -325,9 +325,9 @@ namespace internal * dimensional element by a tensor product and that the zeroth shape * function in zero evaluates to one. */ - template + template void - reinit(const Quadrature<1> & quad, + reinit(const Quadrature & quad, const FiniteElement &fe_dim, const unsigned int base_element = 0); @@ -519,8 +519,8 @@ namespace internal // ------------------------------------------ inline functions template - template - inline ShapeInfo::ShapeInfo(const Quadrature<1> & quad, + template + inline ShapeInfo::ShapeInfo(const Quadrature & quad, const FiniteElement &fe_in, const unsigned int base_element_number) : element_type(tensor_general) diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 83e7362098..873300c9e5 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -128,12 +128,15 @@ namespace internal template - template + template void - ShapeInfo::reinit(const Quadrature<1> & quad, + ShapeInfo::reinit(const Quadrature & quad_in, const FiniteElement &fe_in, const unsigned int base_element_number) { + Assert(quad_in.is_tensor_product(), ExcNotImplemented()); + const auto quad = quad_in.get_tensor_basis()[0]; + const FiniteElement *fe = &fe_in.base_element(base_element_number); n_dimensions = dim; n_components = fe_in.n_components(); diff --git a/source/matrix_free/matrix_free.inst.in b/source/matrix_free/matrix_free.inst.in index 5fae952f73..b256093dcb 100644 --- a/source/matrix_free/matrix_free.inst.in +++ b/source/matrix_free/matrix_free.inst.in @@ -40,7 +40,7 @@ for (deal_II_dimension : DIMENSIONS; const std::vector< const AffineConstraints *> &, const std::vector &, - const std::vector> &, + const std::vector> &, const AdditionalData &); template const DoFHandler & @@ -62,7 +62,7 @@ for (deal_II_dimension : DIMENSIONS; const std::vector *> &, const std::vector *> &, const std::vector &, - const std::vector> &, + const std::vector> &, const AdditionalData &); } diff --git a/source/matrix_free/shape_info.inst.in b/source/matrix_free/shape_info.inst.in index 293099d241..30be5ee7aa 100644 --- a/source/matrix_free/shape_info.inst.in +++ b/source/matrix_free/shape_info.inst.in @@ -16,11 +16,19 @@ for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS) { + template void internal::MatrixFreeFunctions::ShapeInfo:: + reinit( + const Quadrature &, + const FiniteElement &, + const unsigned int); + +#if deal_II_dimension > 1 template void internal::MatrixFreeFunctions::ShapeInfo:: reinit( const Quadrature<1> &, const FiniteElement &, const unsigned int); +#endif } for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) @@ -37,11 +45,19 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) for (deal_II_dimension : DIMENSIONS; deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED) { + template void internal::MatrixFreeFunctions:: + ShapeInfo::reinit( + const Quadrature &, + const FiniteElement &, + const unsigned int); + +#if deal_II_dimension > 1 template void internal::MatrixFreeFunctions:: ShapeInfo::reinit( const Quadrature<1> &, const FiniteElement &, const unsigned int); +#endif } for (deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED) -- 2.39.5