Implemented the inverse of rank-2 SymmetricTensors.

Includes a test case, gratis.

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 4aeeeed31f826658d833691c1507831ffab57cb8..001c32791da847283cee017551821bcdf523f6cc 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -223,6 +223,13 @@ inconvenience this causes.
 
 <ol>
 
+<li> New: The inverse of a rank-2 SymmetricTensor can now be directly computed
+     with SymmetricTensor::invert() instead of having to use the
+     Tensor::invert() function.
+     <br>
+     (Jean-Paul Pelteret, 2016/12/07)
+</li>
+
 <li> Improved: The step-37 tutorial program now shows the matrix-free multigrid
      solver based on MPI parallelization rather than only a serial version.
      Moreover, support for adaptively refined meshes has been added.

diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h
index dd754ec942476aeacd6ca06f73fcceb5c356a4ab..f5f9b69505c8783ef111e4f2bf6bf4a734e3a3b2 100644 (file)
--- a/include/deal.II/base/symmetric_tensor.h
+++ b/include/deal.II/base/symmetric_tensor.h
@@ -32,6 +32,8 @@ template <int dim, typename Number> SymmetricTensor<4,dim,Number>
 deviator_tensor ();
 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
 identity_tensor ();
+template <int dim, typename Number> SymmetricTensor<2,dim,Number>
+invert (const SymmetricTensor<2,dim,Number> &);
 template <int dim, typename Number> SymmetricTensor<4,dim,Number>
 invert (const SymmetricTensor<4,dim,Number> &);
 template <int dim2, typename Number> Number
@@ -823,6 +825,9 @@ private:
   template <int dim2, typename Number2>
   friend SymmetricTensor<4,dim2,Number2> identity_tensor ();
 
+  template <int dim2, typename Number2>
+  friend SymmetricTensor<2,dim2,Number2> invert (const SymmetricTensor<2,dim2,Number2> &);
+
   template <int dim2, typename Number2>
   friend SymmetricTensor<4,dim2,Number2> invert (const SymmetricTensor<4,dim2,Number2> &);
 };
@@ -2551,6 +2556,117 @@ identity_tensor ()
 
 
 
