From: Maximilian Bergbauer Date: Tue, 26 Jul 2022 08:08:20 +0000 (+0200) Subject: Fix Tensor assignment of nested Tensors and provide a function overload of apply_tran... X-Git-Tag: v9.5.0-rc1~1028^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=001d126a99f0068e411c52eeca79d0b7f8d1e6f3;p=dealii.git Fix Tensor assignment of nested Tensors and provide a function overload of apply_transformation to get rid of casts in FEPointEvaluation --- diff --git a/include/deal.II/base/derivative_form.h b/include/deal.II/base/derivative_form.h index bb4fa64e3c..99ffd053a5 100644 --- a/include/deal.II/base/derivative_form.h +++ b/include/deal.II/base/derivative_form.h @@ -108,14 +108,12 @@ public: * Number>. In particular, if order == 1 and the derivative is the Jacobian of * $\mathbf F(\mathbf x)$, then Tensor[i] = $\nabla F_i(\mathbf x)$. */ - template - operator Tensor() const; + operator Tensor() const; /** * Converts a DerivativeForm<1, dim, 1, Number> to Tensor<1, dim, Number>. */ - template - operator Tensor<1, dim, Number2>() const; + operator Tensor<1, dim, Number>() const; /** * Return the transpose of a rectangular DerivativeForm, @@ -275,9 +273,8 @@ DerivativeForm::operator[]( template -template inline DerivativeForm:: -operator Tensor<1, dim, Number2>() const +operator Tensor<1, dim, Number>() const { Assert((1 == spacedim) && (order == 1), ExcMessage("Only allowed for spacedim==1.")); @@ -288,13 +285,12 @@ operator Tensor<1, dim, Number2>() const template -template inline DerivativeForm:: -operator Tensor() const +operator Tensor() const { Assert((dim == spacedim), ExcMessage("Only allowed when dim==spacedim.")); - Tensor t; + Tensor t; if (dim == spacedim) for (unsigned int j = 0; j < dim; ++j) @@ -470,6 +466,31 @@ apply_transformation(const DerivativeForm<1, dim, spacedim, Number1> &grad_F, +/** + * Similar to the previous apply_transformation(), specialized for the case `dim + * == spacedim` where we can return a rank-2 tensor instead of the more general + * `DerivativeForm`. + * Each row of the result corresponds to one of the rows of @p D_X transformed + * by @p grad_F, equivalent to $\mathrm{D\_X} \, \mathrm{grad\_F}^T$ in matrix + * notation. + * + * @relatesalso DerivativeForm + */ +// rank=2 +template +inline Tensor<2, dim, typename ProductType::type> +apply_transformation(const DerivativeForm<1, dim, dim, Number1> &grad_F, + const Tensor<2, dim, Number2> & D_X) +{ + Tensor<2, dim, typename ProductType::type> dest; + for (unsigned int i = 0; i < dim; ++i) + dest[i] = apply_transformation(grad_F, D_X[i]); + + return dest; +} + + + /** * Similar to the previous apply_transformation(). * Each row of the result corresponds to one of the rows of @p D_X transformed diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index c9798ba383..d14ab9bf5d 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -1069,7 +1069,7 @@ constexpr inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> & Tensor<0, dim, Number>::operator=(const Tensor<0, dim, OtherNumber> &p) { - value = internal::NumberType::value(p); + value = internal::NumberType::value(p.value); return *this; } @@ -1361,7 +1361,7 @@ template constexpr DEAL_II_ALWAYS_INLINE Tensor:: operator Tensor<1, dim, Tensor>() const { - return Tensor<1, dim, Tensor>(values); + return Tensor<1, dim, Tensor>(values); } diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index 4e64cfc2df..21969e79c1 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -947,11 +947,11 @@ FEPointEvaluation::evaluate( Number>::set_gradient(val_and_grad.second, j, unit_gradients[i + j]); - gradients[i + j] = static_cast( + gradients[i + j] = apply_transformation(mapping_info->get_mapping_data() .inverse_jacobians[i + j] .transpose(), - unit_gradients[i + j])); + unit_gradients[i + j]); } } } @@ -1087,10 +1087,9 @@ FEPointEvaluation::integrate( if (integration_flags & EvaluationFlags::gradients) for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j) { - gradients[i + j] = - static_cast(apply_transformation( - mapping_info->get_mapping_data().inverse_jacobians[i + j], - gradients[i + j])); + gradients[i + j] = apply_transformation( + mapping_info->get_mapping_data().inverse_jacobians[i + j], + gradients[i + j]); internal::FEPointEvaluation:: EvaluatorTypeTraits::get_gradient( gradient, j, gradients[i + j]); diff --git a/tests/tensors/tensor_mixing_types_03.cc b/tests/tensors/tensor_mixing_types_03.cc new file mode 100644 index 0000000000..69ce6535f1 --- /dev/null +++ b/tests/tensors/tensor_mixing_types_03.cc @@ -0,0 +1,40 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Test for mixed Number type assignments and casts of nested Tensors. + +#include + +#include + +#include "../tests.h" + + +int +main() +{ + initlog(); + + Tensor<1, 1, Tensor<1, 2, float>> nested_f; + Tensor<1, 1, Tensor<1, 2, double>> nested_d; + + nested_d = nested_f; + nested_f = nested_d; + + nested_d = static_cast>>(nested_f); + nested_f = static_cast>>(nested_d); + + deallog << "OK!" << std::endl; +} diff --git a/tests/tensors/tensor_mixing_types_03.output b/tests/tensors/tensor_mixing_types_03.output new file mode 100644 index 0000000000..5cfb783b8f --- /dev/null +++ b/tests/tensors/tensor_mixing_types_03.output @@ -0,0 +1,2 @@ + +DEAL::OK!