From 14d85530c1b0aa28043071b994b3e4fda3a41c5e Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Sun, 23 Jul 2017 10:33:15 +0200 Subject: [PATCH] Add computation of eigenvalues of rank-2 symmetric tensor. This commit adds functions that compute the eignevalues of a rank-2 symmetric tensor by solving a characteristic polynomial equation. --- .../20170727Jean-PaulPelteretEsterComellas | 4 + include/deal.II/base/symmetric_tensor.h | 159 ++++++++++ tests/base/symmetric_tensor_40.cc | 275 ++++++++++++++++++ tests/base/symmetric_tensor_40.output | 2 + 4 files changed, 440 insertions(+) create mode 100644 doc/news/changes/major/20170727Jean-PaulPelteretEsterComellas create mode 100644 tests/base/symmetric_tensor_40.cc create mode 100644 tests/base/symmetric_tensor_40.output diff --git a/doc/news/changes/major/20170727Jean-PaulPelteretEsterComellas b/doc/news/changes/major/20170727Jean-PaulPelteretEsterComellas new file mode 100644 index 0000000000..39a48001c0 --- /dev/null +++ b/doc/news/changes/major/20170727Jean-PaulPelteretEsterComellas @@ -0,0 +1,4 @@ +New: The eigenvalues of a rank-2 symmetric tensor can now be computed using an +analytical approach via the eigenvalues() function. +
+(Jean-Paul Pelteret, Ester Comellas, 2017/07/27) diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index a4d618af71..d45427aba5 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -22,6 +22,10 @@ #include #include +#include +#include +#include + DEAL_II_NAMESPACE_OPEN template class SymmetricTensor; @@ -2320,6 +2324,161 @@ Number second_invariant (const SymmetricTensor<2,3,Number> &t) +/** + * Return the eigenvalues of a symmetric tensor of rank 2. + * + * @relates SymmetricTensor + * @author Jean-Paul Pelteret, 2017 + */ +template +std::array +eigenvalues (const SymmetricTensor<2,1,Number> &T) +{ + return { {T[0][0]} }; +} + + + +/** + * Return the eigenvalues of a symmetric tensor of rank 2. + * The array of eigenvalues is sorted in descending order. + * + * For 2x2 tensors, the eigenvalues of tensor $T$ are the roots of + * the characteristic polynomial + * $0 = \lambda^{2} - \lambda*tr(T) + det(T)$ + * as given by + * $\lambda = \frac{tr(T) \pm \sqrt{tr^{2}(T) - 4*det(T)}}{2}$. + * + * @warning The algorithm employed here determines the eigenvalues by + * computing the roots of the characteristic polynomial. In the case that there + * exists a common root (the eigenvalues are equal), the computation is + * subject to round-off errors + * of order $\sqrt{\epsilon}$. + * + * @relates SymmetricTensor + * @author Jean-Paul Pelteret, 2017 + */ +template +std::array +eigenvalues (const SymmetricTensor<2,2,Number> &T) +{ + const Number upp_tri_sq = T[0][1]*T[0][1]; + if (upp_tri_sq == Number(0.0)) + { + // The tensor is diagonal + std::array eig_vals = + { + {T[0][0], T[1][1]} + }; + + // Sort from largest to smallest. + std::sort(eig_vals.begin(), eig_vals.end(), std::greater()); + return eig_vals; + } + else + { + const Number tr_T = trace(T); + const Number det_T = determinant(T); + const Number descrim = tr_T*tr_T - 4.0*det_T; + Assert(descrim > Number(0.0), ExcMessage("The roots of the characteristic polynomial are complex valued.")); + const Number sqrt_desc = std::sqrt(descrim); + + std::array eig_vals = + { + { + 0.5*(tr_T + sqrt_desc), + 0.5*(tr_T - sqrt_desc) + } + }; + Assert(eig_vals[0] >= eig_vals[1], ExcMessage("The eigenvalue ordering is incorrect.")); + return eig_vals; + } +} + + + +/** + * Return the eigenvalues of a symmetric tensor of rank 2. + * The array of eigenvalues is sorted in descending order. + * + * For 3x3 tensors, the eigenvalues of tensor $T$ are the roots of + * the characteristic polynomial + * $0 = \lambda^{3} - \lambda^{2}*tr(T) - \frac{1}{2} \lambda (tr(T^{2}) - tr^{2}(T)) - det(T)$. + * + * @warning The algorithm employed here determines the eigenvalues by + * computing the roots of the characteristic polynomial. In the case that there + * exists a common root (the eigenvalues are equal), the computation is + * subject to round-off errors + * of order $\sqrt{\epsilon}$. + * + * @relates SymmetricTensor + * @author Jean-Paul Pelteret, 2017 + */ +template +std::array +eigenvalues (const SymmetricTensor<2,3,Number> &T) +{ + const Number upp_tri_sq = T[0][1]*T[0][1] + T[0][2]*T[0][2] + T[1][2]*T[1][2]; + if (upp_tri_sq == Number(0.0)) + { + // The tensor is diagonal + std::array eig_vals + = { {T[0][0], T[1][1], T[2][2]} }; + + // Sort from largest to smallest. + std::sort(eig_vals.begin(), eig_vals.end(), std::greater()); + return eig_vals; + } + else + { + // Perform an affine change to T, and solve a different + // characteristic equation that has a trigonometric solution. + // Decompose T = p*B + q*I , and set q = tr(T)/3 + // and p = (tr((T - q.I)^{2})/6)^{1/2} . Then solve the equation + // 0 = det(\lambda*I - B) = \lambda^{3} - 3*\lambda - det(B) + // which has the solution + // \lambda = 2*cos(1/3 * acos(det(B)/2) +2/3*pi*k ) ; k = 0,1,2 + // when substituting \lambda = 2.cos(theta) and using trig identities. + const Number tr_T = trace(T); + const Number q = tr_T/3.0; + const Number tmp1 = ( T[0][0] - q)*(T[0][0] - q) + + (T[1][1] - q)*(T[1][1] - q) + + (T[2][2] - q)*(T[2][2] - q) + + 2.0 * upp_tri_sq; + const Number p = std::sqrt(tmp1/6.0); + const SymmetricTensor<2,3,Number> B = (1.0/p)*(T - q*unit_symmetric_tensor<3,Number>()); + const Number tmp_2 = determinant(B)/2.0; + + // The value of tmp_2 should be within [-1,1], however + // floating point errors might place it slightly outside + // this range. It is therefore necessary to correct for it + const Number phi = + (tmp_2 <= -1.0 ? M_PI/3.0 : + (tmp_2 >= 1.0 ? 0.0 : + std::acos(tmp_2)/3.0)); + + // Due to the trigonometric solution, the computed eigenvalues + // should be predictably in the order eig1 >= eig2 >= eig3... + std::array eig_vals + = { { + q + 2.0*p *std::cos(phi), + 0.0, + q + 2.0*p *std::cos(phi + (2.0/3.0*M_PI)) + } + }; + // Use the identity tr(T) = eig1 + eig2 + eig3 + eig_vals[1] = tr_T - eig_vals[0] - eig_vals[2]; + + // ... however, when equal roots exist then floating point + // errors may make this no longer be the case. + // Sort from largest to smallest. + std::sort(eig_vals.begin(), eig_vals.end(), std::greater()); + + return eig_vals; + } +} + + /** * Return the transpose of the given symmetric tensor. Since we are working diff --git a/tests/base/symmetric_tensor_40.cc b/tests/base/symmetric_tensor_40.cc new file mode 100644 index 0000000000..f000fcc42a --- /dev/null +++ b/tests/base/symmetric_tensor_40.cc @@ -0,0 +1,275 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test that SymmetricTensor eigenvalue calculations are correct + +#include "../tests.h" + +#include +#include + +void +check_value (const int dim, const int index, + const double expected, const double actual, + const double tol = 1e-12) +{ + const double rel_error = std::abs(expected - actual)/std::abs(actual); + if (rel_error > tol) + { + deallog + << "Incorrect value calculated: " + << "Dim " << dim + << ", eigenvalue " << index + << ". Expected " << expected + << ", actual: " << actual + << ", relative error: " << rel_error + << std::endl; + } +}; + +template +bool +is_unit_vector(const Tensor<1,dim> &v) +{ + return std::abs(v.norm() - 1.0) < 1e-9; +} + +template +bool +check_orientation(Tensor<1,dim> v1, + Tensor<1,dim> v2) +{ + v1 /= v1.norm(); + v2 /= v2.norm(); + return std::abs(std::abs(v1*v2) - 1.0) < 1e-9; +} + +void +test_dim_1 (const double e1, const double tol = 1e-12) +{ + const unsigned int dim = 1; + SymmetricTensor<2,dim> T; + T[0][0] = e1; + const std::array eig_vals = eigenvalues(T); + + check_value(dim, 0, e1, eig_vals[0], tol); +} + +void +test_dim_2 (const double e1, Tensor<1,2> v1, + const double e2, const double tol = 1e-12) +{ + const unsigned int dim = 2; + v1 /= v1.norm(); + const Tensor<1,dim> v2 = cross_product_2d(v1); + + Assert(is_unit_vector(v1), ExcMessage("Vector is not of unit length.")); + Assert(is_unit_vector(v2), ExcMessage("Vector is not of unit length.")); + Assert(e1 >= e2, ExcMessage("Input eigenvalue ordering is not correct.")); + + const SymmetricTensor<2,dim> T + = e1*symmetrize(outer_product(v1,v1)) + + e2*symmetrize(outer_product(v2,v2)); + + const std::array eig_vals = eigenvalues(T); + + check_value(dim, 0, e1, eig_vals[0], tol); + check_value(dim, 1, e2, eig_vals[1], tol); +} + +void +test_dim_3 (const double e1, Tensor<1,3> v1, + const double e2, Tensor<1,3> v2, + const double e3, const double tol = 1e-12) +{ + const unsigned int dim = 3; + + // Note: We do not necessarily expect the seed directors v1,v2 to be orthogonal + v1 /= v1.norm(); + v2 /= v2.norm(); + Tensor<1,dim> v3 = cross_product_3d(v1,v2); + v3 /= v3.norm(); + v1 = cross_product_3d(v2,v3); + v3 = cross_product_3d(v1,v2); // Ensure that system is right-handed + + Assert(is_unit_vector(v1), ExcMessage("Vector is not of unit length.")); + Assert(is_unit_vector(v2), ExcMessage("Vector is not of unit length.")); + Assert(is_unit_vector(v3), ExcMessage("Vector is not of unit length.")); + Assert(check_orientation(v2, cross_product_3d(v3,v1)), ExcMessage("Vectors are not orthogonal.")); + Assert(e1 >= e2, ExcMessage("Input eigenvalue ordering is not correct.")); + Assert(e2 >= e3, ExcMessage("Input eigenvalue ordering is not correct.")); + + + const SymmetricTensor<2,dim> T + = e1*symmetrize(outer_product(v1,v1)) + + e2*symmetrize(outer_product(v2,v2)) + + e3*symmetrize(outer_product(v3,v3)); + + const std::array eig_vals = eigenvalues(T); + + check_value(dim, 0, e1, eig_vals[0], tol); + check_value(dim, 1, e2, eig_vals[1], tol); + check_value(dim, 2, e3, eig_vals[2], tol); +} + +int main() +{ + initlog(); + + const double e2 = 2.8; + const double e3 = 1.2; + + // Dim = 1 + { + deallog.push("Test 1"); + { + test_dim_1(3.6); + } + deallog.pop(); + } + + // Dim = 2 + { + // Diagonal + deallog.push("Test 2a"); + { + test_dim_2( + 3.6, Tensor<1,2>({1,0}), + 2.4 ); + } + deallog.pop(); + + // Diagonal (large difference) + deallog.push("Test 2b"); + { + test_dim_2( + 1.2e7, Tensor<1,2>({1,0}), + -0.2e-8 ); + } + deallog.pop(); + + // Diagonal (equal) + deallog.push("Test 2c"); + { + test_dim_2( + 16.7, Tensor<1,2>({1,0}), + 16.7 ); + } + deallog.pop(); + + // Non-diagonal + deallog.push("Test 2d"); + { + test_dim_2( + 115.7, Tensor<1,2>({1,1}), + 13.6 ); + } + deallog.pop(); + + // Non-diagonal (large difference) + deallog.push("Test 2e"); + { + test_dim_2( + 7.2956e8, Tensor<1,2>({3,2}), + -5.284e3 ); + } + deallog.pop(); + + // Non-diagonal (equal) + deallog.push("Test 2e"); + { + test_dim_2( + -43.2, Tensor<1,2>({1.5,-0.55}), + -43.2 ); + } + deallog.pop(); + } + + // Dim = 3 + { + // Diagonal + deallog.push("Test 3a"); + { + test_dim_3( + 3.6, Tensor<1,3>({1,0,0}), + 2.4, Tensor<1,3>({0,1,0}), + 1.2); + } + deallog.pop(); + + // Diagonal (large difference) + deallog.push("Test 3b"); + { + test_dim_3( + 1.2e7, Tensor<1,3>({1,0,0}), + -0.2e-8, Tensor<1,3>({0,1,0}), + -6.5e8); + } + deallog.pop(); + + // Diagonal (2 equal) + deallog.push("Test 3c"); + { + test_dim_3( + 16.7, Tensor<1,3>({1,0,0}), + 16.7, Tensor<1,3>({0,1,0}), + 1e-6); + } + deallog.pop(); + + // Diagonal (3 equal) + deallog.push("Test 3d"); + { + test_dim_3( + 4.2, Tensor<1,3>({1,0,0}), + 4.2, Tensor<1,3>({0,1,0}), + 4.2); + } + deallog.pop(); + + // Non-diagonal + deallog.push("Test 3e"); + { + test_dim_3( + 115.7, Tensor<1,3>({1,1,1}), + 13.6, Tensor<1,3>({-1,1,-1}), + -45.2); + } + deallog.pop(); + + // Non-diagonal (1 large difference) + deallog.push("Test 3f"); + { + test_dim_3( + 7.2956e8, Tensor<1,3>({3,2,5}), + -4.856e3, Tensor<1,3>({-0.2,3,1}), + -5.284e3, 5e-6); + } + deallog.pop(); + + // Non-diagonal (2 large difference) + deallog.push("Test 3g"); + { + const double tol = 1e-7; + test_dim_3( + 9.274e7, Tensor<1,3>({2,-0.7,1.4}), + 2.59343, Tensor<1,3>({0.5,-0.22,-1.42}), + -5.292e8, tol); + } + deallog.pop(); + } + + deallog << "OK" << std::endl; +} diff --git a/tests/base/symmetric_tensor_40.output b/tests/base/symmetric_tensor_40.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/base/symmetric_tensor_40.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5