From 3afb70dfc68c38559197d8691eb0ffd225e22d7a Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Wed, 24 Nov 2021 14:28:11 -0500 Subject: [PATCH] Separate mapping related computations from FEPointEvaluation --- .../deal.II/matrix_free/fe_point_evaluation.h | 105 +++++++++--------- source/matrix_free/fe_point_evaluation.cc | 21 ++++ .../matrix_free/fe_point_evaluation.inst.in | 2 + 3 files changed, 75 insertions(+), 53 deletions(-) diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index 0ab23f57f8..ddd9007e0e 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -357,13 +357,24 @@ namespace internal is_fast_path_supported(const FiniteElement &fe, const unsigned int base_element_number); + template + bool + is_fast_path_supported(const Mapping &mapping); + template std::vector> get_polynomial_space(const FiniteElement &fe); } // namespace FEPointEvaluation } // namespace internal - +template +void +compute_mapping_data( + const Mapping & mapping, + const typename Triangulation::cell_iterator &cell, + const ArrayView> & unit_points, + dealii::internal::FEValuesImplementation::MappingRelatedData + &mapping_data); /** * This class provides an interface to the evaluation of interpolated solution @@ -615,18 +626,6 @@ private: */ SmartPointer> mapping; - /** - * Pointer to the function of the mapping that computes the necessary - * mapping quantities like quadrature points and Jacobians. - */ - std::function::cell_iterator & /*cell*/, - const ArrayView> & /*unit_points*/, - const UpdateFlags /*update_flags*/, - dealii::internal::FEValuesImplementation::MappingRelatedData - & /*output_data*/)> - fill_mapping_data_for_generic_points; - /** * Pointer to the FiniteElement object passed to the constructor. */ @@ -747,41 +746,6 @@ FEPointEvaluation::FEPointEvaluation( AssertIndexRange(first_selected_component + n_components, fe.n_components() + 1); - if (const MappingQ *mapping_q = - dynamic_cast *>(&mapping)) - { - fill_mapping_data_for_generic_points = - [mapping_q]( - const typename Triangulation::cell_iterator &cell, - const ArrayView> &unit_points, - const UpdateFlags update_flags, - dealii::internal::FEValuesImplementation::MappingRelatedData - &output_data) -> void { - mapping_q->fill_mapping_data_for_generic_points(cell, - unit_points, - update_flags, - output_data); - }; - } - else if (const MappingCartesian *mapping_cartesian = - dynamic_cast *>(&mapping)) - { - fill_mapping_data_for_generic_points = - [mapping_cartesian]( - const typename Triangulation::cell_iterator &cell, - const ArrayView> &unit_points, - const UpdateFlags update_flags, - dealii::internal::FEValuesImplementation::MappingRelatedData - &output_data) -> void { - mapping_cartesian->fill_mapping_data_for_generic_points(cell, - unit_points, - update_flags, - output_data); - }; - } - bool same_base_element = true; unsigned int base_element_number = 0; component_in_base_element = 0; @@ -799,7 +763,7 @@ FEPointEvaluation::FEPointEvaluation( else component += fe.element_multiplicity(base_element_number); - if (fill_mapping_data_for_generic_points && + if (internal::FEPointEvaluation::is_fast_path_supported(mapping) && internal::FEPointEvaluation::is_fast_path_supported( fe, base_element_number) && same_base_element) @@ -847,6 +811,43 @@ FEPointEvaluation::FEPointEvaluation( } +template +void +compute_mapping_data_for_generic_points( + const Mapping & mapping, + const typename Triangulation::cell_iterator &cell, + const ArrayView> & unit_points, + const UpdateFlags update_flags, + dealii::internal::FEValuesImplementation::MappingRelatedData + &mapping_data) +{ + UpdateFlags update_flags_mapping = update_default; + // translate update flags + if (update_flags & update_jacobians) + update_flags_mapping |= update_jacobians; + if (update_flags & update_gradients || + update_flags & update_inverse_jacobians) + update_flags_mapping |= update_inverse_jacobians; + if (update_flags & update_quadrature_points) + update_flags_mapping |= update_quadrature_points; + + if (const MappingQ *mapping_q = + dynamic_cast *>(&mapping)) + { + mapping_q->fill_mapping_data_for_generic_points(cell, + unit_points, + update_flags_mapping, + mapping_data); + } + else if (const MappingCartesian *mapping_cartesian = + dynamic_cast *>(&mapping)) + { + mapping_cartesian->fill_mapping_data_for_generic_points( + cell, unit_points, update_flags_mapping, mapping_data); + } +} + + template void @@ -856,10 +857,8 @@ FEPointEvaluation::reinit( { // If using the fast path, we need to precompute the mapping data. if (!poly.empty()) - fill_mapping_data_for_generic_points(cell, - unit_points, - update_flags_mapping, - mapping_data); + compute_mapping_data_for_generic_points( + *mapping, cell, unit_points, update_flags_mapping, mapping_data); // then call the other version of this function with the precomputed data reinit(cell, unit_points, mapping_data); diff --git a/source/matrix_free/fe_point_evaluation.cc b/source/matrix_free/fe_point_evaluation.cc index 7ecde1907a..3770f18435 100644 --- a/source/matrix_free/fe_point_evaluation.cc +++ b/source/matrix_free/fe_point_evaluation.cc @@ -78,6 +78,27 @@ namespace internal } + + template + bool + is_fast_path_supported(const Mapping &mapping) + { + if (const MappingQ *mapping_q = + dynamic_cast *>(&mapping)) + { + return true; + } + else if (const MappingCartesian *mapping_cartesian = + dynamic_cast *>( + &mapping)) + { + return true; + } + return false; + } + + + template std::vector> get_polynomial_space(const FiniteElement &fe) diff --git a/source/matrix_free/fe_point_evaluation.inst.in b/source/matrix_free/fe_point_evaluation.inst.in index 74736be9ff..6a50a3fc68 100644 --- a/source/matrix_free/fe_point_evaluation.inst.in +++ b/source/matrix_free/fe_point_evaluation.inst.in @@ -21,6 +21,8 @@ for (deal_II_dimension, deal_II_space_dimension : DIMENSIONS) template bool internal::FEPointEvaluation::is_fast_path_supported( const FiniteElement &, const unsigned int); + template bool internal::FEPointEvaluation::is_fast_path_supported( + const Mapping &); template std::vector> internal::FEPointEvaluation::get_polynomial_space( const FiniteElement &); -- 2.39.5