From 4919e33659cf047d4910383362ae808c96325e2a Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Sat, 8 Jan 2022 21:32:55 +0100 Subject: [PATCH] Add some hp-specialized reinit() functions to MeshWorker::ScratchData --- include/deal.II/meshworker/scratch_data.h | 101 ++++++++ source/meshworker/scratch_data.cc | 282 +++++++++++++++++----- 2 files changed, 323 insertions(+), 60 deletions(-) diff --git a/include/deal.II/meshworker/scratch_data.h b/include/deal.II/meshworker/scratch_data.h index 1b7971987c..0c64192be4 100644 --- a/include/deal.II/meshworker/scratch_data.h +++ b/include/deal.II/meshworker/scratch_data.h @@ -622,6 +622,107 @@ namespace MeshWorker /** @} */ // NeighborCellMethods + /** + * @name hp-compatible methods to work on cells and neighbor cells + */ + /** @{ */ // hpCellMethods + + /** + * Initialize the internal FEFaceValues to use the given @p face_no on the given + * @p cell, and return a reference to it. + * + * This variant of the reinit() function compares the active finite element + * on each cell, and chooses the dominated finite element's index to select + * the quadrature rule and mapping for the returned FEFaceValues object. + * This is useful in instances where the shape functions of the elements on + * either side of a face are being evaluated at the same time (such as is + * done in DG methods). See FECollection::find_dominated_fe() for more + * information on the selection process. + * + * After calling this function, get_current_fe_values() will return the + * same object of this method, as an FEValuesBase reference. + */ + const FEFaceValues & + reinit(const typename DoFHandler::active_cell_iterator &cell, + const typename DoFHandler::active_cell_iterator + & neighbor_cell, + const unsigned int face_no); + + /** + * Initialize the internal FESubfaceValues to use the given @p subface_no, + * on @p face_no, on the given @p cell, and return a reference to it. + * + * This variant of the reinit() function compares the active finite element + * on each cell, and chooses the dominated finite element's index to select + * the quadrature rule and mapping for the returned FEFaceValues object. + * This is useful in instances where the shape functions of the elements on + * either side of a face are being evaluated at the same time (such as is + * done in DG methods). See FECollection::find_dominated_fe() for more + * information on the selection process. + * + * After calling this function, get_current_fe_values() will return the + * same object of this method, as an FEValuesBase reference. + * + * If @p subface_no is numbers::invalid_unsigned_int, the reinit() function + * that takes only the @p cell and the @p face_no is called. + */ + const FEFaceValuesBase & + reinit(const typename DoFHandler::active_cell_iterator &cell, + const typename DoFHandler::active_cell_iterator + & neighbor_cell, + const unsigned int face_no, + const unsigned int subface_no); + + /** + * Initialize the internal FEFaceValues to use the given @p face_no on the + * given @p neighbor_cell, and return a reference to it. + * + * This variant of the reinit() function compares the active finite element + * on each cell, and chooses the dominated finite element's index to select + * the quadrature rule and mapping for the returned FEFaceValues object. + * This is useful in instances where the shape functions of the elements on + * either side of a face are being evaluated at the same time (such as is + * done in DG methods). See FECollection::find_dominated_fe() for more + * information on the selection process. + * + * After calling this function, get_current_neighbor_fe_values() will return + * the same object of this method, as an FEValuesBase reference. + */ + const FEFaceValues & + reinit_neighbor( + const typename DoFHandler::active_cell_iterator &cell, + const typename DoFHandler::active_cell_iterator + & neighbor_cell, + const unsigned int face_no); + + /** + * Initialize the internal FESubfaceValues to use the given @p subface_no, + * on @p face_no, on the given @p neighbor_cell, and return a reference to it. + * + * This variant of the reinit() function compares the active finite element + * on each cell, and chooses the dominated finite element's index to select + * the quadrature rule and mapping for the returned FEFaceValues object. + * This is useful in instances where the shape functions of the elements on + * either side of a face are being evaluated at the same time (such as is + * done in DG methods). See FECollection::find_dominated_fe() for more + * information on the selection process. + * + * After calling this function, get_current_neighbor_fe_values() will return + * the same object of this method, as an FEValuesBase reference. + * + * If @p subface_no is numbers::invalid_unsigned_int, the reinit() function + * that takes only the @p neighbor_cell and the @p face_no is called. + */ + const FEFaceValuesBase & + reinit_neighbor( + const typename DoFHandler::active_cell_iterator &cell, + const typename DoFHandler::active_cell_iterator + & neighbor_cell, + const unsigned int face_no, + const unsigned int subface_no); + + /** @} */ // hpCellMethods + /** * Return a GeneralDataStorage object that can be used to store any amount * of data, of any type, which is then made accessible by an identifier diff --git a/source/meshworker/scratch_data.cc b/source/meshworker/scratch_data.cc index 6d5394f341..6ac70b24f2 100644 --- a/source/meshworker/scratch_data.cc +++ b/source/meshworker/scratch_data.cc @@ -283,23 +283,58 @@ namespace MeshWorker } else { - if (!hp_fe_face_values) - hp_fe_face_values = std::make_unique>( - *mapping_collection, - *fe_collection, - face_quadrature_collection, - face_update_flags); + return reinit(cell, cell, face_no); + } + } - hp_fe_face_values->reinit(cell, face_no); - const auto &fe_face_values = hp_fe_face_values->get_present_fe_values(); - const auto &fe = (*fe_collection)[cell->active_fe_index()]; - local_dof_indices.resize(fe.n_dofs_per_cell()); - cell->get_dof_indices(local_dof_indices); - current_fe_values = &fe_face_values; - return fe_face_values; + template + const FEFaceValues & + ScratchData::reinit( + const typename DoFHandler::active_cell_iterator &cell, + const typename DoFHandler::active_cell_iterator + & neighbor_cell, + const unsigned int face_no) + { + Assert(hp_capability_enabled, ExcOnlyAvailableWithHP()); + + if (!hp_fe_face_values) + hp_fe_face_values = std::make_unique>( + *mapping_collection, + *fe_collection, + face_quadrature_collection, + face_update_flags); + + if (neighbor_cell == cell) + { + hp_fe_face_values->reinit(cell, face_no); } + else + { + // 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(), neighbor_cell->active_fe_index()}); + + hp_fe_face_values->reinit(cell, + face_no, + dominated_fe_index, + dominated_fe_index); + } + + const auto &fe_face_values = hp_fe_face_values->get_present_fe_values(); + const auto &fe = (*fe_collection)[cell->active_fe_index()]; + + local_dof_indices.resize(fe.n_dofs_per_cell()); + cell->get_dof_indices(local_dof_indices); + + current_fe_values = &fe_face_values; + return fe_face_values; } @@ -328,28 +363,72 @@ namespace MeshWorker } else { - if (!hp_fe_subface_values) - hp_fe_subface_values = - std::make_unique>( - *mapping_collection, - *fe_collection, - face_quadrature_collection, - face_update_flags); + return reinit(cell, cell, face_no, subface_no); + } + } + else + return reinit(cell, face_no); + } - hp_fe_subface_values->reinit(cell, face_no, subface_no); - const auto &fe_subface_values = - hp_fe_subface_values->get_present_fe_values(); - const auto &fe = (*fe_collection)[cell->active_fe_index()]; - local_dof_indices.resize(fe.n_dofs_per_cell()); - cell->get_dof_indices(local_dof_indices); - current_fe_values = &fe_subface_values; - return fe_subface_values; + template + const FEFaceValuesBase & + ScratchData::reinit( + const typename DoFHandler::active_cell_iterator &cell, + const typename DoFHandler::active_cell_iterator + & neighbor_cell, + const unsigned int face_no, + const unsigned int subface_no) + { + Assert(hp_capability_enabled, ExcOnlyAvailableWithHP()); + + if (subface_no != numbers::invalid_unsigned_int) + { + if (!hp_fe_subface_values) + hp_fe_subface_values = + std::make_unique>( + *mapping_collection, + *fe_collection, + face_quadrature_collection, + face_update_flags); + + if (neighbor_cell == cell) + { + hp_fe_subface_values->reinit(cell, face_no, subface_no); + } + else + { + // 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(), neighbor_cell->active_fe_index()}); + + hp_fe_subface_values->reinit(cell, + face_no, + subface_no, + dominated_fe_index, + dominated_fe_index); } + + const auto &fe_subface_values = + hp_fe_subface_values->get_present_fe_values(); + const auto &fe = (*fe_collection)[cell->active_fe_index()]; + + local_dof_indices.resize(fe.n_dofs_per_cell()); + cell->get_dof_indices(local_dof_indices); + + current_fe_values = &fe_subface_values; + return fe_subface_values; } else - return reinit(cell, face_no); + { + return reinit(cell, neighbor_cell, face_no); + } } @@ -456,25 +535,61 @@ namespace MeshWorker } else { - if (!neighbor_hp_fe_face_values) - neighbor_hp_fe_face_values = - std::make_unique>( - *mapping_collection, - *fe_collection, - face_quadrature_collection, - neighbor_face_update_flags); + return reinit_neighbor(cell, cell, face_no); + } + } - neighbor_hp_fe_face_values->reinit(cell, face_no); - const auto &neighbor_fe_face_values = - neighbor_hp_fe_face_values->get_present_fe_values(); - const auto &neighbor_fe = (*fe_collection)[cell->active_fe_index()]; - neighbor_dof_indices.resize(neighbor_fe.n_dofs_per_cell()); - cell->get_dof_indices(neighbor_dof_indices); - current_neighbor_fe_values = &neighbor_fe_face_values; - return neighbor_fe_face_values; + template + const FEFaceValues & + ScratchData::reinit_neighbor( + const typename DoFHandler::active_cell_iterator &cell, + const typename DoFHandler::active_cell_iterator + & neighbor_cell, + const unsigned int face_no) + { + Assert(hp_capability_enabled, ExcOnlyAvailableWithHP()); + + if (!neighbor_hp_fe_face_values) + neighbor_hp_fe_face_values = + std::make_unique>( + *mapping_collection, + *fe_collection, + face_quadrature_collection, + neighbor_face_update_flags); + + if (neighbor_cell == cell) + { + neighbor_hp_fe_face_values->reinit(neighbor_cell, face_no); + } + else + { + // 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(), neighbor_cell->active_fe_index()}); + + neighbor_hp_fe_face_values->reinit(neighbor_cell, + face_no, + dominated_fe_index, + dominated_fe_index); } + + const auto &neighbor_fe_face_values = + neighbor_hp_fe_face_values->get_present_fe_values(); + const auto &neighbor_fe = + (*fe_collection)[neighbor_cell->active_fe_index()]; + + neighbor_dof_indices.resize(neighbor_fe.n_dofs_per_cell()); + neighbor_cell->get_dof_indices(neighbor_dof_indices); + + current_neighbor_fe_values = &neighbor_fe_face_values; + return neighbor_fe_face_values; } @@ -501,28 +616,75 @@ namespace MeshWorker } else { - if (!neighbor_hp_fe_subface_values) - neighbor_hp_fe_subface_values = - std::make_unique>( - *mapping_collection, - *fe_collection, - face_quadrature_collection, - neighbor_face_update_flags); + return reinit_neighbor(cell, cell, face_no, subface_no); + } + } + else + return reinit_neighbor(cell, face_no); + } - neighbor_hp_fe_subface_values->reinit(cell, face_no, subface_no); - const auto &neighbor_fe_subface_values = - neighbor_hp_fe_subface_values->get_present_fe_values(); - const auto &neighbor_fe = (*fe_collection)[cell->active_fe_index()]; - neighbor_dof_indices.resize(neighbor_fe.n_dofs_per_cell()); - cell->get_dof_indices(neighbor_dof_indices); - current_neighbor_fe_values = &neighbor_fe_subface_values; - return neighbor_fe_subface_values; + template + const FEFaceValuesBase & + ScratchData::reinit_neighbor( + const typename DoFHandler::active_cell_iterator &cell, + const typename DoFHandler::active_cell_iterator + & neighbor_cell, + const unsigned int face_no, + const unsigned int subface_no) + { + Assert(hp_capability_enabled, ExcOnlyAvailableWithHP()); + + if (subface_no != numbers::invalid_unsigned_int) + { + if (!neighbor_hp_fe_subface_values) + neighbor_hp_fe_subface_values = + std::make_unique>( + *mapping_collection, + *fe_collection, + face_quadrature_collection, + neighbor_face_update_flags); + + if (neighbor_cell == cell) + { + neighbor_hp_fe_subface_values->reinit(neighbor_cell, + face_no, + subface_no); + } + else + { + // 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(), neighbor_cell->active_fe_index()}); + + neighbor_hp_fe_subface_values->reinit(neighbor_cell, + face_no, + subface_no, + dominated_fe_index, + dominated_fe_index); } + + const auto &neighbor_fe_subface_values = + neighbor_hp_fe_subface_values->get_present_fe_values(); + const auto &neighbor_fe = + (*fe_collection)[neighbor_cell->active_fe_index()]; + + neighbor_dof_indices.resize(neighbor_fe.n_dofs_per_cell()); + neighbor_cell->get_dof_indices(neighbor_dof_indices); + + current_neighbor_fe_values = &neighbor_fe_subface_values; + return neighbor_fe_subface_values; } else - return reinit_neighbor(cell, face_no); + { + return reinit_neighbor(cell, neighbor_cell, face_no); + } } -- 2.39.5