From: Martin Kronbichler Date: Tue, 25 Feb 2020 18:51:29 +0000 (+0100) Subject: Enable refreshing MatrixFree geometry data for changed mapping X-Git-Tag: v9.2.0-rc1~432^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4f0523453f55ceb3028e8b9a6760c19ea2d270f4;p=dealii.git Enable refreshing MatrixFree geometry data for changed mapping --- diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h index 759b5870df..6d43632194 100644 --- a/include/deal.II/matrix_free/mapping_info.h +++ b/include/deal.II/matrix_free/mapping_info.h @@ -145,6 +145,12 @@ namespace internal */ unsigned int n_q_points; + /** + * Original one-dimensional quadrature formula applied on the given + * cell or face. + */ + Quadrature<1> quadrature_1d; + /** * Quadrature formula applied on the given cell or face. */ @@ -265,6 +271,12 @@ namespace internal */ AlignedVector> quadrature_points; + /** + * Clears all data fields except the descriptor vector. + */ + void + clear_data_fields(); + /** * Returns the quadrature index for a given number of quadrature * points. If not in hp mode or if the index is not found, this @@ -281,7 +293,7 @@ namespace internal template void print_memory_consumption(StreamType & out, - const SizeInfo &task_info) const; + const TaskInfo &task_info) const; /** * Returns the memory consumption in bytes. @@ -324,6 +336,20 @@ namespace internal const UpdateFlags update_flags_inner_faces, const UpdateFlags update_flags_faces_by_cells); + /** + * Update the information in the given cells and faces that is the + * result of a change in the given `mapping` class, keeping the cells, + * quadrature formulas and other unknowns unchanged. This call is only + * valid if MappingInfo::initialize() has been called before. + */ + void + update_mapping( + const dealii::Triangulation & tria, + const std::vector> &cells, + const FaceInfo & faces, + const std::vector &active_fe_index, + const Mapping & mapping); + /** * Return the type of a given cell as detected during initialization. */ @@ -351,6 +377,29 @@ namespace internal print_memory_consumption(StreamType & out, const TaskInfo &task_info) const; + /** + * The given update flags for computing the geometry on the cells. + */ + UpdateFlags update_flags_cells; + + /** + * The given update flags for computing the geometry on the boundary + * faces. + */ + UpdateFlags update_flags_boundary_faces; + + /** + * The given update flags for computing the geometry on the interior + * faces. + */ + UpdateFlags update_flags_inner_faces; + + /** + * The given update flags for computing the geometry on the faces for + * cell-centric loops. + */ + UpdateFlags update_flags_faces_by_cells; + /** * Stores whether a cell is Cartesian (cell type 0), has constant * transform data (Jacobians) (cell type 1), or is general (cell type @@ -400,10 +449,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> &quad, - const UpdateFlags update_flags_cells); + const std::vector &active_fe_index, + const Mapping & mapping); /** * Computes the information in the given faces, called within @@ -415,10 +462,7 @@ namespace internal const std::vector> &cells, const std::vector< FaceToCellTopology> &faces, - const Mapping & mapping, - const std::vector> & quad, - const UpdateFlags update_flags_boundary_faces, - const UpdateFlags update_flags_inner_faces); + const Mapping & mapping); /** * Computes the information in the given faces, called within @@ -428,9 +472,7 @@ namespace internal initialize_faces_by_cells( const dealii::Triangulation & tria, const std::vector> &cells, - const Mapping & mapping, - const std::vector> & quad, - const UpdateFlags update_flags_faces_by_cells); + const Mapping & 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 4074240891..6b192976e4 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -63,8 +63,9 @@ namespace internal Assert(structdim + 1 <= spacedim || update_flags_inner_faces == update_default, ExcMessage("Volume cells do not allow for setting inner faces")); - quadrature = Quadrature(quadrature_1d); - n_q_points = quadrature.size(); + this->quadrature_1d = quadrature_1d; + quadrature = Quadrature(quadrature_1d); + n_q_points = quadrature.size(); quadrature_weights.resize(n_q_points); for (unsigned int i = 0; i < n_q_points; ++i) quadrature_weights[i] = quadrature.weight(i); @@ -125,6 +126,29 @@ namespace internal + template + void + MappingInfoStorage:: + clear_data_fields() + { + data_index_offsets.clear(); + JxW_values.clear(); + normal_vectors.clear(); + for (unsigned int i = 0; i < 2; ++i) + { + jacobians[i].clear(); + jacobian_gradients[i].clear(); + normals_times_jacobians[i].clear(); + } + quadrature_point_offsets.clear(); + quadrature_points.clear(); + } + + + template void MappingInfoStorage:: - print_memory_consumption(StreamType &out, const SizeInfo &task_info) const + print_memory_consumption(StreamType &out, const TaskInfo &task_info) const { // print_memory_statistics involves global communication, so we can // disable the check here only if no processor has any such data @@ -287,19 +311,79 @@ namespace internal clear(); this->mapping = &mapping; + cell_data.resize(quad.size()); + face_data.resize(quad.size()); + face_data_by_cells.resize(quad.size()); + + // dummy FE that is used to set up an FEValues object. Do not need the + // actual finite element because we will only evaluate quantities for + // the mapping that are independent of the FE + this->update_flags_cells = compute_update_flags(update_flags_cells, quad); + + this->update_flags_boundary_faces = + ((update_flags_inner_faces | update_flags_boundary_faces) & + update_quadrature_points ? + update_quadrature_points : + update_default) | + update_normal_vectors | update_JxW_values | update_jacobians; + this->update_flags_inner_faces = this->update_flags_boundary_faces; + this->update_flags_faces_by_cells = update_flags_faces_by_cells; + + for (unsigned int my_q = 0; my_q < quad.size(); ++my_q) + { + const unsigned int n_hp_quads = quad[my_q].size(); + 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); + + 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); + + 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); + } + // 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, quad, update_flags_cells); - initialize_faces(tria, - cells, - face_info.faces, - mapping, - quad, - update_flags_boundary_faces, - update_flags_inner_faces); - initialize_faces_by_cells( - tria, cells, mapping, quad, update_flags_faces_by_cells); + initialize_cells(tria, cells, active_fe_index, mapping); + initialize_faces(tria, cells, face_info.faces, mapping); + initialize_faces_by_cells(tria, cells, mapping); + } + + + + template + void + MappingInfo::update_mapping( + const dealii::Triangulation & tria, + const std::vector> &cells, + const FaceInfo & face_info, + const std::vector & active_fe_index, + const Mapping & mapping) + { + AssertDimension(cells.size() / VectorizedArrayType::n_array_elements, + cell_type.size()); + + for (auto &data : cell_data) + data.clear_data_fields(); + for (auto &data : face_data) + data.clear_data_fields(); + for (auto &data : face_data_by_cells) + data.clear_data_fields(); + + this->mapping = &mapping; + + // 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); } @@ -583,8 +667,6 @@ namespace internal const std::vector> &cells, const std::vector & active_fe_index, const Mapping & mapping, - const std::vector> &quad, - const UpdateFlags update_flags, MappingInfo &mapping_info, std::pair>, @@ -626,7 +708,8 @@ namespace internal fe_values(mapping_info.cell_data.size()); for (unsigned int i = 0; i < fe_values.size(); ++i) fe_values[i].resize(mapping_info.cell_data[i].descriptor.size()); - UpdateFlags update_flags_feval = + const UpdateFlags update_flags = mapping_info.update_flags_cells; + const UpdateFlags update_flags_feval = (update_flags & update_jacobians ? update_jacobians : update_default) | (update_flags & update_jacobian_grads ? update_jacobian_grads : @@ -634,15 +717,18 @@ namespace internal (update_flags & update_quadrature_points ? update_quadrature_points : update_default); - std::vector> n_q_points_1d(quad.size()), - step_size_cartesian(quad.size()); - for (unsigned int my_q = 0; my_q < quad.size(); ++my_q) + std::vector> n_q_points_1d(fe_values.size()), + step_size_cartesian(fe_values.size()); + for (unsigned int my_q = 0; my_q < fe_values.size(); ++my_q) { - n_q_points_1d[my_q].resize(quad[my_q].size()); - step_size_cartesian[my_q].resize(quad[my_q].size()); - for (unsigned int hpq = 0; hpq < quad[my_q].size(); ++hpq) + n_q_points_1d[my_q].resize( + mapping_info.cell_data[my_q].descriptor.size()); + step_size_cartesian[my_q].resize(n_q_points_1d[my_q].size()); + for (unsigned int hpq = 0; hpq < n_q_points_1d[my_q].size(); ++hpq) { - n_q_points_1d[my_q][hpq] = quad[my_q][hpq].size(); + n_q_points_1d[my_q][hpq] = mapping_info.cell_data[my_q] + .descriptor[hpq] + .quadrature_1d.size(); // To walk on the diagonal for lexicographic ordering, we have // to jump one index ahead in each direction. For direction 0, @@ -1024,34 +1110,15 @@ namespace internal const dealii::Triangulation & tria, const std::vector> &cells, const std::vector & active_fe_index, - const Mapping & mapping, - const std::vector> & quad, - const UpdateFlags update_flags_input) + const Mapping & mapping) { - const unsigned int n_quads = quad.size(); const unsigned int n_cells = cells.size(); const unsigned int vectorization_width = VectorizedArrayType::n_array_elements; Assert(n_cells % vectorization_width == 0, ExcInternalError()); const unsigned int n_macro_cells = n_cells / vectorization_width; - cell_data.resize(n_quads); cell_type.resize(n_macro_cells); - // dummy FE that is used to set up an FEValues object. Do not need the - // actual finite element because we will only evaluate quantities for - // the mapping that are independent of the FE - UpdateFlags update_flags = compute_update_flags(update_flags_input, quad); - - for (unsigned int my_q = 0; my_q < n_quads; ++my_q) - { - const unsigned int n_hp_quads = quad[my_q].size(); - 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); - } - if (n_macro_cells == 0) return; @@ -1078,7 +1145,7 @@ namespace internal data_cells_local.push_back(std::make_pair( std::vector< MappingInfoStorage>( - n_quads), + cell_data.size()), ExtractCellHelper:: CompressedCellData( ExtractCellHelper::get_jacobian_size(tria)))); @@ -1090,8 +1157,6 @@ namespace internal cells, active_fe_index, mapping, - quad, - update_flags, *this, data_cells_local.back()); cell_range.first = cell_range.second; @@ -1135,10 +1200,10 @@ namespace internal data_cells_local.back().first[my_q].JxW_values.size()); cell_data[my_q].jacobians[0].resize_fast( cell_data[my_q].JxW_values.size()); - if (update_flags & update_jacobian_grads) + if (update_flags_cells & update_jacobian_grads) cell_data[my_q].jacobian_gradients[0].resize_fast( cell_data[my_q].JxW_values.size()); - if (update_flags & update_quadrature_points) + if (update_flags_cells & update_quadrature_points) { cell_data[my_q].quadrature_point_offsets.resize(cell_type.size()); cell_data[my_q].quadrature_points.resize_fast( @@ -1276,8 +1341,6 @@ namespace internal const std::vector< FaceToCellTopology> &faces, const Mapping & mapping, - const UpdateFlags update_flags_boundary, - const UpdateFlags update_flags_inner, MappingInfo &mapping_info, std::pair< std::vector< @@ -1328,16 +1391,18 @@ namespace internal if (is_boundary_face && fe_boundary_face_values_container[my_q][0] == nullptr) fe_boundary_face_values_container[my_q][0] = - std::make_shared>(mapping, - dummy_fe, - quadrature, - update_flags_boundary); + std::make_shared>( + mapping, + dummy_fe, + quadrature, + mapping_info.update_flags_boundary_faces); else if (fe_face_values_container[my_q][0] == nullptr) fe_face_values_container[my_q][0] = - std::make_shared>(mapping, - dummy_fe, - quadrature, - update_flags_inner); + std::make_shared>( + mapping, + dummy_fe, + quadrature, + mapping_info.update_flags_inner_faces); FEFaceValues &fe_face_values = is_boundary_face ? *fe_boundary_face_values_container[my_q][0] : @@ -1489,7 +1554,7 @@ namespace internal mapping, dummy_fe, quadrature, - update_flags_inner); + mapping_info.update_flags_inner_faces); fe_subface_values_container[my_q][0]->reinit( cell_it, faces[face].exterior_face_no, @@ -1700,39 +1765,9 @@ namespace internal const std::vector> &cells, const std::vector< FaceToCellTopology> &faces, - const Mapping & mapping, - const std::vector> & quad, - const UpdateFlags update_flags_boundary_faces, - const UpdateFlags update_flags_inner_faces) + const Mapping & mapping) { face_type.resize(faces.size(), general); - face_data.resize(quad.size()); - - // We currently always set the same flags on both inner and boundary - // faces to the same value. At some point, we might want to separate the - // two. - UpdateFlags update_flags_compute_boundary = - ((update_flags_inner_faces | update_flags_boundary_faces) & - update_quadrature_points ? - update_quadrature_points : - update_default) | - update_normal_vectors | update_JxW_values | update_jacobians; - UpdateFlags update_flags_compute_inner = - ((update_flags_inner_faces | update_flags_boundary_faces) & - update_quadrature_points ? - update_quadrature_points : - update_default) | - update_normal_vectors | update_JxW_values | update_jacobians; - UpdateFlags update_flags_common = - update_flags_inner_faces | update_flags_boundary_faces; - - for (unsigned int my_q = 0; my_q < quad.size(); ++my_q) - { - face_data[my_q].descriptor.resize(quad[my_q].size()); - for (unsigned int hpq = 0; hpq < quad[my_q].size(); ++hpq) - face_data[my_q].descriptor[hpq].initialize( - quad[my_q][hpq], update_flags_compute_inner); - } if (faces.size() == 0) return; @@ -1761,7 +1796,7 @@ namespace internal data_faces_local.push_back(std::make_pair( std::vector< MappingInfoStorage>( - quad.size()), + face_data.size()), ExtractFaceHelper:: CompressedFaceData( ExtractCellHelper::get_jacobian_size(tria)))); @@ -1773,8 +1808,6 @@ namespace internal cells, faces, mapping, - update_flags_compute_boundary, - update_flags_compute_inner, *this, data_faces_local.back()); face_range.first = face_range.second; @@ -1794,6 +1827,9 @@ namespace internal data_faces_local[0].second.data, indices_compressed[i]); + const UpdateFlags update_flags_common = + update_flags_boundary_faces | update_flags_inner_faces; + // Collect all data in the final data fields. // First allocate the memory const unsigned int n_constant_data = @@ -1926,15 +1962,12 @@ namespace internal MappingInfo::initialize_faces_by_cells( const dealii::Triangulation & tria, const std::vector> &cells, - const Mapping & mapping, - const std::vector> & quad, - const UpdateFlags update_flags_faces_by_cells) + const Mapping & mapping) { if (update_flags_faces_by_cells == update_default) return; - face_data_by_cells.resize(quad.size()); - const unsigned int n_quads = quad.size(); + const unsigned int n_quads = face_data_by_cells.size(); const unsigned int vectorization_width = VectorizedArrayType::n_array_elements; UpdateFlags update_flags = @@ -1945,13 +1978,6 @@ namespace internal for (unsigned int my_q = 0; my_q < n_quads; ++my_q) { - const unsigned int n_hp_quads = quad[my_q].size(); - AssertIndexRange(0, n_hp_quads); - face_data_by_cells[my_q].descriptor.resize(n_hp_quads); - for (unsigned int q = 0; q < n_hp_quads; ++q) - face_data_by_cells[my_q].descriptor[q].initialize(quad[my_q][q], - update_default); - // since we already know the cell type, we can pre-allocate the right // amount of data straight away and we just need to do some basic // counting @@ -2217,7 +2243,7 @@ namespace internal void MappingInfo::print_memory_consumption( StreamType & out, - const SizeInfo &task_info) const + const TaskInfo &task_info) const { out << " Cell types: "; task_info.print_memory_statistics(out, diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 107f680371..32b8ca42b4 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -684,6 +684,18 @@ public: copy_from( const MatrixFree &matrix_free_base); + /** + * Refreshes the geometry data stored in the MappingInfo fields when the + * underlying geometry has changed (e.g. by a mapping that can deform + * through a change in the spatial configuration like MappingFEField) + * whereas the topology of the mesh and unknowns have remained the + * same. Compared to reinit(), this operation only has to re-generate the + * geometry arrays and can thus be significantly cheaper (depending on the + * cost to evaluate the geometry). + */ + void + update_mapping(const Mapping &mapping); + /** * Clear all data fields and brings the class into a condition similar to * after having called the default constructor. diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 7cd7dd18c0..21b77dbe1a 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -681,6 +681,25 @@ MatrixFree::internal_reinit( } + +template +void +MatrixFree::update_mapping( + const Mapping &mapping) +{ + AssertDimension(shape_info.size(1), mapping_info.cell_data.size()); + mapping_info.update_mapping( + dof_handlers.active_dof_handler == DoFHandlers::hp ? + dof_handlers.hp_dof_handler[0]->get_triangulation() : + dof_handlers.dof_handler[0]->get_triangulation(), + cell_level_index, + face_info, + dof_info[0].cell_active_fe_index, + mapping); +} + + + template template bool