From a3b7cea558105a1d139b606573468e2089d36ea2 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 4 Mar 2021 18:02:53 +0100 Subject: [PATCH] MatrixFree: introduce FaceToCellTopology::face_type --- .../deal.II/matrix_free/evaluation_kernels.h | 4 +- include/deal.II/matrix_free/face_info.h | 6 + .../deal.II/matrix_free/face_setup_internal.h | 15 ++ include/deal.II/matrix_free/fe_evaluation.h | 74 +++++--- .../matrix_free/mapping_info.templates.h | 2 - .../matrix_free/matrix_free.templates.h | 171 ++++++++++++------ .../matrix_free/shape_info.templates.h | 7 +- 7 files changed, 189 insertions(+), 90 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 12abbb4ccd..a9ea830a82 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -2549,7 +2549,7 @@ namespace internal if (data.element_type == MatrixFreeFunctions::tensor_none) { const unsigned int n_dofs = data.dofs_per_component_on_cell; - const unsigned int n_q_points = data.n_q_points_face; + const unsigned int n_q_points = data.n_q_points_faces[face_no]; const auto shape_info = data.data.front(); using Eval = EvaluatorTensorProductface(f)->boundary_id(); + info.face_type = + dcell->face(f)->reference_cell() != + dealii::ReferenceCells::get_hypercube(); info.subface_index = GeometryInfo::max_children_per_cell; info.face_orientation = 0; @@ -959,6 +962,9 @@ namespace internal else info.exterior_face_no = cell->neighbor_face_no(face_no); + info.face_type = cell->face(face_no)->reference_cell() != + dealii::ReferenceCells::get_hypercube(); + info.subface_index = GeometryInfo::max_children_per_cell; Assert(neighbor->level() <= cell->level(), ExcInternalError()); if (cell->level() > neighbor->level()) @@ -1037,6 +1043,8 @@ namespace internal return false; if (face1.face_orientation != face2.face_orientation) return false; + if (face1.face_type != face2.face_type) + return false; if (active_fe_indices.size() > 0) { @@ -1072,6 +1080,12 @@ namespace internal operator()(const FaceToCellTopology &face1, const FaceToCellTopology &face2) const { + // check if active FE indices match + if (face1.face_type < face2.face_type) + return true; + else if (face1.face_type > face2.face_type) + return false; + // check if active FE indices match if (active_fe_indices.size() > 0) { @@ -1187,6 +1201,7 @@ namespace internal new_faces(face_comparator); for (const auto &face_type : faces_type) { + macro_face.face_type = faces_in[face_type[0]].face_type; macro_face.interior_face_no = faces_in[face_type[0]].interior_face_no; macro_face.exterior_face_no = diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 45de596e71..1dd0f581d6 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -298,7 +298,8 @@ protected: const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, - const unsigned int active_quad_index); + const unsigned int active_quad_index, + const unsigned int face_type); /** * Constructor that comes with reduced functionality and works similar as @@ -1143,7 +1144,8 @@ protected: const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, - const unsigned int active_quad_index); + const unsigned int active_quad_index, + const unsigned int face_type); /** * Constructor that comes with reduced functionality and works similar as @@ -1428,7 +1430,8 @@ protected: const unsigned int n_q_points, const bool is_interior_face = true, const unsigned int active_fe_index = numbers::invalid_unsigned_int, - const unsigned int active_quad_index = numbers::invalid_unsigned_int); + const unsigned int active_quad_index = numbers::invalid_unsigned_int, + const unsigned int face_type = numbers::invalid_unsigned_int); /** * Constructor with reduced functionality for similar usage of FEEvaluation @@ -1565,7 +1568,8 @@ protected: const unsigned int n_q_points, const bool is_interior_face = true, const unsigned int active_fe_index = numbers::invalid_unsigned_int, - const unsigned int active_quad_index = numbers::invalid_unsigned_int); + const unsigned int active_quad_index = numbers::invalid_unsigned_int, + const unsigned int face_type = numbers::invalid_unsigned_int); /** * Constructor with reduced functionality for similar usage of FEEvaluation @@ -1726,7 +1730,8 @@ protected: const unsigned int n_q_points, const bool is_interior_face = true, const unsigned int active_fe_index = numbers::invalid_unsigned_int, - const unsigned int active_quad_index = numbers::invalid_unsigned_int); + const unsigned int active_quad_index = numbers::invalid_unsigned_int, + const unsigned int face_type = numbers::invalid_unsigned_int); /** * Constructor with reduced functionality for similar usage of FEEvaluation @@ -1872,7 +1877,8 @@ protected: const unsigned int n_q_points, const bool is_interior_face = true, const unsigned int active_fe_index = numbers::invalid_unsigned_int, - const unsigned int active_quad_index = numbers::invalid_unsigned_int); + const unsigned int active_quad_index = numbers::invalid_unsigned_int, + const unsigned int face_type = numbers::invalid_unsigned_int); /** * Constructor with reduced functionality for similar usage of FEEvaluation @@ -3069,6 +3075,9 @@ public: * DoFHandler/AffineConstraints pair the given evaluator should be attached * to. * + * @param face_type In the case of a face, indicate its reference-cell type + * (0 for line or quadrilateral 1 for triangle). + * * @param active_quad_index If matrix_free was set up with hp::Collection * objects, this parameter selects the appropriate number of the quadrature * formula. @@ -3080,7 +3089,8 @@ public: const unsigned int quad_no = 0, const unsigned int first_selected_component = 0, const unsigned int active_fe_index = numbers::invalid_unsigned_int, - const unsigned int active_quad_index = numbers::invalid_unsigned_int); + const unsigned int active_quad_index = numbers::invalid_unsigned_int, + const unsigned int face_type = numbers::invalid_unsigned_int); /** * Constructor. Takes all data stored in MatrixFree for a given face range, @@ -3342,7 +3352,8 @@ inline FEEvaluationBaseData:: const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index_in, - const unsigned int active_quad_index_in) + const unsigned int active_quad_index_in, + const unsigned int face_type) : scratch_data_array(data_in.acquire_scratch_data()) , quad_no(quad_no_in) , matrix_info(&data_in) @@ -3366,7 +3377,12 @@ inline FEEvaluationBaseData:: active_quad_index_in : std::min(active_fe_index, mapping_data->descriptor.size() - 1))) - , descriptor(&mapping_data->descriptor[active_quad_index]) + , descriptor( + &mapping_data->descriptor + [is_face ? + (active_quad_index * std::max(1, dim - 1) + + (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) : + active_quad_index]) , n_quadrature_points(descriptor->n_q_points) , data(&data_in.get_shape_info( dof_no, @@ -4019,7 +4035,8 @@ inline FEEvaluationBase( data_in, dof_no, @@ -4029,7 +4046,8 @@ inline FEEvaluationBase( data_in, dof_no, @@ -6169,7 +6188,8 @@ inline FEEvaluationAccess:: const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, - const unsigned int active_quad_index) + const unsigned int active_quad_index, + const unsigned int face_type) : FEEvaluationBase( data_in, dof_no, @@ -6275,7 +6296,8 @@ inline FEEvaluationAccess:: n_q_points, is_interior_face, active_fe_index, - active_quad_index) + active_quad_index, + face_type) {} @@ -6576,7 +6598,8 @@ inline FEEvaluationAccess:: const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, - const unsigned int active_quad_index) + const unsigned int active_quad_index, + const unsigned int face_type) : FEEvaluationBase( data_in, dof_no, @@ -6586,7 +6609,8 @@ inline FEEvaluationAccess:: n_q_points, is_interior_face, active_fe_index, - active_quad_index) + active_quad_index, + face_type) {} @@ -6985,7 +7009,8 @@ inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>:: const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, - const unsigned int active_quad_index) + const unsigned int active_quad_index, + const unsigned int face_type) : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>( data_in, dof_no, @@ -6995,7 +7020,8 @@ inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>:: n_q_points, is_interior_face, active_fe_index, - active_quad_index) + active_quad_index, + face_type) {} @@ -8412,7 +8438,8 @@ inline FEFaceEvaluationdata->dofs_per_component_on_cell) , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_) , n_q_points(this->n_quadrature_points) @@ -8454,7 +8482,9 @@ inline FEFaceEvaluation &fe_indices, - const bool include) + const bool include, + const bool only_face_type) : fe_indices(fe_indices) , include(include) + , only_face_type(only_face_type) {} template bool operator()(const internal::MatrixFreeFunctions::FaceToCellTopology< - vectorization_width> &face, - const unsigned int & fe_index) + vectorization_width> & face, + const std::array &fe_index) { - const unsigned int fe_index_face = - fe_indices[face.cells_interior[0] / vectorization_width]; + const unsigned int face_type = face.face_type; - return fe_index_face < fe_index || - (include ? (fe_index_face == fe_index) : false); + if (only_face_type) + return include ? (face_type <= fe_index[0]) : (face_type < fe_index[0]); + + const std::array fe_index_face = { + {face_type, fe_indices[face.cells_interior[0] / vectorization_width]}}; + + return include ? (fe_index_face <= fe_index) : (fe_index_face < fe_index); } template bool operator()(const internal::MatrixFreeFunctions::FaceToCellTopology< - vectorization_width> & face, - const std::pair &fe_index) + vectorization_width> & face, + const std::array &fe_index) { - const std::pair fe_index_face = { - fe_indices[face.cells_interior[0] / vectorization_width], - fe_indices[face.cells_exterior[0] / vectorization_width]}; + const unsigned int face_type = face.face_type; + + if (only_face_type) + return include ? (face_type <= fe_index[0]) : (face_type < fe_index[0]); + + const std::array fe_index_face = { + {face_type, + fe_indices[face.cells_interior[0] / vectorization_width], + fe_indices[face.cells_exterior[0] / vectorization_width]}}; return include ? (fe_index_face <= fe_index) : (fe_index_face < fe_index); } @@ -168,6 +180,7 @@ namespace private: const std::vector &fe_indices; const bool include; + const bool only_face_type; }; } // namespace @@ -540,13 +553,18 @@ MatrixFree::internal_reinit( const auto create_inner_face_subrange_hp_by_index = [&](const std::pair &range, + const unsigned int face_type, const unsigned int fe_index_interior, const unsigned int fe_index_exterior, const unsigned int dof_handler_index = 0) -> std::pair { - if (dof_info[dof_handler_index].max_fe_index == 0 || - dof_handlers[dof_handler_index]->get_fe_collection().size() == - 1) + const unsigned int n_face_types = + std::max(dim - 1, 1); + + if ((dof_info[dof_handler_index].max_fe_index == 0 || + dof_handlers[dof_handler_index]->get_fe_collection().size() == + 1) && + n_face_types == 1) return range; AssertIndexRange(fe_index_interior, @@ -555,24 +573,32 @@ MatrixFree::internal_reinit( dof_info[dof_handler_index].max_fe_index); const std::vector &fe_indices = dof_info[dof_handler_index].cell_active_fe_index; - if (fe_indices.empty() == true) + if (fe_indices.empty() == true && n_face_types == 1) return range; else { + const bool only_face_type = + dof_info[dof_handler_index].max_fe_index == 0 || + dof_handlers[dof_handler_index]->get_fe_collection().size() == + 1 || + fe_indices.empty() == true; + std::pair return_range; return_range.first = - std::lower_bound(face_info.faces.begin() + range.first, - face_info.faces.begin() + range.second, - std::pair{ - fe_index_interior, fe_index_exterior}, - FaceRangeCompartor(fe_indices, false)) - + std::lower_bound( + face_info.faces.begin() + range.first, + face_info.faces.begin() + range.second, + std::array{ + {face_type, fe_index_interior, fe_index_exterior}}, + FaceRangeCompartor(fe_indices, false, only_face_type)) - face_info.faces.begin(); return_range.second = - std::lower_bound(face_info.faces.begin() + return_range.first, - face_info.faces.begin() + range.second, - std::pair{ - fe_index_interior, fe_index_exterior}, - FaceRangeCompartor(fe_indices, true)) - + std::lower_bound( + face_info.faces.begin() + return_range.first, + face_info.faces.begin() + range.second, + std::array{ + {face_type, fe_index_interior, fe_index_exterior}}, + FaceRangeCompartor(fe_indices, true, only_face_type)) - face_info.faces.begin(); Assert(return_range.first >= range.first && return_range.second <= range.second, @@ -592,20 +618,27 @@ MatrixFree::internal_reinit( task_info.face_partition_data[i + 1]}; if (range.second > range.first) - for (unsigned int i = 0; i < this->n_active_fe_indices(); - ++i) - for (unsigned int j = 0; j < this->n_active_fe_indices(); - ++j) - { - const auto subrange = - create_inner_face_subrange_hp_by_index(range, i, j); + for (unsigned int t = 0; + t < std::max(dim - 1, 1); + ++t) + for (unsigned int i = 0; i < this->n_active_fe_indices(); + ++i) + for (unsigned int j = 0; + j < this->n_active_fe_indices(); + ++j) + { + const auto subrange = + create_inner_face_subrange_hp_by_index(range, + t, + i, + j); - if (subrange.second <= subrange.first) - continue; + if (subrange.second <= subrange.first) + continue; - data.push_back(subrange.first); - data.push_back(subrange.second); - } + data.push_back(subrange.first); + data.push_back(subrange.second); + } } ptr.push_back(data.size() / 2); @@ -624,34 +657,47 @@ MatrixFree::internal_reinit( const auto create_boundary_face_subrange_hp_by_index = [&](const std::pair &range, + const unsigned int face_type, const unsigned int fe_index, const unsigned int dof_handler_index = 0) -> std::pair { - if (dof_info[dof_handler_index].max_fe_index == 0 || - dof_handlers[dof_handler_index]->get_fe_collection().size() == - 1) + const unsigned int n_face_types = + std::max(dim - 1, 1); + + if ((dof_info[dof_handler_index].max_fe_index == 0 || + dof_handlers[dof_handler_index]->get_fe_collection().size() == + 1) && + n_face_types == 1) return range; AssertIndexRange(fe_index, dof_info[dof_handler_index].max_fe_index); const std::vector &fe_indices = dof_info[dof_handler_index].cell_active_fe_index; - if (fe_indices.empty() == true) + if (fe_indices.empty() == true && n_face_types == 1) return range; else { + const bool only_face_type = + dof_info[dof_handler_index].max_fe_index == 0 || + dof_handlers[dof_handler_index]->get_fe_collection().size() == + 1 || + fe_indices.empty() == true; + std::pair return_range; return_range.first = - std::lower_bound(face_info.faces.begin() + range.first, - face_info.faces.begin() + range.second, - fe_index, - FaceRangeCompartor(fe_indices, false)) - + std::lower_bound( + face_info.faces.begin() + range.first, + face_info.faces.begin() + range.second, + std::array{{face_type, fe_index}}, + FaceRangeCompartor(fe_indices, false, only_face_type)) - face_info.faces.begin(); return_range.second = - std::lower_bound(face_info.faces.begin() + return_range.first, - face_info.faces.begin() + range.second, - fe_index, - FaceRangeCompartor(fe_indices, true)) - + std::lower_bound( + face_info.faces.begin() + return_range.first, + face_info.faces.begin() + range.second, + std::array{{face_type, fe_index}}, + FaceRangeCompartor(fe_indices, true, only_face_type)) - face_info.faces.begin(); Assert(return_range.first >= range.first && return_range.second <= range.second, @@ -672,18 +718,23 @@ MatrixFree::internal_reinit( task_info.boundary_partition_data[i + 1]}; if (range.second > range.first) - for (unsigned int i = 0; i < this->n_active_fe_indices(); - ++i) - { - const auto cell_subrange = - create_boundary_face_subrange_hp_by_index(range, i); + for (unsigned int t = 0; + t < std::max(dim - 1, 1); + ++t) + for (unsigned int i = 0; i < this->n_active_fe_indices(); + ++i) + { + const auto cell_subrange = + create_boundary_face_subrange_hp_by_index(range, + t, + i); - if (cell_subrange.second <= cell_subrange.first) - continue; + if (cell_subrange.second <= cell_subrange.first) + continue; - data.push_back(cell_subrange.first); - data.push_back(cell_subrange.second); - } + data.push_back(cell_subrange.first); + data.push_back(cell_subrange.second); + } } ptr.push_back(data.size() / 2); diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index c3293a0f12..c489a86041 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -93,10 +93,9 @@ namespace internal const auto fe_poly = dynamic_cast *>(&fe); - if (dynamic_cast *>(&fe) != nullptr || - dynamic_cast *>(&fe) != nullptr || - dynamic_cast *>(&fe) != nullptr || - dynamic_cast *>(&fe) != nullptr) + if (dynamic_cast *>(&fe) != nullptr || + dynamic_cast *>(&fe) != nullptr || + dynamic_cast *>(&fe) != nullptr) { scalar_lexicographic.resize(fe.n_dofs_per_cell()); for (unsigned int i = 0; i < scalar_lexicographic.size(); ++i) -- 2.39.5