From: Martin Kronbichler Date: Fri, 21 Apr 2023 16:54:19 +0000 (+0200) Subject: Remove tabulated shape functions from MappingQ X-Git-Tag: v9.5.0-rc1~297^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=939973e9d3571b4ae5afdce0b2a9ae6038bdef69;p=dealii.git Remove tabulated shape functions from MappingQ --- diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h index 8b5e7affff..b3d385f635 100644 --- a/include/deal.II/fe/fe_values.h +++ b/include/deal.II/fe/fe_values.h @@ -3954,7 +3954,6 @@ protected: internal::FEValuesImplementation::MappingRelatedData mapping_output; - /** * A pointer to the finite element object associated with this FEValues * object. diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 28f34de1bf..0643b17e24 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -334,88 +334,6 @@ public: const Quadrature &quadrature, const unsigned int n_original_q_points); - /** - * Compute the values and/or derivatives of the shape functions used for - * the mapping. - * - * Which values, derivatives, or higher order derivatives are computed is - * determined by which of the member arrays have nonzero sizes. They are - * typically set to their appropriate sizes by the initialize() and - * initialize_face() functions, which indeed call this function - * internally. However, it is possible (and at times useful) to do the - * resizing by hand and then call this function directly. An example is in - * a Newton iteration where we update the location of a quadrature point - * (e.g., in MappingQ::transform_real_to_uni_cell()) and need to re- - * compute the mapping and its derivatives at this location, but have - * already sized all internal arrays correctly. - */ - void - compute_shape_function_values(const std::vector> &unit_points); - - /** - * Shape function at quadrature point. Shape functions are in tensor - * product order, so vertices must be reordered to obtain transformation. - */ - const double & - shape(const unsigned int qpoint, const unsigned int shape_nr) const; - - /** - * Shape function at quadrature point. See above. - */ - double & - shape(const unsigned int qpoint, const unsigned int shape_nr); - - /** - * Gradient of shape function in quadrature point. See above. - */ - const Tensor<1, dim> & - derivative(const unsigned int qpoint, const unsigned int shape_nr) const; - - /** - * Gradient of shape function in quadrature point. See above. - */ - Tensor<1, dim> & - derivative(const unsigned int qpoint, const unsigned int shape_nr); - - /** - * Second derivative of shape function in quadrature point. See above. - */ - const Tensor<2, dim> & - second_derivative(const unsigned int qpoint, - const unsigned int shape_nr) const; - - /** - * Second derivative of shape function in quadrature point. See above. - */ - Tensor<2, dim> & - second_derivative(const unsigned int qpoint, const unsigned int shape_nr); - - /** - * third derivative of shape function in quadrature point. See above. - */ - const Tensor<3, dim> & - third_derivative(const unsigned int qpoint, - const unsigned int shape_nr) const; - - /** - * third derivative of shape function in quadrature point. See above. - */ - Tensor<3, dim> & - third_derivative(const unsigned int qpoint, const unsigned int shape_nr); - - /** - * fourth derivative of shape function in quadrature point. See above. - */ - const Tensor<4, dim> & - fourth_derivative(const unsigned int qpoint, - const unsigned int shape_nr) const; - - /** - * fourth derivative of shape function in quadrature point. See above. - */ - Tensor<4, dim> & - fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr); - /** * Return an estimate (in bytes) for the memory consumption of this object. */ @@ -427,45 +345,7 @@ public: * possible orientations. Can be accessed with the correct offset provided * via QProjector::DataSetDescriptor. Not needed/used for cells. */ - AlignedVector> quadrature_points; - - /** - * Values of shape functions. Access by function @p shape. - * - * Computed once. - */ - AlignedVector shape_values; - - /** - * Values of shape function derivatives. Access by function @p derivative. - * - * Computed once. - */ - AlignedVector> shape_derivatives; - - /** - * Values of shape function second derivatives. Access by function @p - * second_derivative. - * - * Computed once. - */ - AlignedVector> shape_second_derivatives; - - /** - * Values of shape function third derivatives. Access by function @p - * second_derivative. - * - * Computed once. - */ - AlignedVector> shape_third_derivatives; - - /** - * Values of shape function fourth derivatives. Access by function @p - * second_derivative. - * - * Computed once. - */ - AlignedVector> shape_fourth_derivatives; + std::vector> quadrature_points; /** * Unit tangential vectors. Used for the computation of boundary forms and @@ -548,26 +428,6 @@ public: */ bool tensor_product_quadrature; - /** - * Tensors of covariant transformation at each of the quadrature points. - * The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} * - * Jacobian, is the first fundamental form of the map; if dim=spacedim - * then it reduces to the transpose of the inverse of the Jacobian matrix, - * which itself is stored in the @p contravariant field of this structure. - * - * Computed on each cell. - */ - mutable AlignedVector> covariant; - - /** - * Tensors of contravariant transformation at each of the quadrature - * points. The contravariant matrix is the Jacobian of the transformation, - * i.e. $J_{ij}=dx_i/d\hat x_j$. - * - * Computed on each cell. - */ - mutable AlignedVector> contravariant; - /** * Auxiliary vectors for internal use. */ @@ -590,6 +450,14 @@ public: * #update_volume_elements. */ mutable AlignedVector volume_elements; + + /** + * Pointer to the mapping output data that holds most of the arrays, + * including the Jacobians representing the covariant and contravariant + * transformations. + */ + mutable internal::FEValuesImplementation::MappingRelatedData + *output_data; }; protected: @@ -840,123 +708,6 @@ using MappingQGeneric = MappingQ; #ifndef DOXYGEN -template -inline const double & -MappingQ::InternalData::shape(const unsigned int qpoint, - const unsigned int shape_nr) const -{ - AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size()); - return shape_values[qpoint * n_shape_functions + shape_nr]; -} - - - -template -inline double & -MappingQ::InternalData::shape(const unsigned int qpoint, - const unsigned int shape_nr) -{ - AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size()); - return shape_values[qpoint * n_shape_functions + shape_nr]; -} - - -template -inline const Tensor<1, dim> & -MappingQ::InternalData::derivative( - const unsigned int qpoint, - const unsigned int shape_nr) const -{ - AssertIndexRange(qpoint * n_shape_functions + shape_nr, - shape_derivatives.size()); - return shape_derivatives[qpoint * n_shape_functions + shape_nr]; -} - - - -template -inline Tensor<1, dim> & -MappingQ::InternalData::derivative(const unsigned int qpoint, - const unsigned int shape_nr) -{ - AssertIndexRange(qpoint * n_shape_functions + shape_nr, - shape_derivatives.size()); - return shape_derivatives[qpoint * n_shape_functions + shape_nr]; -} - - -template -inline const Tensor<2, dim> & -MappingQ::InternalData::second_derivative( - const unsigned int qpoint, - const unsigned int shape_nr) const -{ - AssertIndexRange(qpoint * n_shape_functions + shape_nr, - shape_second_derivatives.size()); - return shape_second_derivatives[qpoint * n_shape_functions + shape_nr]; -} - - -template -inline Tensor<2, dim> & -MappingQ::InternalData::second_derivative( - const unsigned int qpoint, - const unsigned int shape_nr) -{ - AssertIndexRange(qpoint * n_shape_functions + shape_nr, - shape_second_derivatives.size()); - return shape_second_derivatives[qpoint * n_shape_functions + shape_nr]; -} - -template -inline const Tensor<3, dim> & -MappingQ::InternalData::third_derivative( - const unsigned int qpoint, - const unsigned int shape_nr) const -{ - AssertIndexRange(qpoint * n_shape_functions + shape_nr, - shape_third_derivatives.size()); - return shape_third_derivatives[qpoint * n_shape_functions + shape_nr]; -} - - -template -inline Tensor<3, dim> & -MappingQ::InternalData::third_derivative( - const unsigned int qpoint, - const unsigned int shape_nr) -{ - AssertIndexRange(qpoint * n_shape_functions + shape_nr, - shape_third_derivatives.size()); - return shape_third_derivatives[qpoint * n_shape_functions + shape_nr]; -} - - -template -inline const Tensor<4, dim> & -MappingQ::InternalData::fourth_derivative( - const unsigned int qpoint, - const unsigned int shape_nr) const -{ - AssertIndexRange(qpoint * n_shape_functions + shape_nr, - shape_fourth_derivatives.size()); - return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr]; -} - - -template -inline Tensor<4, dim> & -MappingQ::InternalData::fourth_derivative( - const unsigned int qpoint, - const unsigned int shape_nr) -{ - AssertIndexRange(qpoint * n_shape_functions + shape_nr, - shape_fourth_derivatives.size()); - return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr]; -} - - - template inline bool MappingQ::preserves_vertex_locations() const diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index e144e33096..e39202634c 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -733,111 +733,53 @@ namespace internal template inline Point do_transform_real_to_unit_cell_internal_codim1( - const typename dealii::Triangulation::cell_iterator &cell, - const Point & p, - const Point & initial_p_unit, - typename dealii::MappingQ::InternalData &mdata) + const Point & p, + const Point & initial_p_unit, + const std::vector> & points, + const std::vector> &polynomials_1d, + const std::vector & renumber) { const unsigned int spacedim = dim + 1; - const unsigned int n_shapes = mdata.shape_values.size(); - (void)n_shapes; - Assert(n_shapes != 0, ExcInternalError()); - Assert(mdata.shape_derivatives.size() == n_shapes, ExcInternalError()); - Assert(mdata.shape_second_derivatives.size() == n_shapes, - ExcInternalError()); - - std::vector> &points = mdata.mapping_support_points; - Assert(points.size() == n_shapes, ExcInternalError()); - - Point p_minus_F; + AssertDimension(points.size(), + Utilities::pow(polynomials_1d.size(), dim)); - Tensor<1, spacedim> DF[dim]; - Tensor<1, spacedim> D2F[dim][dim]; + Point p_unit = initial_p_unit; - Point p_unit = initial_p_unit; - Point f; - Tensor<2, dim> df; + const double eps = 1.e-12; + const unsigned int loop_limit = 10; - // Evaluate first and second derivatives - mdata.compute_shape_function_values(std::vector>(1, p_unit)); + unsigned int loop = 0; + double f_weighted_norm_square = 1.; - for (unsigned int k = 0; k < mdata.n_shape_functions; ++k) + while (f_weighted_norm_square > eps * eps && loop++ < loop_limit) { - const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k); - const Tensor<2, dim> & hessian_k = mdata.second_derivative(0, k); - const Point &point_k = points[k]; - + const auto p_real = + internal::evaluate_tensor_product_value_and_gradient( + polynomials_1d, + points, + p_unit, + polynomials_1d.size() == 2, + renumber); + Tensor<1, spacedim> p_minus_F = p - p_real.first; + const DerivativeForm<1, spacedim, dim> DF = p_real.second; + + const auto hessian = internal::evaluate_tensor_product_hessian( + polynomials_1d, points, p_unit, renumber); + Point f; + Tensor<2, dim> df; for (unsigned int j = 0; j < dim; ++j) { - DF[j] += grad_phi_k[j] * point_k; + f[j] = DF[j] * p_minus_F; for (unsigned int l = 0; l < dim; ++l) - D2F[j][l] += hessian_k[j][l] * point_k; + df[j][l] = -DF[j] * DF[l] + hessian[j][l] * p_minus_F; } - } - - p_minus_F = p; - p_minus_F -= compute_mapped_location_of_point(mdata); - - - for (unsigned int j = 0; j < dim; ++j) - f[j] = DF[j] * p_minus_F; - - for (unsigned int j = 0; j < dim; ++j) - { - f[j] = DF[j] * p_minus_F; - for (unsigned int l = 0; l < dim; ++l) - df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F; - } - - const double eps = 1.e-12 * cell->diameter(); - const unsigned int loop_limit = 10; - - unsigned int loop = 0; - - while (f.norm() > eps && loop++ < loop_limit) - { // Solve [df(x)]d=f(x) const Tensor<1, dim> d = invert(df) * static_cast &>(f); + f_weighted_norm_square = d.norm_square(); p_unit -= d; - - for (unsigned int j = 0; j < dim; ++j) - { - DF[j].clear(); - for (unsigned int l = 0; l < dim; ++l) - D2F[j][l].clear(); - } - - mdata.compute_shape_function_values( - std::vector>(1, p_unit)); - - for (unsigned int k = 0; k < mdata.n_shape_functions; ++k) - { - const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k); - const Tensor<2, dim> & hessian_k = mdata.second_derivative(0, k); - const Point &point_k = points[k]; - - for (unsigned int j = 0; j < dim; ++j) - { - DF[j] += grad_phi_k[j] * point_k; - for (unsigned int l = 0; l < dim; ++l) - D2F[j][l] += hessian_k[j][l] * point_k; - } - } - - // TODO: implement a line search here in much the same way as for - // the corresponding function above that does so for dim==spacedim - p_minus_F = p; - p_minus_F -= compute_mapped_location_of_point(mdata); - - for (unsigned int j = 0; j < dim; ++j) - { - f[j] = DF[j] * p_minus_F; - for (unsigned int l = 0; l < dim; ++l) - df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F; - } } @@ -1131,6 +1073,15 @@ namespace internal constexpr unsigned int n_comp = 1 + (spacedim - 1) / n_lanes; constexpr unsigned int n_hessians = (dim * (dim + 1)) / 2; + // Since MappingQ::InternalData does not have separate arrays for the + // covariant and contravariant transformations, but uses the arrays in + // the `MappingRelatedData`, it can happen that vectors do not have the + // right size + if (update_flags & update_contravariant_transformation) + jacobians.resize(n_q_points); + if (update_flags & update_covariant_transformation) + inverse_jacobians.resize(n_q_points); + EvaluationFlags::EvaluationFlags evaluation_flag = (update_flags & update_quadrature_points ? EvaluationFlags::values : EvaluationFlags::nothing) | @@ -1151,9 +1102,6 @@ namespace internal Assert(!(evaluation_flag & EvaluationFlags::gradients) || data.n_shape_functions > 0, ExcInternalError()); - Assert(!(evaluation_flag & EvaluationFlags::gradients) || - n_q_points == data.contravariant.size(), - ExcDimensionMismatch(n_q_points, data.contravariant.size())); Assert(!(evaluation_flag & EvaluationFlags::hessians) || n_q_points == jacobian_grads.size(), ExcDimensionMismatch(n_q_points, jacobian_grads.size())); @@ -1214,9 +1162,6 @@ namespace internal if (evaluation_flag & EvaluationFlags::gradients) { - std::fill(data.contravariant.begin(), - data.contravariant.end(), - DerivativeForm<1, dim, spacedim>()); // We need to reinterpret the data after evaluate has been applied. for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp) for (unsigned int point = 0; point < n_q_points; ++point) @@ -1229,43 +1174,28 @@ namespace internal const unsigned int total_number = point * dim + j; const unsigned int new_comp = total_number / n_q_points; const unsigned int new_point = total_number % n_q_points; - data.contravariant[new_point][out_comp * n_lanes + in_comp] - [new_comp] = + jacobians[new_point][out_comp * n_lanes + + in_comp][new_comp] = eval.begin_gradients()[(out_comp * n_q_points + point) * dim + j][in_comp]; } } - if (update_flags & update_covariant_transformation) - if (cell_similarity != CellSimilarity::translation) - for (unsigned int point = 0; point < n_q_points; ++point) - data.covariant[point] = - (data.contravariant[point]).covariant_form(); if (update_flags & update_volume_elements) if (cell_similarity != CellSimilarity::translation) for (unsigned int point = 0; point < n_q_points; ++point) - data.volume_elements[point] = - data.contravariant[point].determinant(); + data.volume_elements[point] = jacobians[point].determinant(); // copy values from InternalData to vector given by reference - if (update_flags & update_jacobians) + if (update_flags & update_covariant_transformation) { - const unsigned int n_q_points = data.contravariant.size(); AssertDimension(jacobians.size(), n_q_points); - if (cell_similarity != CellSimilarity::translation) - for (unsigned int point = 0; point < n_q_points; ++point) - jacobians[point] = data.contravariant[point]; - } - - // copy values from InternalData to vector given by reference - if (update_flags & update_inverse_jacobians) - { - const unsigned int n_q_points = data.contravariant.size(); AssertDimension(inverse_jacobians.size(), n_q_points); if (cell_similarity != CellSimilarity::translation) for (unsigned int point = 0; point < n_q_points; ++point) - inverse_jacobians[point] = data.covariant[point].transpose(); + inverse_jacobians[point] = + jacobians[point].covariant_form().transpose(); } if (evaluation_flag & EvaluationFlags::hessians) @@ -1309,6 +1239,40 @@ namespace internal + template + inline DEAL_II_ALWAYS_INLINE void + store_vectorized_tensor( + const unsigned int n_points, + const unsigned int cur_index, + const DerivativeForm<1, dim, spacedim, VectorizedArray> + & derivative, + std::vector> &result_array) + { + AssertDimension(result_array.size(), n_points); + constexpr unsigned int n_lanes = VectorizedArray::size(); + if (cur_index + n_lanes <= n_points) + { + std::array indices; + for (unsigned int j = 0; j < n_lanes; ++j) + indices[j] = j * dim * spacedim; + const unsigned int even = (dim * spacedim) / 4 * 4; + double * result_ptr = &result_array[cur_index][0][0]; + const VectorizedArray *derivative_ptr = &derivative[0][0]; + vectorized_transpose_and_store( + false, even, derivative_ptr, indices.data(), result_ptr); + for (unsigned int d = even; d < dim * spacedim; ++d) + for (unsigned int j = 0; j < n_lanes; ++j) + result_ptr[j * dim * spacedim + d] = derivative_ptr[d][j]; + } + else + for (unsigned int j = 0; j < n_lanes && cur_index + j < n_points; ++j) + for (unsigned int d = 0; d < spacedim; ++d) + for (unsigned int e = 0; e < dim; ++e) + result_array[cur_index + j][d][e] = derivative[d][e][j]; + } + + + template inline void maybe_update_q_points_Jacobians_generic( @@ -1316,7 +1280,6 @@ namespace internal const typename dealii::MappingQ::InternalData &data, const ArrayView> & unit_points, const std::vector> &polynomials_1d, - const unsigned int polynomial_degree, const std::vector &renumber_lexicographic_to_hierarchic, std::vector> & quadrature_points, std::vector> &jacobians, @@ -1329,6 +1292,15 @@ namespace internal const unsigned int n_points = unit_points.size(); const unsigned int n_lanes = VectorizedArray::size(); + // Since MappingQ::InternalData does not have separate arrays for the + // covariant and contravariant transformations, but uses the arrays in + // the `MappingRelatedData`, it can happen that vectors do not have the + // right size + if (update_flags & update_contravariant_transformation) + jacobians.resize(n_points); + if (update_flags & update_covariant_transformation) + inverse_jacobians.resize(n_points); + // Use the more heavy VectorizedArray code path if there is more than // one point left to compute for (unsigned int i = 0; i < n_points; i += n_lanes) @@ -1348,8 +1320,13 @@ namespace internal polynomials_1d, support_points, p_vec, - polynomial_degree == 1, + polynomials_1d.size() == 2, renumber_lexicographic_to_hierarchic); + DerivativeForm<1, dim, spacedim, VectorizedArray> + derivative; + for (unsigned int d = 0; d < spacedim; ++d) + for (unsigned int e = 0; e < dim; ++e) + derivative[d][e] = result.second[e][d]; if (update_flags & update_quadrature_points) for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) @@ -1360,28 +1337,23 @@ namespace internal continue; if (update_flags & update_contravariant_transformation) - for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int e = 0; e < dim; ++e) - data.contravariant[i + j][d][e] = result.second[e][d][j]; + store_vectorized_tensor(n_points, i, derivative, jacobians); if (update_flags & update_volume_elements) - for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) - data.volume_elements[i + j] = - data.contravariant[i + j].determinant(); - - if (update_flags & update_jacobians) - for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) - jacobians[i + j] = data.contravariant[i + j]; + { + const auto determinant = derivative.determinant(); + for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) + data.volume_elements[i + j] = determinant[j]; + } if (update_flags & update_covariant_transformation) - for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) - data.covariant[i + j] = - data.contravariant[i + j].covariant_form(); - - if (update_flags & update_inverse_jacobians) - for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) - inverse_jacobians[i + j] = data.covariant[i + j].transpose(); + { + const auto covariant = derivative.covariant_form(); + store_vectorized_tensor(n_points, + i, + covariant.transpose(), + inverse_jacobians); + } } else { @@ -1390,7 +1362,7 @@ namespace internal polynomials_1d, support_points, unit_points[i], - polynomial_degree == 1, + polynomials_1d.size() == 2, renumber_lexicographic_to_hierarchic); if (update_flags & update_quadrature_points) @@ -1399,23 +1371,16 @@ namespace internal if (cell_similarity == CellSimilarity::translation) continue; - if (update_flags & update_contravariant_transformation) - { - DerivativeForm<1, spacedim, dim> jac_transposed = result.second; - data.contravariant[i] = jac_transposed.transpose(); - } - - if (update_flags & update_volume_elements) - data.volume_elements[i] = data.contravariant[i].determinant(); + if (dim == spacedim && update_flags & update_volume_elements) + data.volume_elements[i] = + DerivativeForm<1, spacedim, dim>(result.second).determinant(); - if (update_flags & update_jacobians) - jacobians[i] = data.contravariant[i]; + if (update_flags & update_contravariant_transformation) + jacobians[i] = + DerivativeForm<1, spacedim, dim>(result.second).transpose(); if (update_flags & update_covariant_transformation) - data.covariant[i] = data.contravariant[i].covariant_form(); - - if (update_flags & update_inverse_jacobians) - inverse_jacobians[i] = data.covariant[i].transpose(); + inverse_jacobians[i] = jacobians[i].covariant_form().transpose(); } } @@ -1430,38 +1395,33 @@ namespace internal template inline void maybe_update_jacobian_grads( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, + const CellSimilarity::Similarity cell_similarity, const typename dealii::MappingQ::InternalData &data, + const ArrayView> & unit_points, + const std::vector> &polynomials_1d, + const std::vector &renumber_lexicographic_to_hierarchic, std::vector> &jacobian_grads) { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_grads) + if (data.update_each & update_jacobian_grads) { + const std::vector> &support_points = + data.mapping_support_points; const unsigned int n_q_points = jacobian_grads.size(); if (cell_similarity != CellSimilarity::translation) for (unsigned int point = 0; point < n_q_points; ++point) { - const Tensor<2, dim> *second = - &data.second_derivative(point + data_set, 0); - double result[spacedim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - result[i][j][l] = - (second[0][j][l] * data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - result[i][j][l] += - (second[k][j][l] * data.mapping_support_points[k][i]); + const SymmetricTensor<2, dim, Tensor<1, spacedim>> second = + internal::evaluate_tensor_product_hessian( + polynomials_1d, + support_points, + unit_points[point], + renumber_lexicographic_to_hierarchic); for (unsigned int i = 0; i < spacedim; ++i) for (unsigned int j = 0; j < dim; ++j) for (unsigned int l = 0; l < dim; ++l) - jacobian_grads[point][i][j][l] = result[i][j][l]; + jacobian_grads[point][i][j][l] = second[j][l][i]; } } } @@ -1477,14 +1437,17 @@ namespace internal template inline void maybe_update_jacobian_pushed_forward_grads( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, + const CellSimilarity::Similarity cell_similarity, const typename dealii::MappingQ::InternalData &data, + const ArrayView> & unit_points, + const std::vector> &polynomials_1d, + const std::vector & renumber_lexicographic_to_hierarchic, std::vector> &jacobian_pushed_forward_grads) { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_grads) + if (data.update_each & update_jacobian_pushed_forward_grads) { + const std::vector> &support_points = + data.mapping_support_points; const unsigned int n_q_points = jacobian_pushed_forward_grads.size(); if (cell_similarity != CellSimilarity::translation) @@ -1492,33 +1455,25 @@ namespace internal double tmp[spacedim][spacedim][spacedim]; for (unsigned int point = 0; point < n_q_points; ++point) { - const Tensor<2, dim> *second = - &data.second_derivative(point + data_set, 0); - double result[spacedim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - result[i][j][l] = - (second[0][j][l] * data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - result[i][j][l] += - (second[k][j][l] * - data.mapping_support_points[k][i]); + const SymmetricTensor<2, dim, Tensor<1, spacedim>> second = + internal::evaluate_tensor_product_hessian( + polynomials_1d, + support_points, + unit_points[point], + renumber_lexicographic_to_hierarchic); + const DerivativeForm<1, dim, spacedim> covariant = + data.output_data->inverse_jacobians[point].transpose(); // first push forward the j-components for (unsigned int i = 0; i < spacedim; ++i) for (unsigned int j = 0; j < spacedim; ++j) for (unsigned int l = 0; l < dim; ++l) { - tmp[i][j][l] = - result[i][0][l] * data.covariant[point][j][0]; + tmp[i][j][l] = second[0][l][i] * covariant[j][0]; for (unsigned int jr = 1; jr < dim; ++jr) { tmp[i][j][l] += - result[i][jr][l] * data.covariant[point][j][jr]; + second[jr][l][i] * covariant[j][jr]; } } @@ -1528,11 +1483,11 @@ namespace internal for (unsigned int l = 0; l < spacedim; ++l) { jacobian_pushed_forward_grads[point][i][j][l] = - tmp[i][j][0] * data.covariant[point][l][0]; + tmp[i][j][0] * covariant[l][0]; for (unsigned int lr = 1; lr < dim; ++lr) { jacobian_pushed_forward_grads[point][i][j][l] += - tmp[i][j][lr] * data.covariant[point][l][lr]; + tmp[i][j][lr] * covariant[l][lr]; } } } @@ -1542,6 +1497,52 @@ namespace internal + template + inline DerivativeForm<3, dim, spacedim> + expand_3rd_derivatives( + const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed) + { + Assert(dim >= 1 && dim <= 3, ExcNotImplemented()); + DerivativeForm<3, dim, spacedim> result; + for (unsigned int i = 0; i < spacedim; ++i) + { + if (dim == 1) + result[i][0][0][0] = compressed[0][i]; + else if (dim >= 2) + { + for (unsigned int d = 0; d < 2; ++d) + for (unsigned int e = 0; e < 2; ++e) + for (unsigned int f = 0; f < 2; ++f) + result[i][d][e][f] = compressed[d + e + f][i]; + if (dim == 3) + { + AssertDimension(length_tensor, 10); + + // the derivatives are ordered in rows by increasing z + // derivative, and in each row we have x^{(n-j)} y^{(j)} as + // j walks through the columns + for (unsigned int d = 0; d < 2; ++d) + for (unsigned int e = 0; e < 2; ++e) + { + result[i][d][e][2] = compressed[4 + d + e][i]; + result[i][d][2][e] = compressed[4 + d + e][i]; + result[i][2][d][e] = compressed[4 + d + e][i]; + } + for (unsigned int d = 0; d < 2; ++d) + { + result[i][d][2][2] = compressed[7 + d][i]; + result[i][2][d][2] = compressed[7 + d][i]; + result[i][2][2][d] = compressed[7 + d][i]; + } + result[i][2][2][2] = compressed[9][i]; + } + } + } + return result; + } + + + /** * Update the third derivatives of the transformation from unit to real * cell, the Jacobian hessians. @@ -1551,45 +1552,29 @@ namespace internal template inline void maybe_update_jacobian_2nd_derivatives( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, + const CellSimilarity::Similarity cell_similarity, const typename dealii::MappingQ::InternalData &data, + const ArrayView> & unit_points, + const std::vector> &polynomials_1d, + const std::vector &renumber_lexicographic_to_hierarchic, std::vector> &jacobian_2nd_derivatives) { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_2nd_derivatives) + if (data.update_each & update_jacobian_2nd_derivatives) { + const std::vector> &support_points = + data.mapping_support_points; const unsigned int n_q_points = jacobian_2nd_derivatives.size(); if (cell_similarity != CellSimilarity::translation) { for (unsigned int point = 0; point < n_q_points; ++point) { - const Tensor<3, dim> *third = - &data.third_derivative(point + data_set, 0); - double result[spacedim][dim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - result[i][j][l][m] = - (third[0][j][l][m] * - data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - result[i][j][l][m] += - (third[k][j][l][m] * - data.mapping_support_points[k][i]); - - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - jacobian_2nd_derivatives[point][i][j][l][m] = - result[i][j][l][m]; + jacobian_2nd_derivatives[point] = expand_3rd_derivatives( + internal::evaluate_tensor_product_higher_derivatives<3>( + polynomials_1d, + support_points, + unit_points[point], + renumber_lexicographic_to_hierarchic)); } } } @@ -1607,40 +1592,35 @@ namespace internal template inline void maybe_update_jacobian_pushed_forward_2nd_derivatives( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, + const CellSimilarity::Similarity cell_similarity, const typename dealii::MappingQ::InternalData &data, + const ArrayView> & unit_points, + const std::vector> &polynomials_1d, + const std::vector & renumber_lexicographic_to_hierarchic, std::vector> &jacobian_pushed_forward_2nd_derivatives) { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) + if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives) { + const std::vector> &support_points = + data.mapping_support_points; const unsigned int n_q_points = jacobian_pushed_forward_2nd_derivatives.size(); if (cell_similarity != CellSimilarity::translation) { - double tmp[spacedim][spacedim][spacedim][spacedim]; + dealii::ndarray tmp; + dealii::ndarray tmp2; for (unsigned int point = 0; point < n_q_points; ++point) { - const Tensor<3, dim> *third = - &data.third_derivative(point + data_set, 0); - double result[spacedim][dim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - result[i][j][l][m] = - (third[0][j][l][m] * - data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - result[i][j][l][m] += - (third[k][j][l][m] * - data.mapping_support_points[k][i]); + const DerivativeForm<3, dim, spacedim> third = + expand_3rd_derivatives( + internal::evaluate_tensor_product_higher_derivatives<3>( + polynomials_1d, + support_points, + unit_points[point], + renumber_lexicographic_to_hierarchic)); + const DerivativeForm<1, dim, spacedim> covariant = + data.output_data->inverse_jacobians[point].transpose(); // push forward the j-coordinate for (unsigned int i = 0; i < spacedim; ++i) @@ -1648,15 +1628,11 @@ namespace internal for (unsigned int l = 0; l < dim; ++l) for (unsigned int m = 0; m < dim; ++m) { - jacobian_pushed_forward_2nd_derivatives - [point][i][j][l][m] = result[i][0][l][m] * - data.covariant[point][j][0]; + tmp[i][j][l][m] = + third[i][0][l][m] * covariant[j][0]; for (unsigned int jr = 1; jr < dim; ++jr) - jacobian_pushed_forward_2nd_derivatives[point][i] - [j][l] - [m] += - result[i][jr][l][m] * - data.covariant[point][j][jr]; + tmp[i][j][l][m] += + third[i][jr][l][m] * covariant[j][jr]; } // push forward the l-coordinate @@ -1665,16 +1641,11 @@ namespace internal for (unsigned int l = 0; l < spacedim; ++l) for (unsigned int m = 0; m < dim; ++m) { - tmp[i][j][l][m] = - jacobian_pushed_forward_2nd_derivatives[point][i] - [j][0][m] * - data.covariant[point][l][0]; + tmp2[i][j][l][m] = + tmp[i][j][0][m] * covariant[l][0]; for (unsigned int lr = 1; lr < dim; ++lr) - tmp[i][j][l][m] += - jacobian_pushed_forward_2nd_derivatives[point] - [i][j] - [lr][m] * - data.covariant[point][l][lr]; + tmp2[i][j][l][m] += + tmp[i][j][lr][m] * covariant[l][lr]; } // push forward the m-coordinate @@ -1685,12 +1656,12 @@ namespace internal { jacobian_pushed_forward_2nd_derivatives [point][i][j][l][m] = - tmp[i][j][l][0] * data.covariant[point][m][0]; + tmp2[i][j][l][0] * covariant[m][0]; for (unsigned int mr = 1; mr < dim; ++mr) jacobian_pushed_forward_2nd_derivatives[point][i] [j][l] [m] += - tmp[i][j][l][mr] * data.covariant[point][m][mr]; + tmp2[i][j][l][mr] * covariant[m][mr]; } } } @@ -1699,6 +1670,66 @@ namespace internal + template + inline DerivativeForm<4, dim, spacedim> + expand_4th_derivatives( + const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed) + { + Assert(dim >= 1 && dim <= 3, ExcNotImplemented()); + DerivativeForm<4, dim, spacedim> result; + for (unsigned int i = 0; i < spacedim; ++i) + { + if (dim == 1) + result[i][0][0][0][0] = compressed[0][i]; + else if (dim >= 2) + { + for (unsigned int d = 0; d < 2; ++d) + for (unsigned int e = 0; e < 2; ++e) + for (unsigned int f = 0; f < 2; ++f) + for (unsigned int g = 0; g < 2; ++g) + result[i][d][e][f][g] = compressed[d + e + f + g][i]; + if (dim == 3) + { + AssertDimension(length_tensor, 15); + + // the derivatives are ordered in rows by increasing z + // derivative, and in each row we have x^{(n-j)} y^{(j)} as + // j walks through the columns + for (unsigned int d = 0; d < 2; ++d) + for (unsigned int e = 0; e < 2; ++e) + for (unsigned int f = 0; f < 2; ++f) + { + result[i][d][e][f][2] = compressed[5 + d + e + f][i]; + result[i][d][e][2][f] = compressed[5 + d + e + f][i]; + result[i][d][2][e][f] = compressed[5 + d + e + f][i]; + result[i][2][d][e][f] = compressed[5 + d + e + f][i]; + } + for (unsigned int d = 0; d < 2; ++d) + for (unsigned int e = 0; e < 2; ++e) + { + result[i][d][e][2][2] = compressed[9 + d + e][i]; + result[i][d][2][e][2] = compressed[9 + d + e][i]; + result[i][d][2][2][e] = compressed[9 + d + e][i]; + result[i][2][d][e][2] = compressed[9 + d + e][i]; + result[i][2][d][2][e] = compressed[9 + d + e][i]; + result[i][2][2][d][e] = compressed[9 + d + e][i]; + } + for (unsigned int d = 0; d < 2; ++d) + { + result[i][d][2][2][2] = compressed[12 + d][i]; + result[i][2][d][2][2] = compressed[12 + d][i]; + result[i][2][2][d][2] = compressed[12 + d][i]; + result[i][2][2][2][d] = compressed[12 + d][i]; + } + result[i][2][2][2][2] = compressed[14][i]; + } + } + } + return result; + } + + + /** * Update the fourth derivatives of the transformation from unit to real * cell, the Jacobian hessian gradients. @@ -1708,48 +1739,29 @@ namespace internal template inline void maybe_update_jacobian_3rd_derivatives( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, + const CellSimilarity::Similarity cell_similarity, const typename dealii::MappingQ::InternalData &data, + const ArrayView> & unit_points, + const std::vector> &polynomials_1d, + const std::vector &renumber_lexicographic_to_hierarchic, std::vector> &jacobian_3rd_derivatives) { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_3rd_derivatives) + if (data.update_each & update_jacobian_3rd_derivatives) { + const std::vector> &support_points = + data.mapping_support_points; const unsigned int n_q_points = jacobian_3rd_derivatives.size(); if (cell_similarity != CellSimilarity::translation) { for (unsigned int point = 0; point < n_q_points; ++point) { - const Tensor<4, dim> *fourth = - &data.fourth_derivative(point + data_set, 0); - double result[spacedim][dim][dim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - result[i][j][l][m][n] = - (fourth[0][j][l][m][n] * - data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - result[i][j][l][m][n] += - (fourth[k][j][l][m][n] * - data.mapping_support_points[k][i]); - - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - jacobian_3rd_derivatives[point][i][j][l][m][n] = - result[i][j][l][m][n]; + jacobian_3rd_derivatives[point] = expand_4th_derivatives( + internal::evaluate_tensor_product_higher_derivatives<4>( + polynomials_1d, + support_points, + unit_points[point], + renumber_lexicographic_to_hierarchic)); } } } @@ -1767,42 +1779,39 @@ namespace internal template inline void maybe_update_jacobian_pushed_forward_3rd_derivatives( - const CellSimilarity::Similarity cell_similarity, - const typename QProjector::DataSetDescriptor data_set, + const CellSimilarity::Similarity cell_similarity, const typename dealii::MappingQ::InternalData &data, + const ArrayView> & unit_points, + const std::vector> &polynomials_1d, + const std::vector & renumber_lexicographic_to_hierarchic, std::vector> &jacobian_pushed_forward_3rd_derivatives) { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) + if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives) { + const std::vector> &support_points = + data.mapping_support_points; const unsigned int n_q_points = jacobian_pushed_forward_3rd_derivatives.size(); if (cell_similarity != CellSimilarity::translation) { - double tmp[spacedim][spacedim][spacedim][spacedim][spacedim]; + dealii:: + ndarray + tmp; + dealii::ndarray + tmp2; for (unsigned int point = 0; point < n_q_points; ++point) { - const Tensor<4, dim> *fourth = - &data.fourth_derivative(point + data_set, 0); - double result[spacedim][dim][dim][dim][dim]; - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - result[i][j][l][m][n] = - (fourth[0][j][l][m][n] * - data.mapping_support_points[0][i]); - for (unsigned int k = 1; k < data.n_shape_functions; ++k) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) - for (unsigned int n = 0; n < dim; ++n) - result[i][j][l][m][n] += - (fourth[k][j][l][m][n] * - data.mapping_support_points[k][i]); + const DerivativeForm<4, dim, spacedim> fourth = + expand_4th_derivatives( + internal::evaluate_tensor_product_higher_derivatives<4>( + polynomials_1d, + support_points, + unit_points[point], + renumber_lexicographic_to_hierarchic)); + + const DerivativeForm<1, dim, spacedim> covariant = + data.output_data->inverse_jacobians[point].transpose(); // push-forward the j-coordinate for (unsigned int i = 0; i < spacedim; ++i) @@ -1811,12 +1820,11 @@ namespace internal for (unsigned int m = 0; m < dim; ++m) for (unsigned int n = 0; n < dim; ++n) { - tmp[i][j][l][m][n] = result[i][0][l][m][n] * - data.covariant[point][j][0]; + tmp[i][j][l][m][n] = + fourth[i][0][l][m][n] * covariant[j][0]; for (unsigned int jr = 1; jr < dim; ++jr) tmp[i][j][l][m][n] += - result[i][jr][l][m][n] * - data.covariant[point][j][jr]; + fourth[i][jr][l][m][n] * covariant[j][jr]; } // push-forward the l-coordinate @@ -1826,16 +1834,11 @@ namespace internal for (unsigned int m = 0; m < dim; ++m) for (unsigned int n = 0; n < dim; ++n) { - jacobian_pushed_forward_3rd_derivatives - [point][i][j][l][m][n] = - tmp[i][j][0][m][n] * - data.covariant[point][l][0]; + tmp2[i][j][l][m][n] = + tmp[i][j][0][m][n] * covariant[l][0]; for (unsigned int lr = 1; lr < dim; ++lr) - jacobian_pushed_forward_3rd_derivatives[point] - [i][j][l] - [m][n] += - tmp[i][j][lr][m][n] * - data.covariant[point][l][lr]; + tmp2[i][j][l][m][n] += + tmp[i][j][lr][m][n] * covariant[l][lr]; } // push-forward the m-coordinate @@ -1846,15 +1849,10 @@ namespace internal for (unsigned int n = 0; n < dim; ++n) { tmp[i][j][l][m][n] = - jacobian_pushed_forward_3rd_derivatives[point] - [i][j][l] - [0][n] * - data.covariant[point][m][0]; + tmp2[i][j][l][0][n] * covariant[m][0]; for (unsigned int mr = 1; mr < dim; ++mr) tmp[i][j][l][m][n] += - jacobian_pushed_forward_3rd_derivatives - [point][i][j][l][mr][n] * - data.covariant[point][m][mr]; + tmp2[i][j][l][mr][n] * covariant[m][mr]; } // push-forward the n-coordinate @@ -1866,14 +1864,12 @@ namespace internal { jacobian_pushed_forward_3rd_derivatives [point][i][j][l][m][n] = - tmp[i][j][l][m][0] * - data.covariant[point][n][0]; + tmp[i][j][l][m][0] * covariant[n][0]; for (unsigned int nr = 1; nr < dim; ++nr) jacobian_pushed_forward_3rd_derivatives[point] [i][j][l] [m][n] += - tmp[i][j][l][m][nr] * - data.covariant[point][n][nr]; + tmp[i][j][l][m][nr] * covariant[n][nr]; } } } @@ -1907,8 +1903,7 @@ namespace internal const UpdateFlags update_flags = data.update_each; if (update_flags & - (update_boundary_forms | update_normal_vectors | update_jacobians | - update_JxW_values | update_inverse_jacobians)) + (update_boundary_forms | update_normal_vectors | update_JxW_values)) { if (update_flags & update_boundary_forms) AssertDimension(output_data.boundary_forms.size(), n_q_points); @@ -1927,24 +1922,18 @@ namespace internal // 0. for (unsigned int d = 0; d != dim - 1; ++d) { - Assert(face_no + GeometryInfo::faces_per_cell * d < - data.unit_tangentials.size(), + const unsigned int vector_index = + face_no + GeometryInfo::faces_per_cell * d; + Assert(vector_index < data.unit_tangentials.size(), ExcInternalError()); - Assert( - data.aux[d].size() <= - data - .unit_tangentials[face_no + - GeometryInfo::faces_per_cell * d] - .size(), - ExcInternalError()); - - mapping.transform( - make_array_view( - data.unit_tangentials[face_no + - GeometryInfo::faces_per_cell * d]), - mapping_contravariant, - data, - make_array_view(data.aux[d].begin(), data.aux[d].end())); + Assert(data.aux[d].size() <= + data.unit_tangentials[vector_index].size(), + ExcInternalError()); + mapping.transform(make_array_view( + data.unit_tangentials[vector_index]), + mapping_contravariant, + data, + make_array_view(data.aux[d])); } if (update_flags & update_boundary_forms) @@ -1985,15 +1974,18 @@ namespace internal // // to compute the cell normal, use the same method used in // fill_fe_values for cells above - AssertDimension(data.contravariant.size(), n_q_points); + AssertDimension(data.output_data->jacobians.size(), + n_q_points); for (unsigned int point = 0; point < n_q_points; ++point) { + const DerivativeForm<1, dim, spacedim> contravariant = + data.output_data->jacobians[point]; if (dim == 1) { // J is a tangent vector output_data.boundary_forms[point] = - data.contravariant[point].transpose()[0]; + contravariant.transpose()[0]; output_data.boundary_forms[point] /= (face_no == 0 ? -1. : +1.) * output_data.boundary_forms[point].norm(); @@ -2002,7 +1994,7 @@ namespace internal if (dim == 2) { const DerivativeForm<1, spacedim, dim> DX_t = - data.contravariant[point].transpose(); + contravariant.transpose(); Tensor<1, spacedim> cell_normal = cross_product_3d(DX_t[0], DX_t[1]); @@ -2058,11 +2050,12 @@ namespace internal const Quadrature & quadrature, const typename dealii::MappingQ::InternalData &data, const std::vector> &polynomials_1d, - const unsigned int polynomial_degree, const std::vector &renumber_lexicographic_to_hierarchic, internal::FEValuesImplementation::MappingRelatedData &output_data) { + const ArrayView> quadrature_points( + &data.quadrature_points[data_set], quadrature.size()); if (dim > 1 && data.tensor_product_quadrature) { maybe_update_q_points_Jacobians_and_grads_tensor( @@ -2079,42 +2072,54 @@ namespace internal maybe_update_q_points_Jacobians_generic( CellSimilarity::none, data, - make_array_view(&data.quadrature_points[data_set], - &data.quadrature_points[data_set] + - quadrature.size()), + quadrature_points, polynomials_1d, - polynomial_degree, renumber_lexicographic_to_hierarchic, output_data.quadrature_points, output_data.jacobians, output_data.inverse_jacobians); maybe_update_jacobian_grads( - CellSimilarity::none, data_set, data, output_data.jacobian_grads); + CellSimilarity::none, + data, + quadrature_points, + polynomials_1d, + renumber_lexicographic_to_hierarchic, + output_data.jacobian_grads); } maybe_update_jacobian_pushed_forward_grads( CellSimilarity::none, - data_set, data, + quadrature_points, + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_pushed_forward_grads); maybe_update_jacobian_2nd_derivatives( CellSimilarity::none, - data_set, data, + quadrature_points, + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_2nd_derivatives); maybe_update_jacobian_pushed_forward_2nd_derivatives( CellSimilarity::none, - data_set, data, + quadrature_points, + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_pushed_forward_2nd_derivatives); maybe_update_jacobian_3rd_derivatives( CellSimilarity::none, - data_set, data, + quadrature_points, + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_3rd_derivatives); maybe_update_jacobian_pushed_forward_3rd_derivatives( CellSimilarity::none, - data_set, data, + quadrature_points, + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_pushed_forward_3rd_derivatives); maybe_compute_face_data(mapping, @@ -2159,8 +2164,8 @@ namespace internal "update_contravariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) - output[i] = - apply_transformation(data.contravariant[i], input[i]); + output[i] = apply_transformation(data.output_data->jacobians[i], + input[i]); return; } @@ -2180,7 +2185,8 @@ namespace internal for (unsigned int i = 0; i < output.size(); ++i) { output[i] = - apply_transformation(data.contravariant[i], input[i]); + apply_transformation(data.output_data->jacobians[i], + input[i]); output[i] /= data.volume_elements[i]; } return; @@ -2195,7 +2201,8 @@ namespace internal "update_covariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) - output[i] = apply_transformation(data.covariant[i], input[i]); + output[i] = apply_transformation( + data.output_data->inverse_jacobians[i].transpose(), input[i]); return; } @@ -2243,10 +2250,11 @@ namespace internal for (unsigned int i = 0; i < output.size(); ++i) { const DerivativeForm<1, spacedim, dim> A = - apply_transformation(data.contravariant[i], + apply_transformation(data.output_data->jacobians[i], transpose(input[i])); - output[i] = - apply_transformation(data.covariant[i], A.transpose()); + output[i] = apply_transformation( + data.output_data->inverse_jacobians[i].transpose(), + A.transpose()); } return; @@ -2261,11 +2269,11 @@ namespace internal for (unsigned int i = 0; i < output.size(); ++i) { + const DerivativeForm<1, dim, spacedim> covariant = + data.output_data->inverse_jacobians[i].transpose(); const DerivativeForm<1, spacedim, dim> A = - apply_transformation(data.covariant[i], - transpose(input[i])); - output[i] = - apply_transformation(data.covariant[i], A.transpose()); + apply_transformation(covariant, transpose(input[i])); + output[i] = apply_transformation(covariant, A.transpose()); } return; @@ -2286,10 +2294,13 @@ namespace internal for (unsigned int i = 0; i < output.size(); ++i) { + const DerivativeForm<1, dim, spacedim> covariant = + data.output_data->inverse_jacobians[i].transpose(); const DerivativeForm<1, spacedim, dim> A = - apply_transformation(data.covariant[i], input[i]); + apply_transformation(covariant, input[i]); const Tensor<2, spacedim> T = - apply_transformation(data.contravariant[i], A.transpose()); + apply_transformation(data.output_data->jacobians[i], + A.transpose()); output[i] = transpose(T); output[i] /= data.volume_elements[i]; @@ -2338,37 +2349,42 @@ namespace internal "update_contravariant_transformation")); for (unsigned int q = 0; q < output.size(); ++q) - for (unsigned int i = 0; i < spacedim; ++i) - { - double tmp1[dim][dim]; - for (unsigned int J = 0; J < dim; ++J) - for (unsigned int K = 0; K < dim; ++K) - { - tmp1[J][K] = - data.contravariant[q][i][0] * input[q][0][J][K]; - for (unsigned int I = 1; I < dim; ++I) - tmp1[J][K] += - data.contravariant[q][i][I] * input[q][I][J][K]; - } - for (unsigned int j = 0; j < spacedim; ++j) - { - double tmp2[dim]; + { + const DerivativeForm<1, dim, spacedim> covariant = + data.output_data->inverse_jacobians[q].transpose(); + const DerivativeForm<1, dim, spacedim> contravariant = + data.output_data->jacobians[q]; + + for (unsigned int i = 0; i < spacedim; ++i) + { + double tmp1[dim][dim]; + for (unsigned int J = 0; J < dim; ++J) for (unsigned int K = 0; K < dim; ++K) { - tmp2[K] = data.covariant[q][j][0] * tmp1[0][K]; - for (unsigned int J = 1; J < dim; ++J) - tmp2[K] += data.covariant[q][j][J] * tmp1[J][K]; + tmp1[J][K] = + contravariant[i][0] * input[q][0][J][K]; + for (unsigned int I = 1; I < dim; ++I) + tmp1[J][K] += + contravariant[i][I] * input[q][I][J][K]; } - for (unsigned int k = 0; k < spacedim; ++k) - { - output[q][i][j][k] = - data.covariant[q][k][0] * tmp2[0]; - for (unsigned int K = 1; K < dim; ++K) - output[q][i][j][k] += - data.covariant[q][k][K] * tmp2[K]; - } - } - } + for (unsigned int j = 0; j < spacedim; ++j) + { + double tmp2[dim]; + for (unsigned int K = 0; K < dim; ++K) + { + tmp2[K] = covariant[j][0] * tmp1[0][K]; + for (unsigned int J = 1; J < dim; ++J) + tmp2[K] += covariant[j][J] * tmp1[J][K]; + } + for (unsigned int k = 0; k < spacedim; ++k) + { + output[q][i][j][k] = covariant[k][0] * tmp2[0]; + for (unsigned int K = 1; K < dim; ++K) + output[q][i][j][k] += covariant[k][K] * tmp2[K]; + } + } + } + } return; } @@ -2379,37 +2395,38 @@ namespace internal "update_covariant_transformation")); for (unsigned int q = 0; q < output.size(); ++q) - for (unsigned int i = 0; i < spacedim; ++i) - { - double tmp1[dim][dim]; - for (unsigned int J = 0; J < dim; ++J) - for (unsigned int K = 0; K < dim; ++K) - { - tmp1[J][K] = - data.covariant[q][i][0] * input[q][0][J][K]; - for (unsigned int I = 1; I < dim; ++I) - tmp1[J][K] += - data.covariant[q][i][I] * input[q][I][J][K]; - } - for (unsigned int j = 0; j < spacedim; ++j) - { - double tmp2[dim]; + { + const DerivativeForm<1, dim, spacedim> covariant = + data.output_data->inverse_jacobians[q].transpose(); + + for (unsigned int i = 0; i < spacedim; ++i) + { + double tmp1[dim][dim]; + for (unsigned int J = 0; J < dim; ++J) for (unsigned int K = 0; K < dim; ++K) { - tmp2[K] = data.covariant[q][j][0] * tmp1[0][K]; - for (unsigned int J = 1; J < dim; ++J) - tmp2[K] += data.covariant[q][j][J] * tmp1[J][K]; - } - for (unsigned int k = 0; k < spacedim; ++k) - { - output[q][i][j][k] = - data.covariant[q][k][0] * tmp2[0]; - for (unsigned int K = 1; K < dim; ++K) - output[q][i][j][k] += - data.covariant[q][k][K] * tmp2[K]; + tmp1[J][K] = covariant[i][0] * input[q][0][J][K]; + for (unsigned int I = 1; I < dim; ++I) + tmp1[J][K] += covariant[i][I] * input[q][I][J][K]; } - } - } + for (unsigned int j = 0; j < spacedim; ++j) + { + double tmp2[dim]; + for (unsigned int K = 0; K < dim; ++K) + { + tmp2[K] = covariant[j][0] * tmp1[0][K]; + for (unsigned int J = 1; J < dim; ++J) + tmp2[K] += covariant[j][J] * tmp1[J][K]; + } + for (unsigned int k = 0; k < spacedim; ++k) + { + output[q][i][j][k] = covariant[k][0] * tmp2[0]; + for (unsigned int K = 1; K < dim; ++K) + output[q][i][j][k] += covariant[k][K] * tmp2[K]; + } + } + } + } return; } @@ -2427,39 +2444,43 @@ namespace internal "update_volume_elements")); for (unsigned int q = 0; q < output.size(); ++q) - for (unsigned int i = 0; i < spacedim; ++i) - { - double factor[dim]; - for (unsigned int I = 0; I < dim; ++I) - factor[I] = - data.contravariant[q][i][I] / data.volume_elements[q]; - double tmp1[dim][dim]; - for (unsigned int J = 0; J < dim; ++J) - for (unsigned int K = 0; K < dim; ++K) - { - tmp1[J][K] = factor[0] * input[q][0][J][K]; - for (unsigned int I = 1; I < dim; ++I) - tmp1[J][K] += factor[I] * input[q][I][J][K]; - } - for (unsigned int j = 0; j < spacedim; ++j) - { - double tmp2[dim]; + { + const DerivativeForm<1, dim, spacedim> covariant = + data.output_data->inverse_jacobians[q].transpose(); + const DerivativeForm<1, dim, spacedim> contravariant = + data.output_data->jacobians[q]; + for (unsigned int i = 0; i < spacedim; ++i) + { + double factor[dim]; + for (unsigned int I = 0; I < dim; ++I) + factor[I] = + contravariant[i][I] * (1. / data.volume_elements[q]); + double tmp1[dim][dim]; + for (unsigned int J = 0; J < dim; ++J) for (unsigned int K = 0; K < dim; ++K) { - tmp2[K] = data.covariant[q][j][0] * tmp1[0][K]; - for (unsigned int J = 1; J < dim; ++J) - tmp2[K] += data.covariant[q][j][J] * tmp1[J][K]; + tmp1[J][K] = factor[0] * input[q][0][J][K]; + for (unsigned int I = 1; I < dim; ++I) + tmp1[J][K] += factor[I] * input[q][I][J][K]; } - for (unsigned int k = 0; k < spacedim; ++k) - { - output[q][i][j][k] = - data.covariant[q][k][0] * tmp2[0]; - for (unsigned int K = 1; K < dim; ++K) - output[q][i][j][k] += - data.covariant[q][k][K] * tmp2[K]; - } - } - } + for (unsigned int j = 0; j < spacedim; ++j) + { + double tmp2[dim]; + for (unsigned int K = 0; K < dim; ++K) + { + tmp2[K] = covariant[j][0] * tmp1[0][K]; + for (unsigned int J = 1; J < dim; ++J) + tmp2[K] += covariant[j][J] * tmp1[J][K]; + } + for (unsigned int k = 0; k < spacedim; ++k) + { + output[q][i][j][k] = covariant[k][0] * tmp2[0]; + for (unsigned int K = 1; K < dim; ++K) + output[q][i][j][k] += covariant[k][K] * tmp2[K]; + } + } + } + } return; } @@ -2502,7 +2523,8 @@ namespace internal "update_covariant_transformation")); for (unsigned int i = 0; i < output.size(); ++i) - output[i] = apply_transformation(data.covariant[i], input[i]); + output[i] = apply_transformation( + data.output_data->inverse_jacobians[i].transpose(), input[i]); return; } diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 2961f73026..b1e8eb1676 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -54,6 +54,7 @@ MappingQ::InternalData::InternalData( , n_shape_functions(Utilities::fixed_power(polynomial_degree + 1)) , line_support_points(QGaussLobatto<1>(polynomial_degree + 1)) , tensor_product_quadrature(false) + , output_data(nullptr) {} @@ -64,10 +65,7 @@ MappingQ::InternalData::memory_consumption() const { return ( Mapping::InternalDataBase::memory_consumption() + - MemoryConsumption::memory_consumption(shape_values) + - MemoryConsumption::memory_consumption(shape_derivatives) + - MemoryConsumption::memory_consumption(covariant) + - MemoryConsumption::memory_consumption(contravariant) + + MemoryConsumption::memory_consumption(quadrature_points) + MemoryConsumption::memory_consumption(unit_tangentials) + MemoryConsumption::memory_consumption(aux) + MemoryConsumption::memory_consumption(mapping_support_points) + @@ -92,12 +90,6 @@ MappingQ::InternalData::initialize( const unsigned int n_q_points = quadrature.size(); - if (this->update_each & update_covariant_transformation) - covariant.resize(n_original_q_points); - - if (this->update_each & update_contravariant_transformation) - contravariant.resize(n_original_q_points); - if (this->update_each & update_volume_elements) volume_elements.resize(n_original_q_points); @@ -161,49 +153,6 @@ MappingQ::InternalData::initialize( } } } - - const bool needs_higher_order_terms = - this->update_each & - (update_jacobian_pushed_forward_grads | update_jacobian_2nd_derivatives | - update_jacobian_pushed_forward_2nd_derivatives | - update_jacobian_3rd_derivatives | - update_jacobian_pushed_forward_3rd_derivatives); - - const bool needs_higher_order_terms_generic = - !tensor_product_quadrature && - (needs_higher_order_terms || this->update_each & update_jacobian_grads); - - // Only fill the big arrays on demand in case we cannot use the tensor - // product quadrature code path - if (dim == 1 || needs_higher_order_terms_generic || needs_higher_order_terms) - { - // compute shapes and derivatives for codim1 (for - // do_transform_real_to_unit_cell_internal_codim1) - if (dim == (spacedim - 1)) - { - // see if we need the (transformation) shape function values - // and/or gradients and resize the necessary arrays - if (this->update_each & update_quadrature_points) - shape_values.resize(n_shape_functions * n_q_points); - if (this->update_each & update_jacobians) - shape_derivatives.resize(n_shape_functions * n_q_points); - } - - if (this->update_each & - (update_jacobian_grads | update_jacobian_pushed_forward_grads)) - shape_second_derivatives.resize(n_shape_functions * n_q_points); - - if (this->update_each & (update_jacobian_2nd_derivatives | - update_jacobian_pushed_forward_2nd_derivatives)) - shape_third_derivatives.resize(n_shape_functions * n_q_points); - - if (this->update_each & (update_jacobian_3rd_derivatives | - update_jacobian_pushed_forward_3rd_derivatives)) - shape_fourth_derivatives.resize(n_shape_functions * n_q_points); - - // now also fill the various fields with their correct values - compute_shape_function_values(quadrature.get_points()); - } } @@ -217,10 +166,7 @@ MappingQ::InternalData::initialize_face( { initialize(update_flags, quadrature, n_original_q_points); - const unsigned int n_q_points = quadrature.size(); - quadrature_points.resize(n_q_points); - for (unsigned int q = 0; q < n_q_points; ++q) - quadrature_points[q] = quadrature.get_points()[q]; + quadrature_points = quadrature.get_points(); if (dim > 1 && tensor_product_quadrature) { @@ -238,8 +184,7 @@ MappingQ::InternalData::initialize_face( if (dim > 1) { if (this->update_each & - (update_boundary_forms | update_normal_vectors | update_jacobians | - update_JxW_values | update_inverse_jacobians)) + (update_boundary_forms | update_normal_vectors | update_JxW_values)) { aux.resize(dim - 1, AlignedVector>(n_original_q_points)); @@ -269,98 +214,6 @@ MappingQ::InternalData::initialize_face( -template -void -MappingQ::InternalData::compute_shape_function_values( - const std::vector> &unit_points) -{ - const unsigned int n_points = unit_points.size(); - - // Construct the tensor product polynomials used as shape functions for - // the Qp mapping of cells at the boundary. - const TensorProductPolynomials tensor_pols( - Polynomials::generate_complete_Lagrange_basis( - line_support_points.get_points())); - Assert(n_shape_functions == tensor_pols.n(), ExcInternalError()); - - // then also construct the mapping from lexicographic to the Qp shape - // function numbering - const std::vector renumber = - FETools::hierarchic_to_lexicographic_numbering(polynomial_degree); - - std::vector values; - std::vector> grads; - if (shape_values.size() != 0) - { - Assert(shape_values.size() == n_shape_functions * n_points, - ExcInternalError()); - values.resize(n_shape_functions); - } - if (shape_derivatives.size() != 0) - { - Assert(shape_derivatives.size() == n_shape_functions * n_points, - ExcInternalError()); - grads.resize(n_shape_functions); - } - - std::vector> grad2; - if (shape_second_derivatives.size() != 0) - { - Assert(shape_second_derivatives.size() == n_shape_functions * n_points, - ExcInternalError()); - grad2.resize(n_shape_functions); - } - - std::vector> grad3; - if (shape_third_derivatives.size() != 0) - { - Assert(shape_third_derivatives.size() == n_shape_functions * n_points, - ExcInternalError()); - grad3.resize(n_shape_functions); - } - - std::vector> grad4; - if (shape_fourth_derivatives.size() != 0) - { - Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points, - ExcInternalError()); - grad4.resize(n_shape_functions); - } - - - if (shape_values.size() != 0 || shape_derivatives.size() != 0 || - shape_second_derivatives.size() != 0 || - shape_third_derivatives.size() != 0 || - shape_fourth_derivatives.size() != 0) - for (unsigned int point = 0; point < n_points; ++point) - { - tensor_pols.evaluate( - unit_points[point], values, grads, grad2, grad3, grad4); - - if (shape_values.size() != 0) - for (unsigned int i = 0; i < n_shape_functions; ++i) - shape(point, i) = values[renumber[i]]; - - if (shape_derivatives.size() != 0) - for (unsigned int i = 0; i < n_shape_functions; ++i) - derivative(point, i) = grads[renumber[i]]; - - if (shape_second_derivatives.size() != 0) - for (unsigned int i = 0; i < n_shape_functions; ++i) - second_derivative(point, i) = grad2[renumber[i]]; - - if (shape_third_derivatives.size() != 0) - for (unsigned int i = 0; i < n_shape_functions; ++i) - third_derivative(point, i) = grad3[renumber[i]]; - - if (shape_fourth_derivatives.size() != 0) - for (unsigned int i = 0; i < n_shape_functions; ++i) - fourth_derivative(point, i) = grad4[renumber[i]]; - } -} - - - template MappingQ::MappingQ(const unsigned int p) : polynomial_degree(p) @@ -392,30 +245,8 @@ MappingQ::MappingQ(const unsigned int p) template MappingQ::MappingQ(const unsigned int p, const bool) - : polynomial_degree(p) - , line_support_points( - QGaussLobatto<1>(this->polynomial_degree + 1).get_points()) - , polynomials_1d( - Polynomials::generate_complete_Lagrange_basis(line_support_points)) - , renumber_lexicographic_to_hierarchic( - FETools::lexicographic_to_hierarchic_numbering(p)) - , unit_cell_support_points( - internal::MappingQImplementation::unit_support_points( - line_support_points, - renumber_lexicographic_to_hierarchic)) - , support_point_weights_perimeter_to_interior( - internal::MappingQImplementation:: - compute_support_point_weights_perimeter_to_interior( - this->polynomial_degree, - dim)) - , support_point_weights_cell( - internal::MappingQImplementation::compute_support_point_weights_cell( - this->polynomial_degree)) -{ - Assert(p >= 1, - ExcMessage("It only makes sense to create polynomial mappings " - "with a polynomial degree greater or equal to one.")); -} + : MappingQ(p) +{} @@ -579,10 +410,12 @@ MappingQ<1, 2>::transform_real_to_unit_cell_internal( // dispatch to the various specializations for spacedim=dim, // spacedim=dim+1, etc return internal::MappingQImplementation:: - do_transform_real_to_unit_cell_internal_codim1<1>(cell, - p, - initial_p_unit, - *mdata); + do_transform_real_to_unit_cell_internal_codim1<1>( + p, + initial_p_unit, + mdata->mapping_support_points, + polynomials_1d, + renumber_lexicographic_to_hierarchic); } @@ -610,12 +443,16 @@ MappingQ<2, 3>::transform_real_to_unit_cell_internal( // dispatch to the various specializations for spacedim=dim, // spacedim=dim+1, etc return internal::MappingQImplementation:: - do_transform_real_to_unit_cell_internal_codim1<2>(cell, - p, - initial_p_unit, - *mdata); + do_transform_real_to_unit_cell_internal_codim1<2>( + p, + initial_p_unit, + mdata->mapping_support_points, + polynomials_1d, + renumber_lexicographic_to_hierarchic); } + + template <> Point<1> MappingQ<1, 3>::transform_real_to_unit_cell_internal( @@ -857,9 +694,9 @@ MappingQ::requires_update_flags(const UpdateFlags in) const if (out & (update_JxW_values | update_normal_vectors)) out |= update_boundary_forms; - if (out & (update_covariant_transformation | update_JxW_values | - update_jacobians | update_jacobian_grads | - update_boundary_forms | update_normal_vectors)) + if (out & + (update_covariant_transformation | update_jacobian_grads | + update_jacobians | update_boundary_forms | update_normal_vectors)) out |= update_contravariant_transformation; if (out & @@ -957,6 +794,7 @@ MappingQ::fill_fe_values( Assert(dynamic_cast(&internal_data) != nullptr, ExcInternalError()); const InternalData &data = static_cast(internal_data); + data.output_data = &output_data; const unsigned int n_q_points = quadrature.size(); @@ -1000,7 +838,6 @@ MappingQ::fill_fe_values( data, make_array_view(quadrature.get_points()), polynomials_1d, - polynomial_degree, renumber_lexicographic_to_hierarchic, output_data.quadrature_points, output_data.jacobians, @@ -1009,44 +846,56 @@ MappingQ::fill_fe_values( internal::MappingQImplementation::maybe_update_jacobian_grads( computed_cell_similarity, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_grads); } internal::MappingQImplementation::maybe_update_jacobian_pushed_forward_grads< dim, spacedim>(computed_cell_similarity, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_pushed_forward_grads); internal::MappingQImplementation::maybe_update_jacobian_2nd_derivatives< dim, spacedim>(computed_cell_similarity, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_2nd_derivatives); internal::MappingQImplementation:: maybe_update_jacobian_pushed_forward_2nd_derivatives( computed_cell_similarity, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_pushed_forward_2nd_derivatives); internal::MappingQImplementation::maybe_update_jacobian_3rd_derivatives< dim, spacedim>(computed_cell_similarity, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_3rd_derivatives); internal::MappingQImplementation:: maybe_update_jacobian_pushed_forward_3rd_derivatives( computed_cell_similarity, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_pushed_forward_3rd_derivatives); const UpdateFlags update_flags = data.update_each; @@ -1070,7 +919,7 @@ MappingQ::fill_fe_values( { if (dim == spacedim) { - const double det = data.contravariant[point].determinant(); + const double det = data.volume_elements[point]; // check for distorted cells. @@ -1093,7 +942,7 @@ MappingQ::fill_fe_values( Tensor<1, spacedim> DX_t[dim]; for (unsigned int i = 0; i < spacedim; ++i) for (unsigned int j = 0; j < dim; ++j) - DX_t[j][i] = data.contravariant[point][i][j]; + DX_t[j][i] = output_data.jacobians[point][i][j]; Tensor<2, dim> G; // First fundamental form for (unsigned int i = 0; i < dim; ++i) @@ -1163,11 +1012,12 @@ MappingQ::fill_fe_face_values( Assert((dynamic_cast(&internal_data) != nullptr), ExcInternalError()); const InternalData &data = static_cast(internal_data); + data.output_data = &output_data; // if necessary, recompute the support points of the transformation of this // cell (note that we need to first check the triangulation pointer, since - // otherwise the second test might trigger an exception if the triangulations - // are not the same) + // otherwise the second test might trigger an exception if the + // triangulations are not the same) if ((data.mapping_support_points.size() == 0) || (&cell->get_triangulation() != &data.cell_of_current_support_points->get_triangulation()) || @@ -1192,7 +1042,6 @@ MappingQ::fill_fe_face_values( quadrature[0], data, polynomials_1d, - polynomial_degree, renumber_lexicographic_to_hierarchic, output_data); } @@ -1214,11 +1063,12 @@ MappingQ::fill_fe_subface_values( Assert((dynamic_cast(&internal_data) != nullptr), ExcInternalError()); const InternalData &data = static_cast(internal_data); + data.output_data = &output_data; // if necessary, recompute the support points of the transformation of this // cell (note that we need to first check the triangulation pointer, since - // otherwise the second test might trigger an exception if the triangulations - // are not the same) + // otherwise the second test might trigger an exception if the + // triangulations are not the same) if ((data.mapping_support_points.size() == 0) || (&cell->get_triangulation() != &data.cell_of_current_support_points->get_triangulation()) || @@ -1245,7 +1095,6 @@ MappingQ::fill_fe_subface_values( quadrature, data, polynomials_1d, - polynomial_degree, renumber_lexicographic_to_hierarchic, output_data); } @@ -1267,6 +1116,7 @@ MappingQ::fill_fe_immersed_surface_values( Assert(dynamic_cast(&internal_data) != nullptr, ExcInternalError()); const InternalData &data = static_cast(internal_data); + data.output_data = &output_data; const unsigned int n_q_points = quadrature.size(); @@ -1278,7 +1128,6 @@ MappingQ::fill_fe_immersed_surface_values( data, make_array_view(quadrature.get_points()), polynomials_1d, - polynomial_degree, renumber_lexicographic_to_hierarchic, output_data.quadrature_points, output_data.jacobians, @@ -1286,43 +1135,55 @@ MappingQ::fill_fe_immersed_surface_values( internal::MappingQImplementation::maybe_update_jacobian_grads( CellSimilarity::none, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_grads); internal::MappingQImplementation::maybe_update_jacobian_pushed_forward_grads< dim, spacedim>(CellSimilarity::none, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_pushed_forward_grads); internal::MappingQImplementation::maybe_update_jacobian_2nd_derivatives< dim, spacedim>(CellSimilarity::none, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_2nd_derivatives); internal::MappingQImplementation:: maybe_update_jacobian_pushed_forward_2nd_derivatives( CellSimilarity::none, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_pushed_forward_2nd_derivatives); internal::MappingQImplementation::maybe_update_jacobian_3rd_derivatives< dim, spacedim>(CellSimilarity::none, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_3rd_derivatives); internal::MappingQImplementation:: maybe_update_jacobian_pushed_forward_3rd_derivatives( CellSimilarity::none, - QProjector::DataSetDescriptor::cell(), data, + make_array_view(quadrature.get_points()), + polynomials_1d, + renumber_lexicographic_to_hierarchic, output_data.jacobian_pushed_forward_3rd_derivatives); const UpdateFlags update_flags = data.update_each; @@ -1340,7 +1201,7 @@ MappingQ::fill_fe_immersed_surface_values( for (unsigned int point = 0; point < n_q_points; ++point) { - const double det = data.contravariant[point].determinant(); + const double det = data.volume_elements[point]; // check for distorted cells. @@ -1355,8 +1216,8 @@ MappingQ::fill_fe_immersed_surface_values( // The normals are n = J^{-T} * \hat{n} before normalizing. Tensor<1, spacedim> normal; for (unsigned int d = 0; d < spacedim; d++) - normal[d] = - data.covariant[point][d] * quadrature.normal_vector(point); + normal[d] = output_data.inverse_jacobians[point].transpose()[d] * + quadrature.normal_vector(point); output_data.JxW_values[point] = weights[point] * det * normal.norm(); @@ -1395,6 +1256,7 @@ MappingQ::fill_mapping_data_for_generic_points( Quadrature(std::vector>(unit_points.begin(), unit_points.end()))); const InternalData &data = static_cast(*internal_data); + data.output_data = &output_data; data.mapping_support_points = this->compute_mapping_support_points(cell); internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic( @@ -1402,7 +1264,6 @@ MappingQ::fill_mapping_data_for_generic_points( data, unit_points, polynomials_1d, - polynomial_degree, renumber_lexicographic_to_hierarchic, output_data.quadrature_points, output_data.jacobians, @@ -1430,6 +1291,7 @@ MappingQ::fill_mapping_data_for_face_quadrature( const InternalData &data = static_cast(internal_data); data.mapping_support_points = this->compute_mapping_support_points(cell); + data.output_data = &output_data; internal::MappingQImplementation::do_fill_fe_face_values( *this, @@ -1440,7 +1302,6 @@ MappingQ::fill_mapping_data_for_face_quadrature( face_quadrature, data, polynomials_1d, - polynomial_degree, renumber_lexicographic_to_hierarchic, output_data); } @@ -1522,13 +1383,14 @@ MappingQ::transform( AssertDimension(input.size(), output.size()); Assert(dynamic_cast(&mapping_data) != nullptr, ExcInternalError()); - const InternalData &data = static_cast(mapping_data); + const internal::FEValuesImplementation::MappingRelatedData + &data = *static_cast(mapping_data).output_data; switch (mapping_kind) { case mapping_covariant_gradient: { - Assert(data.update_each & update_contravariant_transformation, + Assert(!data.inverse_jacobians.empty(), typename FEValuesBase::ExcAccessToUninitializedField( "update_covariant_transformation")); @@ -1536,18 +1398,20 @@ MappingQ::transform( for (unsigned int i = 0; i < spacedim; ++i) for (unsigned int j = 0; j < spacedim; ++j) { - double tmp[dim]; + double tmp[dim]; + const DerivativeForm<1, dim, spacedim> covariant = + data.inverse_jacobians[q].transpose(); for (unsigned int K = 0; K < dim; ++K) { - tmp[K] = data.covariant[q][j][0] * input[q][i][0][K]; + tmp[K] = covariant[j][0] * input[q][i][0][K]; for (unsigned int J = 1; J < dim; ++J) - tmp[K] += data.covariant[q][j][J] * input[q][i][J][K]; + tmp[K] += covariant[j][J] * input[q][i][J][K]; } for (unsigned int k = 0; k < spacedim; ++k) { - output[q][i][j][k] = data.covariant[q][k][0] * tmp[0]; + output[q][i][j][k] = covariant[k][0] * tmp[0]; for (unsigned int K = 1; K < dim; ++K) - output[q][i][j][k] += data.covariant[q][k][K] * tmp[K]; + output[q][i][j][k] += covariant[k][K] * tmp[K]; } } return;