From 71206dc35d3c2801fa4d06c6257432640cc25492 Mon Sep 17 00:00:00 2001 From: Vinayak Vijay Date: Wed, 31 Jul 2024 19:34:54 +0530 Subject: [PATCH] symmetrize fourth order tensors --- doc/news/changes/minor/20240731Vinayak | 4 + include/deal.II/base/symmetric_tensor.h | 65 +++++++++ tests/tensors/symmetric_tensor_47.cc | 95 +++++++++++++ tests/tensors/symmetric_tensor_47.output | 166 +++++++++++++++++++++++ 4 files changed, 330 insertions(+) create mode 100644 doc/news/changes/minor/20240731Vinayak create mode 100644 tests/tensors/symmetric_tensor_47.cc create mode 100644 tests/tensors/symmetric_tensor_47.output diff --git a/doc/news/changes/minor/20240731Vinayak b/doc/news/changes/minor/20240731Vinayak new file mode 100644 index 0000000000..d3be688fad --- /dev/null +++ b/doc/news/changes/minor/20240731Vinayak @@ -0,0 +1,4 @@ +New: Symmetrize fourth order tensors based +on required symmetry type - major or minor +
+(Vinayak, 2024/07/31) diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index 3fcccf0132..15545ccaf5 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -3546,6 +3546,71 @@ DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE +/** + * Return the symmetrized version of a full rank-4 tensor, i.e. + * as a symmetric rank-4 tensor. The symmetry could be of the following types: + * Only minor: $A_{ijkl}=A_{jikl}=A_{ijlk}=A_{jilk}$. + * Both minor and major: $A_{ijkl}=A_{jikl}=A_{ijlk}=A_{jilk}$ and + * $A_{ijkl}=A_{klij}$. This is the version for general dimensions. + * @param t The tensor to be symmetrized. + * @param major_symmetry This argument decides the presence of major symmetry. + * + * @relatesalso SymmetricTensor + */ +template +DEAL_II_HOST constexpr inline DEAL_II_ALWAYS_INLINE + SymmetricTensor<4, dim, Number> + symmetrize(const Tensor<4, dim, Number> &t, const bool major_symmetry) +{ + SymmetricTensor<4, dim, Number> result; + + const Number half = internal::NumberType::value(0.5); + + // minor symmetry - A_{ijkl}=A_{jikl}=A_{ijlk}=A_{jilk} + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) + { + if (i != j && k == l) + { + // A_{ijkk}=A_{jikk} + result[i][j][k][k] = (t[i][j][k][k] + t[j][i][k][k]) * half; + } + else if (i == j && k != l) + { + // A_{iikl}=A_{iilk} + result[i][i][k][l] = (t[i][i][k][l] + t[i][i][l][k]) * half; + } + else if (i != j && k != l) + { + // A_{ijkl}=A_{jilk} + result[i][j][k][l] = (t[i][j][k][l] + t[j][i][k][l] + + t[i][j][l][k] + t[j][i][l][k]) * + half * half; + } + else + { + // A_{iijj} and A_{iiii} unchanged + result[i][j][k][l] = t[i][j][k][l]; + } + } + + // in case major symmetry is also required + if (major_symmetry) + { + // major symmetry - A_{ijkl}=A_{klij} + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = i; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = k; l < dim; ++l) + result[i][j][k][l] = (t[i][j][k][l] + t[k][l][i][j]) * half; + } + return result; +} + + + /** * Multiplication of a symmetric tensor of general rank with a scalar from the * right. This version of the operator is used if the scalar has the same data diff --git a/tests/tensors/symmetric_tensor_47.cc b/tests/tensors/symmetric_tensor_47.cc new file mode 100644 index 0000000000..3cf5b901b6 --- /dev/null +++ b/tests/tensors/symmetric_tensor_47.cc @@ -0,0 +1,95 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2023 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. +// +// ------------------------------------------------------------------------ + + +// check SymmetricTensor::symmetrize on fourth order Tensors + + +#include +#include +#include + +#include "../tests.h" + + +template +void +make_fourth_order_tensor(Tensor<4, dim, Number> &t) +{ + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) + t[i][j][k][l] = (i + 1) * (j + 2) * (k + 3) * (l + 4); +} + +template +void +print_fourth_order_tensor(const SymmetricTensor<4, dim, Number> &t) +{ + const Number tolerance = 100. * std::numeric_limits::epsilon(); + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) + deallog << i << ' ' << j << ' ' << k << ' ' << l << ' ' + << filter_out_small_numbers(t[i][j][k][l], tolerance) + << std::endl; +} + + +int +main() +{ + initlog(); + const unsigned int dim = 3; + Tensor<4, dim, double> t; + make_fourth_order_tensor(t); + + SymmetricTensor<4, dim, double> st_minor; + + st_minor = symmetrize(t, false); + + // test minor symmetry + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) + AssertThrow(st_minor[i][j][k][l] == st_minor[i][j][l][k] && + st_minor[j][i][k][l] == st_minor[i][j][k][l] && + st_minor[i][j][k][l] == st_minor[j][i][l][k], + ExcInternalError()); + + print_fourth_order_tensor(st_minor); + + deallog << "MINOR SYMMETRY OK" << std::endl; + + + SymmetricTensor<4, dim, double> st_major; + + st_major = symmetrize(t, true); + + // test major symmetry + for (unsigned int i = 0; i < dim; ++i) + for (unsigned int j = 0; j < dim; ++j) + for (unsigned int k = 0; k < dim; ++k) + for (unsigned int l = 0; l < dim; ++l) + AssertThrow(st_major[i][j][k][l] == st_major[k][l][i][j], + ExcInternalError()); + + print_fourth_order_tensor(st_major); + + deallog << "MAJOR SYMMETRY OK" << std::endl; + deallog << "OK" << std::endl; +} diff --git a/tests/tensors/symmetric_tensor_47.output b/tests/tensors/symmetric_tensor_47.output new file mode 100644 index 0000000000..360f4ea1ac --- /dev/null +++ b/tests/tensors/symmetric_tensor_47.output @@ -0,0 +1,166 @@ + +DEAL::0 0 0 0 24.0000 +DEAL::0 0 0 1 31.0000 +DEAL::0 0 0 2 38.0000 +DEAL::0 0 1 0 31.0000 +DEAL::0 0 1 1 40.0000 +DEAL::0 0 1 2 49.0000 +DEAL::0 0 2 0 38.0000 +DEAL::0 0 2 1 49.0000 +DEAL::0 0 2 2 60.0000 +DEAL::0 1 0 0 42.0000 +DEAL::0 1 0 1 54.2500 +DEAL::0 1 0 2 66.5000 +DEAL::0 1 1 0 54.2500 +DEAL::0 1 1 1 70.0000 +DEAL::0 1 1 2 85.7500 +DEAL::0 1 2 0 66.5000 +DEAL::0 1 2 1 85.7500 +DEAL::0 1 2 2 105.000 +DEAL::0 2 0 0 60.0000 +DEAL::0 2 0 1 77.5000 +DEAL::0 2 0 2 95.0000 +DEAL::0 2 1 0 77.5000 +DEAL::0 2 1 1 100.000 +DEAL::0 2 1 2 122.500 +DEAL::0 2 2 0 95.0000 +DEAL::0 2 2 1 122.500 +DEAL::0 2 2 2 150.000 +DEAL::1 0 0 0 42.0000 +DEAL::1 0 0 1 54.2500 +DEAL::1 0 0 2 66.5000 +DEAL::1 0 1 0 54.2500 +DEAL::1 0 1 1 70.0000 +DEAL::1 0 1 2 85.7500 +DEAL::1 0 2 0 66.5000 +DEAL::1 0 2 1 85.7500 +DEAL::1 0 2 2 105.000 +DEAL::1 1 0 0 72.0000 +DEAL::1 1 0 1 93.0000 +DEAL::1 1 0 2 114.000 +DEAL::1 1 1 0 93.0000 +DEAL::1 1 1 1 120.000 +DEAL::1 1 1 2 147.000 +DEAL::1 1 2 0 114.000 +DEAL::1 1 2 1 147.000 +DEAL::1 1 2 2 180.000 +DEAL::1 2 0 0 102.000 +DEAL::1 2 0 1 131.750 +DEAL::1 2 0 2 161.500 +DEAL::1 2 1 0 131.750 +DEAL::1 2 1 1 170.000 +DEAL::1 2 1 2 208.250 +DEAL::1 2 2 0 161.500 +DEAL::1 2 2 1 208.250 +DEAL::1 2 2 2 255.000 +DEAL::2 0 0 0 60.0000 +DEAL::2 0 0 1 77.5000 +DEAL::2 0 0 2 95.0000 +DEAL::2 0 1 0 77.5000 +DEAL::2 0 1 1 100.000 +DEAL::2 0 1 2 122.500 +DEAL::2 0 2 0 95.0000 +DEAL::2 0 2 1 122.500 +DEAL::2 0 2 2 150.000 +DEAL::2 1 0 0 102.000 +DEAL::2 1 0 1 131.750 +DEAL::2 1 0 2 161.500 +DEAL::2 1 1 0 131.750 +DEAL::2 1 1 1 170.000 +DEAL::2 1 1 2 208.250 +DEAL::2 1 2 0 161.500 +DEAL::2 1 2 1 208.250 +DEAL::2 1 2 2 255.000 +DEAL::2 2 0 0 144.000 +DEAL::2 2 0 1 186.000 +DEAL::2 2 0 2 228.000 +DEAL::2 2 1 0 186.000 +DEAL::2 2 1 1 240.000 +DEAL::2 2 1 2 294.000 +DEAL::2 2 2 0 228.000 +DEAL::2 2 2 1 294.000 +DEAL::2 2 2 2 360.000 +DEAL::MINOR SYMMETRY OK +DEAL::0 0 0 0 24.0000 +DEAL::0 0 0 1 33.0000 +DEAL::0 0 0 2 42.0000 +DEAL::0 0 1 0 33.0000 +DEAL::0 0 1 1 56.0000 +DEAL::0 0 1 2 72.0000 +DEAL::0 0 2 0 42.0000 +DEAL::0 0 2 1 72.0000 +DEAL::0 0 2 2 102.000 +DEAL::0 1 0 0 33.0000 +DEAL::0 1 0 1 45.0000 +DEAL::0 1 0 2 57.0000 +DEAL::0 1 1 0 45.0000 +DEAL::0 1 1 1 75.0000 +DEAL::0 1 1 2 96.0000 +DEAL::0 1 2 0 57.0000 +DEAL::0 1 2 1 96.0000 +DEAL::0 1 2 2 135.000 +DEAL::0 2 0 0 42.0000 +DEAL::0 2 0 1 57.0000 +DEAL::0 2 0 2 72.0000 +DEAL::0 2 1 0 57.0000 +DEAL::0 2 1 1 94.0000 +DEAL::0 2 1 2 120.000 +DEAL::0 2 2 0 72.0000 +DEAL::0 2 2 1 120.000 +DEAL::0 2 2 2 168.000 +DEAL::1 0 0 0 33.0000 +DEAL::1 0 0 1 45.0000 +DEAL::1 0 0 2 57.0000 +DEAL::1 0 1 0 45.0000 +DEAL::1 0 1 1 75.0000 +DEAL::1 0 1 2 96.0000 +DEAL::1 0 2 0 57.0000 +DEAL::1 0 2 1 96.0000 +DEAL::1 0 2 2 135.000 +DEAL::1 1 0 0 56.0000 +DEAL::1 1 0 1 75.0000 +DEAL::1 1 0 2 94.0000 +DEAL::1 1 1 0 75.0000 +DEAL::1 1 1 1 120.000 +DEAL::1 1 1 2 152.000 +DEAL::1 1 2 0 94.0000 +DEAL::1 1 2 1 152.000 +DEAL::1 1 2 2 210.000 +DEAL::1 2 0 0 72.0000 +DEAL::1 2 0 1 96.0000 +DEAL::1 2 0 2 120.000 +DEAL::1 2 1 0 96.0000 +DEAL::1 2 1 1 152.000 +DEAL::1 2 1 2 192.000 +DEAL::1 2 2 0 120.000 +DEAL::1 2 2 1 192.000 +DEAL::1 2 2 2 264.000 +DEAL::2 0 0 0 42.0000 +DEAL::2 0 0 1 57.0000 +DEAL::2 0 0 2 72.0000 +DEAL::2 0 1 0 57.0000 +DEAL::2 0 1 1 94.0000 +DEAL::2 0 1 2 120.000 +DEAL::2 0 2 0 72.0000 +DEAL::2 0 2 1 120.000 +DEAL::2 0 2 2 168.000 +DEAL::2 1 0 0 72.0000 +DEAL::2 1 0 1 96.0000 +DEAL::2 1 0 2 120.000 +DEAL::2 1 1 0 96.0000 +DEAL::2 1 1 1 152.000 +DEAL::2 1 1 2 192.000 +DEAL::2 1 2 0 120.000 +DEAL::2 1 2 1 192.000 +DEAL::2 1 2 2 264.000 +DEAL::2 2 0 0 102.000 +DEAL::2 2 0 1 135.000 +DEAL::2 2 0 2 168.000 +DEAL::2 2 1 0 135.000 +DEAL::2 2 1 1 210.000 +DEAL::2 2 1 2 264.000 +DEAL::2 2 2 0 168.000 +DEAL::2 2 2 1 264.000 +DEAL::2 2 2 2 360.000 +DEAL::MAJOR SYMMETRY OK +DEAL::OK -- 2.39.5