From: Peter Munch Date: Fri, 19 Jul 2024 13:11:12 +0000 (+0200) Subject: FEInterfaceValues: fill indices for TriaIterator X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=174cd32d21cf100b6c9f2fc9686359975d22f457;p=dealii.git FEInterfaceValues: fill indices for TriaIterator --- diff --git a/include/deal.II/fe/fe_interface_values.h b/include/deal.II/fe/fe_interface_values.h index 9093cb00d8..3b0b2371f7 100644 --- a/include/deal.II/fe/fe_interface_values.h +++ b/include/deal.II/fe/fe_interface_values.h @@ -2339,6 +2339,9 @@ FEInterfaceValues::reinit( std::is_same_v, typename CellNeighborIteratorType::AccessorType>; + unsigned int active_fe_index = 0; + unsigned int active_fe_index_neighbor = 0; + if (internal_fe_face_values) { if (sub_face_no == numbers::invalid_unsigned_int) @@ -2373,8 +2376,8 @@ FEInterfaceValues::reinit( } else if (internal_hp_fe_face_values) { - unsigned int active_fe_index = fe_index_in; - unsigned int active_fe_index_neighbor = + active_fe_index = fe_index_in; + active_fe_index_neighbor = (fe_index_neighbor_in != numbers::invalid_unsigned_int) ? fe_index_neighbor_in : fe_index_in; @@ -2421,37 +2424,42 @@ FEInterfaceValues::reinit( // Third 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 = - ((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( - {active_fe_index, active_fe_index_neighbor}) : - 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) + if ((used_q_index == numbers::invalid_unsigned_int) || + (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; + const unsigned int dominated_fe_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( + {active_fe_index, active_fe_index_neighbor}) : + 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. @@ -2550,6 +2558,28 @@ FEInterfaceValues::reinit( ++idx; } } + else + { + const unsigned int n_dofs_per_cell_1 = fe_face_values->dofs_per_cell; + const unsigned int n_dofs_per_cell_2 = + fe_face_values_neighbor->dofs_per_cell; + + interface_dof_indices.resize(n_dofs_per_cell_1 + n_dofs_per_cell_2); + dofmap.resize(n_dofs_per_cell_1 + n_dofs_per_cell_2); + + for (unsigned int i = 0; i < n_dofs_per_cell_1; ++i) + { + interface_dof_indices[i] = numbers::invalid_dof_index; + dofmap[i] = {i, numbers::invalid_unsigned_int}; + } + + for (unsigned int i = 0; i < n_dofs_per_cell_2; ++i) + { + interface_dof_indices[i + n_dofs_per_cell_1] = + numbers::invalid_dof_index; + dofmap[i + n_dofs_per_cell_1] = {numbers::invalid_unsigned_int, i}; + } + } } @@ -2589,30 +2619,24 @@ FEInterfaceValues::reinit(const CellIteratorType &cell, fe_face_values_neighbor = nullptr; } + interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell()); + if constexpr (std::is_same_v> || std::is_same_v>) { - if (internal_fe_face_values) - { - interface_dof_indices.resize( - fe_face_values->get_fe().n_dofs_per_cell()); - cell->get_active_or_mg_dof_indices(interface_dof_indices); - } - else if (internal_hp_fe_face_values) - { - interface_dof_indices.resize( - fe_face_values->get_fe().n_dofs_per_cell()); - cell->get_active_or_mg_dof_indices(interface_dof_indices); - } - - dofmap.resize(interface_dof_indices.size()); - for (unsigned int i = 0; i < interface_dof_indices.size(); ++i) - { - dofmap[i] = {{i, numbers::invalid_unsigned_int}}; - } + cell->get_active_or_mg_dof_indices(interface_dof_indices); } + else + { + for (auto &i : interface_dof_indices) + i = numbers::invalid_dof_index; + } + + dofmap.resize(interface_dof_indices.size()); + for (unsigned int i = 0; i < interface_dof_indices.size(); ++i) + dofmap[i] = {{i, numbers::invalid_unsigned_int}}; } diff --git a/tests/feinterface/fe_interface_values_02.cc b/tests/feinterface/fe_interface_values_02.cc index 8d62a08313..4712f3e8b9 100644 --- a/tests/feinterface/fe_interface_values_02.cc +++ b/tests/feinterface/fe_interface_values_02.cc @@ -32,31 +32,10 @@ #include "../test_grids.h" - -template +template void -test(unsigned int fe_degree) +test(const IteratorType &cell, FEInterfaceValues &fiv) { - Triangulation tria; - TestGrids::hyper_line(tria, 2); - - DoFHandler dofh(tria); - FE_DGQ fe(fe_degree); - deallog << fe.get_name() << std::endl; - dofh.distribute_dofs(fe); - - MappingQ mapping(1); - UpdateFlags update_flags = update_values | update_gradients | - update_quadrature_points | update_JxW_values; - - FEInterfaceValues fiv(mapping, - fe, - QGauss(fe.degree + 1), - update_flags); - - - auto cell = dofh.begin(); - for (const unsigned int f : GeometryInfo::face_indices()) if (!cell->at_boundary(f)) { @@ -103,6 +82,32 @@ test(unsigned int fe_degree) } +template +void +test(unsigned int fe_degree) +{ + Triangulation tria; + TestGrids::hyper_line(tria, 2); + + DoFHandler dofh(tria); + FE_DGQ fe(fe_degree); + deallog << fe.get_name() << std::endl; + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + FEInterfaceValues fiv(mapping, + fe, + QGauss(fe.degree + 1), + update_flags); + + test(tria.begin(), fiv); + test(dofh.begin(), fiv); +} + + int main() diff --git a/tests/feinterface/fe_interface_values_02.output b/tests/feinterface/fe_interface_values_02.output index fe3840c855..7360480b12 100644 --- a/tests/feinterface/fe_interface_values_02.output +++ b/tests/feinterface/fe_interface_values_02.output @@ -4,18 +4,34 @@ DEAL::shape_value(true): 1.00000 0.00000 DEAL::shape_value(false): 0.00000 1.00000 DEAL::jump_in_shape_values(): 1.00000 -1.00000 DEAL::average_of_shape_values(): 0.500000 0.500000 +DEAL::shape_value(true): 1.00000 0.00000 +DEAL::shape_value(false): 0.00000 1.00000 +DEAL::jump_in_shape_values(): 1.00000 -1.00000 +DEAL::average_of_shape_values(): 0.500000 0.500000 DEAL::FE_DGQ<2>(1) DEAL::shape_value(true): 0.00000 0.500000 0.00000 0.500000 0.00000 0.00000 0.00000 0.00000 DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 0.500000 0.00000 0.500000 0.00000 DEAL::jump_in_shape_values(): 0.00000 0.500000 0.00000 0.500000 -0.500000 0.00000 -0.500000 0.00000 DEAL::average_of_shape_values(): 0.00000 0.250000 0.00000 0.250000 0.250000 0.00000 0.250000 0.00000 +DEAL::shape_value(true): 0.00000 0.500000 0.00000 0.500000 0.00000 0.00000 0.00000 0.00000 +DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 0.500000 0.00000 0.500000 0.00000 +DEAL::jump_in_shape_values(): 0.00000 0.500000 0.00000 0.500000 -0.500000 0.00000 -0.500000 0.00000 +DEAL::average_of_shape_values(): 0.00000 0.250000 0.00000 0.250000 0.250000 0.00000 0.250000 0.00000 DEAL::FE_DGQ<3>(0) DEAL::shape_value(true): 1.00000 0.00000 DEAL::shape_value(false): 0.00000 1.00000 DEAL::jump_in_shape_values(): 1.00000 -1.00000 DEAL::average_of_shape_values(): 0.500000 0.500000 +DEAL::shape_value(true): 1.00000 0.00000 +DEAL::shape_value(false): 0.00000 1.00000 +DEAL::jump_in_shape_values(): 1.00000 -1.00000 +DEAL::average_of_shape_values(): 0.500000 0.500000 DEAL::FE_DGQ<3>(1) DEAL::shape_value(true): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 DEAL::jump_in_shape_values(): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 -0.250000 0.00000 -0.250000 0.00000 -0.250000 0.00000 -0.250000 0.00000 DEAL::average_of_shape_values(): 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 +DEAL::shape_value(true): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL::shape_value(false): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 +DEAL::jump_in_shape_values(): 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 0.00000 0.250000 -0.250000 0.00000 -0.250000 0.00000 -0.250000 0.00000 -0.250000 0.00000 +DEAL::average_of_shape_values(): 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000 0.125000 0.00000