From 4b19d7fbf4fadc082a3fb55f2ae8e8269a3a1222 Mon Sep 17 00:00:00 2001 From: David Wells Date: Tue, 17 Jun 2025 11:19:43 -0400 Subject: [PATCH] Mapping: refactor the Piola Hessian mappings. --- include/deal.II/fe/mapping_internal.h | 58 +++++++++++++++++++++++++ include/deal.II/fe/mapping_q_internal.h | 36 +++------------ source/fe/mapping_fe.cc | 38 +++------------- source/fe/mapping_manifold.cc | 38 +++------------- 4 files changed, 73 insertions(+), 97 deletions(-) diff --git a/include/deal.II/fe/mapping_internal.h b/include/deal.II/fe/mapping_internal.h index 7fd4592795..8a88e7f8c8 100644 --- a/include/deal.II/fe/mapping_internal.h +++ b/include/deal.II/fe/mapping_internal.h @@ -49,6 +49,19 @@ namespace internal apply_covariant_hessian( const DerivativeForm<1, dim, spacedim, Number> &covariant, const Tensor<3, dim, Number> &input); + + /** + * Map the Hessian of a Piola vector field. For more information see the + * overload of Mapping::transform() which maps 3-differential forms from the + * reference cell to the physical cell. + */ + template + Tensor<3, spacedim, Number> + apply_piola_hessian( + const DerivativeForm<1, dim, spacedim, Number> &covariant, + const DerivativeForm<1, dim, spacedim, Number> &contravariant, + const Number &volume_element, + const Tensor<3, dim, Number> &input); } // namespace internal namespace internal @@ -131,6 +144,51 @@ namespace internal return output; } + + + + template + inline Tensor<3, spacedim, Number> + apply_piola_hessian( + const DerivativeForm<1, dim, spacedim, Number> &covariant, + const DerivativeForm<1, dim, spacedim, Number> &contravariant, + const Number &volume_element, + const Tensor<3, dim, Number> &input) + { + Tensor<3, spacedim, Number> output; + for (unsigned int i = 0; i < spacedim; ++i) + { + Number factor[dim]; + for (unsigned int I = 0; I < dim; ++I) + factor[I] = contravariant[i][I] * (1. / volume_element); + Number 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[0][J][K]; + for (unsigned int I = 1; I < dim; ++I) + tmp1[J][K] += factor[I] * input[I][J][K]; + } + for (unsigned int j = 0; j < spacedim; ++j) + { + Number 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[i][j][k] = covariant[k][0] * tmp2[0]; + for (unsigned int K = 1; K < dim; ++K) + output[i][j][k] += covariant[k][K] * tmp2[K]; + } + } + } + + return output; + } } // namespace internal DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index 1d46c444b8..38e5fb02af 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -2408,37 +2408,11 @@ namespace internal 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) - { - 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]; - 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]; - } - } - } + output[q] = + internal::apply_piola_hessian(covariant, + contravariant, + data.volume_elements[q], + input[q]); } return; diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index 2b0fb25f5d..367c9b287a 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -1940,39 +1940,11 @@ 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]; - 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]; - } - } - } + output[q] = + internal::apply_piola_hessian(data.covariant[q], + data.contravariant[q], + data.volume_elements[q], + input[q]); return; } diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 20df91ffd7..26298d42d5 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -1099,39 +1099,11 @@ 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]; - 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]; - } - } - } + output[q] = + internal::apply_piola_hessian(data.covariant[q], + data.contravariant[q], + data.volume_elements[q], + input[q]); return; } -- 2.39.5