From: David Wells Date: Tue, 17 Jun 2025 14:55:57 +0000 (-0400) Subject: Mapping: refactor the covariant Hessian mappings. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=735825a7d889d0977b19318642d056a7eab514a4;p=dealii.git Mapping: refactor the covariant Hessian mappings. --- diff --git a/include/deal.II/fe/mapping_internal.h b/include/deal.II/fe/mapping_internal.h new file mode 100644 index 0000000000..2b4bf35cff --- /dev/null +++ b/include/deal.II/fe/mapping_internal.h @@ -0,0 +1,85 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2025 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + +#ifndef dealii_mapping_internal_h +#define dealii_mapping_internal_h + +#include + +#include +#include + +// Implementations of transformations used by several Mapping classes (such as +// MappingFE, MappingQ, and MappingFEField, and MappingManifold) + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + /** + * Map the Hessian of a covariant 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_covariant_hessian( + const DerivativeForm<1, dim, spacedim, Number> &covariant, + const Tensor<3, dim, Number> &input); +} // namespace internal + +namespace internal +{ + template + inline Tensor<3, spacedim, Number> + apply_covariant_hessian( + const DerivativeForm<1, dim, spacedim, Number> &covariant, + 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] = covariant[i][0] * input[0][J][K]; + for (unsigned int I = 1; I < dim; ++I) + tmp1[J][K] += covariant[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; + } +} // namespace internal + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h index db35dba380..518aad29e8 100644 --- a/include/deal.II/fe/mapping_q_internal.h +++ b/include/deal.II/fe/mapping_q_internal.h @@ -28,6 +28,7 @@ #include #include #include +#include #include #include @@ -2407,34 +2408,8 @@ namespace internal { 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) - { - 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]; - } - } - } + output[q] = + internal::apply_covariant_hessian(covariant, input[q]); } return; diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index cbb2fa2b65..f335034851 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -1942,37 +1943,9 @@ 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]; - 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_covariant_hessian(data.covariant[q], + input[q]); return; } diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 42af2fbd59..92055faecf 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -26,6 +26,7 @@ #include #include #include +#include #include #include @@ -1101,37 +1102,9 @@ 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]; - 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_covariant_hessian(data.covariant[q], + input[q]); return; }