From 53437d7dd1bb4f6877a19a219b6ad2f33c0652f3 Mon Sep 17 00:00:00 2001 From: David Wells Date: Thu, 19 Jun 2025 10:30:02 -0400 Subject: [PATCH] Mappings: refactor the covariant gradient. --- .../deal.II/fe/mapping_fe_field.templates.h | 22 +++------- include/deal.II/fe/mapping_internal.h | 41 +++++++++++++++++++ source/fe/mapping_fe.cc | 22 ++-------- source/fe/mapping_manifold.cc | 22 ++-------- source/fe/mapping_q.cc | 25 +++-------- 5 files changed, 60 insertions(+), 72 deletions(-) diff --git a/include/deal.II/fe/mapping_fe_field.templates.h b/include/deal.II/fe/mapping_fe_field.templates.h index c132cced57..202aba571c 100644 --- a/include/deal.II/fe/mapping_fe_field.templates.h +++ b/include/deal.II/fe/mapping_fe_field.templates.h @@ -32,6 +32,7 @@ #include #include #include +#include #include @@ -2123,27 +2124,14 @@ MappingFEField::transform( { case mapping_covariant_gradient: { - Assert(data.update_each & update_contravariant_transformation, + Assert(data.update_each & update_covariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField( "update_covariant_transformation")); for (unsigned int q = 0; q < output.size(); ++q) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - for (unsigned int k = 0; k < spacedim; ++k) - { - output[q][i][j][k] = data.covariant[q][j][0] * - data.covariant[q][k][0] * - input[q][i][0][0]; - for (unsigned int J = 0; J < dim; ++J) - { - const unsigned int K0 = (0 == J) ? 1 : 0; - for (unsigned int K = K0; K < dim; ++K) - output[q][i][j][k] += data.covariant[q][j][J] * - data.covariant[q][k][K] * - input[q][i][J][K]; - } - } + output[q] = + internal::apply_covariant_gradient(data.covariant[q], input[q]); + return; } diff --git a/include/deal.II/fe/mapping_internal.h b/include/deal.II/fe/mapping_internal.h index 8a88e7f8c8..f234bd02de 100644 --- a/include/deal.II/fe/mapping_internal.h +++ b/include/deal.II/fe/mapping_internal.h @@ -27,6 +27,17 @@ DEAL_II_NAMESPACE_OPEN namespace internal { + /** + * Map the gradient of a covariant vector field. For more information see the + * overload of Mapping::transform() which maps 2-differential forms from the + * reference cell to the physical cell. + */ + template + Tensor<3, spacedim, Number> + apply_covariant_gradient( + const DerivativeForm<1, dim, spacedim, Number> &covariant, + const DerivativeForm<2, dim, spacedim, Number> &input); + /** * Map the Hessian of a contravariant vector field. For more information see * the overload of Mapping::transform() which maps 3-differential forms from @@ -66,6 +77,36 @@ namespace internal namespace internal { + template + Tensor<3, spacedim, Number> + apply_covariant_gradient( + const DerivativeForm<1, dim, spacedim, Number> &covariant, + const DerivativeForm<2, dim, spacedim, Number> &input) + { + Tensor<3, spacedim, Number> output; + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + { + double tmp[dim]; + for (unsigned int K = 0; K < dim; ++K) + { + tmp[K] = covariant[j][0] * input[i][0][K]; + for (unsigned int J = 1; J < dim; ++J) + tmp[K] += covariant[j][J] * input[i][J][K]; + } + for (unsigned int k = 0; k < spacedim; ++k) + { + output[i][j][k] = covariant[k][0] * tmp[0]; + for (unsigned int K = 1; K < dim; ++K) + output[i][j][k] += covariant[k][K] * tmp[K]; + } + } + + return output; + } + + + template inline Tensor<3, spacedim, Number> apply_contravariant_hessian( diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index 367c9b287a..85108cb2d8 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -2080,28 +2080,14 @@ MappingFE::transform( { case mapping_covariant_gradient: { - Assert(data.update_each & update_contravariant_transformation, + Assert(data.update_each & update_covariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField( "update_covariant_transformation")); for (unsigned int q = 0; q < output.size(); ++q) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - { - double tmp[dim]; - for (unsigned int K = 0; K < dim; ++K) - { - tmp[K] = data.covariant[q][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]; - } - for (unsigned int k = 0; k < spacedim; ++k) - { - output[q][i][j][k] = data.covariant[q][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] = + internal::apply_covariant_gradient(data.covariant[q], input[q]); + return; } diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 26298d42d5..54b6672b6e 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -1303,28 +1303,14 @@ MappingManifold::transform( { case mapping_covariant_gradient: { - Assert(data.update_each & update_contravariant_transformation, + Assert(data.update_each & update_covariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField( "update_covariant_transformation")); for (unsigned int q = 0; q < output.size(); ++q) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - { - double tmp[dim]; - for (unsigned int K = 0; K < dim; ++K) - { - tmp[K] = data.covariant[q][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]; - } - for (unsigned int k = 0; k < spacedim; ++k) - { - output[q][i][j][k] = data.covariant[q][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] = + internal::apply_covariant_gradient(data.covariant[q], input[q]); + return; } diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 5b4d7faa13..78be82a085 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -1497,25 +1497,12 @@ MappingQ::transform( "update_covariant_transformation")); for (unsigned int q = 0; q < output.size(); ++q) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - { - double tmp[dim]; - const DerivativeForm<1, dim, spacedim> covariant = - data.inverse_jacobians[q].transpose(); - for (unsigned int K = 0; K < dim; ++K) - { - tmp[K] = covariant[j][0] * input[q][i][0][K]; - for (unsigned int J = 1; J < dim; ++J) - tmp[K] += covariant[j][J] * input[q][i][J][K]; - } - for (unsigned int k = 0; k < spacedim; ++k) - { - output[q][i][j][k] = covariant[k][0] * tmp[0]; - for (unsigned int K = 1; K < dim; ++K) - output[q][i][j][k] += covariant[k][K] * tmp[K]; - } - } + { + const DerivativeForm<1, dim, spacedim> covariant = + data.inverse_jacobians[q].transpose(); + output[q] = + internal::apply_covariant_gradient(covariant, input[q]); + } return; } -- 2.39.5