+/**
+ * Invert a symmetric rank-2 tensor.
+ *
+ * @note If a tensor is not invertible, then the result is unspecified, but will
+ * likely contain the results of a division by zero or a very small number at
+ * the very least.
+ *
+ * @relates SymmetricTensor
+ * @author Jean-Paul Pelteret, 2016
+ */
+template <int dim, typename Number>
+inline
+SymmetricTensor<2,dim,Number>
+invert (const SymmetricTensor<2,dim,Number> &t)
+{
+  // if desired, take over the
+  // inversion of a 4x4 tensor
+  // from the FullMatrix
+  AssertThrow (false, ExcNotImplemented());
+
+  return SymmetricTensor<2,dim,Number>();
+}
+
+
+
+#ifndef DOXYGEN
+
+template <typename Number>
+inline
+SymmetricTensor<2,1,Number>
+invert (const SymmetricTensor<2,1,Number> &t)
+{
+  SymmetricTensor<2,1,Number> tmp;
+
+  tmp[0][0] = 1.0/t[0][0];
+
+  return tmp;
+}
+
+
+
+template <typename Number>
+inline
+SymmetricTensor<2,2,Number>
+invert (const SymmetricTensor<2,2,Number> &t)
+{
+  SymmetricTensor<2,2,Number> tmp;
+
+  // Sympy result: ([
+  // [ t11/(t00*t11 - t01**2), -t01/(t00*t11 - t01**2)],
+  // [-t01/(t00*t11 - t01**2),  t00/(t00*t11 - t01**2)]  ])
+  const TableIndices<2> idx_00 (0,0);
+  const TableIndices<2> idx_01 (0,1);
+  const TableIndices<2> idx_11 (1,1);
+  const Number inv_det_t
+    = 1.0/(t[idx_00]*t[idx_11]
+           - t[idx_01]*t[idx_01]);
+  tmp[idx_00] = t[idx_11];
+  tmp[idx_01] = -t[idx_01];
+  tmp[idx_11] = t[idx_00];
+  tmp *= inv_det_t;
+
+  return tmp;
+}
+
+
+
+template <typename Number>
+inline
+SymmetricTensor<2,3,Number>
+invert (const SymmetricTensor<2,3,Number> &t)
+{
+  SymmetricTensor<2,3,Number> tmp;
+
+  // Sympy result: ([
+  // [  (t11*t22 - t12**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
+  //    (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
+  //    (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11)],
+  // [  (-t01*t22 + t02*t12)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
+  //    (t00*t22 - t02**2)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
+  //    (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)],
+  // [  (t01*t12 - t02*t11)/(t00*t11*t22 - t00*t12**2 - t01**2*t22 + 2*t01*t02*t12 - t02**2*t11),
+  //    (t00*t12 - t01*t02)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11),
+  //    (-t00*t11 + t01**2)/(-t00*t11*t22 + t00*t12**2 + t01**2*t22 - 2*t01*t02*t12 + t02**2*t11)]  ])
+  const TableIndices<2> idx_00 (0,0);
+  const TableIndices<2> idx_01 (0,1);
+  const TableIndices<2> idx_02 (0,2);
+  const TableIndices<2> idx_11 (1,1);
+  const TableIndices<2> idx_12 (1,2);
+  const TableIndices<2> idx_22 (2,2);
+  const Number inv_det_t
+    = 1.0/(t[idx_00]*t[idx_11]*t[idx_22]
+           - t[idx_00]*t[idx_12]*t[idx_12]
+           - t[idx_01]*t[idx_01]*t[idx_22]
+           + 2.0*t[idx_01]*t[idx_02]*t[idx_12]
+           - t[idx_02]*t[idx_02]*t[idx_11]);
+  tmp[idx_00] = t[idx_11]*t[idx_22] - t[idx_12]*t[idx_12];
+  tmp[idx_01] = -t[idx_01]*t[idx_22] + t[idx_02]*t[idx_12];
+  tmp[idx_02] = t[idx_01]*t[idx_12] - t[idx_02]*t[idx_11];
+  tmp[idx_11] = t[idx_00]*t[idx_22] - t[idx_02]*t[idx_02];
+  tmp[idx_12] = -t[idx_00]*t[idx_12] + t[idx_01]*t[idx_02];
+  tmp[idx_22] = t[idx_00]*t[idx_11] - t[idx_01]*t[idx_01];
+  tmp *= inv_det_t;
+
+  return tmp;
+}
+
+#endif /* DOXYGEN */
+
+
+
 /**
  * Invert a symmetric rank-4 tensor. Since symmetric rank-4 tensors are
  * mappings from and to symmetric rank-2 tensors, they can have an inverse.

diff --git a/tests/base/symmetric_tensor_35.cc b/tests/base/symmetric_tensor_35.cc
new file mode 100644 (file)
index 0000000..21ec1f2
--- /dev/null
+++ b/tests/base/symmetric_tensor_35.cc
@@ -0,0 +1,106 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Check the inverse of a rank-2 Tensor
+
+// Equivalent matlab script:
+/*
+printf ("Symmetric Tensor dim 1\n")
+t0 = [2]
+inv(t0)
+
+printf ("Symmetric Tensor dim 2\n")
+t1 = [2 -1; -1 1.5]
+inv(t1)
+
+printf ("Symmetric Tensor dim 3\n")
+t2 = [2   1    1.5;
+      1   1.5  0.25;
+      1.5 0.25 1.25]
+inv(t2)
+*/
+
+// Symmetric Tensor dim 1
+// t0 =  2
+// ans =  0.50000
+// Symmetric Tensor dim 2
+// t1 =
+//
+//    2.0000  -1.0000
+//   -1.0000   1.5000
+//
+// ans =
+//
+//    0.75000   0.50000
+//    0.50000   1.00000
+//
+// Symmetric Tensor dim 3
+// t2 =
+//
+//    2.00000   1.00000   1.50000
+//    1.00000   1.50000   0.25000
+//    1.50000   0.25000   1.25000
+//
+// ans =
+//
+//   -7.2500   3.5000   8.0000
+//    3.5000  -1.0000  -4.0000
+//    8.0000  -4.0000  -8.0000
+
+#include "../tests.h"
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/logstream.h>
+#include <fstream>
+#include <iomanip>
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(5);
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  deallog << "Symmetric Tensor dim 1" << std::endl;
+  SymmetricTensor<2,1> t1;
+  t1[0][0] = 2.0;
+  deallog << invert(t1) << std::endl;
+  Assert((static_cast<Tensor<2,1> >(invert(t1))*static_cast<Tensor<2,1> >(t1) - unit_symmetric_tensor<1>()).norm() < 1e-12,
+         ExcMessage("Dim 1 inverse symmetric tensor definition is incorrect"));
+
+  deallog << "Symmetric Tensor dim 2" << std::endl;
+  SymmetricTensor<2,2> t2;
+  t2[0][0] = 2.0;
+  t2[0][1] = 1.0;
+  t2[1][1] = 1.5;
+  deallog << invert(t2) << std::endl;
+  Assert((static_cast<Tensor<2,2> >(invert(t2))*static_cast<Tensor<2,2> >(t2) - unit_symmetric_tensor<2>()).norm() < 1e-12,
+         ExcMessage("Dim 2 inverse symmetric tensor definition is incorrect"));
+
+  deallog << "Symmetric Tensor dim 3" << std::endl;
+  SymmetricTensor<2,3> t3;
+  t3[0][0] = 2.0;
+  t3[0][1] = 1.0;
+  t3[0][2] = 1.5;
+  t3[1][1] = 1.5;
+  t3[1][2] = 0.25;
+  t3[2][2] = 1.25;
+  deallog << invert(t3) << std::endl;
+  Assert((static_cast<Tensor<2,3> >(invert(t3))*static_cast<Tensor<2,3> >(t3) - unit_symmetric_tensor<3>()).norm() < 1e-12,
+         ExcMessage("Dim 3 inverse symmetric tensor definition is incorrect"));
+
+  deallog << "OK" << std::endl;
+}

diff --git a/tests/base/symmetric_tensor_35.output b/tests/base/symmetric_tensor_35.output
new file mode 100644 (file)
index 0000000..dcd51a7
--- /dev/null
+++ b/tests/base/symmetric_tensor_35.output
@@ -0,0 +1,8 @@
+
+DEAL::Symmetric Tensor dim 1
+DEAL::0.50000
+DEAL::Symmetric Tensor dim 2
+DEAL::0.75000 -0.50000 -0.50000 1.0000
+DEAL::Symmetric Tensor dim 3
+DEAL::-7.2500 3.5000 8.0000 3.5000 -1.0000 -4.0000 8.0000 -4.0000 -8.0000
+DEAL::OK