From 629fdc839d5a1e7204f96fc365637d9c229efd66 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Tue, 20 Dec 2022 21:29:40 +0100 Subject: [PATCH] Add hp support for interfaces to ScratchData --- include/deal.II/meshworker/scratch_data.h | 2 ++ source/meshworker/scratch_data.cc | 42 +++++++++++++++++------ 2 files changed, 34 insertions(+), 10 deletions(-) diff --git a/include/deal.II/meshworker/scratch_data.h b/include/deal.II/meshworker/scratch_data.h index 345fff7d6f..ee2d98b17d 100644 --- a/include/deal.II/meshworker/scratch_data.h +++ b/include/deal.II/meshworker/scratch_data.h @@ -1405,6 +1405,8 @@ namespace MeshWorker /** * Interface values on facets. */ + // The FEInterfaceValues class supports initialization with hp objects + // as well. std::unique_ptr> interface_fe_values; /** @} */ // non-hp data diff --git a/source/meshworker/scratch_data.cc b/source/meshworker/scratch_data.cc index 2935608346..434307edad 100644 --- a/source/meshworker/scratch_data.cc +++ b/source/meshworker/scratch_data.cc @@ -450,26 +450,48 @@ namespace MeshWorker interface_fe_values = std::make_unique>( *mapping, *fe, face_quadrature, face_update_flags); + interface_fe_values->reinit(cell, face_no, sub_face_no, cell_neighbor, face_no_neighbor, sub_face_no_neighbor); - - current_fe_values = &interface_fe_values->get_fe_face_values(0); - current_neighbor_fe_values = - &interface_fe_values->get_fe_face_values(1); - - cell_neighbor->get_dof_indices(neighbor_dof_indices); - local_dof_indices = interface_fe_values->get_interface_dof_indices(); - return *interface_fe_values; } else { - AssertThrow(false, ExcOnlyAvailableWithoutHP()); - return *interface_fe_values; + if (!interface_fe_values) + interface_fe_values = + std::make_unique>( + *mapping_collection, + *fe_collection, + face_quadrature_collection, + face_update_flags); + + // When we want to ensure some agreement between the cell face and + // its neighbor on the quadrature order and mapping to use on this + // face, then we defer to the dominance of one FE over another. This + // should ensure that the optimal integration order and mapping + // order are selected for this situation. + const unsigned int dominated_fe_index = + fe_collection->find_dominated_fe( + {cell->active_fe_index(), cell_neighbor->active_fe_index()}); + + interface_fe_values->reinit(cell, + face_no, + sub_face_no, + cell_neighbor, + face_no_neighbor, + sub_face_no_neighbor, + dominated_fe_index, + dominated_fe_index); } + + current_fe_values = &interface_fe_values->get_fe_face_values(0); + current_neighbor_fe_values = &interface_fe_values->get_fe_face_values(1); + + local_dof_indices = interface_fe_values->get_interface_dof_indices(); + return *interface_fe_values; } -- 2.39.5