From 29043152f07918e173cf6c31ec458ef83da64834 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 25 Apr 2023 13:32:50 -0600 Subject: [PATCH] Implement a second way to select a quadrature and mapping object. --- include/deal.II/fe/fe_interface_values.h | 84 +++++++++++++++--------- 1 file changed, 54 insertions(+), 30 deletions(-) diff --git a/include/deal.II/fe/fe_interface_values.h b/include/deal.II/fe/fe_interface_values.h index 805d84f48d..835f9cb895 100644 --- a/include/deal.II/fe/fe_interface_values.h +++ b/include/deal.II/fe/fe_interface_values.h @@ -1483,8 +1483,13 @@ public: * -# If one of these arguments is left at its default value, then the * function will need to choose a quadrature and/or mapping that is * appropriate for the two finite element spaces used on the two cells - * adjacent to the current interface. To this end, we try to find out - * which of the two spaces is "larger" (say, if you had used $Q_2$ + * adjacent to the current interface. As the first choice, if the + * quadrature or mapping collection we are considering has only one + * element, then that is clearly the one that should be used. + * -# If the quadrature or mapping collection have multiple elements, + * then we need to dig further. To this end, we try to find out + * which of the two finite element spaces on the two adjacent cells + * is "larger" (say, if you had used $Q_2$ * and $Q_4$ elements on the two adjacent cells, then the $Q_4$ * element is the larger one); the determination of which space * is "larger" is made using the hp::FECollection::find_dominated_fe() @@ -1509,7 +1514,7 @@ public: * hp::MappingCollection to retrieve quadrature and mapping objects * appropriate for the current face. * - There are cases where neither element dominates the other. For - * example, if one uses $Q_2\times Q_1$ and $Q_1\times Q_2$ elements on + * example, if one uses $Q_2\times Q_1$ and $Q_1\times Q_2$ elements on * neighboring cells, neither of the two spaces dominates the other -- * or, in the context of the current function, neither space is "larger" * than the other. In that case, there is no way for the current function @@ -2543,34 +2548,53 @@ FEInterfaceValues::reinit( } else if (internal_hp_fe_face_values) { + unsigned int used_q_index = q_index; + unsigned int used_mapping_index = mapping_index; + + // First check. If there is only one element in a collection, and if none + // had been specified explicitly, then that's clearly the one to take: + if (used_q_index == numbers::invalid_unsigned_int) + if (internal_hp_fe_face_values->get_quadrature_collection().size() == 1) + used_q_index = 0; + + if (used_mapping_index == numbers::invalid_unsigned_int) + if (internal_hp_fe_face_values->get_mapping_collection().size() == 1) + used_mapping_index = 0; + + // Second check, if the above did not already suffice. We see if we + // can get somewhere via the dominated's finite element index. const unsigned int dominated_fe_index = - internal_hp_fe_face_values->get_fe_collection().find_dominated_fe( - {cell->active_fe_index(), cell_neighbor->active_fe_index()}); - - if (q_index == numbers::invalid_unsigned_int) - Assert(dominated_fe_index != numbers::invalid_fe_index, - ExcMessage("You called this function with 'q_index' left at its " - "default value, but this can only work if one of " - "the two finite elements adjacent to this face " - "dominates the other. See the documentation " - "of this function for more information of how " - "to deal with this situation.")); - if (mapping_index == numbers::invalid_unsigned_int) - Assert(dominated_fe_index != numbers::invalid_fe_index, - ExcMessage("You called this function with 'mapping_index' left " - "at its default value, but this can only work if one " - "of the two finite elements adjacent to this face " - "dominates the other. See the documentation " - "of this function for more information of how " - "to deal with this situation.")); - - - const unsigned int used_q_index = - (q_index == numbers::invalid_unsigned_int ? dominated_fe_index : - q_index); - const unsigned int used_mapping_index = - (mapping_index == numbers::invalid_unsigned_int ? dominated_fe_index : - mapping_index); + ((used_q_index == numbers::invalid_unsigned_int) || + (used_mapping_index == numbers::invalid_unsigned_int) ? + internal_hp_fe_face_values->get_fe_collection().find_dominated_fe( + {cell->active_fe_index(), cell_neighbor->active_fe_index()}) : + numbers::invalid_unsigned_int); + + if (used_q_index == numbers::invalid_unsigned_int) + { + Assert(dominated_fe_index != numbers::invalid_fe_index, + ExcMessage( + "You called this function with 'q_index' left at its " + "default value, but this can only work if one of " + "the two finite elements adjacent to this face " + "dominates the other. See the documentation " + "of this function for more information of how " + "to deal with this situation.")); + used_q_index = dominated_fe_index; + } + + if (used_mapping_index == numbers::invalid_unsigned_int) + { + Assert(dominated_fe_index != numbers::invalid_fe_index, + ExcMessage( + "You called this function with 'mapping_index' left " + "at its default value, but this can only work if one " + "of the two finite elements adjacent to this face " + "dominates the other. See the documentation " + "of this function for more information of how " + "to deal with this situation.")); + used_mapping_index = dominated_fe_index; + } // Same as if above, but when hp is enabled. if (sub_face_no == numbers::invalid_unsigned_int) -- 2.39.5