From a3a619e5ae1f7855a1b4d25a6cb8af31c2add833 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 26 May 2025 17:08:11 -0600 Subject: [PATCH] Avoid the use of anonymous namespaces in a header file. --- .../deal.II/fe/mapping_fe_field.templates.h | 1690 ++++++++--------- 1 file changed, 830 insertions(+), 860 deletions(-) diff --git a/include/deal.II/fe/mapping_fe_field.templates.h b/include/deal.II/fe/mapping_fe_field.templates.h index ea7b935647..c132cced57 100644 --- a/include/deal.II/fe/mapping_fe_field.templates.h +++ b/include/deal.II/fe/mapping_fe_field.templates.h @@ -687,819 +687,797 @@ namespace internal { namespace MappingFEFieldImplementation { - namespace + /** + * Compute the locations of quadrature points on the object described by + * the first argument (and the cell for which the mapping support points + * have already been set), but only if the update_flags of the @p data + * argument indicate so. + */ + template + void + maybe_compute_q_points( + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField:: + InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector> &quadrature_points) { - /** - * Compute the locations of quadrature points on the object described by - * the first argument (and the cell for which the mapping support points - * have already been set), but only if the update_flags of the @p data - * argument indicate so. - */ - template - void - maybe_compute_q_points( - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField:: - InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector> &quadrature_points) - { - const UpdateFlags update_flags = data.update_each; - - if (update_flags & update_quadrature_points) - { - for (unsigned int point = 0; point < quadrature_points.size(); - ++point) - { - Point result; - const double *shape = &data.shape(point + data_set, 0); + const UpdateFlags update_flags = data.update_each; - for (unsigned int k = 0; k < data.n_shape_functions; ++k) - { - const unsigned int comp_k = - fe.system_to_component_index(k).first; - if (fe_mask[comp_k]) - result[fe_to_real[comp_k]] += - data.local_dof_values[k] * shape[k]; - } + if (update_flags & update_quadrature_points) + { + for (unsigned int point = 0; point < quadrature_points.size(); + ++point) + { + Point result; + const double *shape = &data.shape(point + data_set, 0); - quadrature_points[point] = result; - } - } - } - - /** - * Update the co- and contravariant matrices as well as their determinant, - * for the cell described stored in the data object, but only if the - * update_flags of the @p data argument indicate so. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - void - maybe_update_Jacobians( - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField:: - InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real) - { - const UpdateFlags update_flags = data.update_each; - - // then Jacobians - if (update_flags & update_contravariant_transformation) - { - const unsigned int n_q_points = data.contravariant.size(); - - Assert(data.n_shape_functions > 0, ExcInternalError()); - - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<1, dim> *data_derv = - &data.derivative(point + data_set, 0); + for (unsigned int k = 0; k < data.n_shape_functions; ++k) + { + const unsigned int comp_k = + fe.system_to_component_index(k).first; + if (fe_mask[comp_k]) + result[fe_to_real[comp_k]] += + data.local_dof_values[k] * shape[k]; + } - Tensor<1, dim> result[spacedim]; + quadrature_points[point] = result; + } + } + } - for (unsigned int k = 0; k < data.n_shape_functions; ++k) - { - const unsigned int comp_k = - fe.system_to_component_index(k).first; - if (fe_mask[comp_k]) - result[fe_to_real[comp_k]] += - data.local_dof_values[k] * data_derv[k]; - } + /** + * Update the co- and contravariant matrices as well as their determinant, + * for the cell described stored in the data object, but only if the + * update_flags of the @p data argument indicate so. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + void + maybe_update_Jacobians( + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField:: + InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real) + { + const UpdateFlags update_flags = data.update_each; - // write result into contravariant data - for (unsigned int i = 0; i < spacedim; ++i) - { - data.contravariant[point][i] = result[i]; - } - } - } - - if (update_flags & update_covariant_transformation) - { - AssertDimension(data.covariant.size(), data.contravariant.size()); - for (unsigned int point = 0; point < data.contravariant.size(); - ++point) - data.covariant[point] = - (data.contravariant[point]).covariant_form(); - } - - if (update_flags & update_volume_elements) - { - AssertDimension(data.contravariant.size(), - data.volume_elements.size()); - for (unsigned int point = 0; point < data.contravariant.size(); - ++point) - data.volume_elements[point] = - data.contravariant[point].determinant(); - } - } - - /** - * Update the Hessian of the transformation from unit to real cell, the - * Jacobian gradients. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - void - maybe_update_jacobian_grads( - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField:: - InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector> &jacobian_grads) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_grads) - { - const unsigned int n_q_points = jacobian_grads.size(); - - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<2, dim> *second = - &data.second_derivative(point + data_set, 0); + // then Jacobians + if (update_flags & update_contravariant_transformation) + { + const unsigned int n_q_points = data.contravariant.size(); - DerivativeForm<2, dim, spacedim> result; + Assert(data.n_shape_functions > 0, ExcInternalError()); - for (unsigned int k = 0; k < data.n_shape_functions; ++k) - { - const unsigned int comp_k = - fe.system_to_component_index(k).first; - if (fe_mask[comp_k]) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - result[fe_to_real[comp_k]][j][l] += - (second[k][j][l] * data.local_dof_values[k]); - } + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<1, dim> *data_derv = + &data.derivative(point + data_set, 0); - // never touch any data for j=dim in case dim - void - maybe_update_jacobian_pushed_forward_grads( - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField:: - InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector> &jacobian_pushed_forward_grads) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_grads) - { - const unsigned int n_q_points = - jacobian_pushed_forward_grads.size(); - - 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); + Tensor<1, dim> result[spacedim]; - DerivativeForm<2, dim, spacedim> result; + for (unsigned int k = 0; k < data.n_shape_functions; ++k) + { + const unsigned int comp_k = + fe.system_to_component_index(k).first; + if (fe_mask[comp_k]) + result[fe_to_real[comp_k]] += + data.local_dof_values[k] * data_derv[k]; + } - for (unsigned int k = 0; k < data.n_shape_functions; ++k) - { - const unsigned int comp_k = - fe.system_to_component_index(k).first; - if (fe_mask[comp_k]) - for (unsigned int j = 0; j < dim; ++j) - for (unsigned int l = 0; l < dim; ++l) - result[fe_to_real[comp_k]][j][l] += - (second[k][j][l] * data.local_dof_values[k]); - } + // write result into contravariant data + for (unsigned int i = 0; i < spacedim; ++i) + { + data.contravariant[point][i] = result[i]; + } + } + } - // 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) + if (update_flags & update_covariant_transformation) + { + AssertDimension(data.covariant.size(), data.contravariant.size()); + for (unsigned int point = 0; point < data.contravariant.size(); + ++point) + data.covariant[point] = + (data.contravariant[point]).covariant_form(); + } + + if (update_flags & update_volume_elements) + { + AssertDimension(data.contravariant.size(), + data.volume_elements.size()); + for (unsigned int point = 0; point < data.contravariant.size(); + ++point) + data.volume_elements[point] = + data.contravariant[point].determinant(); + } + } + + /** + * Update the Hessian of the transformation from unit to real cell, the + * Jacobian gradients. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + void + maybe_update_jacobian_grads( + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField:: + InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector> &jacobian_grads) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_grads) + { + const unsigned int n_q_points = jacobian_grads.size(); + + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<2, dim> *second = + &data.second_derivative(point + data_set, 0); + + DerivativeForm<2, dim, spacedim> result; + + for (unsigned int k = 0; k < data.n_shape_functions; ++k) + { + const unsigned int comp_k = + fe.system_to_component_index(k).first; + if (fe_mask[comp_k]) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + result[fe_to_real[comp_k]][j][l] += + (second[k][j][l] * data.local_dof_values[k]); + } + + // never touch any data for j=dim in case dim + void + maybe_update_jacobian_pushed_forward_grads( + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField:: + InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector> &jacobian_pushed_forward_grads) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_grads) + { + const unsigned int n_q_points = jacobian_pushed_forward_grads.size(); + + 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); + + DerivativeForm<2, dim, spacedim> result; + + for (unsigned int k = 0; k < data.n_shape_functions; ++k) + { + const unsigned int comp_k = + fe.system_to_component_index(k).first; + if (fe_mask[comp_k]) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int l = 0; l < dim; ++l) + result[fe_to_real[comp_k]][j][l] += + (second[k][j][l] * data.local_dof_values[k]); + } + + // 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]; + for (unsigned int jr = 1; jr < dim; ++jr) + { + tmp[i][j][l] += + result[i][jr][l] * data.covariant[point][j][jr]; + } + } + + // now, pushing forward the l-components + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + 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]; + 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]; + } + } + } + } + } + + /** + * Update the third derivative of the transformation from unit to real + * cell, the Jacobian hessians. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + void + maybe_update_jacobian_2nd_derivatives( + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField:: + InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector> &jacobian_2nd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_2nd_derivatives) + { + const unsigned int n_q_points = jacobian_2nd_derivatives.size(); + + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<3, dim> *third = + &data.third_derivative(point + data_set, 0); + + DerivativeForm<3, dim, spacedim> result; + + for (unsigned int k = 0; k < data.n_shape_functions; ++k) + { + const unsigned int comp_k = + fe.system_to_component_index(k).first; + if (fe_mask[comp_k]) + 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[fe_to_real[comp_k]][j][l][m] += + (third[k][j][l][m] * data.local_dof_values[k]); + } + + // never touch any data for j=dim in case dim + void + maybe_update_jacobian_pushed_forward_2nd_derivatives( + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField:: + InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector> &jacobian_pushed_forward_2nd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) + { + const unsigned int n_q_points = + jacobian_pushed_forward_2nd_derivatives.size(); + + double tmp[spacedim][spacedim][spacedim][spacedim]; + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<3, dim> *third = + &data.third_derivative(point + data_set, 0); + + DerivativeForm<3, dim, spacedim> result; + + for (unsigned int k = 0; k < data.n_shape_functions; ++k) + { + const unsigned int comp_k = + fe.system_to_component_index(k).first; + if (fe_mask[comp_k]) + 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[fe_to_real[comp_k]][j][l][m] += + (third[k][j][l][m] * data.local_dof_values[k]); + } + + // push forward the j-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < dim; ++l) + for (unsigned int m = 0; m < dim; ++m) { - tmp[i][j][l] = - result[i][0][l] * data.covariant[point][j][0]; + jacobian_pushed_forward_2nd_derivatives + [point][i][j][l][m] = + result[i][0][l][m] * data.covariant[point][j][0]; for (unsigned int jr = 1; jr < dim; ++jr) - { - tmp[i][j][l] += - result[i][jr][l] * data.covariant[point][j][jr]; - } + jacobian_pushed_forward_2nd_derivatives[point][i][j] + [l][m] += + result[i][jr][l][m] * data.covariant[point][j][jr]; } - // now, pushing forward the l-components - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) + // push forward the l-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < spacedim; ++l) + for (unsigned int m = 0; m < dim; ++m) { - jacobian_pushed_forward_grads[point][i][j][l] = - tmp[i][j][0] * data.covariant[point][l][0]; + tmp[i][j][l][m] = + jacobian_pushed_forward_2nd_derivatives[point][i][j] + [0][m] * + data.covariant[point][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][l][m] += + jacobian_pushed_forward_2nd_derivatives[point][i][j] + [lr][m] * + data.covariant[point][l][lr]; } - } - } - } - - /** - * Update the third derivative of the transformation from unit to real - * cell, the Jacobian hessians. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - void - maybe_update_jacobian_2nd_derivatives( - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField:: - InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector> &jacobian_2nd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_2nd_derivatives) - { - const unsigned int n_q_points = jacobian_2nd_derivatives.size(); - - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<3, dim> *third = - &data.third_derivative(point + data_set, 0); - DerivativeForm<3, dim, spacedim> result; + // push forward the m-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < spacedim; ++l) + for (unsigned int m = 0; m < spacedim; ++m) + { + jacobian_pushed_forward_2nd_derivatives + [point][i][j][l][m] = + tmp[i][j][l][0] * data.covariant[point][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]; + } + } + } + } - for (unsigned int k = 0; k < data.n_shape_functions; ++k) - { - const unsigned int comp_k = - fe.system_to_component_index(k).first; - if (fe_mask[comp_k]) - 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[fe_to_real[comp_k]][j][l][m] += - (third[k][j][l][m] * data.local_dof_values[k]); - } + /** + * Update the fourth derivative of the transformation from unit to real + * cell, the Jacobian hessian gradients. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + void + maybe_update_jacobian_3rd_derivatives( + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField:: + InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector> &jacobian_3rd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_3rd_derivatives) + { + const unsigned int n_q_points = jacobian_3rd_derivatives.size(); - // never touch any data for j=dim in case dim - void - maybe_update_jacobian_pushed_forward_2nd_derivatives( - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField:: - InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector> - &jacobian_pushed_forward_2nd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_2nd_derivatives) - { - const unsigned int n_q_points = - jacobian_pushed_forward_2nd_derivatives.size(); - - double tmp[spacedim][spacedim][spacedim][spacedim]; - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<3, dim> *third = - &data.third_derivative(point + data_set, 0); + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<4, dim> *fourth = + &data.fourth_derivative(point + data_set, 0); - DerivativeForm<3, dim, spacedim> result; + DerivativeForm<4, dim, spacedim> result; - for (unsigned int k = 0; k < data.n_shape_functions; ++k) - { - const unsigned int comp_k = - fe.system_to_component_index(k).first; - if (fe_mask[comp_k]) - 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[fe_to_real[comp_k]][j][l][m] += - (third[k][j][l][m] * data.local_dof_values[k]); - } + for (unsigned int k = 0; k < data.n_shape_functions; ++k) + { + const unsigned int comp_k = + fe.system_to_component_index(k).first; + if (fe_mask[comp_k]) + 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[fe_to_real[comp_k]][j][l][m][n] += + (fourth[k][j][l][m][n] * + data.local_dof_values[k]); + } + + // never touch any data for j,l,m,n=dim in case + // dim + void + maybe_update_jacobian_pushed_forward_3rd_derivatives( + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField:: + InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + std::vector> &jacobian_pushed_forward_3rd_derivatives) + { + const UpdateFlags update_flags = data.update_each; + if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) + { + const unsigned int n_q_points = + jacobian_pushed_forward_3rd_derivatives.size(); + + double tmp[spacedim][spacedim][spacedim][spacedim][spacedim]; + for (unsigned int point = 0; point < n_q_points; ++point) + { + const Tensor<4, dim> *fourth = + &data.fourth_derivative(point + data_set, 0); + + DerivativeForm<4, dim, spacedim> result; - // push forward the j-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < dim; ++l) - for (unsigned int m = 0; m < dim; ++m) + for (unsigned int k = 0; k < data.n_shape_functions; ++k) + { + const unsigned int comp_k = + fe.system_to_component_index(k).first; + if (fe_mask[comp_k]) + 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[fe_to_real[comp_k]][j][l][m][n] += + (fourth[k][j][l][m][n] * + data.local_dof_values[k]); + } + + // push-forward the j-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++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_pushed_forward_2nd_derivatives - [point][i][j][l][m] = - result[i][0][l][m] * data.covariant[point][j][0]; + tmp[i][j][l][m][n] = + result[i][0][l][m][n] * data.covariant[point][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][n] += result[i][jr][l][m][n] * + data.covariant[point][j][jr]; } - // push forward the l-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) - for (unsigned int m = 0; m < dim; ++m) + // push-forward the l-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < spacedim; ++l) + for (unsigned int m = 0; m < dim; ++m) + for (unsigned int n = 0; n < dim; ++n) { - tmp[i][j][l][m] = - jacobian_pushed_forward_2nd_derivatives[point][i][j] - [0][m] * - data.covariant[point][l][0]; + jacobian_pushed_forward_3rd_derivatives + [point][i][j][l][m][n] = + tmp[i][j][0][m][n] * data.covariant[point][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] * + jacobian_pushed_forward_3rd_derivatives[point][i][j] + [l][m][n] += + tmp[i][j][lr][m][n] * data.covariant[point][l][lr]; } - // push forward the m-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) - for (unsigned int m = 0; m < spacedim; ++m) + // push-forward the m-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < spacedim; ++l) + for (unsigned int m = 0; m < spacedim; ++m) + for (unsigned int n = 0; n < dim; ++n) { - jacobian_pushed_forward_2nd_derivatives - [point][i][j][l][m] = - tmp[i][j][l][0] * data.covariant[point][m][0]; + tmp[i][j][l][m][n] = + jacobian_pushed_forward_3rd_derivatives[point][i][j] + [l][0][n] * + data.covariant[point][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]; + tmp[i][j][l][m][n] += + jacobian_pushed_forward_3rd_derivatives[point][i] + [j][l][mr] + [n] * + data.covariant[point][m][mr]; } - } - } - } - - /** - * Update the fourth derivative of the transformation from unit to real - * cell, the Jacobian hessian gradients. - * - * Skip the computation if possible as indicated by the first argument. - */ - template - void - maybe_update_jacobian_3rd_derivatives( - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField:: - InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector> &jacobian_3rd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_3rd_derivatives) - { - const unsigned int n_q_points = jacobian_3rd_derivatives.size(); - - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<4, dim> *fourth = - &data.fourth_derivative(point + data_set, 0); - DerivativeForm<4, dim, spacedim> result; + // push-forward the n-coordinate + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int l = 0; l < spacedim; ++l) + for (unsigned int m = 0; m < spacedim; ++m) + for (unsigned int n = 0; n < spacedim; ++n) + { + jacobian_pushed_forward_3rd_derivatives + [point][i][j][l][m][n] = + tmp[i][j][l][m][0] * data.covariant[point][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]; + } + } + } + } - for (unsigned int k = 0; k < data.n_shape_functions; ++k) - { - const unsigned int comp_k = - fe.system_to_component_index(k).first; - if (fe_mask[comp_k]) - 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[fe_to_real[comp_k]][j][l][m][n] += - (fourth[k][j][l][m][n] * - data.local_dof_values[k]); - } - // never touch any data for j,l,m,n=dim in case - // dim - void - maybe_update_jacobian_pushed_forward_3rd_derivatives( - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField:: - InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - std::vector> - &jacobian_pushed_forward_3rd_derivatives) - { - const UpdateFlags update_flags = data.update_each; - if (update_flags & update_jacobian_pushed_forward_3rd_derivatives) - { - const unsigned int n_q_points = - jacobian_pushed_forward_3rd_derivatives.size(); - - double tmp[spacedim][spacedim][spacedim][spacedim][spacedim]; - for (unsigned int point = 0; point < n_q_points; ++point) - { - const Tensor<4, dim> *fourth = - &data.fourth_derivative(point + data_set, 0); + /** + * Depending on what information is called for in the update flags of the + * @p data object, compute the various pieces of information that is + * required by the fill_fe_face_values() and fill_fe_subface_values() + * functions. This function simply unifies the work that would be done by + * those two functions. + * + * The resulting data is put into the @p output_data argument. + */ + template + void + maybe_compute_face_data( + const dealii::Mapping &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField:: + InternalData &data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) + { + const UpdateFlags update_flags = data.update_each; - DerivativeForm<4, dim, spacedim> result; + if (update_flags & update_boundary_forms) + { + const unsigned int n_q_points = output_data.boundary_forms.size(); + if (update_flags & update_normal_vectors) + AssertDimension(output_data.normal_vectors.size(), n_q_points); + if (update_flags & update_JxW_values) + AssertDimension(output_data.JxW_values.size(), n_q_points); + + // map the unit tangentials to the real cell. checking for d!=dim-1 + // eliminates compiler warnings regarding unsigned int expressions < + // 0. + for (unsigned int d = 0; d != dim - 1; ++d) + { + Assert(face_no + cell->n_faces() * d < + data.unit_tangentials.size(), + ExcInternalError()); + Assert( + data.aux[d].size() <= + data.unit_tangentials[face_no + cell->n_faces() * d].size(), + ExcInternalError()); + + mapping.transform( + make_array_view( + data.unit_tangentials[face_no + cell->n_faces() * d]), + mapping_contravariant, + data, + make_array_view(data.aux[d])); + } - for (unsigned int k = 0; k < data.n_shape_functions; ++k) + // if dim==spacedim, we can use the unit tangentials to compute the + // boundary form by simply taking the cross product + if (dim == spacedim) + { + for (unsigned int i = 0; i < n_q_points; ++i) + switch (dim) { - const unsigned int comp_k = - fe.system_to_component_index(k).first; - if (fe_mask[comp_k]) - 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[fe_to_real[comp_k]][j][l][m][n] += - (fourth[k][j][l][m][n] * - data.local_dof_values[k]); + case 1: + // in 1d, we don't have access to any of the data.aux + // fields (because it has only dim-1 components), but we + // can still compute the boundary form by simply looking + // at the number of the face + output_data.boundary_forms[i][0] = + (face_no == 0 ? -1 : +1); + break; + case 2: + output_data.boundary_forms[i] = + cross_product_2d(data.aux[0][i]); + break; + case 3: + output_data.boundary_forms[i] = + cross_product_3d(data.aux[0][i], data.aux[1][i]); + break; + default: + DEAL_II_NOT_IMPLEMENTED(); } - - // push-forward the j-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++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) - { - tmp[i][j][l][m][n] = result[i][0][l][m][n] * - data.covariant[point][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]; - } - - // push-forward the l-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) - 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]; - 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]; - } - - // push-forward the m-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) - for (unsigned int m = 0; m < spacedim; ++m) - 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]; - 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]; - } - - // push-forward the n-coordinate - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int l = 0; l < spacedim; ++l) - for (unsigned int m = 0; m < spacedim; ++m) - for (unsigned int n = 0; n < spacedim; ++n) - { - jacobian_pushed_forward_3rd_derivatives - [point][i][j][l][m][n] = - tmp[i][j][l][m][0] * - data.covariant[point][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]; - } - } - } - } - - - /** - * Depending on what information is called for in the update flags of the - * @p data object, compute the various pieces of information that is - * required by the fill_fe_face_values() and fill_fe_subface_values() - * functions. This function simply unifies the work that would be done by - * those two functions. - * - * The resulting data is put into the @p output_data argument. - */ - template - void - maybe_compute_face_data( - const dealii::Mapping &mapping, - const typename dealii::Triangulation::cell_iterator - &cell, - const unsigned int face_no, - const unsigned int subface_no, - const typename QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField:: - InternalData &data, - internal::FEValuesImplementation::MappingRelatedData - &output_data) - { - const UpdateFlags update_flags = data.update_each; - - if (update_flags & update_boundary_forms) - { - const unsigned int n_q_points = output_data.boundary_forms.size(); - if (update_flags & update_normal_vectors) - AssertDimension(output_data.normal_vectors.size(), n_q_points); - if (update_flags & update_JxW_values) - AssertDimension(output_data.JxW_values.size(), n_q_points); - - // map the unit tangentials to the real cell. checking for d!=dim-1 - // eliminates compiler warnings regarding unsigned int expressions < - // 0. - for (unsigned int d = 0; d != dim - 1; ++d) - { - Assert(face_no + cell->n_faces() * d < - data.unit_tangentials.size(), - ExcInternalError()); - Assert( - data.aux[d].size() <= - data.unit_tangentials[face_no + cell->n_faces() * d].size(), - ExcInternalError()); - - mapping.transform( - make_array_view( - data.unit_tangentials[face_no + cell->n_faces() * d]), - mapping_contravariant, - data, - make_array_view(data.aux[d])); - } - - // if dim==spacedim, we can use the unit tangentials to compute the - // boundary form by simply taking the cross product - if (dim == spacedim) - { - for (unsigned int i = 0; i < n_q_points; ++i) - switch (dim) + } + else //(dim < spacedim) + { + // in the codim-one case, the boundary form results from the + // cross product of all the face tangential vectors and the cell + // normal vector + // + // to compute the cell normal, use the same method used in + // fill_fe_values for cells above + AssertDimension(data.contravariant.size(), n_q_points); + + for (unsigned int point = 0; point < n_q_points; ++point) + { + if (dim == 1) { - case 1: - // in 1d, we don't have access to any of the data.aux - // fields (because it has only dim-1 components), but we - // can still compute the boundary form by simply looking - // at the number of the face - output_data.boundary_forms[i][0] = - (face_no == 0 ? -1 : +1); - break; - case 2: - output_data.boundary_forms[i] = - cross_product_2d(data.aux[0][i]); - break; - case 3: - output_data.boundary_forms[i] = - cross_product_3d(data.aux[0][i], data.aux[1][i]); - break; - default: - DEAL_II_NOT_IMPLEMENTED(); + // J is a tangent vector + output_data.boundary_forms[point] = + data.contravariant[point].transpose()[0]; + output_data.boundary_forms[point] /= + (face_no == 0 ? -1. : +1.) * + output_data.boundary_forms[point].norm(); } - } - else //(dim < spacedim) - { - // in the codim-one case, the boundary form results from the - // cross product of all the face tangential vectors and the cell - // normal vector - // - // to compute the cell normal, use the same method used in - // fill_fe_values for cells above - AssertDimension(data.contravariant.size(), n_q_points); - - for (unsigned int point = 0; point < n_q_points; ++point) - { - if (dim == 1) - { - // J is a tangent vector - output_data.boundary_forms[point] = - data.contravariant[point].transpose()[0]; - output_data.boundary_forms[point] /= - (face_no == 0 ? -1. : +1.) * - output_data.boundary_forms[point].norm(); - } - - if (dim == 2) - { - const DerivativeForm<1, spacedim, dim> DX_t = - data.contravariant[point].transpose(); - - Tensor<1, spacedim> cell_normal = - cross_product_3d(DX_t[0], DX_t[1]); - cell_normal /= cell_normal.norm(); - - // then compute the face normal from the face tangent - // and the cell normal: - output_data.boundary_forms[point] = - cross_product_3d(data.aux[0][point], cell_normal); - } - } - } - if (update_flags & (update_normal_vectors | update_JxW_values)) - for (unsigned int i = 0; i < output_data.boundary_forms.size(); - ++i) - { - if (update_flags & update_JxW_values) + if (dim == 2) { - output_data.JxW_values[i] = - output_data.boundary_forms[i].norm() * - data.quadrature_weights[i + data_set]; + const DerivativeForm<1, spacedim, dim> DX_t = + data.contravariant[point].transpose(); - if (subface_no != numbers::invalid_unsigned_int) - { - // TODO - const double area_ratio = - GeometryInfo::subface_ratio( - cell->subface_case(face_no), subface_no); - output_data.JxW_values[i] *= area_ratio; - } - } + Tensor<1, spacedim> cell_normal = + cross_product_3d(DX_t[0], DX_t[1]); + cell_normal /= cell_normal.norm(); - if (update_flags & update_normal_vectors) - output_data.normal_vectors[i] = - Point(output_data.boundary_forms[i] / - output_data.boundary_forms[i].norm()); + // then compute the face normal from the face tangent + // and the cell normal: + output_data.boundary_forms[point] = + cross_product_3d(data.aux[0][point], cell_normal); + } } - } - } - - /** - * Do the work of MappingFEField::fill_fe_face_values() and - * MappingFEField::fill_fe_subface_values() in a generic way, using the - * 'data_set' to differentiate whether we will work on a face (and if so, - * which one) or subface. - */ - template - void - do_fill_fe_face_values( - const dealii::Mapping &mapping, - const typename dealii::Triangulation::cell_iterator - &cell, - const unsigned int face_no, - const unsigned int subface_no, - const typename dealii::QProjector::DataSetDescriptor data_set, - const typename dealii::MappingFEField:: - InternalData &data, - const FiniteElement &fe, - const ComponentMask &fe_mask, - const std::vector &fe_to_real, - internal::FEValuesImplementation::MappingRelatedData - &output_data) - { - maybe_compute_q_points( - data_set, - data, - fe, - fe_mask, - fe_to_real, - output_data.quadrature_points); - - maybe_update_Jacobians( - data_set, data, fe, fe_mask, fe_to_real); - - const UpdateFlags update_flags = data.update_each; - const unsigned int n_q_points = data.contravariant.size(); - - if (update_flags & update_jacobians) - for (unsigned int point = 0; point < n_q_points; ++point) - output_data.jacobians[point] = data.contravariant[point]; + } - if (update_flags & update_inverse_jacobians) - for (unsigned int point = 0; point < n_q_points; ++point) - output_data.inverse_jacobians[point] = - data.covariant[point].transpose(); + if (update_flags & (update_normal_vectors | update_JxW_values)) + for (unsigned int i = 0; i < output_data.boundary_forms.size(); ++i) + { + if (update_flags & update_JxW_values) + { + output_data.JxW_values[i] = + output_data.boundary_forms[i].norm() * + data.quadrature_weights[i + data_set]; - maybe_update_jacobian_grads( - data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads); + if (subface_no != numbers::invalid_unsigned_int) + { + // TODO + const double area_ratio = + GeometryInfo::subface_ratio( + cell->subface_case(face_no), subface_no); + output_data.JxW_values[i] *= area_ratio; + } + } - maybe_update_jacobian_pushed_forward_grads( - data_set, - data, - fe, - fe_mask, - fe_to_real, - output_data.jacobian_pushed_forward_grads); + if (update_flags & update_normal_vectors) + output_data.normal_vectors[i] = + Point(output_data.boundary_forms[i] / + output_data.boundary_forms[i].norm()); + } + } + } - maybe_update_jacobian_2nd_derivatives( - data_set, - data, - fe, - fe_mask, - fe_to_real, - output_data.jacobian_2nd_derivatives); + /** + * Do the work of MappingFEField::fill_fe_face_values() and + * MappingFEField::fill_fe_subface_values() in a generic way, using the + * 'data_set' to differentiate whether we will work on a face (and if so, + * which one) or subface. + */ + template + void + do_fill_fe_face_values( + const dealii::Mapping &mapping, + const typename dealii::Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const typename dealii::QProjector::DataSetDescriptor data_set, + const typename dealii::MappingFEField:: + InternalData &data, + const FiniteElement &fe, + const ComponentMask &fe_mask, + const std::vector &fe_to_real, + internal::FEValuesImplementation::MappingRelatedData + &output_data) + { + maybe_compute_q_points( + data_set, data, fe, fe_mask, fe_to_real, output_data.quadrature_points); - maybe_update_jacobian_pushed_forward_2nd_derivatives( - data_set, - data, - fe, - fe_mask, - fe_to_real, - output_data.jacobian_pushed_forward_2nd_derivatives); + maybe_update_Jacobians( + data_set, data, fe, fe_mask, fe_to_real); - maybe_update_jacobian_3rd_derivatives( - data_set, - data, - fe, - fe_mask, - fe_to_real, - output_data.jacobian_3rd_derivatives); + const UpdateFlags update_flags = data.update_each; + const unsigned int n_q_points = data.contravariant.size(); - maybe_update_jacobian_pushed_forward_3rd_derivatives( - data_set, - data, - fe, - fe_mask, - fe_to_real, - output_data.jacobian_pushed_forward_3rd_derivatives); + if (update_flags & update_jacobians) + for (unsigned int point = 0; point < n_q_points; ++point) + output_data.jacobians[point] = data.contravariant[point]; - maybe_compute_face_data( - mapping, cell, face_no, subface_no, data_set, data, output_data); - } - } // namespace - } // namespace MappingFEFieldImplementation + if (update_flags & update_inverse_jacobians) + for (unsigned int point = 0; point < n_q_points; ++point) + output_data.inverse_jacobians[point] = + data.covariant[point].transpose(); + + maybe_update_jacobian_grads( + data_set, data, fe, fe_mask, fe_to_real, output_data.jacobian_grads); + + maybe_update_jacobian_pushed_forward_grads( + data_set, + data, + fe, + fe_mask, + fe_to_real, + output_data.jacobian_pushed_forward_grads); + + maybe_update_jacobian_2nd_derivatives( + data_set, + data, + fe, + fe_mask, + fe_to_real, + output_data.jacobian_2nd_derivatives); + + maybe_update_jacobian_pushed_forward_2nd_derivatives( + data_set, + data, + fe, + fe_mask, + fe_to_real, + output_data.jacobian_pushed_forward_2nd_derivatives); + + maybe_update_jacobian_3rd_derivatives( + data_set, + data, + fe, + fe_mask, + fe_to_real, + output_data.jacobian_3rd_derivatives); + + maybe_update_jacobian_pushed_forward_3rd_derivatives( + data_set, + data, + fe, + fe_mask, + fe_to_real, + output_data.jacobian_pushed_forward_3rd_derivatives); + + maybe_compute_face_data( + mapping, cell, face_no, subface_no, data_set, data, output_data); + } + } // namespace MappingFEFieldImplementation } // namespace internal @@ -1957,125 +1935,117 @@ namespace internal { namespace MappingFEFieldImplementation { - namespace + template + void + transform_fields( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) { - template - void - transform_fields( - const ArrayView> &input, - const MappingKind mapping_kind, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView> &output) - { - AssertDimension(input.size(), output.size()); - Assert((dynamic_cast< - const typename dealii:: - MappingFEField::InternalData *>( - &mapping_data) != nullptr), - ExcInternalError()); - const typename dealii::MappingFEField:: - InternalData &data = static_cast< - const typename dealii::MappingFEField:: - InternalData &>(mapping_data); - - switch (mapping_kind) - { - case mapping_contravariant: - { - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_contravariant_transformation")); + AssertDimension(input.size(), output.size()); + Assert((dynamic_cast< + const typename dealii:: + MappingFEField::InternalData *>( + &mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingFEField:: + InternalData &data = static_cast< + const typename dealii::MappingFEField:: + InternalData &>(mapping_data); + + switch (mapping_kind) + { + case mapping_contravariant: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + output[i] = + apply_transformation(data.contravariant[i], input[i]); - for (unsigned int i = 0; i < output.size(); ++i) + return; + } + + case mapping_piola: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_volume_elements")); + Assert(rank == 1, ExcMessage("Only for rank 1")); + for (unsigned int i = 0; i < output.size(); ++i) + { output[i] = apply_transformation(data.contravariant[i], input[i]); + output[i] /= data.volume_elements[i]; + } + return; + } - return; - } - case mapping_piola: - { - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_contravariant_transformation")); - Assert( - data.update_each & update_volume_elements, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_volume_elements")); - Assert(rank == 1, ExcMessage("Only for rank 1")); - for (unsigned int i = 0; i < output.size(); ++i) - { - output[i] = - apply_transformation(data.contravariant[i], input[i]); - output[i] /= data.volume_elements[i]; - } - return; - } + // We still allow this operation as in the + // reference cell Derivatives are Tensor + // rather than DerivativeForm + case mapping_covariant: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + for (unsigned int i = 0; i < output.size(); ++i) + output[i] = apply_transformation(data.covariant[i], input[i]); - // We still allow this operation as in the - // reference cell Derivatives are Tensor - // rather than DerivativeForm - case mapping_covariant: - { - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_contravariant_transformation")); + return; + } - for (unsigned int i = 0; i < output.size(); ++i) - output[i] = apply_transformation(data.covariant[i], input[i]); + default: + DEAL_II_NOT_IMPLEMENTED(); + } + } - return; - } - default: - DEAL_II_NOT_IMPLEMENTED(); - } - } - - - template - void - transform_differential_forms( - const ArrayView> &input, - const MappingKind mapping_kind, - const typename Mapping::InternalDataBase &mapping_data, - const ArrayView> &output) - { - AssertDimension(input.size(), output.size()); - Assert((dynamic_cast< - const typename dealii:: - MappingFEField::InternalData *>( - &mapping_data) != nullptr), - ExcInternalError()); - const typename dealii::MappingFEField:: - InternalData &data = static_cast< - const typename dealii::MappingFEField:: - InternalData &>(mapping_data); - - switch (mapping_kind) - { - case mapping_covariant: - { - Assert( - data.update_each & update_contravariant_transformation, - typename FEValuesBase::ExcAccessToUninitializedField( - "update_contravariant_transformation")); + template + void + transform_differential_forms( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) + { + AssertDimension(input.size(), output.size()); + Assert((dynamic_cast< + const typename dealii:: + MappingFEField::InternalData *>( + &mapping_data) != nullptr), + ExcInternalError()); + const typename dealii::MappingFEField:: + InternalData &data = static_cast< + const typename dealii::MappingFEField:: + InternalData &>(mapping_data); + + switch (mapping_kind) + { + case mapping_covariant: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); - for (unsigned int i = 0; i < output.size(); ++i) - output[i] = apply_transformation(data.covariant[i], input[i]); + for (unsigned int i = 0; i < output.size(); ++i) + output[i] = apply_transformation(data.covariant[i], input[i]); - return; - } - default: - DEAL_II_NOT_IMPLEMENTED(); - } - } - } // namespace - } // namespace MappingFEFieldImplementation + return; + } + default: + DEAL_II_NOT_IMPLEMENTED(); + } + } + } // namespace MappingFEFieldImplementation } // namespace internal -- 2.39.5