From 57da06d4aef4cb8551f1f8f2f42b810d75fdffde Mon Sep 17 00:00:00 2001 From: Maximilian Bergbauer Date: Fri, 24 Mar 2023 23:06:41 +0100 Subject: [PATCH] Compute values_of_array only once --- .../deal.II/matrix_free/fe_point_evaluation.h | 191 ++++++--- .../matrix_free/tensor_product_kernels.h | 397 ++++++++++++------ include/deal.II/non_matching/mapping_info.h | 4 +- 3 files changed, 393 insertions(+), 199 deletions(-) diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index 03e7592746..40b65738b1 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -649,6 +649,14 @@ private: void setup(const unsigned int first_selected_component); + /** + * Shared functionality of all @p reinit() functions. Resizes data fields and + * precomputes the @p shapes vector, holding the evaluation of 1D basis + * functions of tensor product polynomials, if necessary. + */ + void + do_reinit(); + /** * Number of quadrature points of the current cell/face. */ @@ -775,6 +783,12 @@ private: * Bool indicating if class is reinitialized and data vectors a resized. */ bool is_reinitialized; + + /** + * Vector containing tensor product shape functions evaluated (during + * reinit()) at the vectorized unit points. + */ + AlignedVector, 2, dim>> shapes; }; // ----------------------- template and inline function ---------------------- @@ -830,6 +844,8 @@ FEPointEvaluation::setup( AssertIndexRange(first_selected_component + n_components, fe->n_components() + 1); + shapes.reserve(100); + bool same_base_element = true; unsigned int base_element_number = 0; component_in_base_element = 0; @@ -914,14 +930,7 @@ FEPointEvaluation::reinit( fe_values->reinit(cell); } - const_cast(n_q_points) = unit_points.size(); - - 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()); - - is_reinitialized = true; + do_reinit(); } @@ -934,16 +943,7 @@ FEPointEvaluation::reinit( current_cell_index = cell_index; current_face_number = numbers::invalid_unsigned_int; - const_cast(n_q_points) = - mapping_info->get_unit_points(current_cell_index, current_face_number) - .size(); - - 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()); - - is_reinitialized = true; + do_reinit(); } @@ -957,15 +957,49 @@ FEPointEvaluation::reinit( current_cell_index = cell_index; current_face_number = face_number; - const_cast(n_q_points) = - mapping_info->get_unit_points(current_cell_index, current_face_number) - .size(); + do_reinit(); +} + + + +template +void +FEPointEvaluation::do_reinit() +{ + const auto unit_points = + mapping_info->get_unit_points(current_cell_index, current_face_number); + + const_cast(n_q_points) = unit_points.size(); 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) + { + const std::size_t n_points = unit_points.size(); + const std::size_t n_lanes = VectorizedArray::size(); + const std::size_t n_batches = + n_points / n_lanes + (n_points % n_lanes > 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); + } + } + is_reinitialized = true; } @@ -977,6 +1011,9 @@ FEPointEvaluation::evaluate( const ArrayView & solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flag) { + if (!is_reinitialized) + reinit(numbers::invalid_unsigned_int); + if (n_q_points == 0) return; @@ -986,14 +1023,6 @@ FEPointEvaluation::evaluate( fast_path) { // fast path with tensor product evaluation - const auto unit_points = - mapping_info->get_unit_points(current_cell_index, current_face_number); - - // we need to call reinit() here if we reuse the same MappingInfo object - // for several FEPointEvaluation objects to resize the data fields - if (!is_reinitialized || n_q_points != unit_points.size()) - reinit(numbers::invalid_unsigned_int); - if (solution_renumbered.size() != dofs_per_component) solution_renumbered.resize(dofs_per_component); for (unsigned int comp = 0; comp < n_components; ++comp) @@ -1011,23 +1040,41 @@ FEPointEvaluation::evaluate( 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 = VectorizedArray::size(); - for (unsigned int i = 0; i < n_points; i += n_lanes) + 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]; - // compute - const auto val_and_grad = - internal::evaluate_tensor_product_value_and_gradient( - poly, - solution_renumbered, - vectorized_points, - polynomials_are_hat_functions); + const unsigned int n_shapes = poly.size(); + const auto val_and_grad = [&]() { + if (polynomials_are_hat_functions) + { + // 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]; + + return internal:: + evaluate_tensor_product_value_and_gradient_linear( + poly, solution_renumbered, vectorized_points); + } + else + return internal:: + evaluate_tensor_product_value_and_gradient_shapes< + dim, + value_type, + VectorizedArray>( + make_array_view(shapes.begin() + qb * n_shapes, + shapes.begin() + (qb * n_shapes + n_shapes)), + poly.size(), + solution_renumbered); + }(); // convert back to standard format if (evaluation_flag & EvaluationFlags::values) @@ -1048,9 +1095,6 @@ FEPointEvaluation::evaluate( Number>::set_gradient(val_and_grad.second, j, unit_gradients[i + j]); - const auto &mapping_data = - mapping_info->get_mapping_data(current_cell_index, - current_face_number); gradients[i + j] = apply_transformation( mapping_data.inverse_jacobians[i + j].transpose(), unit_gradients[i + j]); @@ -1134,13 +1178,6 @@ FEPointEvaluation::integrate( fast_path) { // fast path with tensor product integration - const auto unit_points = - mapping_info->get_unit_points(current_cell_index, current_face_number); - - // we need to call reinit() here if we reuse the same MappingInfo object - // for several FEPointEvaluation objects to resize the data fields - if (!is_reinitialized || n_q_points != unit_points.size()) - reinit(numbers::invalid_unsigned_int); if (integration_flags & EvaluationFlags::values) AssertIndexRange(n_q_points, values.size() + 1); @@ -1156,16 +1193,15 @@ FEPointEvaluation::integrate( n_components, VectorizedArray>::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 = VectorizedArray::size(); - for (unsigned int i = 0; i < n_points; i += n_lanes) + 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]; - typename internal::ProductTypeNoPoint>::type value = {}; @@ -1184,9 +1220,6 @@ FEPointEvaluation::integrate( if (integration_flags & EvaluationFlags::gradients) for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) { - const auto &mapping_data = - mapping_info->get_mapping_data(current_cell_index, - current_face_number); gradients[i + j] = apply_transformation(mapping_data.inverse_jacobians[i + j], gradients[i + j]); @@ -1196,12 +1229,34 @@ FEPointEvaluation::integrate( } // compute - internal::integrate_add_tensor_product_value_and_gradient( - poly, - value, - gradient, - vectorized_points, - solution_renumbered_vectorized); + 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 && i + j < n_points; ++j) + for (unsigned int d = 0; d < dim; ++d) + vectorized_points[d][j] = unit_points[i + j][d]; + + internal::integrate_add_tensor_product_value_and_gradient_linear( + poly, + value, + gradient, + solution_renumbered_vectorized, + vectorized_points); + } + else + internal::integrate_add_tensor_product_value_and_gradient_shapes< + dim, + VectorizedArray, + 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); } // add between the lanes and write into the result diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 18f0faddee..74b7279fe2 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -2982,14 +2982,37 @@ namespace internal + /** + * Computes the values and derivatives of the 1d polynomials @p poly at the + * specified point @p p and stores it in @p shapes. + */ + template + inline void + compute_values_of_array( + ArrayView> & shapes, + const std::vector> &poly, + const Point & p) + { + const int n_shapes = poly.size(); + + // Evaluate 1d polynomials and their derivatives + std::array point; + for (unsigned int d = 0; d < dim; ++d) + point[d] = p[d]; + for (int i = 0; i < n_shapes; ++i) + poly[i].values_of_array(point, 1, &shapes[i][0]); + } + + + /** * Interpolate inner dimensions of tensor product shape functions. */ template inline std::array::type, 3> - do_interpolate_xy(const std::vector & values, - const std::vector & renumber, - const dealii::ndarray &shapes, + do_interpolate_xy(const std::vector & values, + const std::vector & renumber, + const ArrayView> &shapes, const int n_shapes_runtime, int & i) { @@ -3035,122 +3058,27 @@ namespace internal /** - * Compute the polynomial interpolation of a tensor product shape function - * $\varphi_i$ given a vector of coefficients $u_i$ in the form - * $u_h(\mathbf{x}) = \sum_{i=1}^{k^d} \varphi_i(\mathbf{x}) u_i$. The shape - * functions $\varphi_i(\mathbf{x}) = - * \prod_{d=1}^{\text{dim}}\varphi_{i_d}^\text{1d}(x_d)$ represent a tensor - * product. The function returns a pair with the value of the interpolation - * as the first component and the gradient in reference coordinates as the - * second component. Note that for compound types (e.g. the `values` field - * begin a Point argument), the components of the gradient are - * sorted as Tensor<1, dim, Tensor<1, spacedim>> with the derivatives - * as the first index; this is a consequence of the generic arguments in the - * function. - * - * @param poly The underlying one-dimensional polynomial basis - * $\{\varphi^{1d}_{i_1}\}$ given as a vector of polynomials. - * - * @param values The expansion coefficients $u_i$ of type `Number` in - * the polynomial interpolation. The coefficients can be simply `double` - * variables but e.g. also Point in case they define arithmetic - * operations with the type `Number2`. - * - * @param p The position in reference coordinates where the interpolation - * should be evaluated. - * - * @param d_linear Flag to specify whether a d-linear (linear in 1d, - * bi-linear in 2d, tri-linear in 3d) interpolation should be made, which - * allows to unroll loops and considerably speed up evaluation. - * - * @param renumber Optional parameter to specify a renumbering in the - * coefficient vector, assuming that `values[renumber[i]]` returns - * the lexicographic (tensor product) entry of the coefficients. If the - * vector is entry, the values are assumed to be sorted lexicographically. + * Interpolates the values and gradients into the points specified in + * @p compute_values_of_array() with help of the precomputed @p shapes. */ template inline std::pair< typename ProductTypeNoPoint::type, Tensor<1, dim, typename ProductTypeNoPoint::type>> - evaluate_tensor_product_value_and_gradient( - const std::vector> &poly, - const std::vector & values, - const Point & p, - const bool d_linear = false, - const std::vector & renumber = {}) + evaluate_tensor_product_value_and_gradient_shapes( + const ArrayView> &shapes, + const int n_shapes, + const std::vector & values, + const std::vector & renumber = {}) { static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented"); using Number3 = typename ProductTypeNoPoint::type; - // use `int` type for this variable and the loops below to inform the - // compiler that the loops below will never overflow, which allows it to - // generate more optimized code for the variable loop bounds in the - // present context - const int n_shapes = poly.size(); AssertDimension(Utilities::pow(n_shapes, dim), values.size()); Assert(renumber.empty() || renumber.size() == values.size(), ExcDimensionMismatch(renumber.size(), values.size())); - // shortcut for linear interpolation to speed up evaluation - if (d_linear) - { - AssertDimension(poly.size(), 2); - for (unsigned int i = 0; i < renumber.size(); ++i) - AssertDimension(renumber[i], i); - - if (dim == 1) - { - Tensor<1, dim, Number3> derivative; - derivative[0] = values[1] - values[0]; - return std::make_pair((1. - p[0]) * values[0] + p[0] * values[1], - derivative); - } - else if (dim == 2) - { - const Number2 x0 = 1. - p[0], x1 = p[0]; - const Number3 tmp0 = x0 * values[0] + x1 * values[1]; - const Number3 tmp1 = x0 * values[2] + x1 * values[3]; - const Number3 mapped = (1. - p[1]) * tmp0 + p[1] * tmp1; - Tensor<1, dim, Number3> derivative; - derivative[0] = (1. - p[1]) * (values[1] - values[0]) + - p[1] * (values[3] - values[2]); - derivative[1] = tmp1 - tmp0; - return std::make_pair(mapped, derivative); - } - else if (dim == 3) - { - const Number2 x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1], - z0 = 1. - p[2], z1 = p[2]; - const Number3 tmp0 = x0 * values[0] + x1 * values[1]; - const Number3 tmp1 = x0 * values[2] + x1 * values[3]; - const Number3 tmpy0 = y0 * tmp0 + y1 * tmp1; - const Number3 tmp2 = x0 * values[4] + x1 * values[5]; - const Number3 tmp3 = x0 * values[6] + x1 * values[7]; - const Number3 tmpy1 = y0 * tmp2 + y1 * tmp3; - const Number3 mapped = z0 * tmpy0 + z1 * tmpy1; - Tensor<1, dim, Number3> derivative; - derivative[2] = tmpy1 - tmpy0; - derivative[1] = z0 * (tmp1 - tmp0) + z1 * (tmp3 - tmp2); - derivative[0] = - z0 * - (y0 * (values[1] - values[0]) + y1 * (values[3] - values[2])) + - z1 * - (y0 * (values[5] - values[4]) + y1 * (values[7] - values[6])); - return std::make_pair(mapped, derivative); - } - } - - AssertIndexRange(n_shapes, 200); - dealii::ndarray shapes; - - // Evaluate 1d polynomials and their derivatives - std::array point; - for (unsigned int d = 0; d < dim; ++d) - point[d] = p[d]; - for (int i = 0; i < n_shapes; ++i) - poly[i].values_of_array(point, 1, &shapes[i][0]); - // Go through the tensor product of shape functions and interpolate // with optimal algorithm std::pair> result = {}; @@ -3199,6 +3127,147 @@ namespace internal + /** + * Specializes @p evaluate_tensor_product_value_and_gradient() for linear + * polynomials which massively reduces the necessary instructions. + */ + template + inline std::pair< + typename ProductTypeNoPoint::type, + Tensor<1, dim, typename ProductTypeNoPoint::type>> + evaluate_tensor_product_value_and_gradient_linear( + const std::vector> &poly, + const std::vector & values, + const Point & p, + const std::vector & renumber = {}) + { + (void)poly; + static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented"); + + using Number3 = typename ProductTypeNoPoint::type; + + AssertDimension(Utilities::pow(poly.size(), dim), values.size()); + Assert(renumber.empty() || renumber.size() == values.size(), + ExcDimensionMismatch(renumber.size(), values.size())); + + AssertDimension(poly.size(), 2); + for (unsigned int i = 0; i < renumber.size(); ++i) + AssertDimension(renumber[i], i); + + if (dim == 1) + { + Tensor<1, dim, Number3> derivative; + derivative[0] = values[1] - values[0]; + return std::make_pair((1. - p[0]) * values[0] + p[0] * values[1], + derivative); + } + else if (dim == 2) + { + const Number2 x0 = 1. - p[0], x1 = p[0]; + const Number3 tmp0 = x0 * values[0] + x1 * values[1]; + const Number3 tmp1 = x0 * values[2] + x1 * values[3]; + const Number3 mapped = (1. - p[1]) * tmp0 + p[1] * tmp1; + Tensor<1, dim, Number3> derivative; + derivative[0] = (1. - p[1]) * (values[1] - values[0]) + + p[1] * (values[3] - values[2]); + derivative[1] = tmp1 - tmp0; + return std::make_pair(mapped, derivative); + } + else if (dim == 3) + { + const Number2 x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1], + z0 = 1. - p[2], z1 = p[2]; + const Number3 tmp0 = x0 * values[0] + x1 * values[1]; + const Number3 tmp1 = x0 * values[2] + x1 * values[3]; + const Number3 tmpy0 = y0 * tmp0 + y1 * tmp1; + const Number3 tmp2 = x0 * values[4] + x1 * values[5]; + const Number3 tmp3 = x0 * values[6] + x1 * values[7]; + const Number3 tmpy1 = y0 * tmp2 + y1 * tmp3; + const Number3 mapped = z0 * tmpy0 + z1 * tmpy1; + Tensor<1, dim, Number3> derivative; + derivative[2] = tmpy1 - tmpy0; + derivative[1] = z0 * (tmp1 - tmp0) + z1 * (tmp3 - tmp2); + derivative[0] = + z0 * (y0 * (values[1] - values[0]) + y1 * (values[3] - values[2])) + + z1 * (y0 * (values[5] - values[4]) + y1 * (values[7] - values[6])); + return std::make_pair(mapped, derivative); + } + + // work around a compile error: missing return statement + return std::make_pair(Number3(), Tensor<1, dim, Number3>()); + } + + + + /** + * Compute the polynomial interpolation of a tensor product shape function + * $\varphi_i$ given a vector of coefficients $u_i$ in the form + * $u_h(\mathbf{x}) = \sum_{i=1}^{k^d} \varphi_i(\mathbf{x}) u_i$. The shape + * functions $\varphi_i(\mathbf{x}) = + * \prod_{d=1}^{\text{dim}}\varphi_{i_d}^\text{1d}(x_d)$ represent a tensor + * product. The function returns a pair with the value of the interpolation + * as the first component and the gradient in reference coordinates as the + * second component. Note that for compound types (e.g. the `values` field + * begin a Point argument), the components of the gradient are + * sorted as Tensor<1, dim, Tensor<1, spacedim>> with the derivatives + * as the first index; this is a consequence of the generic arguments in the + * function. + * + * @param poly The underlying one-dimensional polynomial basis + * $\{\varphi^{1d}_{i_1}\}$ given as a vector of polynomials. + * + * @param values The expansion coefficients $u_i$ of type `Number` in + * the polynomial interpolation. The coefficients can be simply `double` + * variables but e.g. also Point in case they define arithmetic + * operations with the type `Number2`. + * + * @param p The position in reference coordinates where the interpolation + * should be evaluated. + * + * @param d_linear Flag to specify whether a d-linear (linear in 1d, + * bi-linear in 2d, tri-linear in 3d) interpolation should be made, which + * allows to unroll loops and considerably speed up evaluation. + * + * @param renumber Optional parameter to specify a renumbering in the + * coefficient vector, assuming that `values[renumber[i]]` returns + * the lexicographic (tensor product) entry of the coefficients. If the + * vector is entry, the values are assumed to be sorted lexicographically. + */ + template + inline std::pair< + typename ProductTypeNoPoint::type, + Tensor<1, dim, typename ProductTypeNoPoint::type>> + evaluate_tensor_product_value_and_gradient( + const std::vector> &poly, + const std::vector & values, + const Point & p, + const bool d_linear = false, + const std::vector & renumber = {}) + { + if (d_linear) + { + return evaluate_tensor_product_value_and_gradient_linear(poly, + values, + p, + renumber); + } + else + { + std::array, 200> shapes; + + auto view = make_array_view(shapes); + + compute_values_of_array(view, poly, p); + + return evaluate_tensor_product_value_and_gradient_shapes( + view, poly.size(), values, renumber); + } + } + + + template SymmetricTensor<2, dim, typename ProductTypeNoPoint::type> evaluate_tensor_product_hessian( @@ -3305,12 +3374,13 @@ namespace internal */ template inline void - do_apply_test_functions_xy(AlignedVector & values, - const std::vector &renumber, - const dealii::ndarray &shapes, - const std::array &test_grads_value, - const int n_shapes_runtime, - int & i) + do_apply_test_functions_xy( + AlignedVector & values, + const std::vector & renumber, + const ArrayView> &shapes, + const std::array & test_grads_value, + const int n_shapes_runtime, + int & i) { const int n_shapes = length > 0 ? length : n_shapes_runtime; for (int i1 = 0; i1 < (dim > 1 ? n_shapes : 1); ++i1) @@ -3344,32 +3414,21 @@ namespace internal */ template inline void - integrate_add_tensor_product_value_and_gradient( - const std::vector> &poly, - const Number2 & value, - const Tensor<1, dim, Number2> & gradient, - const Point & p, - AlignedVector & values, - const std::vector & renumber = {}) + integrate_add_tensor_product_value_and_gradient_shapes( + const ArrayView> &shapes, + const int n_shapes, + const Number2 & value, + const Tensor<1, dim, Number2> & gradient, + AlignedVector & values, + const std::vector & renumber = {}) { static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented"); // as in evaluate, use `int` type to produce better code in this context - const int n_shapes = poly.size(); AssertDimension(Utilities::pow(n_shapes, dim), values.size()); Assert(renumber.empty() || renumber.size() == values.size(), ExcDimensionMismatch(renumber.size(), values.size())); - AssertIndexRange(n_shapes, 200); - dealii::ndarray shapes; - - // Evaluate 1d polynomials and their derivatives - std::array point; - for (unsigned int d = 0; d < dim; ++d) - point[d] = p[d]; - for (int i = 0; i < n_shapes; ++i) - poly[i].values_of_array(point, 1, &shapes[i][0]); - // Implement the transpose of the function above std::array test_grads_value; for (int i2 = 0, i = 0; i2 < (dim > 2 ? n_shapes : 1); ++i2) @@ -3411,6 +3470,86 @@ namespace internal + /** + * Specializes @p evaluate_tensor_product_value_and_gradient() for linear + * polynomials which massively reduces the necessary instructions. + */ + template + inline void + integrate_add_tensor_product_value_and_gradient_linear( + const std::vector> &poly, + const Number2 & value, + const Tensor<1, dim, Number2> & gradient, + AlignedVector & values, + const Point & p, + const std::vector & renumber = {}) + { + (void)poly; + static_assert(dim >= 1 && dim <= 3, "Only dim=1,2,3 implemented"); + + AssertDimension(Utilities::pow(poly.size(), dim), values.size()); + Assert(renumber.empty() || renumber.size() == values.size(), + ExcDimensionMismatch(renumber.size(), values.size())); + + AssertDimension(poly.size(), 2); + for (unsigned int i = 0; i < renumber.size(); ++i) + AssertDimension(renumber[i], i); + + if (dim == 1) + { + const auto x0 = 1. - p[0], x1 = p[0]; + + values[0] = value * x0 - gradient[0]; + values[1] = value * x1 + gradient[0]; + } + else if (dim == 2) + { + const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1]; + + const auto test_value_y0 = value * y0 - gradient[1]; + const auto test_grad_xy0 = gradient[0] * y0; + const auto test_value_y1 = value * y1 + gradient[1]; + const auto test_grad_xy1 = gradient[0] * y1; + + values[0] += x0 * test_value_y0 - test_grad_xy0; + values[1] += x1 * test_value_y0 + test_grad_xy0; + values[2] += x0 * test_value_y1 - test_grad_xy1; + values[3] += x1 * test_value_y1 + test_grad_xy1; + } + else if (dim == 3) + { + const auto x0 = 1. - p[0], x1 = p[0], y0 = 1. - p[1], y1 = p[1], + z0 = 1. - p[2], z1 = p[2]; + + const auto test_value_z0 = value * z0 - gradient[2]; + const auto test_grad_x0 = gradient[0] * z0; + const auto test_grad_y0 = gradient[1] * z0; + const auto test_value_z1 = value * z1 + gradient[2]; + const auto test_grad_x1 = gradient[0] * z1; + const auto test_grad_y1 = gradient[1] * z1; + + const auto test_value_y00 = test_value_z0 * y0 - test_grad_y0; + const auto test_grad_xy00 = test_grad_x0 * y0; + const auto test_value_y01 = test_value_z0 * y1 + test_grad_y0; + const auto test_grad_xy01 = test_grad_x0 * y1; + const auto test_value_y10 = test_value_z1 * y0 - test_grad_y1; + const auto test_grad_xy10 = test_grad_x1 * y0; + const auto test_value_y11 = test_value_z1 * y1 + test_grad_y1; + const auto test_grad_xy11 = test_grad_x1 * y1; + + values[0] += x0 * test_value_y00 - test_grad_xy00; + values[1] += x1 * test_value_y00 + test_grad_xy00; + values[2] += x0 * test_value_y01 - test_grad_xy01; + values[3] += x1 * test_value_y01 + test_grad_xy01; + values[4] += x0 * test_value_y10 - test_grad_xy10; + values[5] += x1 * test_value_y10 + test_grad_xy10; + values[6] += x0 * test_value_y11 - test_grad_xy11; + values[7] += x1 * test_value_y11 + test_grad_xy11; + } + } + + + template inline void weight_fe_q_dofs_by_entity(const Number * weights, diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h index 21347f0530..63ad40f87d 100644 --- a/include/deal.II/non_matching/mapping_info.h +++ b/include/deal.II/non_matching/mapping_info.h @@ -616,7 +616,7 @@ namespace NonMatching template - const ArrayView> + inline const ArrayView> MappingInfo::get_unit_points( const unsigned int cell_index, const unsigned int face_number) const @@ -703,7 +703,7 @@ namespace NonMatching template - const typename MappingInfo::MappingData & + inline const typename MappingInfo::MappingData & MappingInfo::get_mapping_data( const unsigned int cell_index, const unsigned int face_number) const -- 2.39.5