From: Peter Munch <peterrmuench@gmail.com> Date: Thu, 4 Mar 2021 17:02:53 +0000 (+0100) Subject: MatrixFree: introduce FaceToCellTopology::face_type X-Git-Tag: v9.3.0-rc1~372^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11401%2Fhead;p=dealii.git MatrixFree: introduce FaceToCellTopology::face_type --- 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 = EvaluatorTensorProduct<evaluate_general, @@ -2704,7 +2704,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 = EvaluatorTensorProduct<evaluate_general, diff --git a/include/deal.II/matrix_free/face_info.h b/include/deal.II/matrix_free/face_info.h index b411c27d7e..2e1e01ef5a 100644 --- a/include/deal.II/matrix_free/face_info.h +++ b/include/deal.II/matrix_free/face_info.h @@ -110,6 +110,12 @@ namespace internal */ unsigned char face_orientation; + /** + * Reference-cell type of the given face: 0 for line or quadrilateral, + * 1 for triangle. + */ + unsigned char face_type; + /** * Return the memory consumption of the present data structure. */ diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h index d480a4db25..e3905e6fad 100644 --- a/include/deal.II/matrix_free/face_setup_internal.h +++ b/include/deal.II/matrix_free/face_setup_internal.h @@ -780,6 +780,9 @@ namespace internal info.cells_exterior[0] = numbers::invalid_unsigned_int; info.interior_face_no = f; info.exterior_face_no = dcell->face(f)->boundary_id(); + info.face_type = + dcell->face(f)->reference_cell() != + dealii::ReferenceCells::get_hypercube<dim - 1>(); info.subface_index = GeometryInfo<dim>::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<dim - 1>(); + info.subface_index = GeometryInfo<dim>::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<length> &face1, const FaceToCellTopology<length> &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<dim, Number, is_face, VectorizedArrayType>:: 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<dim, Number, is_face, VectorizedArrayType>:: active_quad_index_in : std::min<unsigned int>(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<unsigned int>(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<dim, 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) : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>( data_in, dof_no, @@ -4029,7 +4046,8 @@ inline FEEvaluationBase<dim, n_q_points, is_interior_face, active_fe_index, - active_quad_index) + active_quad_index, + face_type) , n_fe_components(data_in.get_dof_info(dof_no).start_components.back()) , dof_values_initialized(false) , values_quad_initialized(false) @@ -6159,7 +6177,8 @@ inline FEEvaluationAccess<dim, 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<dim, n_components_, Number, is_face, VectorizedArrayType>( data_in, dof_no, @@ -6169,7 +6188,8 @@ inline FEEvaluationAccess<dim, n_q_points, is_interior_face, active_fe_index, - active_quad_index) + active_quad_index, + face_type) {} @@ -6265,7 +6285,8 @@ inline FEEvaluationAccess<dim, 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<dim, 1, Number, is_face, VectorizedArrayType>( data_in, dof_no, @@ -6275,7 +6296,8 @@ inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>:: n_q_points, is_interior_face, active_fe_index, - active_quad_index) + active_quad_index, + face_type) {} @@ -6576,7 +6598,8 @@ inline FEEvaluationAccess<dim, dim, 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<dim, dim, Number, is_face, VectorizedArrayType>( data_in, dof_no, @@ -6586,7 +6609,8 @@ inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>:: 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 FEFaceEvaluation<dim, const unsigned int quad_no, const unsigned int first_selected_component, const unsigned int active_fe_index, - const unsigned int active_quad_index) + const unsigned int active_quad_index, + const unsigned int face_type) : BaseClass(matrix_free, dof_no, first_selected_component, @@ -8421,7 +8448,8 @@ inline FEFaceEvaluation<dim, static_n_q_points, is_interior_face, active_fe_index, - active_quad_index) + active_quad_index, + face_type) , dofs_per_component(this->data->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<dim, quad_no, first_selected_component, matrix_free.get_face_active_fe_index(range, - is_interior_face)) + is_interior_face), + numbers::invalid_unsigned_int, + matrix_free.get_face_info(range.first).face_type) {} diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index d51a855ba7..445e28dea7 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -1843,8 +1843,6 @@ namespace internal #ifndef DEAL_II_WITH_SIMPLEX_SUPPORT // currently only non-hp-case... 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 diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 4229b31ba3..147b71f572 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -134,33 +134,45 @@ namespace { public: FaceRangeCompartor(const std::vector<unsigned int> &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 <int vectorization_width> bool operator()(const internal::MatrixFreeFunctions::FaceToCellTopology< - vectorization_width> &face, - const unsigned int & fe_index) + vectorization_width> & face, + const std::array<unsigned int, 2> &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<unsigned int, 2> 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 <int vectorization_width> bool operator()(const internal::MatrixFreeFunctions::FaceToCellTopology< - vectorization_width> & face, - const std::pair<unsigned int, unsigned int> &fe_index) + vectorization_width> & face, + const std::array<unsigned int, 3> &fe_index) { - const std::pair<unsigned int, unsigned int> 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<unsigned int, 3> 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<unsigned int> &fe_indices; const bool include; + const bool only_face_type; }; } // namespace @@ -540,13 +553,18 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit( const auto create_inner_face_subrange_hp_by_index = [&](const std::pair<unsigned int, unsigned int> &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<unsigned int, unsigned int> { - 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<unsigned int>(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<dim, Number, VectorizedArrayType>::internal_reinit( dof_info[dof_handler_index].max_fe_index); const std::vector<unsigned int> &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<unsigned int, unsigned int> return_range; return_range.first = - std::lower_bound(face_info.faces.begin() + range.first, - face_info.faces.begin() + range.second, - std::pair<unsigned int, unsigned int>{ - 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<unsigned int, 3>{ + {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<unsigned int, unsigned int>{ - 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<unsigned int, 3>{ + {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<dim, Number, VectorizedArrayType>::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<unsigned int>(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<dim, Number, VectorizedArrayType>::internal_reinit( const auto create_boundary_face_subrange_hp_by_index = [&](const std::pair<unsigned int, unsigned int> &range, + const unsigned int face_type, const unsigned int fe_index, const unsigned int dof_handler_index = 0) -> std::pair<unsigned int, unsigned int> { - 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<unsigned int>(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<unsigned int> &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<unsigned int, unsigned int> 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<unsigned int, 2>{{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<unsigned int, 2>{{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<dim, Number, VectorizedArrayType>::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<unsigned int>(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<const FE_Poly<dim, dim> *>(&fe); - if (dynamic_cast<const FE_SimplexP<dim, dim> *>(&fe) != nullptr || - dynamic_cast<const FE_SimplexDGP<dim, dim> *>(&fe) != nullptr || - dynamic_cast<const FE_WedgeP<dim, dim> *>(&fe) != nullptr || - dynamic_cast<const FE_PyramidP<dim, dim> *>(&fe) != nullptr) + if (dynamic_cast<const FE_SimplexPoly<dim, dim> *>(&fe) != nullptr || + dynamic_cast<const FE_Wedge<dim, dim> *>(&fe) != nullptr || + dynamic_cast<const FE_Pyramid<dim, dim> *>(&fe) != nullptr) { scalar_lexicographic.resize(fe.n_dofs_per_cell()); for (unsigned int i = 0; i < scalar_lexicographic.size(); ++i)