From da5dffb99b959ab994003d5c8f9aef1d0d6caa43 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 25 Nov 2020 15:25:14 +0100 Subject: [PATCH] Make MappingInfo::initialize_faces work for multiple fe indices --- include/deal.II/matrix_free/mapping_info.h | 1 + .../matrix_free/mapping_info.templates.h | 115 ++++++++++++++---- 2 files changed, 90 insertions(+), 26 deletions(-) diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h index 941d530556..dc63926076 100644 --- a/include/deal.II/matrix_free/mapping_info.h +++ b/include/deal.II/matrix_free/mapping_info.h @@ -503,6 +503,7 @@ namespace internal const std::vector> &cells, const std::vector> & faces, + const std::vector & active_fe_index, const dealii::hp::MappingCollection &mapping); /** diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 2a669655d5..6bea2049f1 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -450,7 +450,8 @@ namespace internal // 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( + tria, cells, face_info.faces, active_fe_index, *mapping); initialize_faces_by_cells(tria, cells, *mapping); } } @@ -487,7 +488,8 @@ namespace internal // 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( + tria, cells, face_info.faces, active_fe_index, *mapping); initialize_faces_by_cells(tria, cells, *mapping); } } @@ -889,14 +891,13 @@ namespace internal std::vector>>> fe_values(mapping_info.cell_data.size()); - const unsigned int n_active_fe_indices = + const unsigned int max_active_fe_index = active_fe_index.size() > 0 ? - (*std::max_element(active_fe_index.begin(), active_fe_index.end()) + - 1) : - 1; + *std::max_element(active_fe_index.begin(), active_fe_index.end()) : + 0; for (unsigned int i = 0; i < fe_values.size(); ++i) - fe_values[i].resize(n_active_fe_indices); + fe_values[i].resize(max_active_fe_index + 1); const UpdateFlags update_flags = mapping_info.update_flags_cells; const UpdateFlags update_flags_feval = (update_flags & update_jacobians ? update_jacobians : @@ -1759,6 +1760,7 @@ namespace internal const std::vector> &cells, const std::vector> & faces, + const std::vector & active_fe_index, const dealii::hp::MappingCollection & mapping_in, MappingInfo &mapping_info, std::pair< @@ -1766,24 +1768,25 @@ namespace internal MappingInfoStorage>, CompressedFaceData> &data) { - AssertDimension(mapping_in.size(), 1); - const auto &mapping = mapping_in[0]; - FE_Nothing dummy_fe; + const unsigned int max_active_fe_index = + active_fe_index.size() > 0 ? + *std::max_element(active_fe_index.begin(), active_fe_index.end()) : + 0; + std::vector>>> fe_face_values_container(mapping_info.face_data.size()); for (unsigned int my_q = 0; my_q < mapping_info.face_data.size(); ++my_q) - fe_face_values_container[my_q].resize( - mapping_info.face_data[my_q].descriptor.size()); + fe_face_values_container[my_q].resize(max_active_fe_index + 1); std::vector>>> fe_boundary_face_values_container(mapping_info.face_data.size()); for (unsigned int my_q = 0; my_q < mapping_info.face_data.size(); ++my_q) - fe_boundary_face_values_container[my_q].resize( - mapping_info.face_data[my_q].descriptor.size()); + fe_boundary_face_values_container[my_q].resize(max_active_fe_index + + 1); std::vector>>> fe_subface_values_container(mapping_info.face_data.size()); @@ -1801,25 +1804,47 @@ namespace internal for (unsigned int my_q = 0; my_q < mapping_info.face_data.size(); ++my_q) { +#ifndef DEAL_II_WITH_SIMPLEX_SUPPORT // currently only non-hp case... - Assert(mapping_info.face_data[my_q].descriptor.size() == 1, - ExcNotImplemented()); + AssertDimension(mapping_in.size(), 1); + AssertDimension(mapping_info.face_data[my_q].descriptor.size(), + 1); +#endif + + // We assume that we have the faces sorted by the active fe + // indices so that the active fe index of the interior side of the + // face batch is the same as the fe index of the interior side of + // its first entry. + const unsigned int fe_index = + active_fe_index.size() > 0 ? + active_fe_index[faces[face].cells_interior[0] / + VectorizedArrayType::size()] : + 0; + const unsigned int hp_quad_index = + mapping_info.cell_data[my_q].descriptor.size() == 1 ? 0 : + fe_index; + const unsigned int hp_mapping_index = + mapping_in.size() == 1 ? 0 : fe_index; + + const auto & mapping = mapping_in[hp_mapping_index]; const Quadrature &quadrature = - mapping_info.face_data[my_q].descriptor[0].quadrature; + mapping_info.face_data[my_q] + .descriptor[hp_quad_index] + .quadrature; const bool is_boundary_face = faces[face].cells_exterior[0] == numbers::invalid_unsigned_int; if (is_boundary_face && - fe_boundary_face_values_container[my_q][0] == nullptr) - fe_boundary_face_values_container[my_q][0] = + fe_boundary_face_values_container[my_q][fe_index] == nullptr) + fe_boundary_face_values_container[my_q][fe_index] = 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] = + else if (fe_face_values_container[my_q][fe_index] == nullptr) + fe_face_values_container[my_q][fe_index] = std::make_shared>( mapping, dummy_fe, @@ -1827,8 +1852,9 @@ namespace internal mapping_info.update_flags_inner_faces); FEFaceValues &fe_face_values = - is_boundary_face ? *fe_boundary_face_values_container[my_q][0] : - *fe_face_values_container[my_q][0]; + is_boundary_face ? + *fe_boundary_face_values_container[my_q][fe_index] : + *fe_face_values_container[my_q][fe_index]; const unsigned int n_q_points = fe_face_values.n_quadrature_points; face_data.resize(n_q_points); @@ -1962,9 +1988,44 @@ namespace internal if (faces[face].subface_index >= GeometryInfo::max_children_per_cell) { - fe_face_values.reinit(cell_it, - faces[face].exterior_face_no); - actual_fe_face_values = &fe_face_values; + // We assume that we have the faces sorted by the + // active fe indices so that the active fe index of + // the exterior side of the face batch is the same as + // the fe index of the exterior side of its first + // entry. + const unsigned int fe_index = + active_fe_index.size() > 0 ? + active_fe_index[faces[face].cells_exterior[0] / + VectorizedArrayType::size()] : + 0; + const unsigned int hp_quad_index = + mapping_info.cell_data[my_q].descriptor.size() == + 1 ? + 0 : + fe_index; + const unsigned int hp_mapping_index = + mapping_in.size() == 1 ? 0 : fe_index; + + const auto &mapping = mapping_in[hp_mapping_index]; + const Quadrature &quadrature = + mapping_info.face_data[my_q] + .descriptor[hp_quad_index] + .quadrature; + + if (fe_face_values_container[my_q][fe_index] == + nullptr) + fe_face_values_container[my_q][fe_index] = + std::make_shared>( + mapping, + dummy_fe, + quadrature, + mapping_info.update_flags_boundary_faces); + + fe_face_values_container[my_q][fe_index]->reinit( + cell_it, faces[face].exterior_face_no); + + actual_fe_face_values = + fe_face_values_container[my_q][fe_index].get(); } else { @@ -2406,6 +2467,7 @@ namespace internal const dealii::Triangulation & tria, const std::vector> & cells, const std::vector> &faces, + const std::vector & active_fe_index, const dealii::hp::MappingCollection &mapping) { face_type.resize(faces.size(), general); @@ -2449,6 +2511,7 @@ namespace internal tria, cells, faces, + active_fe_index, mapping, *this, data_faces_local.back()); -- 2.39.5