]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add computation of eigenvalues of rank-2 symmetric tensor.
authorJean-Paul Pelteret <jppelteret@gmail.com>
Sun, 23 Jul 2017 08:33:15 +0000 (10:33 +0200)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Fri, 28 Jul 2017 14:59:52 +0000 (16:59 +0200)
This commit adds functions that compute the eignevalues of a rank-2
symmetric tensor by solving a characteristic polynomial equation.

doc/news/changes/major/20170727Jean-PaulPelteretEsterComellas [new file with mode: 0644]
include/deal.II/base/symmetric_tensor.h
tests/base/symmetric_tensor_40.cc [new file with mode: 0644]
tests/base/symmetric_tensor_40.output [new file with mode: 0644]

diff --git a/doc/news/changes/major/20170727Jean-PaulPelteretEsterComellas b/doc/news/changes/major/20170727Jean-PaulPelteretEsterComellas
new file mode 100644 (file)
index 0000000..39a4800
--- /dev/null
@@ -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.
+<br>
+(Jean-Paul Pelteret, Ester Comellas, 2017/07/27)
index a4d618af71eab2e258371b58e6364d098f9f5709..d45427aba58545b0b3c0de62c066621295c5d7d0 100644 (file)
 #include <deal.II/base/table_indices.h>
 #include <deal.II/base/template_constraints.h>
 
+#include <array>
+#include <algorithm>
+#include <functional>
+
 DEAL_II_NAMESPACE_OPEN
 
 template <int rank, int dim, typename Number=double> 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 <typename Number>
+std::array<Number,1>
+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
+ * <a href="https://en.wikipedia.org/wiki/Eigenvalue_algorithm#2.C3.972_matrices">the characteristic polynomial</a>
+ * $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
+ * <a href="https://scicomp.stackexchange.com/q/23686">subject to round-off errors</a>
+ * of order $\sqrt{\epsilon}$.
+ *
+ * @relates SymmetricTensor
+ * @author Jean-Paul Pelteret, 2017
+ */
+template <typename Number>
+std::array<Number,2>
+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<Number,2> eig_vals =
+      {
+        {T[0][0], T[1][1]}
+      };
+
+      // Sort from largest to smallest.
+      std::sort(eig_vals.begin(), eig_vals.end(), std::greater<Number>());
+      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<Number,2> 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
+ * <a href="https://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices">the characteristic polynomial</a>
+ * $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
+ * <a href="https://scicomp.stackexchange.com/q/23686">subject to round-off errors</a>
+ * of order $\sqrt{\epsilon}$.
+ *
+ * @relates SymmetricTensor
+ * @author Jean-Paul Pelteret, 2017
+ */
+template <typename Number>
+std::array<Number,3>
+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<Number,3> 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<Number>());
+      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<Number,3> 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<Number>());
+
+      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 (file)
index 0000000..f000fcc
--- /dev/null
@@ -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 <deal.II/base/symmetric_tensor.h>
+#include <array>
+
+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<int dim>
+bool
+is_unit_vector(const Tensor<1,dim> &v)
+{
+  return std::abs(v.norm() - 1.0) < 1e-9;
+}
+
+template<int dim>
+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<double,dim> 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<double,dim> 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<double,dim> 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 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.