+/**
+ * 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 <int dim, typename Number>
+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<Number>::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
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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 <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/vectorization.h>
+
+#include "../tests.h"
+
+
+template <int dim, typename Number>
+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 <int dim, typename Number>
+void
+print_fourth_order_tensor(const SymmetricTensor<4, dim, Number> &t)
+{
+ const Number tolerance = 100. * std::numeric_limits<Number>::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;
+}
projects / dealii.git / commitdiff