From 17a4b285d867a9d830dd1c232a8320edb05d7fa5 Mon Sep 17 00:00:00 2001 From: Maximilian Bergbauer Date: Fri, 31 Mar 2023 16:58:24 +0200 Subject: [PATCH] Vectorize MappingInfo --- include/deal.II/base/derivative_form.h | 20 + include/deal.II/base/vectorization.h | 125 +- include/deal.II/fe/mapping.h | 9 +- .../deal.II/matrix_free/fe_point_evaluation.h | 786 +++++------ .../matrix_free/tensor_product_kernels.h | 2 +- include/deal.II/non_matching/mapping_info.h | 1149 ++++++++++++----- source/fe/mapping_q.cc | 6 +- tests/non_matching/mapping_info.cc | 10 +- tests/non_matching/mapping_info.output | 8 + 9 files changed, 1265 insertions(+), 850 deletions(-) diff --git a/include/deal.II/base/derivative_form.h b/include/deal.II/base/derivative_form.h index 3cb74f5e9e..9c36923c66 100644 --- a/include/deal.II/base/derivative_form.h +++ b/include/deal.II/base/derivative_form.h @@ -106,6 +106,13 @@ public: DerivativeForm & operator=(const Tensor<1, dim, Number> &); + /** + * Number conversion operator. + */ + template + DerivativeForm & + operator=(const DerivativeForm &df); + /** * Converts a DerivativeForm to Tensor. In particular, if order == 1 and the derivative is the Jacobian of @@ -252,6 +259,19 @@ DerivativeForm::operator=( +template +template +inline DerivativeForm & +DerivativeForm::operator=( + const DerivativeForm &df) +{ + for (unsigned int j = 0; j < spacedim; ++j) + (*this)[j] = df[j]; + return *this; +} + + + template inline Tensor & DerivativeForm::operator[](const unsigned int i) diff --git a/include/deal.II/base/vectorization.h b/include/deal.II/base/vectorization.h index 9aea1cb6a4..7bac5b49cf 100644 --- a/include/deal.II/base/vectorization.h +++ b/include/deal.II/base/vectorization.h @@ -5478,49 +5478,150 @@ namespace internal template struct VectorizedArrayTrait { - using value_type = T; + /** + * Define scalar value type. + */ + using value_type = T; + + /** + * Define width of template type. + */ static constexpr std::size_t width = 1; - static T & - get(T &value, unsigned int c) + /** + * Define vectorized value type for internal vectorization. + */ + using vectorized_value_type = VectorizedArray; + + /** + * Define a stride which defines how often the template type T fits into the + * vectorized_value_type. This is useful to write vectorized templated code + * where the internal computation is vectorized and the user interface is + * optionally scalar or also vectorized. + */ + static constexpr std::size_t stride = vectorized_value_type::size(); + + /** + * Get a reference to scalar value (on lane 0). + */ + static value_type & + get(value_type &value, unsigned int c) { - AssertDimension(c, 0); + AssertIndexRange(c, 1); (void)c; return value; } - static const T & - get(const T &value, unsigned int c) + /** + * Get a read-only reference to scalar value (on lane 0). + */ + static const value_type & + get(const value_type &value, unsigned int c) { - AssertDimension(c, 0); + AssertIndexRange(c, 1); (void)c; return value; } + + /** + * Get a reference to scalar value on lane c from a vectorized values field. + */ + static value_type & + get_from_vectorized(vectorized_value_type &values, unsigned int c) + { + AssertIndexRange(c, stride); + + return values[c]; + } + + /** + * Get a read-only reference to scalar value on lane c from a vectorized + * values field. + */ + static const value_type & + get_from_vectorized(const vectorized_value_type &values, unsigned int c) + { + AssertIndexRange(c, stride); + + return values[c]; + } }; template struct VectorizedArrayTrait> { - using value_type = T; + /** + * Define scalar value type. + */ + using value_type = T; + + /** + * Define width of template type. + */ static constexpr std::size_t width = width_; - static T & - get(VectorizedArray &values, unsigned int c) + /** + * Define vectorized value type for internal vectorization. + */ + using vectorized_value_type = VectorizedArray; + + /** + * Define a stride which defines how often the template type + * VectorizedArray fits into the vectorized value type. This is + * useful to write vectorized templated code where the internal computation + * is vectorized and the user interface is optionally scalar or also + * vectorized. + */ + static constexpr std::size_t stride = 1; + + /** + * Get a reference to scalar value on lane c. + */ + static value_type & + get(vectorized_value_type &values, unsigned int c) { AssertIndexRange(c, width_); return values[c]; } - static const T & - get(const VectorizedArray &values, unsigned int c) + /** + * Get a read-only reference to scalar value on lane c. + */ + static const value_type & + get(const vectorized_value_type &values, unsigned int c) { AssertIndexRange(c, width_); return values[c]; } + + /** + * Get a reference to vectorized values from a vectorized values field. + */ + static vectorized_value_type & + get_from_vectorized(vectorized_value_type &values, unsigned int c) + { + (void)c; + AssertIndexRange(c, stride); + + return values; + } + + /** + * Get a read-only reference to vectorized values from a vectorized values + * field. + */ + static const vectorized_value_type & + get_from_vectorized(const vectorized_value_type &values, unsigned int c) + { + (void)c; + AssertIndexRange(c, stride); + + return values; + } }; } // namespace internal diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index b426848215..f2af092067 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -55,8 +55,11 @@ namespace NonMatching { template class FEImmersedSurfaceValues; - template - class MappingInfo; + namespace internal + { + template + class ComputeMappingDataHelper; + } } // namespace NonMatching @@ -1320,7 +1323,7 @@ public: friend class FEFaceValues; friend class FESubfaceValues; friend class NonMatching::FEImmersedSurfaceValues; - friend class NonMatching::MappingInfo; + friend class NonMatching::internal::ComputeMappingDataHelper; }; diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index 4e07a9a447..00efffd856 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -40,20 +40,16 @@ namespace internal { namespace FEPointEvaluation { - template - struct VectorizedPointsTypeTraits - { - using value_type = VectorizedArray; - static constexpr std::size_t stride = value_type::size(); - }; - - template - struct VectorizedPointsTypeTraits> - { - using value_type = VectorizedArray; - static constexpr std::size_t stride = 1; - }; - + DeclException1( + ExcFEPointEvaluationAccessToUninitializedMappingField, + std::string, + << "You are requesting information from an FEPointEvaluation " + << "object for which this kind of information has not been computed. " + << "What information these objects compute is determined by the update_* " + << "flags you pass to MappingInfo() in the Constructor. Here, " + << "the operation you are attempting requires the <" << arg1 + << "> flag to be set, but it was apparently not specified " + << "upon initialization."); /** * Struct to distinguish between the value and gradient types of different @@ -65,8 +61,8 @@ namespace internal using ScalarNumber = typename internal::VectorizedArrayTrait::value_type; using VectorizedArrayType = - typename internal::FEPointEvaluation::VectorizedPointsTypeTraits< - Number>::value_type; + typename dealii::internal::VectorizedArrayTrait< + Number>::vectorized_value_type; using value_type = Tensor<1, n_components, Number>; using scalar_value_type = Tensor<1, n_components, ScalarNumber>; using vectorized_value_type = @@ -80,18 +76,18 @@ namespace internal Tensor<1, dim, Tensor<1, n_components, VectorizedArrayType>>; static void - read_value(const Number vector_entry, + read_value(const ScalarNumber vector_entry, const unsigned int component, - value_type & result) + scalar_value_type &result) { AssertIndexRange(component, n_components); result[component] = vector_entry; } static void - write_value(Number & vector_entry, - const unsigned int component, - const value_type & result) + write_value(VectorizedArrayType & vector_entry, + const unsigned int component, + const vectorized_value_type &result) { AssertIndexRange(component, n_components); vector_entry = result[component]; @@ -100,41 +96,24 @@ namespace internal static void set_gradient(const interface_vectorized_gradient_type &value, const unsigned int vector_lane, - scalar_gradient_type & result) + gradient_type & result) { for (unsigned int i = 0; i < n_components; ++i) for (unsigned int d = 0; d < dim; ++d) - result[i][d] = value[d][i][vector_lane]; - } - - static void - set_gradient(const interface_vectorized_gradient_type &value, - const unsigned int, - vectorized_gradient_type &result) - { - for (unsigned int i = 0; i < n_components; ++i) - for (unsigned int d = 0; d < dim; ++d) - result[i][d] = value[d][i]; + result[i][d] = + internal::VectorizedArrayTrait::get_from_vectorized( + value[d][i], vector_lane); } static void get_gradient(interface_vectorized_gradient_type &value, const unsigned int vector_lane, - const scalar_gradient_type & result) + const gradient_type & result) { for (unsigned int i = 0; i < n_components; ++i) for (unsigned int d = 0; d < dim; ++d) - value[d][i][vector_lane] = result[i][d]; - } - - static void - get_gradient(interface_vectorized_gradient_type &value, - const unsigned int, - const vectorized_gradient_type &result) - { - for (unsigned int i = 0; i < n_components; ++i) - for (unsigned int d = 0; d < dim; ++d) - value[d][i] = result[i][d]; + internal::VectorizedArrayTrait::get_from_vectorized( + value[d][i], vector_lane) = result[i][d]; } static void @@ -239,8 +218,8 @@ namespace internal using ScalarNumber = typename internal::VectorizedArrayTrait::value_type; using VectorizedArrayType = - typename internal::FEPointEvaluation::VectorizedPointsTypeTraits< - Number>::value_type; + typename dealii::internal::VectorizedArrayTrait< + Number>::vectorized_value_type; using value_type = Number; using scalar_value_type = ScalarNumber; using vectorized_value_type = VectorizedArrayType; @@ -249,17 +228,17 @@ namespace internal using vectorized_gradient_type = Tensor<1, dim, VectorizedArrayType>; static void - read_value(const Number vector_entry, + read_value(const ScalarNumber vector_entry, const unsigned int, - value_type &result) + scalar_value_type &result) { result = vector_entry; } static void - write_value(Number &vector_entry, + write_value(VectorizedArrayType &vector_entry, const unsigned int, - const value_type &result) + const vectorized_value_type &result) { vector_entry = result; } @@ -392,8 +371,8 @@ namespace internal using ScalarNumber = typename internal::VectorizedArrayTrait::value_type; using VectorizedArrayType = - typename internal::FEPointEvaluation::VectorizedPointsTypeTraits< - Number>::value_type; + typename dealii::internal::VectorizedArrayTrait< + Number>::vectorized_value_type; using value_type = Tensor<1, dim, Number>; using scalar_value_type = Tensor<1, dim, ScalarNumber>; using vectorized_value_type = Tensor<1, dim, VectorizedArrayType>; @@ -404,17 +383,17 @@ namespace internal Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>>; static void - read_value(const Number vector_entry, + read_value(const ScalarNumber vector_entry, const unsigned int component, - value_type & result) + scalar_value_type &result) { result[component] = vector_entry; } static void - write_value(Number & vector_entry, - const unsigned int component, - const value_type & result) + write_value(VectorizedArrayType & vector_entry, + const unsigned int component, + const vectorized_value_type &result) { vector_entry = result[component]; } @@ -422,41 +401,24 @@ namespace internal static void set_gradient(const interface_vectorized_gradient_type &value, const unsigned int vector_lane, - scalar_gradient_type & result) + gradient_type & result) { for (unsigned int i = 0; i < dim; ++i) for (unsigned int d = 0; d < dim; ++d) - result[i][d] = value[d][i][vector_lane]; - } - - static void - set_gradient(const interface_vectorized_gradient_type &value, - const unsigned int, - vectorized_gradient_type &result) - { - for (unsigned int i = 0; i < dim; ++i) - for (unsigned int d = 0; d < dim; ++d) - result[i][d] = value[d][i]; + result[i][d] = + internal::VectorizedArrayTrait::get_from_vectorized( + value[d][i], vector_lane); } static void get_gradient(interface_vectorized_gradient_type &value, const unsigned int vector_lane, - const scalar_gradient_type & result) - { - for (unsigned int i = 0; i < dim; ++i) - for (unsigned int d = 0; d < dim; ++d) - value[d][i][vector_lane] = result[i][d]; - } - - static void - get_gradient(interface_vectorized_gradient_type &value, - const unsigned int, - const vectorized_gradient_type &result) + const gradient_type & result) { for (unsigned int i = 0; i < dim; ++i) for (unsigned int d = 0; d < dim; ++d) - value[d][i] = result[i][d]; + internal::VectorizedArrayTrait::get_from_vectorized( + value[d][i], vector_lane) = result[i][d]; } static void @@ -561,8 +523,8 @@ namespace internal using ScalarNumber = typename internal::VectorizedArrayTrait::value_type; using VectorizedArrayType = - typename internal::FEPointEvaluation::VectorizedPointsTypeTraits< - Number>::value_type; + typename dealii::internal::VectorizedArrayTrait< + Number>::vectorized_value_type; using value_type = Number; using scalar_value_type = ScalarNumber; using vectorized_value_type = VectorizedArrayType; @@ -571,17 +533,17 @@ namespace internal using vectorized_gradient_type = Tensor<1, 1, VectorizedArrayType>; static void - read_value(const Number vector_entry, + read_value(const ScalarNumber vector_entry, const unsigned int, - value_type &result) + scalar_value_type &result) { result = vector_entry; } static void - write_value(Number &vector_entry, + write_value(VectorizedArrayType &vector_entry, const unsigned int, - const value_type &result) + const vectorized_value_type &result) { vector_entry = result; } @@ -763,9 +725,8 @@ public: using ScalarNumber = typename internal::VectorizedArrayTrait::value_type; - using VectorizedArrayType = - typename internal::FEPointEvaluation::VectorizedPointsTypeTraits< - Number>::value_type; + using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait< + Number>::vectorized_value_type; using ETT = typename internal::FEPointEvaluation:: EvaluatorTypeTraits; using value_type = typename ETT::value_type; @@ -812,9 +773,10 @@ public: * objects, this parameter allows to select a range of `n_components` * components starting from this parameter. */ - FEPointEvaluation(NonMatching::MappingInfo &mapping_info, - const FiniteElement & fe, - const unsigned int first_selected_component = 0); + FEPointEvaluation( + NonMatching::MappingInfo &mapping_info, + const FiniteElement & fe, + const unsigned int first_selected_component = 0); /** * Copy constructor. @@ -1021,6 +983,13 @@ public: quadrature_point_indices() const; private: + static constexpr std::size_t n_lanes_user_interface = + internal::VectorizedArrayTrait::width; + static constexpr std::size_t n_lanes_internal = + internal::VectorizedArrayTrait::width; + static constexpr std::size_t stride = + internal::VectorizedArrayTrait::stride; + /** * Common setup function for both constructors. Does the setup for both fast * and slow path. @@ -1068,15 +1037,16 @@ private: integrate_slow(const ArrayView & solution_values, const EvaluationFlags::EvaluationFlags &integration_flags); + /** - * Number of quadrature points of the current cell/face. + * Number of quadrature batches of the current cell/face. */ const unsigned int n_q_points; /** - * Number of active quadrature points of the last quadrature point batch. + * Number of quadrature points of the current cell/face. */ - const unsigned int n_filled_lanes_last_batch; + const unsigned int n_q_points_scalar; /** * Pointer to the Mapping object passed to the constructor. @@ -1137,6 +1107,43 @@ private: */ std::vector gradients; + /** + * Pointer to first unit point batch of current cell/face from MappingInfo, + * set internally during do_reinit(). + */ + const Point *unit_point_ptr; + + /** + * Pointer to real point of first quadrature point of current cell/face from + * MappingInfo, set internally during do_reinit(). + */ + const Point *real_point_ptr; + + /** + * Pointer to Jacobian of first quadrature point of current cell/face from + * MappingInfo, set internally during do_reinit(). + */ + const DerivativeForm<1, dim, spacedim, Number> *jacobian_ptr; + + /** + * Pointer to inverse Jacobian of first quadrature point of current cell/face + * from MappingInfo, set internally during do_reinit(). + */ + const DerivativeForm<1, spacedim, dim, Number> *inverse_jacobian_ptr; + + /** + * Pointer to normal vector of first quadrature point of current cell/face + * from MappingInfo, set internally during do_reinit(). + */ + const Tensor<1, spacedim, Number> *normal_ptr; + + /** + * Pointer to Jacobian determinant times quadrature weight of first quadrature + * point of current cell/face from MappingInfo, set internally during + * do_reinit(). + */ + const Number *JxW_ptr; + /** * Number of unknowns per component, i.e., number of unique basis functions, * for the chosen FiniteElement (or base element). @@ -1168,14 +1175,14 @@ private: /** * Pointer to mapping info on the fly computed during reinit. */ - std::unique_ptr> + std::unique_ptr> mapping_info_on_the_fly; /** * Pointer to currently used mapping info (either on the fly or external * precomputed). */ - SmartPointer> mapping_info; + SmartPointer> mapping_info; /** * The current cell index to access mapping data from mapping info. @@ -1220,13 +1227,14 @@ FEPointEvaluation::FEPointEvaluation( const UpdateFlags update_flags, const unsigned int first_selected_component) : n_q_points(numbers::invalid_unsigned_int) - , n_filled_lanes_last_batch(numbers::invalid_unsigned_int) + , n_q_points_scalar(numbers::invalid_unsigned_int) , mapping(&mapping) , fe(&fe) , update_flags(update_flags) , mapping_info_on_the_fly( - std::make_unique>(mapping, - update_flags)) + std::make_unique>( + mapping, + update_flags)) , mapping_info(mapping_info_on_the_fly.get()) , current_cell_index(numbers::invalid_unsigned_int) , current_face_number(numbers::invalid_unsigned_int) @@ -1239,11 +1247,11 @@ FEPointEvaluation::FEPointEvaluation( template FEPointEvaluation::FEPointEvaluation( - NonMatching::MappingInfo &mapping_info, - const FiniteElement & fe, - const unsigned int first_selected_component) + NonMatching::MappingInfo &mapping_info, + const FiniteElement & fe, + const unsigned int first_selected_component) : n_q_points(numbers::invalid_unsigned_int) - , n_filled_lanes_last_batch(numbers::invalid_unsigned_int) + , n_q_points_scalar(numbers::invalid_unsigned_int) , mapping(&mapping_info.get_mapping()) , fe(&fe) , update_flags(mapping_info.get_update_flags()) @@ -1263,7 +1271,7 @@ template FEPointEvaluation::FEPointEvaluation( FEPointEvaluation &other) noexcept : n_q_points(other.n_q_points) - , n_filled_lanes_last_batch(numbers::invalid_unsigned_int) + , n_q_points_scalar(other.n_q_points_scalar) , mapping(other.mapping) , fe(other.fe) , poly(other.poly) @@ -1281,7 +1289,7 @@ FEPointEvaluation::FEPointEvaluation( , fe_values(other.fe_values) , mapping_info_on_the_fly( other.mapping_info_on_the_fly ? - std::make_unique>( + std::make_unique>( *mapping, update_flags) : nullptr) @@ -1302,7 +1310,7 @@ template FEPointEvaluation::FEPointEvaluation( FEPointEvaluation &&other) noexcept : n_q_points(other.n_q_points) - , n_filled_lanes_last_batch(numbers::invalid_unsigned_int) + , n_q_points_scalar(other.n_q_points_scalar) , mapping(other.mapping) , fe(other.fe) , poly(other.poly) @@ -1482,54 +1490,55 @@ template void FEPointEvaluation::do_reinit() { - const auto unit_points = - mapping_info->get_unit_points(current_cell_index, current_face_number); - - const unsigned int n_q_points_unvectorized = unit_points.size(); - - if (std::is_same::value) - { - const_cast(n_q_points) = unit_points.size(); - } - else - { - const unsigned int n_lanes = - internal::VectorizedArrayTrait::width; - const_cast(n_filled_lanes_last_batch) = - n_q_points_unvectorized % n_lanes; - const_cast(n_q_points) = - n_q_points_unvectorized / n_lanes; - if (n_filled_lanes_last_batch > 0) - ++const_cast(n_q_points); - } + const_cast(n_q_points_scalar) = + mapping_info->get_n_q_points_unvectorized(current_cell_index, + current_face_number); + + const_cast(n_q_points) = + n_q_points_scalar / n_lanes_user_interface + + (n_q_points_scalar % n_lanes_user_interface > 0 ? 1 : 0); + + // set unit point pointer + const unsigned int unit_point_offset = + mapping_info->compute_unit_point_index_offset(current_cell_index, + current_face_number); + unit_point_ptr = mapping_info->get_unit_point(unit_point_offset); + + // set data pointers + const UpdateFlags update_flags_mapping = + mapping_info->get_update_flags_mapping(); + const unsigned int data_offset = + mapping_info->compute_data_index_offset(current_cell_index, + current_face_number); + if (update_flags_mapping & UpdateFlags::update_quadrature_points) + real_point_ptr = mapping_info->get_real_point(data_offset); + if (update_flags_mapping & UpdateFlags::update_jacobians) + jacobian_ptr = mapping_info->get_jacobian(data_offset); + if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) + inverse_jacobian_ptr = mapping_info->get_inverse_jacobian(data_offset); + if (update_flags_mapping & UpdateFlags::update_normal_vectors) + normal_ptr = mapping_info->get_normal_vector(data_offset); + if (update_flags_mapping & UpdateFlags::update_JxW_values) + JxW_ptr = mapping_info->get_JxW(data_offset); if (update_flags & update_values) values.resize(n_q_points, numbers::signaling_nan()); if (update_flags & update_gradients) gradients.resize(n_q_points, numbers::signaling_nan()); - if (!polynomials_are_hat_functions) + if (fast_path && !polynomials_are_hat_functions) { - const std::size_t n_points = unit_points.size(); - const std::size_t n_lanes = VectorizedArrayType::size(); const std::size_t n_batches = - n_points / n_lanes + (n_points % n_lanes > 0 ? 1 : 0); + n_q_points_scalar / n_lanes_internal + + (n_q_points_scalar % n_lanes_internal > 0 ? 1 : 0); const std::size_t n_shapes = poly.size(); shapes.resize_fast(n_batches * n_shapes); - for (unsigned int i = 0, qb = 0; i < n_points; i += n_lanes, ++qb) - { - // convert to vectorized format - Point vectorized_points; - for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) - for (unsigned int d = 0; d < dim; ++d) - vectorized_points[d][j] = unit_points[i + j][d]; - - auto view = - make_array_view(shapes.begin() + qb * n_shapes, - shapes.begin() + (qb * n_shapes + n_shapes)); - - internal::compute_values_of_array(view, poly, vectorized_points); - } + for (unsigned int qb = 0; qb < n_batches; ++qb) + internal::compute_values_of_array(make_array_view(shapes, + qb * n_shapes, + n_shapes), + poly, + unit_point_ptr[qb]); } is_reinitialized = true; @@ -1548,103 +1557,54 @@ FEPointEvaluation::evaluate_fast( solution_renumbered.resize(dofs_per_component); for (unsigned int comp = 0; comp < n_components; ++comp) for (unsigned int i = 0; i < dofs_per_component; ++i) - internal::FEPointEvaluation:: - EvaluatorTypeTraits::read_value( - solution_values[renumber[(component_in_base_element + comp) * - dofs_per_component + - i]], - comp, - solution_renumbered[i]); + ETT::read_value( + solution_values[renumber[(component_in_base_element + comp) * + dofs_per_component + + i]], + comp, + solution_renumbered[i]); // unit gradients are currently only implemented with the fast tensor // path unit_gradients.resize(n_q_points, numbers::signaling_nan()); - const auto unit_points = - mapping_info->get_unit_points(current_cell_index, current_face_number); - const auto &mapping_data = - mapping_info->get_mapping_data(current_cell_index, current_face_number); - - const std::size_t n_points = unit_points.size(); - const std::size_t n_lanes = - internal::VectorizedArrayTrait::width; - const std::size_t stride = - internal::FEPointEvaluation::VectorizedPointsTypeTraits::stride; - // loop over quadrature batches qb / points q - for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes) + const unsigned int n_shapes = poly.size(); + for (unsigned int qb = 0, q = 0; q < n_q_points_scalar; + ++qb, q += n_lanes_internal) { // compute - const unsigned int n_shapes = poly.size(); - const auto val_and_grad = [&]() { - if (polynomials_are_hat_functions) - { - // convert quadrature points to vectorized format - Point vectorized_points; - for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v) - for (unsigned int d = 0; d < dim; ++d) - vectorized_points[d][v] = unit_points[q + v][d]; - - return internal::evaluate_tensor_product_value_and_gradient_linear( - poly, solution_renumbered, vectorized_points); - } - else - return internal::evaluate_tensor_product_value_and_gradient_shapes< + const auto val_and_grad = + polynomials_are_hat_functions ? + internal::evaluate_tensor_product_value_and_gradient_linear( + poly, solution_renumbered, unit_point_ptr[qb]) : + internal::evaluate_tensor_product_value_and_gradient_shapes< dim, scalar_value_type, - VectorizedArrayType>(make_array_view(shapes.begin() + qb * n_shapes, - shapes.begin() + - (qb * n_shapes + n_shapes)), - poly.size(), + VectorizedArrayType>(make_array_view(shapes, + qb * n_shapes, + n_shapes), + n_shapes, solution_renumbered); - }(); if (evaluation_flags & EvaluationFlags::values) { - for (unsigned int v = 0; v < stride && q + v < n_points; ++v) - internal::FEPointEvaluation:: - EvaluatorTypeTraits::set_value( - val_and_grad.first, v, values[qb * stride + v]); + for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v) + ETT::set_value(val_and_grad.first, v, values[qb * stride + v]); } if (evaluation_flags & EvaluationFlags::gradients) { Assert(update_flags & update_gradients || update_flags & update_inverse_jacobians, ExcNotInitialized()); - if (std::is_same::value) - { - // convert back to standard format - for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v) - { - internal::FEPointEvaluation::EvaluatorTypeTraits< - dim, - n_components, - Number>::set_gradient(val_and_grad.second, - v, - unit_gradients[q + v]); - gradients[q + v] = apply_transformation( - mapping_data.inverse_jacobians[q + v].transpose(), - unit_gradients[q + v]); - } - } - else + + for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v) { - // convert to vectorized format - DerivativeForm<1, spacedim, dim, Number> - vectorized_inverse_transposed_jacobians; - for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int s = 0; s < spacedim; ++s) - internal::VectorizedArrayTrait::get( - vectorized_inverse_transposed_jacobians[s][d], v) = - mapping_data.inverse_jacobians[q + v].transpose()[s][d]; - - internal::FEPointEvaluation:: - EvaluatorTypeTraits::set_gradient( - val_and_grad.second, 0, unit_gradients[qb]); - gradients[qb] = - apply_transformation(vectorized_inverse_transposed_jacobians, - unit_gradients[qb]); + const unsigned int offset = qb * stride + v; + ETT::set_gradient(val_and_grad.second, v, unit_gradients[offset]); + gradients[offset] = + apply_transformation(inverse_jacobian_ptr[offset].transpose(), + unit_gradients[offset]); } } } @@ -1664,7 +1624,6 @@ FEPointEvaluation::evaluate_slow( "Not initialized. Please call FEPointEvaluation::reinit()!")); const std::size_t n_points = fe_values->get_quadrature().size(); - const std::size_t n_lanes = internal::VectorizedArrayTrait::width; if (evaluation_flags & EvaluationFlags::values) { @@ -1676,29 +1635,26 @@ FEPointEvaluation::evaluate_slow( for (unsigned int d = 0; d < n_components; ++d) if (nonzero_shape_function_component[i][d] && (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) for (unsigned int v = 0; - v < (q + n_lanes > n_points ? n_filled_lanes_last_batch : - n_lanes); + v < n_lanes_user_interface && q + v < n_points; ++v) - internal::FEPointEvaluation:: - EvaluatorTypeTraits::access( - values[qb], - v, - d, - fe_values->shape_value(i, q + v) * value); + ETT::access(values[qb], + v, + d, + fe_values->shape_value(i, q + v) * value); else if (nonzero_shape_function_component[i][d]) - for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) for (unsigned int v = 0; - v < (q + n_lanes > n_points ? n_filled_lanes_last_batch : - n_lanes); + v < n_lanes_user_interface && q + v < n_points; ++v) - internal::FEPointEvaluation:: - EvaluatorTypeTraits::access( - values[qb], - v, - d, - fe_values->shape_value_component(i, q + v, d) * value); + ETT::access(values[qb], + v, + d, + fe_values->shape_value_component(i, q + v, d) * + value); } } @@ -1712,29 +1668,26 @@ FEPointEvaluation::evaluate_slow( for (unsigned int d = 0; d < n_components; ++d) if (nonzero_shape_function_component[i][d] && (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) for (unsigned int v = 0; - v < (q + n_lanes > n_points ? n_filled_lanes_last_batch : - n_lanes); + v < n_lanes_user_interface && q + v < n_points; ++v) - internal::FEPointEvaluation:: - EvaluatorTypeTraits::access( - gradients[qb], - v, - d, - fe_values->shape_grad(i, q + v) * value); + ETT::access(gradients[qb], + v, + d, + fe_values->shape_grad(i, q + v) * value); else if (nonzero_shape_function_component[i][d]) - for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) for (unsigned int v = 0; - v < (q + n_lanes > n_points ? n_filled_lanes_last_batch : - n_lanes); + v < n_lanes_user_interface && q + v < n_points; ++v) - internal::FEPointEvaluation:: - EvaluatorTypeTraits::access( - gradients[qb], - v, - d, - fe_values->shape_grad_component(i, q + v, d) * value); + ETT::access(gradients[qb], + v, + d, + fe_values->shape_grad_component(i, q + v, d) * + value); } } } @@ -1773,135 +1726,83 @@ inline void FEPointEvaluation::integrate_fast( const ArrayView & solution_values, const EvaluationFlags::EvaluationFlags &integration_flags) - { // fast path with tensor product integration if (solution_renumbered_vectorized.size() != dofs_per_component) solution_renumbered_vectorized.resize(dofs_per_component); // zero content - solution_renumbered_vectorized.fill( - typename internal::FEPointEvaluation::EvaluatorTypeTraits< - dim, - n_components, - VectorizedArrayType>::value_type()); - - const auto unit_points = - mapping_info->get_unit_points(current_cell_index, current_face_number); - const auto &mapping_data = - mapping_info->get_mapping_data(current_cell_index, current_face_number); - - const std::size_t n_points = unit_points.size(); - const std::size_t n_lanes = - internal::VectorizedArrayTrait::width; - const std::size_t stride = - internal::FEPointEvaluation::VectorizedPointsTypeTraits::stride; + solution_renumbered_vectorized.fill(vectorized_value_type()); // loop over quadrature batches qb / points q - for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes) + const unsigned int n_shapes = poly.size(); + for (unsigned int qb = 0, q = 0; q < n_q_points_scalar; + ++qb, q += n_lanes_internal) { const bool incomplete_last_batch = - (qb == (n_q_points - 1)) && (n_filled_lanes_last_batch > 0); + q + n_lanes_user_interface > n_q_points_scalar; - typename internal::ProductTypeNoPoint::type value = - {}; - Tensor<1, - dim, - typename internal::ProductTypeNoPoint::type> - gradient; + vectorized_value_type value = {}; + Tensor<1, dim, vectorized_value_type> gradient; if (integration_flags & EvaluationFlags::values) { // zero out lanes of incomplete last quadrature point batch if (incomplete_last_batch) { - for (unsigned int v = n_filled_lanes_last_batch; v < n_lanes; ++v) - internal::FEPointEvaluation::EvaluatorTypeTraits< - dim, - n_components, - Number>::set_zero_value(values[qb], v); + const unsigned int n_filled_lanes_last_batch = + n_q_points_scalar % n_lanes_internal; + for (unsigned int v = n_filled_lanes_last_batch; + v < n_lanes_internal; + ++v) + ETT::set_zero_value(values[qb], v); } - for (unsigned int v = 0; v < stride && q + v < n_points; ++v) - internal::FEPointEvaluation:: - EvaluatorTypeTraits::get_value( - value, v, values[qb * stride + v]); + + for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v) + ETT::get_value(value, v, values[qb * stride + v]); } if (integration_flags & EvaluationFlags::gradients) { - if (std::is_same::value) + // zero out lanes of incomplete last quadrature point batch + if (incomplete_last_batch) { - for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v) - { - gradients[q + v] = - apply_transformation(mapping_data.inverse_jacobians[q + v], - gradients[q + v]); - internal::FEPointEvaluation::EvaluatorTypeTraits< - dim, - n_components, - Number>::get_gradient(gradient, v, gradients[q + v]); - } + const unsigned int n_filled_lanes_last_batch = + n_q_points_scalar % n_lanes_internal; + for (unsigned int v = n_filled_lanes_last_batch; + v < n_lanes_internal; + ++v) + ETT::set_zero_gradient(gradients[qb], v); } - else + + for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v) { - // convert to vectorized format - DerivativeForm<1, spacedim, dim, Number> - vectorized_inverse_jacobians; - for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int s = 0; s < spacedim; ++s) - internal::VectorizedArrayTrait::get( - vectorized_inverse_jacobians[s][d], v) = - mapping_data.inverse_jacobians[q + v][s][d]; - - // zero out lanes of incomplete last quadrature point batch - if (incomplete_last_batch) - { - for (unsigned int v = n_filled_lanes_last_batch; v < n_lanes; - ++v) - internal::FEPointEvaluation::EvaluatorTypeTraits< - dim, - n_components, - Number>::set_zero_gradient(gradients[qb], v); - } - - gradients[qb] = apply_transformation(vectorized_inverse_jacobians, - gradients[qb]); - internal::FEPointEvaluation:: - EvaluatorTypeTraits::get_gradient( - gradient, 0, gradients[qb]); + const unsigned int offset = qb * stride + v; + ETT::get_gradient( + gradient, + v, + apply_transformation(inverse_jacobian_ptr[offset], + gradients[offset])); } } // compute - const unsigned int n_shapes = poly.size(); if (polynomials_are_hat_functions) - { - // convert to vectorized format - Point vectorized_points; - for (unsigned int j = 0; j < n_lanes && q + j < n_points; ++j) - for (unsigned int d = 0; d < dim; ++d) - vectorized_points[d][j] = unit_points[q + j][d]; - - internal::integrate_add_tensor_product_value_and_gradient_linear( - poly, - value, - gradient, - solution_renumbered_vectorized, - vectorized_points); - } + internal::integrate_add_tensor_product_value_and_gradient_linear( + poly, + value, + gradient, + solution_renumbered_vectorized, + unit_point_ptr[qb]); else internal::integrate_add_tensor_product_value_and_gradient_shapes< dim, VectorizedArrayType, - typename internal::ProductTypeNoPoint::type>( - make_array_view(shapes.begin() + qb * n_shapes, - shapes.begin() + (qb * n_shapes + n_shapes)), - n_shapes, - value, - gradient, - solution_renumbered_vectorized); + vectorized_value_type>(make_array_view(shapes, + qb * n_shapes, + n_shapes), + n_shapes, + value, + gradient, + solution_renumbered_vectorized); } // add between the lanes and write into the result @@ -1910,13 +1811,8 @@ FEPointEvaluation::integrate_fast( for (unsigned int i = 0; i < dofs_per_component; ++i) { VectorizedArrayType result; - internal::FEPointEvaluation::EvaluatorTypeTraits< - dim, - n_components, - VectorizedArrayType>::write_value(result, - comp, - solution_renumbered_vectorized[i]); - for (unsigned int lane = n_lanes / 2; lane > 0; lane /= 2) + ETT::write_value(result, comp, solution_renumbered_vectorized[i]); + for (unsigned int lane = n_lanes_internal / 2; lane > 0; lane /= 2) for (unsigned int j = 0; j < lane; ++j) result[j] += result[lane + j]; solution_values[renumber[comp * dofs_per_component + i]] = result[0]; @@ -1938,7 +1834,6 @@ FEPointEvaluation::integrate_slow( std::fill(solution_values.begin(), solution_values.end(), 0.0); const std::size_t n_points = fe_values->get_quadrature().size(); - const std::size_t n_lanes = internal::VectorizedArrayTrait::width; if (integration_flags & EvaluationFlags::values) { @@ -1948,27 +1843,22 @@ FEPointEvaluation::integrate_slow( for (unsigned int d = 0; d < n_components; ++d) if (nonzero_shape_function_component[i][d] && (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) for (unsigned int v = 0; - v < (q + n_lanes > n_points ? n_filled_lanes_last_batch : - n_lanes); + v < n_lanes_user_interface && q + v < n_points; ++v) - solution_values[i] += - fe_values->shape_value(i, q + v) * - internal::FEPointEvaluation:: - EvaluatorTypeTraits::access( - values[qb], v, d); + solution_values[i] += fe_values->shape_value(i, q + v) * + ETT::access(values[qb], v, d); else if (nonzero_shape_function_component[i][d]) - for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) for (unsigned int v = 0; - v < (q + n_lanes > n_points ? n_filled_lanes_last_batch : - n_lanes); + v < n_lanes_user_interface && q + v < n_points; ++v) solution_values[i] += fe_values->shape_value_component(i, q + v, d) * - internal::FEPointEvaluation:: - EvaluatorTypeTraits::access( - values[qb], v, d); + ETT::access(values[qb], v, d); } } @@ -1980,27 +1870,22 @@ FEPointEvaluation::integrate_slow( for (unsigned int d = 0; d < n_components; ++d) if (nonzero_shape_function_component[i][d] && (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) for (unsigned int v = 0; - v < (q + n_lanes > n_points ? n_filled_lanes_last_batch : - n_lanes); + v < n_lanes_user_interface && q + v < n_points; ++v) - solution_values[i] += - fe_values->shape_grad(i, q + v) * - internal::FEPointEvaluation:: - EvaluatorTypeTraits::access( - gradients[qb], v, d); + solution_values[i] += fe_values->shape_grad(i, q + v) * + ETT::access(gradients[qb], v, d); else if (nonzero_shape_function_component[i][d]) - for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes) + for (unsigned int qb = 0, q = 0; q < n_points; + ++qb, q += n_lanes_user_interface) for (unsigned int v = 0; - v < (q + n_lanes > n_points ? n_filled_lanes_last_batch : - n_lanes); + v < n_lanes_user_interface && q + v < n_points; ++v) solution_values[i] += fe_values->shape_grad_component(i, q + v, d) * - internal::FEPointEvaluation:: - EvaluatorTypeTraits::access( - gradients[qb], v, d); + ETT::access(gradients[qb], v, d); } } } @@ -2110,20 +1995,12 @@ inline DerivativeForm<1, dim, spacedim, Number> FEPointEvaluation::jacobian( const unsigned int point_index) const { - const auto &mapping_data = - mapping_info->get_mapping_data(current_cell_index, current_face_number); - - const unsigned int n_lanes = internal::VectorizedArrayTrait::width; - DerivativeForm<1, dim, spacedim, Number> vectorized_jacobian; - for (unsigned int v = 0; - v < n_lanes && point_index * n_lanes + v < mapping_data.jacobians.size(); - ++v) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int s = 0; s < spacedim; ++s) - internal::VectorizedArrayTrait::get(vectorized_jacobian[d][s], - v) = - mapping_data.jacobians[point_index * n_lanes + v][d][s]; - return vectorized_jacobian; + AssertIndexRange(point_index, n_q_points); + Assert(jacobian_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_jacobians")); + return jacobian_ptr[point_index]; } @@ -2133,21 +2010,12 @@ inline DerivativeForm<1, spacedim, dim, Number> FEPointEvaluation::inverse_jacobian( const unsigned int point_index) const { - const auto &mapping_data = - mapping_info->get_mapping_data(current_cell_index, current_face_number); - - const unsigned int n_lanes = internal::VectorizedArrayTrait::width; - DerivativeForm<1, dim, spacedim, Number> vectorized_inverse_jacobian; - for (unsigned int v = 0; - v < n_lanes && - point_index * n_lanes + v < mapping_data.inverse_jacobians.size(); - ++v) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int s = 0; s < spacedim; ++s) - internal::VectorizedArrayTrait::get( - vectorized_inverse_jacobian[d][s], v) = - mapping_data.inverse_jacobians[point_index * n_lanes + v][d][s]; - return vectorized_inverse_jacobian; + AssertIndexRange(point_index, n_q_points); + Assert(inverse_jacobian_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_inverse_jacobians")); + return inverse_jacobian_ptr[point_index]; } @@ -2157,38 +2025,27 @@ inline Number FEPointEvaluation::JxW( const unsigned int point_index) const { - const auto &mapping_data = - mapping_info->get_mapping_data(current_cell_index, current_face_number); - - const unsigned int n_lanes = internal::VectorizedArrayTrait::width; - Number vectorized_JxW = 0; - for (unsigned int v = 0; v < n_lanes && point_index * n_lanes + v < - mapping_data.JxW_values.size(); - ++v) - internal::VectorizedArrayTrait::get(vectorized_JxW, v) = - mapping_data.JxW_values[point_index * n_lanes + v]; - return vectorized_JxW; + AssertIndexRange(point_index, n_q_points); + Assert(JxW_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_JxW_values")); + return JxW_ptr[point_index]; } + template inline Tensor<1, spacedim, Number> FEPointEvaluation::normal_vector( const unsigned int point_index) const { - const auto &mapping_data = - mapping_info->get_mapping_data(current_cell_index, current_face_number); - - const unsigned int n_lanes = internal::VectorizedArrayTrait::width; - Tensor<1, spacedim, Number> vectorized_normal_vectors; - for (unsigned int v = 0; v < n_lanes && point_index * n_lanes + v < - mapping_data.normal_vectors.size(); - ++v) - for (unsigned int s = 0; s < spacedim; ++s) - internal::VectorizedArrayTrait::get(vectorized_normal_vectors[s], - v) = - mapping_data.normal_vectors[point_index * n_lanes + v][s]; - return vectorized_normal_vectors; + AssertIndexRange(point_index, n_q_points); + Assert(normal_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_normal_vectors")); + return normal_ptr[point_index]; } @@ -2198,19 +2055,12 @@ inline Point FEPointEvaluation::real_point( const unsigned int point_index) const { - const auto &mapping_data = - mapping_info->get_mapping_data(current_cell_index, current_face_number); - - const unsigned int n_lanes = internal::VectorizedArrayTrait::width; - Point vectorized_real_point; - for (unsigned int v = 0; - v < n_lanes && - point_index * n_lanes + v < mapping_data.quadrature_points.size(); - ++v) - for (unsigned int s = 0; s < spacedim; ++s) - internal::VectorizedArrayTrait::get(vectorized_real_point[s], v) = - mapping_data.quadrature_points[point_index * n_lanes + v][s]; - return vectorized_real_point; + AssertIndexRange(point_index, n_q_points); + Assert(real_point_ptr != nullptr, + internal::FEPointEvaluation:: + ExcFEPointEvaluationAccessToUninitializedMappingField( + "update_quadrature_points")); + return real_point_ptr[point_index]; } @@ -2220,17 +2070,13 @@ inline Point FEPointEvaluation::unit_point( const unsigned int point_index) const { - const auto unit_points = - mapping_info->get_unit_points(current_cell_index, current_face_number); - const unsigned int n_lanes = internal::VectorizedArrayTrait::width; - Point vectorized_unit_point; - for (unsigned int v = 0; - v < n_lanes && point_index * n_lanes + v < unit_points.size(); - ++v) - for (unsigned int d = 0; d < dim; ++d) - internal::VectorizedArrayTrait::get(vectorized_unit_point[d], v) = - unit_points[point_index * n_lanes + v][d]; - return vectorized_unit_point; + AssertIndexRange(point_index, n_q_points); + Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!")); + Point unit_point; + for (unsigned int d = 0; d < dim; ++d) + unit_point[d] = internal::VectorizedArrayTrait::get_from_vectorized( + unit_point_ptr[point_index / stride][d], point_index % stride); + return unit_point; } diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 29af641639..99b69eb97f 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -2989,7 +2989,7 @@ namespace internal template inline void compute_values_of_array( - ArrayView> & shapes, + ArrayView> shapes, const std::vector> &poly, const Point & p) { diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h index 7a106c18ac..418cd25f78 100644 --- a/include/deal.II/non_matching/mapping_info.h +++ b/include/deal.II/non_matching/mapping_info.h @@ -38,17 +38,174 @@ DEAL_II_NAMESPACE_OPEN namespace NonMatching { + namespace internal + { + template + class ComputeMappingDataHelper + { + using MappingData = + dealii::internal::FEValuesImplementation::MappingRelatedData; + + public: + static UpdateFlags + required_update_flags( + const SmartPointer> mapping, + const UpdateFlags & update_flags) + { + return mapping->requires_update_flags(update_flags); + } + + static void + compute_mapping_data_for_quadrature( + const SmartPointer> mapping, + const UpdateFlags & update_flags_mapping, + const typename Triangulation::cell_iterator &cell, + CellSimilarity::Similarity &cell_similarity, + const Quadrature & quadrature, + std::shared_ptr::InternalDataBase> + internal_mapping_data, + MappingData &mapping_data) + { + mapping_data.initialize(quadrature.size(), update_flags_mapping); + + // reuse internal_mapping_data for MappingQ to avoid memory allocations + if (const MappingQ *mapping_q = + dynamic_cast *>(&(*mapping))) + { + (void)mapping_q; + auto &data = + dynamic_cast::InternalData &>( + *internal_mapping_data); + data.initialize(update_flags_mapping, + quadrature, + quadrature.size()); + } + else + { + internal_mapping_data = + mapping->get_data(update_flags_mapping, quadrature); + } + + cell_similarity = mapping->fill_fe_values(cell, + cell_similarity, + quadrature, + *internal_mapping_data, + mapping_data); + } + + + + static void + compute_mapping_data_for_immersed_surface_quadrature( + const SmartPointer> mapping, + const UpdateFlags & update_flags_mapping, + const typename Triangulation::cell_iterator &cell, + const ImmersedSurfaceQuadrature & quadrature, + std::shared_ptr::InternalDataBase> + internal_mapping_data, + MappingData &mapping_data) + { + mapping_data.initialize(quadrature.size(), update_flags_mapping); + + // reuse internal_mapping_data for MappingQ to avoid memory allocations + if (dynamic_cast *>(&(*mapping))) + { + auto &data = + dynamic_cast::InternalData &>( + *internal_mapping_data); + data.initialize(update_flags_mapping, + quadrature, + quadrature.size()); + } + else + { + internal_mapping_data = + mapping->get_data(update_flags_mapping, quadrature); + } + + mapping->fill_fe_immersed_surface_values(cell, + quadrature, + *internal_mapping_data, + mapping_data); + } + + + + static void + compute_mapping_data_for_face_quadrature( + const SmartPointer> mapping, + const UpdateFlags & update_flags_mapping, + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature & quadrature, + std::shared_ptr::InternalDataBase> + internal_mapping_data, + MappingData &mapping_data) + { + mapping_data.initialize(quadrature.size(), update_flags_mapping); + + // reuse internal_mapping_data for MappingQ to avoid memory allocations + if (const MappingQ *mapping_q = + dynamic_cast *>(&(*mapping))) + { + auto &data = + dynamic_cast::InternalData &>( + *internal_mapping_data); + data.initialize_face(update_flags_mapping, + QProjector::project_to_oriented_face( + ReferenceCells::get_hypercube(), + quadrature, + face_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no)), + quadrature.size()); + + mapping_q->fill_mapping_data_for_face_quadrature( + cell, face_no, quadrature, *internal_mapping_data, mapping_data); + } + else + { + auto internal_mapping_data = + mapping->get_face_data(update_flags_mapping, + hp::QCollection(quadrature)); + + mapping->fill_fe_face_values(cell, + face_no, + hp::QCollection(quadrature), + *internal_mapping_data, + mapping_data); + } + } + }; + } // namespace internal + /** * This class provides the mapping information computation and mapping data * storage to be used together with FEPointEvaluation. + * + * MappingInfo is vectorized across quadrature points, which means data can be + * provided in vectorized format. + * + * Two different modes are available: partially vectorized and fully + * vectorized across quadrature points. Partially vectorized (scalar Number + * template argument) means the computed mapping data is provided in scalar + * format and only unit points are provided vectorized (for seamless + * interaction with FEPointEvaluation). Fully vectorized (vectorized Number + * template argument, e.g. VectorizedArray) provides both mapping data + * and unit points in vectorized format. The Number template parameter of + * MappingInfo and FEPointEvaluation has to be identical. */ - template + template class MappingInfo : public Subscriptor { public: - using MappingData = - dealii::internal::FEValuesImplementation::MappingRelatedData; + /** + * The VectorizedArray type the unit points are stored in inside this class. + */ + using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait< + Number>::vectorized_value_type; /** * Constructor. @@ -139,37 +296,46 @@ namespace NonMatching const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int); /** - * Getter function for current unit points. - * - * @p cell_index and @p face_number are the indices - * into the compressed, CRS like data storage of unit points. - * - * If you have initialized this object with reinit_cells() you can access - * the stored unit points of the cell with the respective @p cell_index - * (and the default argument for @p face_number). - * - * If you have initialized this object with reinit_faces() you can access - * the stored unit points of the faces on the cell with the respective @p cell_index - * and the respective local @p face_number. - * - * If you have initialized this object with reinit() you can access the - * stored unit points of a single cell with the default arguments. - * - * The correct state of this object is checked in this call (in debug mode). + * Getter function for unit points. The offset can be obtained with + * compute_unit_point_index_offset(). + */ + const Point * + get_unit_point(const unsigned int offset) const; + + /** + * Getter function for Jacobians. The offset can be obtained with + * compute_data_index_offset(). + */ + const DerivativeForm<1, dim, spacedim, Number> * + get_jacobian(const unsigned int offset) const; + + /** + * Getter function for inverse Jacobians. The offset can be obtained with + * compute_data_index_offset(). */ - const ArrayView> - get_unit_points( - const unsigned int cell_index = numbers::invalid_unsigned_int, - const unsigned int face_number = numbers::invalid_unsigned_int) const; + const DerivativeForm<1, spacedim, dim, Number> * + get_inverse_jacobian(const unsigned int offset) const; /** - * Getter function for computed mapping data. This function accesses - * internal data and is therefore not a stable interface. + * Getter function for normal vectors. The offset can be obtained with + * compute_data_index_offset(). */ - const MappingData & - get_mapping_data( - const unsigned int cell_index = numbers::invalid_unsigned_int, - const unsigned int face_number = numbers::invalid_unsigned_int) const; + const Tensor<1, spacedim, Number> * + get_normal_vector(const unsigned int offset) const; + + /** + * Getter function for Jacobian times quadrature weight (JxW). The offset + * can be obtained with compute_data_index_offset(). + */ + const Number * + get_JxW(const unsigned int offset) const; + + /** + * Getter function for real points. The offset can be obtained with + * compute_data_index_offset(). + */ + const Point * + get_real_point(const unsigned int offset) const; /** * Getter function for underlying mapping. @@ -183,79 +349,127 @@ namespace NonMatching UpdateFlags get_update_flags() const; + /** + * Getter function for the mapping update flags. + */ + UpdateFlags + get_update_flags_mapping() const; + /** * Connects to is_reinitialized(). */ boost::signals2::connection connect_is_reinitialized(const std::function &set_is_reinitialized); + /** + * Compute the unit points index offset for the current cell/face. + */ + unsigned int + compute_unit_point_index_offset(const unsigned int cell_index, + const unsigned int face_number) const; + + /** + * Compute the data index offset for the current cell/face. + */ + unsigned int + compute_data_index_offset(const unsigned int cell_index, + const unsigned int face_number) const; + + /** + * Get number of unvectorized quadrature points. + */ + unsigned int + get_n_q_points_unvectorized(const unsigned int cell_index, + const unsigned int face_number) const; + private: + using MappingData = + dealii::internal::FEValuesImplementation::MappingRelatedData; + /** - * Enum class for reinitialized states. + * Compute number of quadrature point batches depending on NumberType. */ - enum class State - { - invalid, - single_cell, - cell_vector, - faces_on_cells_in_vector - }; + template + unsigned int + compute_n_q_points(const unsigned int n_q_points_unvectorized); /** - * Enum class that stores the currently initialized state - * upon the last call of reinit(). + * Resize the unit_point data field. */ - State state; + void + resize_unit_points(const unsigned int n_unit_point_batches); /** - * Compute the mapping related data for the given @p mapping, - * @p cell and @p unit_points that is required by the FEPointEvaluation - * class. + * Resize the mapping data fields. */ void - compute_mapping_data_for_quadrature( - const typename Triangulation::cell_iterator &cell, - CellSimilarity::Similarity &cell_similarity, - const Quadrature & quadrature, - MappingData & mapping_data); + resize_data_fields(const unsigned int n_data_point_batches); /** - * Compute the mapping related data for the given @p mapping, - * @p cell and @p quadrature that is required by the FEPointEvaluation - * class. + * Store the unit points. */ void - compute_mapping_data_for_immersed_surface_quadrature( - const typename Triangulation::cell_iterator &cell, - const ImmersedSurfaceQuadrature & quadrature, - MappingData & mapping_data); + store_unit_points(const unsigned int unit_points_index_offset, + const unsigned int n_q_points, + const unsigned int n_q_points_unvectorized, + const std::vector> &points); /** - * Compute the mapping related data for the given @p mapping, @p cell, - * @p face_no and @p quadrature that is required by the FEPointEvaluation - * class. + * Store the requested mapping data. */ void - compute_mapping_data_for_face_quadrature( - const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const Quadrature & quadrature, - MappingData & mapping_data); + store_mapping_data(const unsigned int unit_points_index_offset, + const unsigned int n_q_points, + const unsigned int n_q_points_unvectorized, + const MappingData &mapping_data); + + /** + * Compute the compressed cell index. + */ + unsigned int + compute_compressed_cell_index(const unsigned int cell_index) const; + + /** + * Compute the geometry index offset of the current cell/face. + */ + unsigned int + compute_geometry_index_offset(const unsigned int cell_index, + const unsigned int face_number) const; + + /** + * Enum class for reinitialized states. + */ + enum class State + { + invalid, + single_cell, + cell_vector, + faces_on_cells_in_vector + }; + + /** + * Enum class that stores the currently initialized state + * upon the last call of reinit(). + */ + State state; /** * The reference points specified at reinit(). + * + * Indexed by @p unit_points_index. */ - std::vector> unit_points; + AlignedVector> unit_points; /** - * Offset to point to the first unit point of a cell/face + * Offset to point to the first unit point of a cell/face. */ - std::vector unit_points_index; + AlignedVector unit_points_index; /** * A pointer to the internal data of the underlying mapping. */ - std::unique_ptr::InternalDataBase> + std::shared_ptr::InternalDataBase> internal_mapping_data; /** @@ -274,10 +488,54 @@ namespace NonMatching UpdateFlags update_flags_mapping; /** - * The internal data container for mapping information. The implementation - * is subject to future changes. + * Stores the index offset into the arrays @p JxW_values, @p jacobians, + * @p inverse_jacobians and @p normal_vectors. */ - std::vector mapping_data; + AlignedVector data_index_offsets; + + /** + * The storage of the Jacobian determinant times the quadrature weight on + * quadrature points. + * + * Indexed by @p data_index_offsets. + */ + AlignedVector JxW_values; + + /** + * Stores the normal vectors. + * + * Indexed by @p data_index_offsets. + */ + AlignedVector> normal_vectors; + + /** + * The storage of contravariant transformation on quadrature points, i.e., + * the Jacobians of the transformation from the unit to the real cell. + * + * Indexed by @p data_index_offsets. + */ + AlignedVector> jacobians; + + /** + * The storage of covariant transformation on quadrature points, i.e., + * the inverse Jacobians of the transformation from the + * unit to the real cell. + * + * Indexed by @p data_index_offsets. + */ + AlignedVector> inverse_jacobians; + + /** + * The mapped real points. + * + * Indexed by @p data_index_offsets. + */ + AlignedVector> real_points; + + /** + * Number of unvectorized unit points per geometric entity (cell/face). + */ + std::vector n_q_points_unvectorized; /** * Offset to point to the first element of a cell in internal data @@ -306,15 +564,16 @@ namespace NonMatching // ----------------------- template functions ---------------------- - template - MappingInfo::MappingInfo(const Mapping &mapping, - const UpdateFlags update_flags) + template + MappingInfo::MappingInfo( + const Mapping &mapping, + const UpdateFlags update_flags) : mapping(&mapping) , update_flags(update_flags) + , update_flags_mapping(update_default) { - update_flags_mapping = update_default; // translate update flags - if (update_flags & update_jacobians || update_flags & update_JxW_values) + if (update_flags & update_jacobians) update_flags_mapping |= update_jacobians; if (update_flags & update_JxW_values) update_flags_mapping |= update_JxW_values; @@ -327,6 +586,10 @@ namespace NonMatching // always save quadrature points for now update_flags_mapping |= update_quadrature_points; + update_flags_mapping = + internal::ComputeMappingDataHelper::required_update_flags( + this->mapping, update_flags_mapping); + // construct internal_mapping_data for MappingQ to be able to reuse it in // reinit() calls to avoid memory allocations if (const MappingQ *mapping_q = @@ -340,9 +603,9 @@ namespace NonMatching - template + template void - MappingInfo::reinit( + MappingInfo::reinit( const typename Triangulation::cell_iterator &cell, const std::vector> & unit_points_in) { @@ -351,9 +614,9 @@ namespace NonMatching - template + template void - MappingInfo::reinit( + MappingInfo::reinit( const typename Triangulation::cell_iterator &cell, const ArrayView> & unit_points_in) { @@ -364,21 +627,48 @@ namespace NonMatching - template + template void - MappingInfo::reinit( + MappingInfo::reinit( const typename Triangulation::cell_iterator &cell, const Quadrature & quadrature) { - unit_points = quadrature.get_points(); - - mapping_data.resize(1); - CellSimilarity::Similarity cell_similarity = - CellSimilarity::Similarity::none; - compute_mapping_data_for_quadrature(cell, - cell_similarity, - quadrature, - mapping_data[0]); + n_q_points_unvectorized.resize(1); + n_q_points_unvectorized[0] = quadrature.size(); + + const unsigned int n_q_points = + compute_n_q_points(n_q_points_unvectorized[0]); + + const unsigned int n_q_points_data = + compute_n_q_points(n_q_points_unvectorized[0]); + + // resize data vectors + resize_unit_points(n_q_points); + resize_data_fields(n_q_points_data); + + // store unit points + store_unit_points(0, + n_q_points, + n_q_points_unvectorized[0], + quadrature.get_points()); + + // compute mapping data + MappingData mapping_data; + CellSimilarity::Similarity cell_similarity = CellSimilarity::none; + internal::ComputeMappingDataHelper:: + compute_mapping_data_for_quadrature(mapping, + update_flags_mapping, + cell, + cell_similarity, + quadrature, + internal_mapping_data, + mapping_data); + + // store mapping data + store_mapping_data(0, + n_q_points_data, + n_q_points_unvectorized[0], + mapping_data); state = State::single_cell; is_reinitialized(); @@ -386,10 +676,10 @@ namespace NonMatching - template + template template void - MappingInfo::reinit_cells( + MappingInfo::reinit_cells( const ContainerType & cell_iterator_range, const std::vector>> &unit_points_vector, const unsigned int n_unfiltered_cells) @@ -409,10 +699,10 @@ namespace NonMatching - template + template template void - MappingInfo::reinit_cells( + MappingInfo::reinit_cells( const ContainerType & cell_iterator_range, const std::vector> &quadrature_vector, const unsigned int n_unfiltered_cells) @@ -425,37 +715,70 @@ namespace NonMatching std::distance(cell_iterator_range.begin(), cell_iterator_range.end())); + n_q_points_unvectorized.reserve(n_cells); + // fill unit points index offset vector unit_points_index.reserve(n_cells + 1); unit_points_index.push_back(0); + data_index_offsets.reserve(n_cells + 1); + data_index_offsets.push_back(0); for (const auto &quadrature : quadrature_vector) - unit_points_index.push_back(unit_points_index.back() + quadrature.size()); + { + const unsigned int n_points = quadrature.size(); + n_q_points_unvectorized.push_back(n_points); + + const unsigned int n_q_points = + compute_n_q_points(n_points); + unit_points_index.push_back(unit_points_index.back() + n_q_points); + + const unsigned int n_q_points_data = + compute_n_q_points(n_points); + data_index_offsets.push_back(data_index_offsets.back() + + n_q_points_data); + } const unsigned int n_unit_points = unit_points_index.back(); + const unsigned int n_data_points = data_index_offsets.back(); - unit_points.resize(n_unit_points); - mapping_data.resize(n_cells); + // resize data vectors + resize_unit_points(n_unit_points); + resize_data_fields(n_data_points); if (do_cell_index_compression) cell_index_to_compressed_cell_index.resize(n_unfiltered_cells, numbers::invalid_unsigned_int); + + MappingData mapping_data; CellSimilarity::Similarity cell_similarity = CellSimilarity::Similarity::none; unsigned int cell_index = 0; for (const auto &cell : cell_iterator_range) { - auto it = unit_points.begin() + unit_points_index[cell_index]; - for (const auto &unit_point : - quadrature_vector[cell_index].get_points()) - { - *it = unit_point; - ++it; - } - - compute_mapping_data_for_quadrature(cell, - cell_similarity, - quadrature_vector[cell_index], - mapping_data[cell_index]); + // store unit points + const unsigned int n_q_points = compute_n_q_points( + n_q_points_unvectorized[cell_index]); + store_unit_points(unit_points_index[cell_index], + n_q_points, + n_q_points_unvectorized[cell_index], + quadrature_vector[cell_index].get_points()); + + // compute mapping data + internal::ComputeMappingDataHelper:: + compute_mapping_data_for_quadrature(mapping, + update_flags_mapping, + cell, + cell_similarity, + quadrature_vector[cell_index], + internal_mapping_data, + mapping_data); + + // store mapping data + const unsigned int n_q_points_data = + compute_n_q_points(n_q_points_unvectorized[cell_index]); + store_mapping_data(data_index_offsets[cell_index], + n_q_points_data, + n_q_points_unvectorized[cell_index], + mapping_data); if (do_cell_index_compression) cell_index_to_compressed_cell_index[cell->active_cell_index()] = @@ -470,10 +793,10 @@ namespace NonMatching - template + template template void - MappingInfo::reinit_surface( + MappingInfo::reinit_surface( const IteratorRange & cell_iterator_range, const std::vector> &quadrature_vector, const unsigned int n_unfiltered_cells) @@ -489,35 +812,70 @@ namespace NonMatching std::distance(cell_iterator_range.begin(), cell_iterator_range.end())); + n_q_points_unvectorized.reserve(n_cells); + // fill unit points index offset vector unit_points_index.reserve(n_cells + 1); unit_points_index.push_back(0); + data_index_offsets.reserve(n_cells + 1); + data_index_offsets.push_back(0); for (const auto &quadrature : quadrature_vector) - unit_points_index.push_back(unit_points_index.back() + - quadrature.get_points().size()); + { + const unsigned int n_points = quadrature.size(); + n_q_points_unvectorized.push_back(n_points); + + const unsigned int n_q_points = + compute_n_q_points(n_points); + unit_points_index.push_back(unit_points_index.back() + n_q_points); + + const unsigned int n_q_points_data = + compute_n_q_points(n_points); + data_index_offsets.push_back(data_index_offsets.back() + + n_q_points_data); + } const unsigned int n_unit_points = unit_points_index.back(); + const unsigned int n_data_points = data_index_offsets.back(); - unit_points.resize(n_unit_points); - mapping_data.resize(n_cells); + // resize data vectors + resize_unit_points(n_unit_points); + resize_data_fields(n_data_points); if (do_cell_index_compression) cell_index_to_compressed_cell_index.resize(n_unfiltered_cells, numbers::invalid_unsigned_int); + + MappingData mapping_data; unsigned int cell_index = 0; for (const auto &cell : cell_iterator_range) { const auto &quadrature = quadrature_vector[cell_index]; - auto it = unit_points.begin() + unit_points_index[cell_index]; - for (const auto &unit_point : quadrature.get_points()) - { - *it = unit_point; - ++it; - } - - compute_mapping_data_for_immersed_surface_quadrature( - cell, quadrature, mapping_data[cell_index]); + // store unit points + const unsigned int n_q_points = compute_n_q_points( + n_q_points_unvectorized[cell_index]); + store_unit_points(unit_points_index[cell_index], + n_q_points, + n_q_points_unvectorized[cell_index], + quadrature_vector[cell_index].get_points()); + + // compute mapping data + internal::ComputeMappingDataHelper:: + compute_mapping_data_for_immersed_surface_quadrature( + mapping, + update_flags_mapping, + cell, + quadrature, + internal_mapping_data, + mapping_data); + + // store mapping data + const unsigned int n_q_points_data = + compute_n_q_points(n_q_points_unvectorized[cell_index]); + store_mapping_data(data_index_offsets[cell_index], + n_q_points_data, + n_q_points_unvectorized[cell_index], + mapping_data); if (do_cell_index_compression) cell_index_to_compressed_cell_index[cell->active_cell_index()] = @@ -532,10 +890,10 @@ namespace NonMatching - template + template template void - MappingInfo::reinit_faces( + MappingInfo::reinit_faces( const IteratorRange & cell_iterator_range, const std::vector>> &quadrature_vector, const unsigned int n_unfiltered_cells) @@ -559,10 +917,14 @@ namespace NonMatching ++cell_index; } + n_q_points_unvectorized.reserve(n_faces); + // fill unit points index offset vector unit_points_index.resize(n_faces + 1); + data_index_offsets.resize(n_faces + 1); cell_index = 0; unsigned int n_unit_points = 0; + unsigned int n_data_points = 0; for (const auto &cell : cell_iterator_range) { for (const auto &f : cell->face_indices()) @@ -570,14 +932,26 @@ namespace NonMatching const unsigned int current_face_index = cell_index_offset[cell_index] + f; - unit_points_index[current_face_index] = n_unit_points; - n_unit_points += - quadrature_vector[cell_index][f].get_points().size(); + unit_points_index[current_face_index] = n_unit_points; + data_index_offsets[current_face_index] = n_data_points; + + const unsigned int n_points = + quadrature_vector[cell_index][f].size(); + n_q_points_unvectorized.push_back(n_points); + + const unsigned int n_q_points = + compute_n_q_points(n_points); + n_unit_points += n_q_points; + + const unsigned int n_q_points_data = + compute_n_q_points(n_points); + n_data_points += n_q_points_data; } ++cell_index; } - unit_points_index[n_faces] = n_unit_points; + unit_points_index[n_faces] = n_unit_points; + data_index_offsets[n_faces] = n_data_points; // compress indices if (do_cell_index_compression) @@ -585,8 +959,11 @@ namespace NonMatching numbers::invalid_unsigned_int); // fill unit points and mapping data for every face of all cells - unit_points.resize(n_unit_points); - mapping_data.resize(n_faces); + // resize data vectors + resize_unit_points(n_unit_points); + resize_data_fields(n_data_points); + + MappingData mapping_data; cell_index = 0; QProjector q_projector; for (const auto &cell : cell_iterator_range) @@ -611,16 +988,30 @@ namespace NonMatching const unsigned int current_face_index = cell_index_offset[cell_index] + f; - auto it = - unit_points.begin() + unit_points_index[current_face_index]; - for (const auto &unit_point : unit_points_on_cell) - { - *it = unit_point; - ++it; - } - - compute_mapping_data_for_face_quadrature( - cell, f, quadrature_on_face, mapping_data[current_face_index]); + // store unit points + const unsigned int n_q_points = + compute_n_q_points( + n_q_points_unvectorized[current_face_index]); + store_unit_points(unit_points_index[current_face_index], + n_q_points, + n_q_points_unvectorized[current_face_index], + quadrature_on_cell.get_points()); + + internal::ComputeMappingDataHelper:: + compute_mapping_data_for_face_quadrature(mapping, + update_flags_mapping, + cell, + f, + quadrature_on_face, + internal_mapping_data, + mapping_data); + + const unsigned int n_q_points_data = compute_n_q_points( + n_q_points_unvectorized[current_face_index]); + store_mapping_data(data_index_offsets[current_face_index], + n_q_points_data, + n_q_points_unvectorized[current_face_index], + mapping_data); } if (do_cell_index_compression) cell_index_to_compressed_cell_index[cell->active_cell_index()] = @@ -635,9 +1026,9 @@ namespace NonMatching - template - inline const ArrayView> - MappingInfo::get_unit_points( + template + unsigned int + MappingInfo::get_n_q_points_unvectorized( const unsigned int cell_index, const unsigned int face_number) const { @@ -647,287 +1038,327 @@ namespace NonMatching Assert(state == State::single_cell, ExcMessage( "This mapping info is not reinitialized for a single cell!")); - return unit_points; + return n_q_points_unvectorized[0]; + } + else + { + return n_q_points_unvectorized[compute_geometry_index_offset( + cell_index, face_number)]; } - else if (face_number == numbers::invalid_unsigned_int) + } + + + + template + template + unsigned int + MappingInfo::compute_n_q_points( + const unsigned int n_q_points_unvectorized) + { + const unsigned int n_lanes = + dealii::internal::VectorizedArrayTrait::width; + const unsigned int n_filled_lanes_last_batch = + n_q_points_unvectorized % n_lanes; + unsigned int n_q_points = n_q_points_unvectorized / n_lanes; + if (n_filled_lanes_last_batch > 0) + ++n_q_points; + return n_q_points; + } + + + + template + unsigned int + MappingInfo::compute_geometry_index_offset( + const unsigned int cell_index, + const unsigned int face_number) const + { + const unsigned int compressed_cell_index = + compute_compressed_cell_index(cell_index); + if (face_number == numbers::invalid_unsigned_int) { Assert(state == State::cell_vector, ExcMessage( "This mapping info is not reinitialized for a cell vector!")); - if (do_cell_index_compression) - { - Assert(cell_index_to_compressed_cell_index[cell_index] != - numbers::invalid_unsigned_int, - ExcMessage("Mapping info object was not initialized for this" - " active cell index!")); - const auto it_begin = - unit_points.begin() + - unit_points_index - [cell_index_to_compressed_cell_index[cell_index]]; - const auto it_end = - unit_points.begin() + - unit_points_index - [cell_index_to_compressed_cell_index[cell_index] + 1]; - return make_array_view(it_begin, it_end); - } - else - { - const auto it_begin = - unit_points.begin() + unit_points_index[cell_index]; - const auto it_end = - unit_points.begin() + unit_points_index[cell_index + 1]; - return make_array_view(it_begin, it_end); - } + return compressed_cell_index; } - else if (cell_index != numbers::invalid_unsigned_int) + else { + Assert(cell_index != numbers::invalid_unsigned_int, + ExcMessage( + "cell_index has to be set if face_number is specified!")); Assert(state == State::faces_on_cells_in_vector, ExcMessage("This mapping info is not reinitialized for faces" " on cells in a vector!")); - if (do_cell_index_compression) - { - Assert( - cell_index_to_compressed_cell_index[cell_index] != - numbers::invalid_unsigned_int, - ExcMessage( - "Mapping info object was not initialized for this active cell index" - " and corresponding face numbers!")); - const unsigned int current_face_index = - cell_index_offset - [cell_index_to_compressed_cell_index[cell_index]] + - face_number; - const auto it_begin = - unit_points.begin() + unit_points_index[current_face_index]; - const auto it_end = - unit_points.begin() + unit_points_index[current_face_index + 1]; - return make_array_view(it_begin, it_end); - } - else - { - const unsigned int current_face_index = - cell_index_offset[cell_index] + face_number; - const auto it_begin = - unit_points.begin() + unit_points_index[current_face_index]; - const auto it_end = - unit_points.begin() + unit_points_index[current_face_index + 1]; - return make_array_view(it_begin, it_end); - } + return cell_index_offset[compressed_cell_index] + face_number; } - else - AssertThrow( - false, - ExcMessage( - "cell_index has to be specified if face_number is specified!")); } - template - inline const typename MappingInfo::MappingData & - MappingInfo::get_mapping_data( - const unsigned int cell_index, - const unsigned int face_number) const + template + unsigned int + MappingInfo::compute_compressed_cell_index( + const unsigned int cell_index) const { - if (cell_index == numbers::invalid_unsigned_int && - face_number == numbers::invalid_unsigned_int) + if (do_cell_index_compression) { - Assert(state == State::single_cell, - ExcMessage( - "This mapping info is not reinitialized for a single cell!")); - return mapping_data[0]; + Assert(cell_index_to_compressed_cell_index[cell_index] != + numbers::invalid_unsigned_int, + ExcMessage("Mapping info object was not initialized for this" + " active cell index!")); + return cell_index_to_compressed_cell_index[cell_index]; } - else if (face_number == numbers::invalid_unsigned_int) + else + return cell_index; + } + + + template + void + MappingInfo::store_unit_points( + const unsigned int unit_points_index_offset, + const unsigned int n_q_points, + const unsigned int n_q_points_unvectorized, + const std::vector> &points) + { + const unsigned int n_lanes = + dealii::internal::VectorizedArrayTrait::width; + + for (unsigned int q = 0; q < n_q_points; ++q) { - Assert(state == State::cell_vector, - ExcMessage( - "This mapping info is not reinitialized for a cell vector!")); - if (do_cell_index_compression) - { - Assert(cell_index_to_compressed_cell_index[cell_index] != - numbers::invalid_unsigned_int, - ExcMessage("Mapping info object was not initialized for this" - " active cell index!")); - return mapping_data - [cell_index_to_compressed_cell_index[cell_index]]; - } - else - return mapping_data[cell_index]; + const unsigned int offset = unit_points_index_offset + q; + for (unsigned int v = 0; + v < n_lanes && q * n_lanes + v < n_q_points_unvectorized; + ++v) + for (unsigned int d = 0; d < dim; ++d) + dealii::internal::VectorizedArrayTrait::get( + unit_points[offset][d], v) = points[q * n_lanes + v][d]; } - else if (cell_index != numbers::invalid_unsigned_int) + } + + + + template + void + MappingInfo::store_mapping_data( + const unsigned int unit_points_index_offset, + const unsigned int n_q_points, + const unsigned int n_q_points_unvectorized, + const MappingInfo::MappingData &mapping_data) + { + const unsigned int n_lanes = + dealii::internal::VectorizedArrayTrait::width; + + for (unsigned int q = 0; q < n_q_points; ++q) { - Assert(state == State::faces_on_cells_in_vector, - ExcMessage("This mapping info is not reinitialized for faces" - " on cells in a vector!")); - if (do_cell_index_compression) + const unsigned int offset = unit_points_index_offset + q; + for (unsigned int v = 0; + v < n_lanes && q * n_lanes + v < n_q_points_unvectorized; + ++v) { - Assert( - cell_index_to_compressed_cell_index[cell_index] != - numbers::invalid_unsigned_int, - ExcMessage( - "Mapping info object was not initialized for this active cell index" - " and corresponding face numbers!")); - return mapping_data - [cell_index_offset - [cell_index_to_compressed_cell_index[cell_index]] + - face_number]; + if (update_flags_mapping & UpdateFlags::update_jacobians) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int s = 0; s < spacedim; ++s) + dealii::internal::VectorizedArrayTrait::get( + jacobians[offset][d][s], v) = + mapping_data.jacobians[q * n_lanes + v][d][s]; + if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int s = 0; s < spacedim; ++s) + dealii::internal::VectorizedArrayTrait::get( + inverse_jacobians[offset][s][d], v) = + mapping_data.inverse_jacobians[q * n_lanes + v][s][d]; + if (update_flags_mapping & UpdateFlags::update_JxW_values) + dealii::internal::VectorizedArrayTrait::get( + JxW_values[offset], v) = + mapping_data.JxW_values[q * n_lanes + v]; + if (update_flags_mapping & UpdateFlags::update_normal_vectors) + for (unsigned int s = 0; s < spacedim; ++s) + dealii::internal::VectorizedArrayTrait::get( + normal_vectors[offset][s], v) = + mapping_data.normal_vectors[q * n_lanes + v][s]; + if (update_flags_mapping & UpdateFlags::update_quadrature_points) + for (unsigned int s = 0; s < spacedim; ++s) + dealii::internal::VectorizedArrayTrait::get( + real_points[offset][s], v) = + mapping_data.quadrature_points[q * n_lanes + v][s]; } - else - return mapping_data[cell_index_offset[cell_index] + face_number]; } - else - AssertThrow( - false, - ExcMessage( - "cell_index has to be specified if face_number is specified!")); } - template - const Mapping & - MappingInfo::get_mapping() const + template + void + MappingInfo::resize_unit_points( + const unsigned int n_unit_point_batches) { - return *mapping; + unit_points.resize(n_unit_point_batches); } - template - UpdateFlags - MappingInfo::get_update_flags() const + template + void + MappingInfo::resize_data_fields( + const unsigned int n_data_point_batches) { - return update_flags; + if (update_flags_mapping & UpdateFlags::update_jacobians) + jacobians.resize(n_data_point_batches); + if (update_flags_mapping & UpdateFlags::update_inverse_jacobians) + inverse_jacobians.resize(n_data_point_batches); + if (update_flags_mapping & UpdateFlags::update_JxW_values) + JxW_values.resize(n_data_point_batches); + if (update_flags_mapping & UpdateFlags::update_normal_vectors) + normal_vectors.resize(n_data_point_batches); + if (update_flags_mapping & UpdateFlags::update_quadrature_points) + real_points.resize(n_data_point_batches); } - template - boost::signals2::connection - MappingInfo::connect_is_reinitialized( - const std::function &set_is_reinitialized) + template + inline const Point< + dim, + typename MappingInfo::VectorizedArrayType> * + MappingInfo::get_unit_point( + const unsigned int offset) const { - return is_reinitialized.connect(set_is_reinitialized); + return &unit_points[offset]; } - template - void - MappingInfo::compute_mapping_data_for_quadrature( - const typename Triangulation::cell_iterator &cell, - CellSimilarity::Similarity & cell_similarity, - const Quadrature & quadrature, - MappingData & mapping_data) + template + inline const Point * + MappingInfo::get_real_point( + const unsigned int offset) const { - update_flags_mapping |= - mapping->requires_update_flags(update_flags_mapping); + return &real_points[offset]; + } - mapping_data.initialize(quadrature.size(), update_flags_mapping); - // reuse internal_mapping_data for MappingQ to avoid memory allocations - if (const MappingQ *mapping_q = - dynamic_cast *>(&(*mapping))) + + template + unsigned int + MappingInfo::compute_unit_point_index_offset( + const unsigned int cell_index, + const unsigned int face_number) const + { + if (cell_index == numbers::invalid_unsigned_int && + face_number == numbers::invalid_unsigned_int) { - (void)mapping_q; - auto &data = - dynamic_cast::InternalData &>( - *internal_mapping_data); - data.initialize(update_flags_mapping, quadrature, quadrature.size()); + Assert(state == State::single_cell, + ExcMessage( + "This mapping info is not reinitialized for a single cell!")); + return 0; } else { - internal_mapping_data = - mapping->get_data(update_flags_mapping, quadrature); + const unsigned int offset = + compute_geometry_index_offset(cell_index, face_number); + return unit_points_index[offset]; } - - cell_similarity = mapping->fill_fe_values( - cell, cell_similarity, quadrature, *internal_mapping_data, mapping_data); } - template - void - MappingInfo:: - compute_mapping_data_for_immersed_surface_quadrature( - const typename Triangulation::cell_iterator &cell, - const ImmersedSurfaceQuadrature & quadrature, - MappingData & mapping_data) + template + unsigned int + MappingInfo::compute_data_index_offset( + const unsigned int cell_index, + const unsigned int face_number) const { - update_flags_mapping |= - mapping->requires_update_flags(update_flags_mapping); - - mapping_data.initialize(quadrature.size(), update_flags_mapping); - - // reuse internal_mapping_data for MappingQ to avoid memory allocations - if (const MappingQ *mapping_q = - dynamic_cast *>(&(*mapping))) + if (cell_index == numbers::invalid_unsigned_int && + face_number == numbers::invalid_unsigned_int) { - (void)mapping_q; - auto &data = - dynamic_cast::InternalData &>( - *internal_mapping_data); - data.initialize(update_flags_mapping, quadrature, quadrature.size()); + Assert(state == State::single_cell, + ExcMessage( + "This mapping info is not reinitialized for a single cell!")); + return 0; } else { - internal_mapping_data = - mapping->get_data(update_flags_mapping, quadrature); + const unsigned int offset = + compute_geometry_index_offset(cell_index, face_number); + return data_index_offsets[offset]; } + } + - mapping->fill_fe_immersed_surface_values(cell, - quadrature, - *internal_mapping_data, - mapping_data); + + template + inline const DerivativeForm<1, dim, spacedim, Number> * + MappingInfo::get_jacobian( + const unsigned int offset) const + { + return &jacobians[offset]; } - template - void - MappingInfo::compute_mapping_data_for_face_quadrature( - const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const Quadrature & quadrature, - MappingData & mapping_data) + template + inline const DerivativeForm<1, spacedim, dim, Number> * + MappingInfo::get_inverse_jacobian( + const unsigned int offset) const { - update_flags_mapping |= - mapping->requires_update_flags(update_flags_mapping); + return &inverse_jacobians[offset]; + } - mapping_data.initialize(quadrature.size(), update_flags_mapping); - // reuse internal_mapping_data for MappingQ to avoid memory allocations - if (const MappingQ *mapping_q = - dynamic_cast *>(&(*mapping))) - { - auto &data = - dynamic_cast::InternalData &>( - *internal_mapping_data); - data.initialize_face(update_flags_mapping, - QProjector::project_to_oriented_face( - ReferenceCells::get_hypercube(), - quadrature, - face_no, - cell->face_orientation(face_no), - cell->face_flip(face_no), - cell->face_rotation(face_no)), - quadrature.size()); - - mapping_q->fill_mapping_data_for_face_quadrature( - cell, face_no, quadrature, *internal_mapping_data, mapping_data); - } - else - { - auto internal_mapping_data = - mapping->get_face_data(update_flags_mapping, - hp::QCollection(quadrature)); - - mapping->fill_fe_face_values(cell, - face_no, - hp::QCollection(quadrature), - *internal_mapping_data, - mapping_data); - } + template + inline const Tensor<1, spacedim, Number> * + MappingInfo::get_normal_vector( + const unsigned int offset) const + { + return &normal_vectors[offset]; + } + + + + template + inline const Number * + MappingInfo::get_JxW(const unsigned int offset) const + { + return &JxW_values[offset]; + } + + + + template + const Mapping & + MappingInfo::get_mapping() const + { + return *mapping; + } + + + + template + UpdateFlags + MappingInfo::get_update_flags() const + { + return update_flags; + } + + + + template + UpdateFlags + MappingInfo::get_update_flags_mapping() const + { + return update_flags_mapping; + } + + + + template + boost::signals2::connection + MappingInfo::connect_is_reinitialized( + const std::function &set_is_reinitialized) + { + return is_reinitialized.connect(set_is_reinitialized); } } // namespace NonMatching diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index b1e8eb1676..bd3d1dcf50 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -186,8 +186,10 @@ MappingQ::InternalData::initialize_face( if (this->update_each & (update_boundary_forms | update_normal_vectors | update_JxW_values)) { - aux.resize(dim - 1, - AlignedVector>(n_original_q_points)); + aux.resize(dim - 1); + aux[0].resize(n_original_q_points); + if (dim > 2) + aux[1].resize(n_original_q_points); // Compute tangentials to the unit cell. for (const unsigned int i : GeometryInfo::face_indices()) diff --git a/tests/non_matching/mapping_info.cc b/tests/non_matching/mapping_info.cc index 612ed9185e..a2b5b98dce 100644 --- a/tests/non_matching/mapping_info.cc +++ b/tests/non_matching/mapping_info.cc @@ -40,10 +40,10 @@ using namespace dealii; +template void test(const bool filtered_compression) { - constexpr unsigned int dim = 2; constexpr unsigned int degree = 1; FE_Q fe_q(degree); @@ -471,7 +471,11 @@ main(int argc, char **argv) initlog(); - test(true); + test<2>(true); deallog << std::endl; - test(false); + test<2>(false); + deallog << std::endl; + test<3>(true); + deallog << std::endl; + test<3>(false); } diff --git a/tests/non_matching/mapping_info.output b/tests/non_matching/mapping_info.output index b4737523f9..5da28926fd 100644 --- a/tests/non_matching/mapping_info.output +++ b/tests/non_matching/mapping_info.output @@ -6,3 +6,11 @@ DEAL:: DEAL::check difference l2 norm cell: 0.00000 DEAL::check difference l2 norm surface: 0.00000 DEAL::check difference l2 norm faces: 0.00000 +DEAL:: +DEAL::check difference l2 norm cell: 0.00000 +DEAL::check difference l2 norm surface: 0.00000 +DEAL::check difference l2 norm faces: 0.00000 +DEAL:: +DEAL::check difference l2 norm cell: 0.00000 +DEAL::check difference l2 norm surface: 0.00000 +DEAL::check difference l2 norm faces: 0.00000 -- 2.39.5