From: Maximilian Bergbauer Date: Mon, 28 Aug 2023 11:20:46 +0000 (+0200) Subject: Optimize FEPointEvaluation/NM::MappingInfo X-Git-Tag: relicensing~489^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15971%2Fhead;p=dealii.git Optimize FEPointEvaluation/NM::MappingInfo --- diff --git a/include/deal.II/base/derivative_form.h b/include/deal.II/base/derivative_form.h index e3a73fa695..8c1e4710f6 100644 --- a/include/deal.II/base/derivative_form.h +++ b/include/deal.II/base/derivative_form.h @@ -591,6 +591,82 @@ transpose(const DerivativeForm<1, dim, spacedim, Number> &DF) } + +/** + * Specialization of apply_transformation() for a diagonal DerivativeForm. + * + * @relatesalso DerivativeForm + */ +template +inline Tensor<1, spacedim, typename ProductType::type> +apply_diagonal_transformation( + const DerivativeForm<1, dim, spacedim, Number1> &grad_F, + const Tensor<1, dim, Number2> &d_x) +{ + Assert(dim == spacedim, + ExcMessage("Only dim = spacedim allowed for diagonal transformation")); + Tensor<1, spacedim, typename ProductType::type> dest; + for (unsigned int i = 0; i < spacedim; ++i) + dest[i] = grad_F[i][i] * d_x[i]; + return dest; +} + + +/** + * Similar to the previous apply_diagonal_transformation(), specialized for the + * case `dim == spacedim` where we can return a rank-2 tensor instead of the + * more general `DerivativeForm`. + * Each row of the result corresponds to one of the rows of @p D_X transformed + * by @p grad_F, equivalent to $\mathrm{D\_X} \, \mathrm{grad\_F}^T$ in matrix + * notation. + * + * @relatesalso DerivativeForm + */ +template +inline Tensor<2, dim, typename ProductType::type> +apply_diagonal_transformation( + const DerivativeForm<1, dim, dim, Number1> &grad_F, + const Tensor<2, dim, Number2> &D_X) +{ + Tensor<2, dim, typename ProductType::type> dest; + for (unsigned int i = 0; i < dim; ++i) + dest[i] = apply_diagonal_transformation(grad_F, D_X[i]); + + return dest; +} + + + +/** + * Similar to the previous apply_diagonal_transformation(). + * Each row of the result corresponds to one of the rows of @p D_X transformed + * by @p grad_F. + * + * @relatesalso DerivativeForm + */ +template +inline Tensor<1, + n_components, + Tensor<1, spacedim, typename ProductType::type>> +apply_diagonal_transformation( + const DerivativeForm<1, dim, spacedim, Number1> &grad_F, + const Tensor<1, n_components, Tensor<1, dim, Number2>> &D_X) +{ + Tensor<1, + n_components, + Tensor<1, spacedim, typename ProductType::type>> + dest; + for (unsigned int i = 0; i < n_components; ++i) + dest[i] = apply_diagonal_transformation(grad_F, D_X[i]); + + return dest; +} + + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/include/deal.II/base/floating_point_comparator.h b/include/deal.II/base/floating_point_comparator.h index 9c3e1214d9..fc76666ab3 100644 --- a/include/deal.II/base/floating_point_comparator.h +++ b/include/deal.II/base/floating_point_comparator.h @@ -18,6 +18,7 @@ #include +#include #include #include #include @@ -105,6 +106,14 @@ struct FloatingPointComparator ComparisonResult compare(const Tensor &t1, const Tensor &t2) const; + /** + * Compare two derivative forms. + */ + template + ComparisonResult + compare(const DerivativeForm &t1, + const DerivativeForm &t2) const; + /** * Compare two tables. */ @@ -217,6 +226,24 @@ FloatingPointComparator::compare(const Tensor &t1, +template +template +typename FloatingPointComparator::ComparisonResult +FloatingPointComparator::compare( + const DerivativeForm &t1, + const DerivativeForm &t2) const +{ + for (unsigned int i = 0; i < dim; ++i) + { + const ComparisonResult result = compare(t1[i], t2[i]); + if (result != ComparisonResult::equal) + return result; + } + return ComparisonResult::equal; +} + + + template template typename FloatingPointComparator::ComparisonResult diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index 9c7a568542..9189fd0b2e 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -30,6 +30,7 @@ #include #include +#include #include #include @@ -85,13 +86,12 @@ namespace internal result[component] = vector_entry; } - static void - write_value(VectorizedArrayType &vector_entry, - const unsigned int component, - const vectorized_value_type &result) + static ScalarNumber + sum_value(const unsigned int component, + const vectorized_value_type &result) { AssertIndexRange(component, n_components); - vector_entry = result[component]; + return result[component].sum(); } static void @@ -237,12 +237,10 @@ namespace internal result = vector_entry; } - static void - write_value(VectorizedArrayType &vector_entry, - const unsigned int, - const vectorized_value_type &result) + static ScalarNumber + sum_value(const unsigned int, const vectorized_value_type &result) { - vector_entry = result; + return result.sum(); } static void @@ -389,15 +387,16 @@ namespace internal const unsigned int component, scalar_value_type &result) { + AssertIndexRange(component, dim); result[component] = vector_entry; } - static void - write_value(VectorizedArrayType &vector_entry, - const unsigned int component, - const vectorized_value_type &result) + static ScalarNumber + sum_value(const unsigned int component, + const vectorized_value_type &result) { - vector_entry = result[component]; + AssertIndexRange(component, dim); + return result[component].sum(); } static void @@ -543,12 +542,10 @@ namespace internal result = vector_entry; } - static void - write_value(VectorizedArrayType &vector_entry, - const unsigned int, - const vectorized_value_type &result) + static ScalarNumber + sum_value(const unsigned int, const vectorized_value_type &result) { - vector_entry = result; + return result.sum(); } static void @@ -834,6 +831,23 @@ public: void reinit(const unsigned int cell_index, const unsigned int face_number); + /** + * This function interpolates the finite element solution, represented by + * `solution_values`, on the cell and `unit_points` passed to reinit(). + * + * @param[in] solution_values This array is supposed to contain the unknown + * values on the element read out by + * `FEEvaluation::read_dof_values(global_vector)`. + * + * @param[in] evaluation_flags Flags specifying which quantities should be + * evaluated at the points. + */ + template + void + evaluate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); + /** * This function interpolates the finite element solution, represented by * `solution_values`, on the cell and `unit_points` passed to reinit(). @@ -847,7 +861,36 @@ public: */ void evaluate(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); + const EvaluationFlags::EvaluationFlags &evaluation_flags) + { + evaluate(StridedArrayView(solution_values.data(), + solution_values.size()), + evaluation_flags); + } + + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points multiplied be + * the Jacobian determinant times the quadrature weight (JxW). + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during + * `FEEvaluation::set_dof_values(global_vector)` or + * `FEEvaluation::distribute_local_to_global(global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + */ + template + void + integrate(const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags); + /** * This function multiplies the quantities passed in by previous @@ -869,7 +912,12 @@ public: */ void integrate(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags); + const EvaluationFlags::EvaluationFlags &integration_flags) + { + integrate(StridedArrayView(solution_values.data(), + solution_values.size()), + integration_flags); + } /** * This function multiplies the quantities passed in by previous @@ -882,7 +930,35 @@ public: * finite element formulation. * * @param[out] solution_values This array will contain the result of the - * integral, which can be used to during + * integral, which can be used during + * `FEEvaluation::set_dof_values(global_vector)` or + * `FEEvaluation::distribute_local_to_global(global_vector)`. Note + * that for multi-component systems where only some of the components are + * selected by the present class, the entries in `solution_values` not touched + * by this class will be set to zero. + * + * @param[in] integration_flags Flags specifying which quantities should be + * integrated at the points. + * + */ + template + void + test_and_sum( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags); + + /** + * This function multiplies the quantities passed in by previous + * submit_value() or submit_gradient() calls by the value or gradient of the + * test functions, and performs summation over all given points. This is + * similar to the integration of a bilinear form in terms of the test + * function, with the difference that this formula does not include a `JxW` + * factor (in contrast to the integrate function of this class). This allows + * the class to naturally embed point information (e.g. particles) into a + * finite element formulation. + * + * @param[out] solution_values This array will contain the result of the + * integral, which can be used during * `cell->set_dof_values(solution_values, global_vector)` or * `cell->distribute_local_to_global(solution_values, global_vector)`. Note * that for multi-component systems where only some of the components are @@ -895,7 +971,12 @@ public: */ void test_and_sum(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags); + const EvaluationFlags::EvaluationFlags &integration_flags) + { + test_and_sum(StridedArrayView(solution_values.data(), + solution_values.size()), + integration_flags); + } /** * Return the value at quadrature point number @p point_index after a call to @@ -930,18 +1011,6 @@ public: const gradient_type & get_gradient(const unsigned int point_index) const; - /** - * Return the gradient in unit coordinates at the point with index - * `point_index` after a call to FEPointEvaluation::evaluate() with - * EvaluationFlags::gradients set, or the gradient that has been stored there - * with a call to FEPointEvaluation::submit_gradient(). If the object is - * vector-valued, a vector-valued return argument is given. Note that when - * vectorization is enabled, values from several points are grouped - * together. - */ - const gradient_type & - get_unit_gradient(const unsigned int point_index) const; - /** * Write a contribution that is tested by the gradient to the field * containing the values on points with the given `point_index`. Access to @@ -1032,6 +1101,7 @@ private: * precomputes the @p shapes vector, holding the evaluation of 1D basis * functions of tensor product polynomials, if necessary. */ + template void do_reinit(); @@ -1039,87 +1109,109 @@ private: * Resizes necessary data fields, reads in and renumbers solution values. * Interpolates onto face if face path is selected. */ + template void prepare_evaluate_fast( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); /** * Evaluates the actual interpolation on the cell or face for a quadrature * batch. */ + template void compute_evaluate_fast( - const EvaluationFlags::EvaluationFlags &evaluation_flags, - const unsigned int n_shapes, - const unsigned int qb, - vectorized_value_type &value, - interface_vectorized_gradient_type &gradient); + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags, + const unsigned int n_shapes, + const unsigned int qb, + vectorized_value_type &value, + interface_vectorized_gradient_type &gradient); /** * Fast path of the evaluate function. */ + template void - evaluate_fast(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); + evaluate_fast( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); /** * Slow path of the evaluate function using FEValues. */ + template void - evaluate_slow(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags); + evaluate_slow( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags); /** * Integrates the product of the data passed in by submit_value() and * submit_gradient() with the values or gradients of test functions on the * cell or face for a given quadrature batch. */ + template void compute_integrate_fast( - const EvaluationFlags::EvaluationFlags &integration_flags, - const unsigned int n_shapes, - const unsigned int qb, - const vectorized_value_type &value, - const interface_vectorized_gradient_type &gradient); + const EvaluationFlags::EvaluationFlags &integration_flags, + const unsigned int n_shapes, + const unsigned int qb, + const vectorized_value_type value, + const interface_vectorized_gradient_type gradient, + vectorized_value_type *solution_values_vectorized_linear); /** * Addition across the lanes of VectorizedArray as accumulated by the * compute_integrate_fast_function(), writing the sum into the result vector. * Applies face contributions to cell contributions for face path. */ + template void finish_integrate_fast( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags); + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + vectorized_value_type *solution_values_vectorized_linear); /** * Fast path of the integrate function. */ - template + template void - integrate_fast(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags); + integrate_fast( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags); /** * Slow path of the integrate function using FEValues. */ - template + template void - integrate_slow(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags); + integrate_slow( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags); /** * Implementation of the integrate/test_and_sum function. */ - template + template void - do_integrate(const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags); + do_integrate( + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags); /** * Number of quadrature batches of the current cell/face. */ + const unsigned int n_q_batches; + + /** + * Number of quadrature points/batches of the current cell/face. + */ const unsigned int n_q_points; /** @@ -1144,9 +1236,9 @@ private: std::vector> poly; /** - * Store whether the polynomials are linear with nodes at 0 and 1. + * Store whether the linear path should be used. */ - bool polynomials_are_hat_functions; + bool use_linear_path; /** * Renumbering between the unknowns of unknowns implied by the FiniteElement @@ -1181,11 +1273,6 @@ private: */ std::vector values; - /** - * Temporary array to store the gradients in unit coordinates at the points. - */ - std::vector unit_gradients; - /** * Temporary array to store the gradients in real coordinates at the points. */ @@ -1234,6 +1321,11 @@ private: */ const Number *JxW_ptr; + /** + * Cell type describing the geometry of the cell and compression of jacobians. + */ + internal::MatrixFreeFunctions::GeometryType cell_type; + /** * Number of unknowns per component, i.e., number of unique basis functions, * for the chosen FiniteElement (or base element). @@ -1340,7 +1432,8 @@ FEPointEvaluation::FEPointEvaluation( const FiniteElement &fe, const UpdateFlags update_flags, const unsigned int first_selected_component) - : n_q_points(numbers::invalid_unsigned_int) + : n_q_batches(numbers::invalid_unsigned_int) + , n_q_points(numbers::invalid_unsigned_int) , n_q_points_scalar(numbers::invalid_unsigned_int) , mapping(&mapping) , fe(&fe) @@ -1366,7 +1459,8 @@ FEPointEvaluation::FEPointEvaluation( NonMatching::MappingInfo &mapping_info, const FiniteElement &fe, const unsigned int first_selected_component) - : n_q_points(numbers::invalid_unsigned_int) + : n_q_batches(numbers::invalid_unsigned_int) + , n_q_points(numbers::invalid_unsigned_int) , n_q_points_scalar(numbers::invalid_unsigned_int) , mapping(&mapping_info.get_mapping()) , fe(&fe) @@ -1388,17 +1482,17 @@ FEPointEvaluation::FEPointEvaluation( template FEPointEvaluation::FEPointEvaluation( FEPointEvaluation &other) noexcept - : n_q_points(other.n_q_points) + : n_q_batches(other.n_q_batches) + , n_q_points(other.n_q_points) , n_q_points_scalar(other.n_q_points_scalar) , mapping(other.mapping) , fe(other.fe) , poly(other.poly) - , polynomials_are_hat_functions(other.polynomials_are_hat_functions) + , use_linear_path(other.use_linear_path) , renumber(other.renumber) , solution_renumbered(other.solution_renumbered) , solution_renumbered_vectorized(other.solution_renumbered_vectorized) , values(other.values) - , unit_gradients(other.unit_gradients) , gradients(other.gradients) , dofs_per_component(other.dofs_per_component) , dofs_per_component_face(other.dofs_per_component_face) @@ -1430,17 +1524,17 @@ FEPointEvaluation::FEPointEvaluation( template FEPointEvaluation::FEPointEvaluation( FEPointEvaluation &&other) noexcept - : n_q_points(other.n_q_points) + : n_q_batches(other.n_q_batches) + , n_q_points(other.n_q_points) , n_q_points_scalar(other.n_q_points_scalar) , mapping(other.mapping) , fe(other.fe) , poly(other.poly) - , polynomials_are_hat_functions(other.polynomials_are_hat_functions) + , use_linear_path(other.use_linear_path) , renumber(other.renumber) , solution_renumbered(other.solution_renumbered) , solution_renumbered_vectorized(other.solution_renumbered_vectorized) , values(other.values) - , unit_gradients(other.unit_gradients) , gradients(other.gradients) , dofs_per_component(other.dofs_per_component) , dofs_per_component_face(other.dofs_per_component_face) @@ -1519,15 +1613,18 @@ FEPointEvaluation::setup( if (is_lexicographic) renumber.clear(); - polynomials_are_hat_functions = - (poly.size() == 2 && poly[0].value(0.) == 1. && - poly[0].value(1.) == 0. && poly[1].value(0.) == 0. && - poly[1].value(1.) == 1.); + use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. && + poly[0].value(1.) == 0. && poly[1].value(0.) == 0. && + poly[1].value(1.) == 1.) && + (fe->n_components() == n_components); const unsigned int size_face = 2 * dofs_per_component_face; const unsigned int size_cell = dofs_per_component; scratch_data_scalar.resize(size_face + size_cell); + solution_renumbered.resize(dofs_per_component); + solution_renumbered_vectorized.resize(dofs_per_component); + fast_path = true; } else @@ -1556,7 +1653,7 @@ FEPointEvaluation::setup( template -void +inline void FEPointEvaluation::reinit( const typename Triangulation::cell_iterator &cell, const ArrayView> &unit_points) @@ -1577,48 +1674,50 @@ FEPointEvaluation::reinit( fe_values->reinit(cell); } - do_reinit(); + if (use_linear_path) + do_reinit(); + else + do_reinit(); } template -void +inline void FEPointEvaluation::reinit() { current_cell_index = numbers::invalid_unsigned_int; current_face_number = numbers::invalid_unsigned_int; - do_reinit(); + if (use_linear_path) + do_reinit(); + else + do_reinit(); } template -void +inline void FEPointEvaluation::reinit( const unsigned int cell_index) { current_cell_index = cell_index; current_face_number = numbers::invalid_unsigned_int; - do_reinit(); + if (use_linear_path) + do_reinit(); + else + do_reinit(); if (!fast_path) { - const auto unit_points_vectorized = mapping_info->get_unit_point( - mapping_info->compute_unit_point_index_offset(current_cell_index, - current_face_number)); - const unsigned int n_q_points_unvectorized = - mapping_info->get_n_q_points_unvectorized(current_cell_index, - current_face_number); - - std::vector> unit_points(n_q_points_unvectorized); + std::vector> unit_points(n_q_points_scalar); - for (unsigned int v = 0; v < n_q_points_unvectorized; ++v) + for (unsigned int v = 0; v < n_q_points_scalar; ++v) for (unsigned int d = 0; d < dim; ++d) - unit_points[v][d] = unit_points_vectorized[v / n_lanes_internal][d] - [v % n_lanes_internal]; + unit_points[v][d] = + unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal]; fe_values = std::make_shared>( *mapping, @@ -1634,7 +1733,7 @@ FEPointEvaluation::reinit( template -void +inline void FEPointEvaluation::reinit( const unsigned int cell_index, const unsigned int face_number) @@ -1642,27 +1741,44 @@ FEPointEvaluation::reinit( current_cell_index = cell_index; current_face_number = face_number; - do_reinit(); + if (use_linear_path) + do_reinit(); + else + do_reinit(); } template +template inline void FEPointEvaluation::do_reinit() { + const unsigned int geometry_index = + mapping_info->template compute_geometry_index_offset( + current_cell_index, current_face_number); + + cell_type = mapping_info->get_cell_type(geometry_index); + const_cast(n_q_points_scalar) = - mapping_info->get_n_q_points_unvectorized(current_cell_index, - current_face_number); + mapping_info->get_n_q_points_unvectorized(geometry_index); + + // round up n_q_points_scalar / n_lanes_internal + const_cast(n_q_batches) = + (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal; + + const unsigned int n_q_points_before = n_q_points; - // round up n_points_scalar / n_lanes_user_interface const_cast(n_q_points) = - (n_q_points_scalar + n_lanes_user_interface - 1) / n_lanes_user_interface; + (stride == 1) ? n_q_batches : n_q_points_scalar; - 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 (n_q_points != n_q_points_before) + { + if (update_flags & update_values) + values.resize(n_q_points); + if (update_flags & update_gradients) + gradients.resize(n_q_points); + } if (n_q_points == 0) { @@ -1672,12 +1788,12 @@ FEPointEvaluation::do_reinit() // use face path if mapping_info in face state and number of quadrature points // is large enough - use_face_path = mapping_info->is_face_state() && n_q_points_scalar >= 6; + use_face_path = + is_face && (mapping_info->is_face_state() && n_q_points_scalar >= 6); // set unit point pointer const unsigned int unit_point_offset = - mapping_info->compute_unit_point_index_offset(current_cell_index, - current_face_number); + mapping_info->compute_unit_point_index_offset(geometry_index); if (use_face_path) unit_point_faces_ptr = @@ -1686,35 +1802,43 @@ FEPointEvaluation::do_reinit() unit_point_ptr = mapping_info->get_unit_point(unit_point_offset); // set data pointers + const unsigned int data_offset = + mapping_info->compute_data_index_offset(geometry_index); + const unsigned int compressed_data_offset = + mapping_info->compute_compressed_data_index_offset(geometry_index); +#ifdef DEBUG 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); + jacobian_ptr = mapping_info->get_jacobian(compressed_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) + inverse_jacobian_ptr = + mapping_info->get_inverse_jacobian(compressed_data_offset); + if (is_face && 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); +#else + real_point_ptr = mapping_info->get_real_point(data_offset); + jacobian_ptr = mapping_info->get_jacobian(compressed_data_offset); + inverse_jacobian_ptr = + mapping_info->get_inverse_jacobian(compressed_data_offset); + normal_ptr = mapping_info->get_normal_vector(data_offset); + JxW_ptr = mapping_info->get_JxW(data_offset); +#endif - if (fast_path && !polynomials_are_hat_functions) + if (!is_linear && fast_path) { - // round up n_q_points_scalar / n_lanes_internal - const std::size_t n_batches = - (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal; const std::size_t n_shapes = poly.size(); - for (unsigned int qb = 0; qb < n_batches; ++qb) + for (unsigned int qb = 0; qb < n_q_batches; ++qb) if (use_face_path) { if (dim > 1) { - shapes_faces.resize_fast(n_batches * n_shapes); + shapes_faces.resize_fast(n_q_batches * n_shapes); internal::compute_values_of_array( shapes_faces.data() + qb * n_shapes, poly, @@ -1724,7 +1848,7 @@ FEPointEvaluation::do_reinit() } else { - shapes.resize_fast(n_batches * n_shapes); + shapes.resize_fast(n_q_batches * n_shapes); internal::compute_values_of_array( shapes.data() + qb * n_shapes, poly, @@ -1739,45 +1863,38 @@ FEPointEvaluation::do_reinit() template +template inline void FEPointEvaluation::prepare_evaluate_fast( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags) + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) { - if (use_face_path) - { - if (solution_renumbered.size() != 2 * dofs_per_component_face) - solution_renumbered.resize(2 * dofs_per_component_face); - } - else - { - if (solution_renumbered.size() != dofs_per_component) - solution_renumbered.resize(dofs_per_component); - } + const unsigned int dofs_per_comp = + is_linear ? Utilities::pow(2, dim) : dofs_per_component; + for (unsigned int comp = 0; comp < n_components; ++comp) { const std::size_t offset = - (component_in_base_element + comp) * dofs_per_component; + (component_in_base_element + comp) * dofs_per_comp; - if (use_face_path) + if (is_face_path) { const ScalarNumber *input; - if (renumber.empty()) + if (is_linear || renumber.empty()) { - for (unsigned int i = 0; i < dofs_per_component; ++i) + for (unsigned int i = 0; i < dofs_per_comp; ++i) scratch_data_scalar[i] = solution_values[i + offset]; input = scratch_data_scalar.data(); } else { const unsigned int *renumber_ptr = renumber.data() + offset; - for (unsigned int i = 0; i < dofs_per_component; ++i) + for (unsigned int i = 0; i < dofs_per_comp; ++i) scratch_data_scalar[i] = solution_values[renumber_ptr[i]]; input = scratch_data_scalar.data(); } - ScalarNumber *output = - scratch_data_scalar.begin() + dofs_per_component; + ScalarNumber *output = scratch_data_scalar.begin() + dofs_per_comp; internal::FEFaceNormalEvaluationImpl:: template interpolate(1, @@ -1787,14 +1904,16 @@ FEPointEvaluation::prepare_evaluate_fast( output, current_face_number); - for (unsigned int i = 0; i < 2 * dofs_per_component_face; ++i) + const unsigned int dofs_per_comp_face = + is_linear ? Utilities::pow(2, dim - 1) : dofs_per_component_face; + for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i) ETT::read_value(output[i], comp, solution_renumbered[i]); } else { - if (renumber.empty()) + if ((is_linear && n_components == 1) || renumber.empty()) { - for (unsigned int i = 0; i < dofs_per_component; ++i) + for (unsigned int i = 0; i < dofs_per_comp; ++i) ETT::read_value(solution_values[i + offset], comp, solution_renumbered[i]); @@ -1802,43 +1921,39 @@ FEPointEvaluation::prepare_evaluate_fast( else { const unsigned int *renumber_ptr = renumber.data() + offset; - for (unsigned int i = 0; i < dofs_per_component; ++i) + for (unsigned int i = 0; i < dofs_per_comp; ++i) ETT::read_value(solution_values[renumber_ptr[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()); } template +template inline void FEPointEvaluation::compute_evaluate_fast( - const EvaluationFlags::EvaluationFlags &evaluation_flags, - const unsigned int n_shapes, - const unsigned int qb, - vectorized_value_type &value, - interface_vectorized_gradient_type &gradient) + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags, + const unsigned int n_shapes, + const unsigned int qb, + vectorized_value_type &value, + interface_vectorized_gradient_type &gradient) { - if (use_face_path) + if (is_face_path) { if (evaluation_flags & EvaluationFlags::gradients) { const std::array interpolated_value = - polynomials_are_hat_functions ? + is_linear ? internal::evaluate_tensor_product_value_and_gradient_linear< dim - 1, scalar_value_type, VectorizedArrayType, - 2>(n_shapes, - solution_renumbered.data(), - unit_point_faces_ptr[qb]) : + 2>(solution_renumbered.data(), unit_point_faces_ptr[qb]) : internal::evaluate_tensor_product_value_and_gradient_shapes< dim - 1, scalar_value_type, @@ -1889,30 +2004,42 @@ FEPointEvaluation::compute_evaluate_fast( } else { - value = polynomials_are_hat_functions ? - internal::evaluate_tensor_product_value_linear< - dim - 1, - scalar_value_type, - VectorizedArrayType>(n_shapes, - solution_renumbered.data(), - unit_point_faces_ptr[qb]) : - internal::evaluate_tensor_product_value_shapes< - dim - 1, - scalar_value_type, - VectorizedArrayType, - false>(shapes_faces.data() + qb * n_shapes, - n_shapes, - solution_renumbered.data()); + value = is_linear ? internal::evaluate_tensor_product_value_linear< + dim - 1, + scalar_value_type, + VectorizedArrayType>(solution_renumbered.data(), + unit_point_faces_ptr[qb]) : + internal::evaluate_tensor_product_value_shapes< + dim - 1, + scalar_value_type, + VectorizedArrayType, + false>(shapes_faces.data() + qb * n_shapes, + n_shapes, + solution_renumbered.data()); } } else { if (evaluation_flags & EvaluationFlags::gradients) { - const std::array result = - polynomials_are_hat_functions ? - internal::evaluate_tensor_product_value_and_gradient_linear( - n_shapes, solution_renumbered.data(), unit_point_ptr[qb]) : + std::array result; + if constexpr (is_linear) + { + if constexpr (n_components == 1) + result = + internal::evaluate_tensor_product_value_and_gradient_linear< + dim, + scalar_value_type, + VectorizedArrayType, + 1, + stride_view>(solution_values.data(), unit_point_ptr[qb]); + else + result = + internal::evaluate_tensor_product_value_and_gradient_linear( + solution_renumbered.data(), unit_point_ptr[qb]); + } + else + result = internal::evaluate_tensor_product_value_and_gradient_shapes< dim, scalar_value_type, @@ -1930,17 +2057,26 @@ FEPointEvaluation::compute_evaluate_fast( } else { - value = - polynomials_are_hat_functions ? - internal::evaluate_tensor_product_value_linear( - n_shapes, solution_renumbered.data(), unit_point_ptr[qb]) : - internal::evaluate_tensor_product_value_shapes< - dim, - scalar_value_type, - VectorizedArrayType, - false>(shapes.data() + qb * n_shapes, - n_shapes, - solution_renumbered.data()); + if constexpr (is_linear) + { + if constexpr (n_components == 1) + value = internal::evaluate_tensor_product_value_linear< + dim, + scalar_value_type, + VectorizedArrayType, + stride_view>(solution_values.data(), unit_point_ptr[qb]); + else + value = internal::evaluate_tensor_product_value_linear( + solution_renumbered.data(), unit_point_ptr[qb]); + } + else + value = internal::evaluate_tensor_product_value_shapes< + dim, + scalar_value_type, + VectorizedArrayType, + false>(shapes.data() + qb * n_shapes, + n_shapes, + solution_renumbered.data()); } } } @@ -1948,28 +2084,33 @@ FEPointEvaluation::compute_evaluate_fast( template +template inline void FEPointEvaluation::evaluate_fast( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags) + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) { - prepare_evaluate_fast(solution_values, evaluation_flags); - - // loop over quadrature batches qb / points q - const unsigned int n_shapes = poly.size(); - vectorized_value_type value; - interface_vectorized_gradient_type gradient; - for (unsigned int qb = 0, q = 0; q < n_q_points_scalar; - ++qb, q += n_lanes_internal) + if (!(is_linear && n_components == 1) || is_face_path) + prepare_evaluate_fast(solution_values, + evaluation_flags); + + // loop over quadrature batches qb + const unsigned int n_shapes = is_linear ? 2 : poly.size(); + + for (unsigned int qb = 0; qb < n_q_batches; ++qb) { - compute_evaluate_fast(evaluation_flags, n_shapes, qb, value, gradient); + vectorized_value_type value; + interface_vectorized_gradient_type gradient; + + compute_evaluate_fast( + solution_values, evaluation_flags, n_shapes, qb, value, gradient); if (evaluation_flags & EvaluationFlags::values) { - for (unsigned int v = 0; - v < stride && (stride == 1 || q + v < n_q_points_scalar); - ++v) - ETT::set_value(value, v, values[qb * stride + v]); + for (unsigned int v = 0, offset = qb * stride; + v < stride && (stride == 1 || offset < n_q_points_scalar); + ++v, ++offset) + ETT::set_value(value, v, values[offset]); } if (evaluation_flags & EvaluationFlags::gradients) { @@ -1977,15 +2118,25 @@ FEPointEvaluation::evaluate_fast( update_flags & update_inverse_jacobians, ExcNotInitialized()); - for (unsigned int v = 0; - v < stride && (stride == 1 || q + v < n_q_points_scalar); - ++v) + for (unsigned int v = 0, offset = qb * stride; + v < stride && (stride == 1 || offset < n_q_points_scalar); + ++v, ++offset) { - const unsigned int offset = qb * stride + v; - ETT::set_gradient(gradient, v, unit_gradients[offset]); + gradient_type unit_gradient; + ETT::set_gradient(gradient, v, unit_gradient); gradients[offset] = - apply_transformation(inverse_jacobian_ptr[offset].transpose(), - unit_gradients[offset]); + cell_type <= + internal::MatrixFreeFunctions::GeometryType::cartesian ? + apply_diagonal_transformation(inverse_jacobian_ptr[0], + unit_gradient) : + apply_transformation( + inverse_jacobian_ptr + [cell_type <= + internal::MatrixFreeFunctions::GeometryType::affine ? + 0 : + offset] + .transpose(), + unit_gradient); } } } @@ -1994,10 +2145,11 @@ FEPointEvaluation::evaluate_fast( template +template inline void FEPointEvaluation::evaluate_slow( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags) + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) { // slow path with FEValues Assert(fe_values.get() != nullptr, @@ -2076,10 +2228,11 @@ FEPointEvaluation::evaluate_slow( template +template void FEPointEvaluation::evaluate( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &evaluation_flags) + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &evaluation_flags) { if (!is_reinitialized) reinit(); @@ -2095,7 +2248,22 @@ FEPointEvaluation::evaluate( AssertDimension(solution_values.size(), fe->dofs_per_cell); if (fast_path) - evaluate_fast(solution_values, evaluation_flags); + { + if (use_face_path) + { + if (use_linear_path) + evaluate_fast(solution_values, evaluation_flags); + else + evaluate_fast(solution_values, evaluation_flags); + } + else + { + if (use_linear_path) + evaluate_fast(solution_values, evaluation_flags); + else + evaluate_fast(solution_values, evaluation_flags); + } + } else evaluate_slow(solution_values, evaluation_flags); } @@ -2103,15 +2271,17 @@ FEPointEvaluation::evaluate( template +template inline void FEPointEvaluation::compute_integrate_fast( - const EvaluationFlags::EvaluationFlags &integration_flags, - const unsigned int n_shapes, - const unsigned int qb, - const vectorized_value_type &value, - const interface_vectorized_gradient_type &gradient) + const EvaluationFlags::EvaluationFlags &integration_flags, + const unsigned int n_shapes, + const unsigned int qb, + const vectorized_value_type value, + const interface_vectorized_gradient_type gradient, + vectorized_value_type *solution_values_vectorized_linear) { - if (use_face_path) + if (is_face_path) { if (integration_flags & EvaluationFlags::gradients) { @@ -2158,6 +2328,7 @@ FEPointEvaluation::compute_integrate_fast( Assert(false, ExcInternalError()); internal::integrate_tensor_product_value_and_gradient< + is_linear, dim - 1, VectorizedArrayType, vectorized_value_type, @@ -2165,47 +2336,51 @@ FEPointEvaluation::compute_integrate_fast( n_shapes, value_face.data(), gradient_in_face, - solution_renumbered_vectorized.data(), + is_linear ? solution_values_vectorized_linear : + solution_renumbered_vectorized.data(), unit_point_faces_ptr[qb], - polynomials_are_hat_functions, qb != 0); } else - internal::integrate_tensor_product_value( shapes_faces.data() + qb * n_shapes, n_shapes, value, - solution_renumbered_vectorized.data(), + is_linear ? solution_values_vectorized_linear : + solution_renumbered_vectorized.data(), unit_point_faces_ptr[qb], - polynomials_are_hat_functions, qb != 0); } else { if (integration_flags & EvaluationFlags::gradients) internal::integrate_tensor_product_value_and_gradient< + is_linear, dim, VectorizedArrayType, vectorized_value_type>(shapes.data() + qb * n_shapes, n_shapes, &value, gradient, - solution_renumbered_vectorized.data(), + is_linear ? + solution_values_vectorized_linear : + solution_renumbered_vectorized.data(), unit_point_ptr[qb], - polynomials_are_hat_functions, qb != 0); else - internal::integrate_tensor_product_value( shapes.data() + qb * n_shapes, n_shapes, value, - solution_renumbered_vectorized.data(), + is_linear ? solution_values_vectorized_linear : + solution_renumbered_vectorized.data(), unit_point_ptr[qb], - polynomials_are_hat_functions, qb != 0); } } @@ -2213,31 +2388,40 @@ FEPointEvaluation::compute_integrate_fast( template +template inline void FEPointEvaluation::finish_integrate_fast( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags, + vectorized_value_type *solution_values_vectorized_linear) { - std::fill(solution_values.begin(), solution_values.end(), ScalarNumber()); + if (fe->n_components() > n_components) + for (unsigned int i = 0; i < solution_values.size(); ++i) + solution_values[i] = 0; + + const unsigned int dofs_per_comp = + is_linear ? Utilities::pow(2, dim) : dofs_per_component; + for (unsigned int comp = 0; comp < n_components; ++comp) { const std::size_t offset = - (component_in_base_element + comp) * dofs_per_component; + (component_in_base_element + comp) * dofs_per_comp; - if (use_face_path) + if (is_face_path) { - const unsigned int size_input = 2 * dofs_per_component_face; + const unsigned int dofs_per_comp_face = + is_linear ? Utilities::pow(2, dim - 1) : dofs_per_component_face; + + const unsigned int size_input = 2 * dofs_per_comp_face; ScalarNumber *input = scratch_data_scalar.begin(); ScalarNumber *output = input + size_input; - for (unsigned int i = 0; i < 2 * dofs_per_component_face; ++i) - { - VectorizedArrayType vectorized_input; - ETT::write_value(vectorized_input, - comp, + for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i) + input[i] = + ETT::sum_value(comp, + is_linear ? + *(solution_values_vectorized_linear + i) : solution_renumbered_vectorized[i]); - input[i] = vectorized_input.sum(); - } internal::FEFaceNormalEvaluationImpl:: template interpolate(1, @@ -2247,33 +2431,36 @@ FEPointEvaluation::finish_integrate_fast( output, current_face_number); - if (renumber.empty()) - for (unsigned int i = 0; i < dofs_per_component; ++i) - solution_values[i + offset] = output[i]; + if (is_linear || renumber.empty()) + { + for (unsigned int i = 0; i < dofs_per_comp; ++i) + solution_values[i + offset] = output[i]; + } else - for (unsigned int i = 0; i < dofs_per_component; ++i) - solution_values[renumber[i + offset]] = output[i]; + { + const unsigned int *renumber_ptr = renumber.data() + offset; + for (unsigned int i = 0; i < dofs_per_comp; ++i) + solution_values[renumber_ptr[i]] = output[i]; + } } else { - if (renumber.empty()) - for (unsigned int i = 0; i < dofs_per_component; ++i) - { - VectorizedArrayType result; - ETT::write_value(result, - comp, - solution_renumbered_vectorized[i]); - solution_values[i + offset] = result.sum(); - } + if (is_linear || renumber.empty()) + { + for (unsigned int i = 0; i < dofs_per_comp; ++i) + solution_values[i + offset] = + ETT::sum_value(comp, + is_linear ? + *(solution_values_vectorized_linear + i) : + solution_renumbered_vectorized[i]); + } else - for (unsigned int i = 0; i < dofs_per_component; ++i) - { - VectorizedArrayType result; - ETT::write_value(result, - comp, - solution_renumbered_vectorized[i]); - solution_values[renumber[i + offset]] = result.sum(); - } + { + const unsigned int *renumber_ptr = renumber.data() + offset; + for (unsigned int i = 0; i < dofs_per_comp; ++i) + solution_values[renumber_ptr[i]] = + ETT::sum_value(comp, solution_renumbered_vectorized[i]); + } } } } @@ -2281,107 +2468,103 @@ FEPointEvaluation::finish_integrate_fast( template -template +template inline void FEPointEvaluation::integrate_fast( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags) { - // fast path with tensor product integration - if (use_face_path) - { - if (solution_renumbered_vectorized.size() != 2 * dofs_per_component_face) - solution_renumbered_vectorized.resize(2 * dofs_per_component_face); - } - else - { - if (solution_renumbered_vectorized.size() != dofs_per_component) - solution_renumbered_vectorized.resize(dofs_per_component); - } + // zero out lanes of incomplete last quadrature point batch + if constexpr (stride == 1) + if (const unsigned int n_filled_lanes = + n_q_points_scalar & (n_lanes_internal - 1); + n_filled_lanes > 0) + { + if (integration_flags & EvaluationFlags::values) + for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v) + ETT::set_zero_value(values.back(), v); + if (integration_flags & EvaluationFlags::gradients) + for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v) + ETT::set_zero_gradient(gradients.back(), v); + } - // loop over quadrature batches qb / points q - 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 = - q + n_lanes_user_interface > n_q_points_scalar; + std::array + solution_values_vectorized_linear; + // loop over quadrature batches qb + const unsigned int n_shapes = is_linear ? 2 : poly.size(); + + const bool cartesian_cell = + cell_type <= internal::MatrixFreeFunctions::GeometryType::cartesian; + const bool affine_cell = + cell_type <= internal::MatrixFreeFunctions::GeometryType::affine; + for (unsigned int qb = 0; qb < n_q_batches; ++qb) + { 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) - { - 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, offset = qb * stride; + v < stride && (stride == 1 || offset < n_q_points_scalar); + ++v, ++offset) + ETT::get_value(value, + v, + do_JxW ? values[offset] * JxW_ptr[offset] : + values[offset]); - for (unsigned int v = 0; - v < stride && (stride == 1 || q + v < n_q_points_scalar); - ++v) - { - const unsigned int offset = qb * stride + v; - if (do_JxW) - values[offset] *= JxW_ptr[offset]; - ETT::get_value(value, v, values[offset]); - } - } if (integration_flags & EvaluationFlags::gradients) - { - // zero out lanes of incomplete last quadrature point batch - if (incomplete_last_batch) - { - 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); - } - - for (unsigned int v = 0; - v < stride && (stride == 1 || q + v < n_q_points_scalar); - ++v) - { - const unsigned int offset = qb * stride + v; - if (do_JxW) - gradients[offset] *= JxW_ptr[offset]; - ETT::get_gradient( - gradient, - v, - apply_transformation(inverse_jacobian_ptr[offset], - gradients[offset])); - } - } + for (unsigned int v = 0, offset = qb * stride; + v < stride && (stride == 1 || offset < n_q_points_scalar); + ++v, ++offset) + { + const auto grad_w = + do_JxW ? gradients[offset] * JxW_ptr[offset] : gradients[offset]; + ETT::get_gradient( + gradient, + v, + cartesian_cell ? + apply_diagonal_transformation(inverse_jacobian_ptr[0], grad_w) : + apply_transformation( + inverse_jacobian_ptr[affine_cell ? 0 : offset], grad_w)); + } - compute_integrate_fast(integration_flags, n_shapes, qb, value, gradient); + compute_integrate_fast( + integration_flags, + n_shapes, + qb, + value, + gradient, + solution_values_vectorized_linear.data()); } // add between the lanes and write into the result - finish_integrate_fast(solution_values, integration_flags); + finish_integrate_fast( + solution_values, + integration_flags, + solution_values_vectorized_linear.data()); } template -template +template inline void FEPointEvaluation::integrate_slow( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags) { // slow path with FEValues Assert(fe_values.get() != nullptr, ExcMessage( "Not initialized. Please call FEPointEvaluation::reinit()!")); - std::fill(solution_values.begin(), solution_values.end(), 0.0); + for (unsigned int i = 0; i < solution_values.size(); ++i) + solution_values[i] = 0; const std::size_t n_points = fe_values->get_quadrature().size(); @@ -2447,28 +2630,24 @@ FEPointEvaluation::integrate_slow( template -template +template void FEPointEvaluation::do_integrate( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags) { if (!is_reinitialized) reinit(); - if (n_q_points == 0) // no evaluation points provided - { - std::fill(solution_values.begin(), solution_values.end(), 0.0); - return; - } - Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented()); - if (!((integration_flags & EvaluationFlags::values) || + if (n_q_points == 0 || // no evaluation points provided + !((integration_flags & EvaluationFlags::values) || (integration_flags & EvaluationFlags::gradients))) // no integration flags { - std::fill(solution_values.begin(), solution_values.end(), 0.0); + for (unsigned int i = 0; i < solution_values.size(); ++i) + solution_values[i] = 0; return; } @@ -2479,7 +2658,26 @@ FEPointEvaluation::do_integrate( AssertDimension(solution_values.size(), fe->dofs_per_cell); if (fast_path) - integrate_fast(solution_values, integration_flags); + { + if (use_face_path) + { + if (use_linear_path) + integrate_fast(solution_values, + integration_flags); + else + integrate_fast(solution_values, + integration_flags); + } + else + { + if (use_linear_path) + integrate_fast(solution_values, + integration_flags); + else + integrate_fast(solution_values, + integration_flags); + } + } else integrate_slow(solution_values, integration_flags); } @@ -2487,10 +2685,11 @@ FEPointEvaluation::do_integrate( template +template void FEPointEvaluation::integrate( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags) { do_integrate(solution_values, integration_flags); } @@ -2498,10 +2697,11 @@ FEPointEvaluation::integrate( template +template void FEPointEvaluation::test_and_sum( - const ArrayView &solution_values, - const EvaluationFlags::EvaluationFlags &integration_flags) + const StridedArrayView &solution_values, + const EvaluationFlags::EvaluationFlags &integration_flags) { do_integrate(solution_values, integration_flags); } @@ -2532,22 +2732,6 @@ inline const typename FEPointEvaluation:: -template -inline const typename FEPointEvaluation:: - gradient_type & - FEPointEvaluation::get_unit_gradient( - const unsigned int point_index) const -{ - Assert(fast_path, - ExcMessage("Unit gradients are currently only implemented for tensor " - "product finite elements combined with MappingQ " - "mappings")); - AssertIndexRange(point_index, unit_gradients.size()); - return unit_gradients[point_index]; -} - - - template inline void FEPointEvaluation::submit_value( @@ -2582,7 +2766,10 @@ FEPointEvaluation::jacobian( internal::FEPointEvaluation:: ExcFEPointEvaluationAccessToUninitializedMappingField( "update_jacobians")); - return jacobian_ptr[point_index]; + return jacobian_ptr[cell_type <= ::dealii::internal::MatrixFreeFunctions:: + GeometryType::affine ? + 0 : + point_index]; } @@ -2597,7 +2784,11 @@ FEPointEvaluation::inverse_jacobian( internal::FEPointEvaluation:: ExcFEPointEvaluationAccessToUninitializedMappingField( "update_inverse_jacobians")); - return inverse_jacobian_ptr[point_index]; + return inverse_jacobian_ptr + [cell_type <= + ::dealii::internal::MatrixFreeFunctions::GeometryType::affine ? + 0 : + point_index]; } diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 837510b986..e6f1fa7b95 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -2272,14 +2272,16 @@ namespace internal * Specializes @p evaluate_tensor_product_value_and_gradient() for linear * polynomials which massively reduces the necessary instructions. */ - template + template inline std::array::type, dim + n_values> evaluate_tensor_product_value_and_gradient_linear( - const unsigned int n_shapes, - const Number *values, - const Point &p, - const std::vector &renumber = {}) + const Number *values, + const Point &p) { static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); static_assert(1 <= n_values && n_values <= 2, @@ -2287,15 +2289,12 @@ namespace internal using Number3 = typename ProductTypeNoPoint::type; + static_assert( + n_values == 1 || stride == 1, + "Either n_values or stride has to be one for correct data access!"); // If n_values > 1, we want to interpolate from a second array, // placed in the same array immediately after the main data. This // is used to interpolate normal derivatives onto faces. - const Number *values_2 = - n_values > 1 ? values + Utilities::fixed_power(n_shapes) : nullptr; - - AssertDimension(n_shapes, 2); - for (unsigned int i = 0; i < renumber.size(); ++i) - AssertDimension(renumber[i], i); std::array result; if (dim == 0) @@ -2303,23 +2302,23 @@ namespace internal // we only need the value on faces of a 1d element result[0] = values[0]; if (n_values > 1) - result[1] = values_2[0]; + result[1] = values[1]; } else if (dim == 1) { // gradient - result[0] = Number3(values[1] - values[0]); + result[0] = Number3(values[stride] - values[0]); // values result[1] = Number3(values[0]) + p[0] * result[0]; if (n_values > 1) - result[2] = Number3(values_2[0]) + p[0] * (values_2[1] - values_2[0]); + result[2] = Number3(values[2]) + p[0] * (values[3] - values[2]); } else if (dim == 2) { - const Number3 val10 = Number3(values[1] - values[0]); - const Number3 val32 = Number3(values[3] - values[2]); + const Number3 val10 = Number3(values[stride] - values[0]); + const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); const Number3 tmp0 = Number3(values[0]) + p[0] * val10; - const Number3 tmp1 = Number3(values[2]) + p[0] * val32; + const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; // gradient result[0] = val10 + p[1] * (val32 - val10); @@ -2331,25 +2330,25 @@ namespace internal if (n_values > 1) { const Number3 tmp0_2 = - Number3(values_2[0]) + p[0] * (values_2[1] - values_2[0]); + Number3(values[4]) + p[0] * (values[5] - values[4]); const Number3 tmp1_2 = - Number3(values_2[2]) + p[0] * (values_2[3] - values_2[0]); + Number3(values[6]) + p[0] * (values[7] - values[6]); result[3] = tmp0_2 + p[1] * (tmp1_2 - tmp0_2); } } else if (dim == 3) { - const Number3 val10 = Number3(values[1] - values[0]); - const Number3 val32 = Number3(values[3] - values[2]); + const Number3 val10 = Number3(values[stride] - values[0]); + const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); const Number3 tmp0 = Number3(values[0]) + p[0] * val10; - const Number3 tmp1 = Number3(values[2]) + p[0] * val32; + const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; const Number3 tmp10 = tmp1 - tmp0; const Number3 tmpy0 = tmp0 + p[1] * tmp10; - const Number3 val54 = Number3(values[5] - values[4]); - const Number3 val76 = Number3(values[7] - values[6]); - const Number3 tmp2 = Number3(values[4]) + p[0] * val54; - const Number3 tmp3 = Number3(values[6]) + p[0] * val76; + const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]); + const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]); + const Number3 tmp2 = Number3(values[4 * stride]) + p[0] * val54; + const Number3 tmp3 = Number3(values[6 * stride]) + p[0] * val76; const Number3 tmp32 = tmp3 - tmp2; const Number3 tmpy1 = tmp2 + p[1] * tmp32; @@ -2419,8 +2418,8 @@ namespace internal std::array result; if (d_linear) { - result = evaluate_tensor_product_value_and_gradient_linear( - poly.size(), values.data(), p, renumber); + result = + evaluate_tensor_product_value_and_gradient_linear(values.data(), p); } else { @@ -2548,23 +2547,15 @@ namespace internal - template + template inline typename ProductTypeNoPoint::type - evaluate_tensor_product_value_linear( - const unsigned int n_shapes, - const Number *values, - const Point &p, - const std::vector &renumber = {}) + evaluate_tensor_product_value_linear(const Number *values, + const Point &p) { - (void)n_shapes; static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented"); using Number3 = typename ProductTypeNoPoint::type; - AssertDimension(n_shapes, 2); - for (unsigned int i = 0; i < renumber.size(); ++i) - AssertDimension(renumber[i], i); - if (dim == 0) { // we only need the value on faces of a 1d element @@ -2572,28 +2563,28 @@ namespace internal } else if (dim == 1) { - return Number3(values[0]) + p[0] * Number3(values[1] - values[0]); + return Number3(values[0]) + p[0] * Number3(values[stride] - values[0]); } else if (dim == 2) { - const Number3 val10 = Number3(values[1] - values[0]); - const Number3 val32 = Number3(values[3] - values[2]); + const Number3 val10 = Number3(values[stride] - values[0]); + const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); const Number3 tmp0 = Number3(values[0]) + p[0] * val10; - const Number3 tmp1 = Number3(values[2]) + p[0] * val32; + const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; return tmp0 + p[1] * (tmp1 - tmp0); } else if (dim == 3) { - const Number3 val10 = Number3(values[1] - values[0]); - const Number3 val32 = Number3(values[3] - values[2]); + const Number3 val10 = Number3(values[stride] - values[0]); + const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]); const Number3 tmp0 = Number3(values[0]) + p[0] * val10; - const Number3 tmp1 = Number3(values[2]) + p[0] * val32; + const Number3 tmp1 = Number3(values[2 * stride]) + p[0] * val32; const Number3 tmpy0 = tmp0 + p[1] * (tmp1 - tmp0); - const Number3 val54 = Number3(values[5] - values[4]); - const Number3 val76 = Number3(values[7] - values[6]); - const Number3 tmp2 = Number3(values[4]) + p[0] * val54; - const Number3 tmp3 = Number3(values[6]) + p[0] * val76; + const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]); + const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]); + const Number3 tmp2 = Number3(values[4 * stride]) + p[0] * val54; + const Number3 tmp3 = Number3(values[6 * stride]) + p[0] * val76; const Number3 tmpy1 = tmp2 + p[1] * (tmp3 - tmp2); return tmpy0 + p[2] * (tmpy1 - tmpy0); @@ -2617,10 +2608,7 @@ namespace internal typename ProductTypeNoPoint::type result; if (d_linear) { - result = evaluate_tensor_product_value_linear(poly.size(), - values.data(), - p, - renumber); + result = evaluate_tensor_product_value_linear(values.data(), p); } else { @@ -3049,19 +3037,15 @@ namespace internal int n_values = 1> inline void integrate_add_tensor_product_value_and_gradient_linear( - const unsigned int n_shapes, const Number2 *value, const Tensor<1, dim, Number2> &gradient, Number2 *values, const Point &p) { - (void)n_shapes; static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented"); static_assert(1 <= n_values && n_values <= 2, "Only n_values=1,2 implemented"); - AssertDimension(n_shapes, 2); - // Note that 'add' is a template argument, so the compiler will remove // these checks if (dim == 0) @@ -3210,7 +3194,11 @@ namespace internal * function depending on if values should be added to or set and if * polynomials are linear. */ - template + template inline void integrate_tensor_product_value_and_gradient( const dealii::ndarray *shapes, @@ -3219,7 +3207,6 @@ namespace internal const Tensor<1, dim, Number2> &gradient, Number2 *values, const Point &p, - const bool is_linear, const bool do_add) { if (do_add) @@ -3230,7 +3217,7 @@ namespace internal Number, Number2, true, - n_values>(n_shapes, value, gradient, values, p); + n_values>(value, gradient, values, p); else internal::integrate_add_tensor_product_value_and_gradient_shapes< dim, @@ -3247,7 +3234,7 @@ namespace internal Number, Number2, false, - n_values>(n_shapes, value, gradient, values, p); + n_values>(value, gradient, values, p); else internal::integrate_add_tensor_product_value_and_gradient_shapes< dim, @@ -3382,16 +3369,12 @@ namespace internal */ template inline void - integrate_add_tensor_product_value_linear(const unsigned int n_shapes, - const Number2 &value, + integrate_add_tensor_product_value_linear(const Number2 &value, Number2 *values, const Point &p) { - (void)n_shapes; static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented"); - AssertDimension(n_shapes, 2); - if (dim == 0) { if (add) @@ -3481,14 +3464,13 @@ namespace internal * function depending on if values should be added to or set and if * polynomials are linear. */ - template + template inline void integrate_tensor_product_value(const dealii::ndarray *shapes, const unsigned int n_shapes, const Number2 &value, Number2 *values, const Point &p, - const bool is_linear, const bool do_add) { if (do_add) @@ -3497,8 +3479,7 @@ namespace internal internal::integrate_add_tensor_product_value_linear(n_shapes, - value, + true>(value, values, p); else @@ -3516,8 +3497,7 @@ namespace internal internal::integrate_add_tensor_product_value_linear(n_shapes, - value, + false>(value, values, p); else diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h index eb094cc432..9e674a4f71 100644 --- a/include/deal.II/non_matching/mapping_info.h +++ b/include/deal.II/non_matching/mapping_info.h @@ -20,6 +20,7 @@ #include +#include "deal.II/base/floating_point_comparator.h" #include #include #include @@ -31,6 +32,8 @@ #include #include +#include + #include @@ -179,6 +182,49 @@ namespace NonMatching } } }; + + template + dealii::internal::MatrixFreeFunctions::GeometryType + compute_geometry_type( + const double diameter, + const std::vector> + &inverse_jacobians) + { + const auto jac_0 = inverse_jacobians[0]; + const double zero_tolerance_double = + 1. / diameter * std::numeric_limits::epsilon() * 1024.; + bool jacobian_constant = true; + for (unsigned int q = 1; q < inverse_jacobians.size(); ++q) + { + const DerivativeForm<1, dim, spacedim> &jac = inverse_jacobians[q]; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < spacedim; ++e) + if (std::fabs(jac_0[d][e] - jac[d][e]) > zero_tolerance_double) + jacobian_constant = false; + if (!jacobian_constant) + break; + } + + // check whether the Jacobian is diagonal to machine + // accuracy + bool cell_cartesian = jacobian_constant; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + if (d != e) + if (std::fabs(jac_0[d][e]) > zero_tolerance_double) + { + cell_cartesian = false; + break; + } + + // return cell type + if (cell_cartesian) + return dealii::internal::MatrixFreeFunctions::GeometryType::cartesian; + else if (jacobian_constant) + return dealii::internal::MatrixFreeFunctions::GeometryType::affine; + else + return dealii::internal::MatrixFreeFunctions::GeometryType::general; + } } // namespace internal /** @@ -412,26 +458,44 @@ namespace NonMatching boost::signals2::connection connect_is_reinitialized(const std::function &set_is_reinitialized); + /** + * Compute the geometry index offset of the current cell/face. + */ + template + unsigned int + compute_geometry_index_offset(const unsigned int cell_index, + const unsigned int face_number) const; + /** * 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_unit_point_index_offset(const unsigned int geometry_index) 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; + compute_data_index_offset(const unsigned int geometry_index) const; + + /** + * Compute the data index offset for the current cell/face. + */ + unsigned int + compute_compressed_data_index_offset( + const unsigned int geometry_index) 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; + get_n_q_points_unvectorized(const unsigned int geometry_index) const; + + /** + * Get cell geometry type. + */ + dealii::internal::MatrixFreeFunctions::GeometryType + get_cell_type(const unsigned int geometry_index) const; /** * Return cell iterator. @@ -503,7 +567,9 @@ namespace NonMatching const unsigned int n_q_points, const unsigned int n_q_points_unvectorized, const MappingData &mapping_data, - const std::vector &weights); + const std::vector &weights, + const unsigned int compressed_unit_point_index_offset, + const bool affine_cell); /** * Compute the compressed cell index. @@ -512,11 +578,18 @@ namespace NonMatching compute_compressed_cell_index(const unsigned int cell_index) const; /** - * Compute the geometry index offset of the current cell/face. + * Compute the mapping information for cells/surface. */ - unsigned int - compute_geometry_index_offset(const unsigned int cell_index, - const unsigned int face_number) const; + template + void + do_reinit_cells( + const ContainerType &cell_iterator_range, + const std::vector &quadrature_vector, + const unsigned int n_unfiltered_cells, + const std::function< + void(const typename Triangulation::cell_iterator &cell, + const QuadratureType &quadrature, + MappingData &mapping_data)> &compute_mapping_data); /** * Enum class for reinitialized states. @@ -552,7 +625,7 @@ namespace NonMatching /** * Offset to point to the first unit point of a cell/face. */ - AlignedVector unit_points_index; + std::vector unit_points_index; /** * A pointer to the internal data of the underlying mapping. @@ -581,10 +654,22 @@ namespace NonMatching const AdditionalData additional_data; /** - * Stores the index offset into the arrays @p JxW_values, @p jacobians, - * @p inverse_jacobians and @p normal_vectors. + * Stores whether a cell is Cartesian (cell type 0), has constant + * transform data (Jacobians) (cell type 1), or is general (cell type + * 3). Type 2 is only used for faces and no cells are assigned this + * value. + */ + std::vector cell_type; + + /** + * Stores the index offset into the arrays @p JxW_values and @p normal_vectors. + */ + std::vector data_index_offsets; + + /** + * Stores the index offset into the arrays @p jacobians and @p inverse_jacobians. */ - AlignedVector data_index_offsets; + std::vector compressed_data_index_offsets; /** * The storage of the Jacobian determinant times the quadrature weight on @@ -605,7 +690,7 @@ namespace NonMatching * 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. + * Indexed by @p compressed_data_index_offsets. */ AlignedVector> jacobians; @@ -614,7 +699,7 @@ namespace NonMatching * the inverse Jacobians of the transformation from the * unit to the real cell. * - * Indexed by @p data_index_offsets. + * Indexed by @p compressed_data_index_offsets. */ AlignedVector> inverse_jacobians; @@ -718,6 +803,8 @@ namespace NonMatching n_q_points_unvectorized.clear(); unit_points_index.clear(); data_index_offsets.clear(); + compressed_data_index_offsets.clear(); + cell_type.clear(); } @@ -783,12 +870,32 @@ namespace NonMatching internal_mapping_data, mapping_data); + // check for cartesian/affine cell + if (!quadrature.empty() && + update_flags_mapping & UpdateFlags::update_inverse_jacobians) + { + cell_type.push_back( + internal::compute_geometry_type(cell->diameter(), + mapping_data.inverse_jacobians)); + } + else + cell_type.push_back( + dealii::internal::MatrixFreeFunctions::GeometryType::general); + // store mapping data - store_mapping_data(0, - n_q_points_data, - n_q_points_unvectorized[0], - mapping_data, - quadrature.get_weights()); + store_mapping_data( + 0, + n_q_points_data, + n_q_points_unvectorized[0], + mapping_data, + quadrature.get_weights(), + 0, + cell_type.back() <= + dealii::internal::MatrixFreeFunctions::GeometryType::affine); + + unit_points_index.push_back(0); + data_index_offsets.push_back(0); + compressed_data_index_offsets.push_back(0); state = State::single_cell; is_reinitialized(); @@ -820,12 +927,16 @@ namespace NonMatching template - template + template void - MappingInfo::reinit_cells( - const ContainerType &cell_iterator_range, - const std::vector> &quadrature_vector, - const unsigned int n_unfiltered_cells) + MappingInfo::do_reinit_cells( + const ContainerType &cell_iterator_range, + const std::vector &quadrature_vector, + const unsigned int n_unfiltered_cells, + const std::function< + void(const typename Triangulation::cell_iterator &cell, + const QuadratureType &quadrature, + MappingData &mapping_data)> &compute_mapping_data) { clear(); @@ -839,6 +950,8 @@ namespace NonMatching n_q_points_unvectorized.reserve(n_cells); + cell_type.reserve(n_cells); + if (additional_data.store_cells) cell_level_and_indices.resize(n_cells); @@ -873,9 +986,8 @@ namespace NonMatching cell_index_to_compressed_cell_index.resize(n_unfiltered_cells, numbers::invalid_unsigned_int); - MappingData mapping_data; - CellSimilarity::Similarity cell_similarity = - CellSimilarity::Similarity::none; + MappingData mapping_data; + MappingData mapping_data_last_cell; unsigned int cell_index = 0; for (const auto &cell : cell_iterator_range) { @@ -885,32 +997,97 @@ namespace NonMatching cell_level_and_indices[cell_index] = {cell->level(), cell->index()}; } + const auto &quadrature = quadrature_vector[cell_index]; + const bool empty = quadrature.empty(); + // 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()); + quadrature.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); + compute_mapping_data(cell, quadrature, mapping_data); // store mapping data const unsigned int n_q_points_data = compute_n_q_points(n_q_points_unvectorized[cell_index]); + + // check for cartesian/affine cell + if (!empty && + update_flags_mapping & UpdateFlags::update_inverse_jacobians) + { + cell_type.push_back( + internal::compute_geometry_type(cell->diameter(), + mapping_data.inverse_jacobians)); + } + else + cell_type.push_back( + dealii::internal::MatrixFreeFunctions::GeometryType::general); + + if (cell_index > 0) + { + // check if current and last cell are affine + const bool affine_cells = + cell_type[cell_index] <= + dealii::internal::MatrixFreeFunctions::affine && + cell_type[cell_index - 1] <= + dealii::internal::MatrixFreeFunctions::affine; + + // create a comparator to compare inverse Jacobian of current + // and last cell + FloatingPointComparator comparator( + 1e4 / cell->diameter() * std::numeric_limits::epsilon() * + 1024.); + + // we can only compare if current and last cell have at least + // one quadrature point and both cells are at least affine + const auto comparison_result = + (!affine_cells || mapping_data.inverse_jacobians.empty() || + mapping_data_last_cell.inverse_jacobians.empty()) ? + FloatingPointComparator::ComparisonResult::less : + comparator.compare(mapping_data.inverse_jacobians[0], + mapping_data_last_cell.inverse_jacobians[0]); + + // we can compress the Jacobians and inverse Jacobians if + // inverse Jacobians are equal and cells are affine + if (affine_cells && + comparison_result == + FloatingPointComparator::ComparisonResult::equal) + { + compressed_data_index_offsets.push_back( + compressed_data_index_offsets.back()); + } + else + { + const unsigned int n_compressed_data_last_cell = + cell_type[cell_index - 1] <= + dealii::internal::MatrixFreeFunctions::affine ? + 1 : + compute_n_q_points( + n_q_points_unvectorized[cell_index - 1]); + + compressed_data_index_offsets.push_back( + compressed_data_index_offsets.back() + + n_compressed_data_last_cell); + } + } + else + compressed_data_index_offsets.push_back(0); + + // cache mapping_data from last cell + mapping_data_last_cell = mapping_data; + store_mapping_data(data_index_offsets[cell_index], n_q_points_data, n_q_points_unvectorized[cell_index], mapping_data, - quadrature_vector[cell_index].get_weights()); + quadrature.get_weights(), + compressed_data_index_offsets[cell_index], + cell_type[cell_index] <= + dealii::internal::MatrixFreeFunctions::affine); if (do_cell_index_compression) cell_index_to_compressed_cell_index[cell->active_cell_index()] = @@ -919,12 +1096,47 @@ namespace NonMatching ++cell_index; } + // TODO: release allocated memory from compressed data vectors + state = State::cell_vector; is_reinitialized(); } + template + template + void + MappingInfo::reinit_cells( + const ContainerType &cell_iterator_range, + const std::vector> &quadrature_vector, + const unsigned int n_unfiltered_cells) + { + auto compute_mapping_data_for_cells = + [&](const typename Triangulation::cell_iterator &cell, + const Quadrature &quadrature, + MappingData &mapping_data) { + CellSimilarity::Similarity cell_similarity = + CellSimilarity::Similarity::none; + internal::ComputeMappingDataHelper:: + compute_mapping_data_for_quadrature(mapping, + update_flags_mapping, + cell, + cell_similarity, + quadrature, + internal_mapping_data, + mapping_data); + }; + + do_reinit_cells>( + cell_iterator_range, + quadrature_vector, + n_unfiltered_cells, + compute_mapping_data_for_cells); + } + + + template template void @@ -933,8 +1145,6 @@ namespace NonMatching const std::vector> &quadrature_vector, const unsigned int n_unfiltered_cells) { - clear(); - Assert( additional_data.use_global_weights == false, ExcMessage( @@ -942,65 +1152,13 @@ namespace NonMatching Assert(additional_data.store_cells == false, ExcNotImplemented()); - do_cell_index_compression = - n_unfiltered_cells != numbers::invalid_unsigned_int; - if (update_flags_mapping & (update_JxW_values | update_normal_vectors)) update_flags_mapping |= update_covariant_transformation; - const unsigned int n_cells = quadrature_vector.size(); - AssertDimension(n_cells, - 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) - { - 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(); - - // 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]; - - // 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 + auto compute_mapping_data_for_surface = + [&](const typename Triangulation::cell_iterator &cell, + const ImmersedSurfaceQuadrature &quadrature, + MappingData &mapping_data) { internal::ComputeMappingDataHelper:: compute_mapping_data_for_immersed_surface_quadrature( mapping, @@ -1009,25 +1167,13 @@ namespace NonMatching 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, - quadrature_vector[cell_index].get_weights()); - - if (do_cell_index_compression) - cell_index_to_compressed_cell_index[cell->active_cell_index()] = - cell_index; - - ++cell_index; - } - - state = State::cell_vector; - is_reinitialized(); + do_reinit_cells, ImmersedSurfaceQuadrature>( + cell_iterator_range, + quadrature_vector, + n_unfiltered_cells, + compute_mapping_data_for_surface); } @@ -1065,6 +1211,8 @@ namespace NonMatching n_q_points_unvectorized.reserve(n_faces); + cell_type.reserve(n_faces); + // fill unit points index offset vector unit_points_index.resize(n_faces + 1); data_index_offsets.resize(n_faces + 1); @@ -1156,13 +1304,21 @@ namespace NonMatching internal_mapping_data, mapping_data); + cell_type.push_back( + dealii::internal::MatrixFreeFunctions::GeometryType::general); + + compressed_data_index_offsets.push_back( + data_index_offsets[current_face_index]); + 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, - quadrature_on_face.get_weights()); + quadrature_on_face.get_weights(), + data_index_offsets[current_face_index], + false); } if (do_cell_index_compression) cell_index_to_compressed_cell_index[cell->active_cell_index()] = @@ -1189,22 +1345,18 @@ namespace NonMatching template unsigned int MappingInfo::get_n_q_points_unvectorized( - const unsigned int cell_index, - const unsigned int face_number) const + const unsigned int geometry_index) const { - if (cell_index == numbers::invalid_unsigned_int && - face_number == numbers::invalid_unsigned_int) - { - Assert(state == State::single_cell, - ExcMessage( - "This mapping info is not reinitialized for a single cell!")); - return n_q_points_unvectorized[0]; - } - else - { - return n_q_points_unvectorized[compute_geometry_index_offset( - cell_index, face_number)]; - } + return n_q_points_unvectorized[geometry_index]; + } + + + template + dealii::internal::MatrixFreeFunctions::GeometryType + MappingInfo::get_cell_type( + const unsigned int geometry_index) const + { + return cell_type[geometry_index]; } @@ -1244,14 +1396,18 @@ namespace NonMatching template + template unsigned int MappingInfo::compute_geometry_index_offset( const unsigned int cell_index, const unsigned int face_number) const { + if (cell_index == numbers::invalid_unsigned_int) + return 0; + const unsigned int compressed_cell_index = compute_compressed_cell_index(cell_index); - if (face_number == numbers::invalid_unsigned_int) + if (!is_face) { Assert(state == State::cell_vector, ExcMessage( @@ -1347,7 +1503,9 @@ namespace NonMatching const unsigned int n_q_points, const unsigned int n_q_points_unvectorized, const MappingInfo::MappingData &mapping_data, - const std::vector &weights) + const std::vector &weights, + const unsigned int compressed_unit_point_index_offset, + const bool affine_cell) { const unsigned int n_lanes = dealii::internal::VectorizedArrayTrait::width(); @@ -1355,22 +1513,28 @@ namespace NonMatching for (unsigned int q = 0; q < n_q_points; ++q) { const unsigned int offset = unit_points_index_offset + q; + const unsigned int compressed_offset = + compressed_unit_point_index_offset + q; for (unsigned int v = 0; v < n_lanes && q * n_lanes + v < n_q_points_unvectorized; ++v) { - 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 (q == 0 || !affine_cell) + { + 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[compressed_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[compressed_offset][s][d], v) = + mapping_data.inverse_jacobians[q * n_lanes + v][s][d]; + } if (update_flags_mapping & UpdateFlags::update_JxW_values) { if (additional_data.use_global_weights) @@ -1447,7 +1611,7 @@ namespace NonMatching MappingInfo::get_unit_point( const unsigned int offset) const { - return &unit_points[offset]; + return unit_points.data() + offset; } @@ -1459,7 +1623,7 @@ namespace NonMatching MappingInfo::get_unit_point_faces( const unsigned int offset) const { - return &unit_points_faces[offset]; + return unit_points_faces.data() + offset; } @@ -1469,7 +1633,7 @@ namespace NonMatching MappingInfo::get_real_point( const unsigned int offset) const { - return &real_points[offset]; + return real_points.data() + offset; } @@ -1477,23 +1641,9 @@ namespace NonMatching template unsigned int MappingInfo::compute_unit_point_index_offset( - const unsigned int cell_index, - const unsigned int face_number) const + const unsigned int geometry_index) const { - if (cell_index == numbers::invalid_unsigned_int && - face_number == numbers::invalid_unsigned_int) - { - Assert(state == State::single_cell, - ExcMessage( - "This mapping info is not reinitialized for a single cell!")); - return 0; - } - else - { - const unsigned int offset = - compute_geometry_index_offset(cell_index, face_number); - return unit_points_index[offset]; - } + return unit_points_index[geometry_index]; } @@ -1501,23 +1651,18 @@ namespace NonMatching template unsigned int MappingInfo::compute_data_index_offset( - const unsigned int cell_index, - const unsigned int face_number) const + const unsigned int geometry_index) const { - if (cell_index == numbers::invalid_unsigned_int && - face_number == numbers::invalid_unsigned_int) - { - Assert(state == State::single_cell, - ExcMessage( - "This mapping info is not reinitialized for a single cell!")); - return 0; - } - else - { - const unsigned int offset = - compute_geometry_index_offset(cell_index, face_number); - return data_index_offsets[offset]; - } + return data_index_offsets[geometry_index]; + } + + + template + unsigned int + MappingInfo::compute_compressed_data_index_offset( + const unsigned int geometry_index) const + { + return compressed_data_index_offsets[geometry_index]; } @@ -1527,7 +1672,7 @@ namespace NonMatching MappingInfo::get_jacobian( const unsigned int offset) const { - return &jacobians[offset]; + return jacobians.data() + offset; } @@ -1537,7 +1682,7 @@ namespace NonMatching MappingInfo::get_inverse_jacobian( const unsigned int offset) const { - return &inverse_jacobians[offset]; + return inverse_jacobians.data() + offset; } @@ -1546,7 +1691,7 @@ namespace NonMatching MappingInfo::get_normal_vector( const unsigned int offset) const { - return &normal_vectors[offset]; + return normal_vectors.data() + offset; } @@ -1555,7 +1700,7 @@ namespace NonMatching inline const Number * MappingInfo::get_JxW(const unsigned int offset) const { - return &JxW_values[offset]; + return JxW_values.data() + offset; }