From d7cbec4bfc5b356a3d19f218894734654064db15 Mon Sep 17 00:00:00 2001 From: David Wells Date: Tue, 17 Jun 2025 11:40:11 -0400 Subject: [PATCH] Mapping: refactor the contravariant Hessian mappings. --- include/deal.II/fe/mapping_internal.h | 53 +++++++++++++++++++++++++ include/deal.II/fe/mapping_q_internal.h | 35 +++------------- source/fe/mapping_fe.cc | 36 +++-------------- source/fe/mapping_manifold.cc | 36 +++-------------- 4 files changed, 68 insertions(+), 92 deletions(-) diff --git a/include/deal.II/fe/mapping_internal.h b/include/deal.II/fe/mapping_internal.h index 2b4bf35cff..7fd4592795 100644 --- a/include/deal.II/fe/mapping_internal.h +++ b/include/deal.II/fe/mapping_internal.h @@ -27,6 +27,18 @@ DEAL_II_NAMESPACE_OPEN namespace internal { + /** + * Map the Hessian of a contravariant 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_contravariant_hessian( + const DerivativeForm<1, dim, spacedim, Number> &covariant, + const DerivativeForm<1, dim, spacedim, Number> &contravariant, + const Tensor<3, dim, Number> &input); + /** * Map the Hessian of a covariant vector field. For more information see the * overload of Mapping::transform() which maps 3-differential forms from the @@ -41,6 +53,47 @@ namespace internal namespace internal { + template + inline Tensor<3, spacedim, Number> + apply_contravariant_hessian( + const DerivativeForm<1, dim, spacedim, Number> &covariant, + const DerivativeForm<1, dim, spacedim, Number> &contravariant, + const Tensor<3, dim, Number> &input) + { + Tensor<3, spacedim, Number> output; + for (unsigned int i = 0; i < spacedim; ++i) + { + Number tmp1[dim][dim]; + for (unsigned int J = 0; J < dim; ++J) + for (unsigned int K = 0; K < dim; ++K) + { + tmp1[J][K] = contravariant[i][0] * input[0][J][K]; + for (unsigned int I = 1; I < dim; ++I) + tmp1[J][K] += contravariant[i][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; + } + + + template inline Tensor<3, spacedim, Number> apply_covariant_hessian( diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index 518aad29e8..1d46c444b8 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -2364,37 +2364,12 @@ 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 tmp1[dim][dim]; - for (unsigned int J = 0; J < dim; ++J) - for (unsigned int K = 0; K < dim; ++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 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_contravariant_hessian(covariant, + contravariant, + input[q]); } + return; } diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index f335034851..2b0fb25f5d 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -1901,37 +1901,11 @@ 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]; - 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_contravariant_hessian(data.covariant[q], + data.contravariant[q], + input[q]); + return; } diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 92055faecf..20df91ffd7 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -1060,37 +1060,11 @@ 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]; - 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_contravariant_hessian(data.covariant[q], + data.contravariant[q], + input[q]); + return; } -- 2.39.5