From 1c87117796b129bb9759a84b9af42425caefaa5d Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 19 Nov 2020 21:39:36 +0100 Subject: [PATCH] Introduce MappingCollection in MatrixFree --- include/deal.II/matrix_free/mapping_info.h | 28 ++++--- .../matrix_free/mapping_info.templates.h | 62 ++++++++------ include/deal.II/matrix_free/matrix_free.h | 80 ++++++++++++------- .../matrix_free/matrix_free.templates.h | 12 ++- source/matrix_free/matrix_free.inst.in | 4 +- 5 files changed, 122 insertions(+), 64 deletions(-) diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h index d717e3fc17..941d530556 100644 --- a/include/deal.II/matrix_free/mapping_info.h +++ b/include/deal.II/matrix_free/mapping_info.h @@ -27,6 +27,7 @@ #include #include +#include #include #include @@ -335,10 +336,10 @@ namespace internal 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 std::shared_ptr> &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); @@ -355,7 +356,7 @@ namespace internal const std::vector> &cells, const FaceInfo & faces, const std::vector &active_fe_index, - const Mapping & mapping); + const std::shared_ptr> &mapping); /** * Return the type of a given cell as detected during initialization. @@ -444,7 +445,12 @@ namespace internal face_data_by_cells; /** - * The pointer to the underlying Mapping object. + * The pointer to the underlying hp::MappingCollection object. + */ + std::shared_ptr> mapping_collection; + + /** + * The pointer to the first entry of mapping_collection. */ SmartPointer> mapping; @@ -484,8 +490,8 @@ namespace internal initialize_cells( const dealii::Triangulation & tria, const std::vector> &cells, - const std::vector &active_fe_index, - const Mapping & mapping); + const std::vector & active_fe_index, + const dealii::hp::MappingCollection &mapping); /** * Computes the information in the given faces, called within @@ -496,8 +502,8 @@ namespace internal const dealii::Triangulation & tria, const std::vector> &cells, const std::vector> - & faces, - const Mapping &mapping); + & faces, + const dealii::hp::MappingCollection &mapping); /** * Computes the information in the given faces, called within @@ -507,7 +513,7 @@ namespace internal initialize_faces_by_cells( const dealii::Triangulation & tria, const std::vector> &cells, - const Mapping & mapping); + const dealii::hp::MappingCollection & mapping); /** * Helper function to determine which update flags must be set in the diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index e81fcb0fe8..4c12701240 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -272,7 +272,8 @@ namespace internal face_data_by_cells.clear(); cell_type.clear(); face_type.clear(); - mapping = nullptr; + mapping_collection = nullptr; + mapping = nullptr; } @@ -331,15 +332,16 @@ namespace internal const std::vector> &cells, const FaceInfo & face_info, const std::vector & active_fe_index, - const Mapping & mapping, - const std::vector> & quad, + const std::shared_ptr> &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) { clear(); - this->mapping = &mapping; + this->mapping_collection = mapping; + this->mapping = &mapping->operator[](0); cell_data.resize(quad.size()); face_data.resize(quad.size()); @@ -435,16 +437,16 @@ namespace internal // In case we have no hp adaptivity (active_fe_index is empty), we have // cells, and the mapping is MappingQGeneric or a derived class, we can // use the fast method. - if (active_fe_index.empty() && !cells.empty() && - dynamic_cast *>(&mapping)) + if (active_fe_index.empty() && !cells.empty() && mapping->size() == 1 && + dynamic_cast *>(&mapping->operator[](0))) compute_mapping_q(tria, cells, face_info.faces); else { // Could call these functions in parallel, but not useful because // the work inside is nicely split up already - initialize_cells(tria, cells, active_fe_index, mapping); - initialize_faces(tria, cells, face_info.faces, mapping); - initialize_faces_by_cells(tria, cells, mapping); + initialize_cells(tria, cells, active_fe_index, *mapping); + initialize_faces(tria, cells, face_info.faces, *mapping); + initialize_faces_by_cells(tria, cells, *mapping); } } @@ -457,7 +459,7 @@ namespace internal const std::vector> &cells, const FaceInfo & face_info, const std::vector & active_fe_index, - const Mapping & mapping) + const std::shared_ptr> &mapping) { AssertDimension(cells.size() / VectorizedArrayType::size(), cell_type.size()); @@ -469,18 +471,19 @@ namespace internal for (auto &data : face_data_by_cells) data.clear_data_fields(); - this->mapping = &mapping; + this->mapping_collection = mapping; + this->mapping = &mapping->operator[](0); - if (active_fe_index.empty() && !cells.empty() && - dynamic_cast *>(&mapping)) + if (active_fe_index.empty() && !cells.empty() && mapping->size() == 1 && + dynamic_cast *>(&mapping->operator[](0))) compute_mapping_q(tria, cells, face_info.faces); else { // Could call these functions in parallel, but not useful because // the work inside is nicely split up already - initialize_cells(tria, cells, active_fe_index, mapping); - initialize_faces(tria, cells, face_info.faces, mapping); - initialize_faces_by_cells(tria, cells, mapping); + initialize_cells(tria, cells, active_fe_index, *mapping); + initialize_faces(tria, cells, face_info.faces, *mapping); + initialize_faces_by_cells(tria, cells, *mapping); } } @@ -843,12 +846,15 @@ namespace internal const dealii::Triangulation & tria, const std::vector> &cells, const std::vector & active_fe_index, - const Mapping & mapping, + const dealii::hp::MappingCollection & mapping_in, MappingInfo &mapping_info, std::pair>, CompressedCellData> &data) { + AssertDimension(mapping_in.size(), 1); + const auto &mapping = mapping_in[0]; + FE_Nothing dummy_fe; // when we make comparisons about the size of Jacobians we need to @@ -1515,7 +1521,7 @@ namespace internal const dealii::Triangulation & tria, const std::vector> &cells, const std::vector & active_fe_index, - const Mapping & mapping) + const dealii::hp::MappingCollection & mapping) { const unsigned int n_cells = cells.size(); const unsigned int n_lanes = VectorizedArrayType::size(); @@ -1748,13 +1754,16 @@ namespace internal const std::vector> &cells, const std::vector> & faces, - const Mapping & mapping, + const dealii::hp::MappingCollection & mapping_in, MappingInfo &mapping_info, std::pair< std::vector< MappingInfoStorage>, CompressedFaceData> &data) { + AssertDimension(mapping_in.size(), 1); + const auto &mapping = mapping_in[0]; + FE_Nothing dummy_fe; std::vector>>> @@ -2392,7 +2401,7 @@ namespace internal const dealii::Triangulation & tria, const std::vector> & cells, const std::vector> &faces, - const Mapping &mapping) + const dealii::hp::MappingCollection &mapping) { face_type.resize(faces.size(), general); @@ -2591,8 +2600,11 @@ namespace internal { // step 1: extract quadrature point data with the data appropriate for // MappingQGeneric + AssertDimension(this->mapping_collection->size(), 1); + const MappingQGeneric *mapping_q = - dynamic_cast *>(&*this->mapping); + dynamic_cast *>( + &this->mapping_collection->operator[](0)); Assert(mapping_q != nullptr, ExcInternalError()); const unsigned int mapping_degree = mapping_q->get_degree(); @@ -2937,7 +2949,7 @@ namespace internal // transitioned to extracting the information from cell quadrature // points but we need to figure out the correct indices of neighbors // within the list of arrays still - initialize_faces_by_cells(tria, cell_array, *this->mapping); + initialize_faces_by_cells(tria, cell_array, *this->mapping_collection); } @@ -2947,11 +2959,15 @@ namespace internal MappingInfo::initialize_faces_by_cells( const dealii::Triangulation & tria, const std::vector> &cells, - const Mapping & mapping) + const dealii::hp::MappingCollection & mapping_in) { if (update_flags_faces_by_cells == update_default) return; + + AssertDimension(mapping_in.size(), 1); + const auto &mapping = mapping_in[0]; + const unsigned int n_quads = face_data_by_cells.size(); const unsigned int n_lanes = VectorizedArrayType::size(); UpdateFlags update_flags = diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 165e5a23ee..43d3f88261 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -35,6 +35,7 @@ #include #include +#include #include #include @@ -580,9 +581,9 @@ public: * not allowed. In that case, use the initialization function with * several DoFHandler arguments. */ - template + template void - reinit(const Mapping & mapping, + reinit(const MappingType & mapping, const DoFHandler & dof_handler, const AffineConstraints &constraint, const QuadratureType & quad, @@ -622,9 +623,9 @@ public: * can be sets independently from the number of DoFHandlers, when several * elements are always integrated with the same quadrature formula. */ - template + template void - reinit(const Mapping & mapping, + reinit(const MappingType & mapping, const std::vector *> & dof_handler, const std::vector *> &constraint, const std::vector & quad, @@ -635,9 +636,12 @@ public: * * @deprecated Use the overload taking a DoFHandler object instead. */ - template + template DEAL_II_DEPRECATED void - reinit(const Mapping & mapping, + reinit(const MappingType & mapping, const std::vector & dof_handler, const std::vector *> &constraint, const std::vector & quad, @@ -675,9 +679,9 @@ public: * as might be necessary when several components in a vector-valued problem * are integrated together based on the same quadrature formula. */ - template + template void - reinit(const Mapping & mapping, + reinit(const MappingType & mapping, const std::vector *> & dof_handler, const std::vector *> &constraint, const QuadratureType & quad, @@ -688,9 +692,12 @@ public: * * @deprecated Use the overload taking a DoFHandler object instead. */ - template + template DEAL_II_DEPRECATED void - reinit(const Mapping & mapping, + reinit(const MappingType & mapping, const std::vector & dof_handler, const std::vector *> &constraint, const QuadratureType & quad, @@ -742,6 +749,12 @@ public: void update_mapping(const Mapping &mapping); + /** + * Same as above but with hp::MappingCollection. + */ + void + update_mapping(const std::shared_ptr> &mapping); + /** * Clear all data fields and brings the class into a condition similar to * after having called the default constructor. @@ -2030,7 +2043,7 @@ private: template void internal_reinit( - const Mapping & mapping, + const std::shared_ptr> & mapping, const std::vector *> & dof_handlers, const std::vector *> &constraint, const std::vector & locally_owned_set, @@ -2853,7 +2866,8 @@ MatrixFree::reinit( std::vector> quad_hp; quad_hp.emplace_back(quad); - internal_reinit(StaticMappingQ1::mapping, + internal_reinit(std::make_shared>( + StaticMappingQ1::mapping), dof_handlers, constraints, locally_owned_sets, @@ -2864,10 +2878,10 @@ MatrixFree::reinit( template -template +template void MatrixFree::reinit( - const Mapping & mapping, + const MappingType & mapping, const DoFHandler & dof_handler, const AffineConstraints &constraints_in, const QuadratureType & quad, @@ -2887,7 +2901,7 @@ MatrixFree::reinit( std::vector> quad_hp; quad_hp.emplace_back(quad); - internal_reinit(mapping, + internal_reinit(std::make_shared>(mapping), dof_handlers, constraints, locally_owned_sets, @@ -2913,7 +2927,9 @@ MatrixFree::reinit( std::vector> quad_hp; for (unsigned int q = 0; q < quad.size(); ++q) quad_hp.emplace_back(quad[q]); - internal_reinit(StaticMappingQ1::mapping, + + internal_reinit(std::make_shared>( + StaticMappingQ1::mapping), dof_handler, constraint, locally_owned_set, @@ -2924,10 +2940,13 @@ MatrixFree::reinit( template -template +template void MatrixFree::reinit( - const Mapping & mapping, + const MappingType & mapping, const std::vector & dof_handler, const std::vector *> &constraint, const std::vector & quad, @@ -2961,7 +2980,9 @@ MatrixFree::reinit( dof_handler, additional_data.mg_level); std::vector> quad_hp; quad_hp.emplace_back(quad); - internal_reinit(StaticMappingQ1::mapping, + + internal_reinit(std::make_shared>( + StaticMappingQ1::mapping), dof_handler, constraint, locally_owned_set, @@ -2994,10 +3015,10 @@ MatrixFree::reinit( template -template +template void MatrixFree::reinit( - const Mapping & mapping, + const MappingType & mapping, const std::vector *> & dof_handler, const std::vector *> &constraint, const QuadratureType & quad, @@ -3009,7 +3030,8 @@ MatrixFree::reinit( dof_handler, additional_data.mg_level); std::vector> quad_hp; quad_hp.emplace_back(quad); - internal_reinit(mapping, + + internal_reinit(std::make_shared>(mapping), dof_handler, constraint, locally_owned_set, @@ -3020,10 +3042,10 @@ MatrixFree::reinit( template -template +template void MatrixFree::reinit( - const Mapping & mapping, + const MappingType & mapping, const std::vector *> & dof_handler, const std::vector *> &constraint, const std::vector & quad, @@ -3036,7 +3058,8 @@ MatrixFree::reinit( std::vector> quad_hp; for (unsigned int q = 0; q < quad.size(); ++q) quad_hp.emplace_back(quad[q]); - internal_reinit(mapping, + + internal_reinit(std::make_shared>(mapping), dof_handler, constraint, locally_owned_set, @@ -3047,10 +3070,13 @@ MatrixFree::reinit( template -template +template void MatrixFree::reinit( - const Mapping & mapping, + const MappingType & mapping, const std::vector & dof_handler, const std::vector *> &constraint, const QuadratureType & quad, diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 3b5643dd84..01809a8ea9 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -290,7 +290,7 @@ template template void MatrixFree::internal_reinit( - const Mapping & mapping, + const std::shared_ptr> & mapping, const std::vector *> & dof_handler, const std::vector *> &constraint, const std::vector & locally_owned_dofs, @@ -466,6 +466,16 @@ template void MatrixFree::update_mapping( const Mapping &mapping) +{ + update_mapping(std::make_shared>(mapping)); +} + + + +template +void +MatrixFree::update_mapping( + const std::shared_ptr> &mapping) { AssertDimension(shape_info.size(1), mapping_info.cell_data.size()); mapping_info.update_mapping(dof_handlers[0]->get_triangulation(), diff --git a/source/matrix_free/matrix_free.inst.in b/source/matrix_free/matrix_free.inst.in index aaacd0001c..0ae3ae5d0d 100644 --- a/source/matrix_free/matrix_free.inst.in +++ b/source/matrix_free/matrix_free.inst.in @@ -35,7 +35,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_scalar_vectorized::value_type, deal_II_scalar_vectorized>:: internal_reinit( - const Mapping &, + const std::shared_ptr> &, const std::vector *> &, const std::vector< const AffineConstraints *> &, @@ -52,7 +52,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_float_vectorized::value_type, deal_II_float_vectorized>:: internal_reinit( - const Mapping &, + const std::shared_ptr> &, const std::vector *> &, const std::vector *> &, const std::vector &, -- 2.39.5