]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Refactor definition of rank-2 Tensor inverse function; add test case
authorJean-Paul Pelteret <jppelteret@gmail.com>
Wed, 7 Dec 2016 09:05:24 +0000 (10:05 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Wed, 7 Dec 2016 14:38:24 +0000 (15:38 +0100)
This silences warnings produced by `-Warray-bounds`.

include/deal.II/base/tensor.h
tests/base/tensor_27.cc [new file with mode: 0644]
tests/base/tensor_27.output [new file with mode: 0644]

index d5f600d615b6d0d6b50fabe29e799f0532dce84a..78f857f80dcc23ecd901823471ce4cf658107c71 100644 (file)
@@ -1842,59 +1842,83 @@ Tensor<2,dim,Number>
 invert (const Tensor<2,dim,Number> &t)
 {
   Number return_tensor [dim][dim];
-  switch (dim)
-    {
-    case 1:
-      return_tensor[0][0] = 1.0/t[0][0];
-      break;
 
-    case 2:
-      // this is Maple output,
-      // thus a bit unstructured
-    {
-      const Number det = t[0][0]*t[1][1]-t[1][0]*t[0][1];
-      const Number t4 = 1.0/det;
-      return_tensor[0][0] = t[1][1]*t4;
-      return_tensor[0][1] = -t[0][1]*t4;
-      return_tensor[1][0] = -t[1][0]*t4;
-      return_tensor[1][1] = t[0][0]*t4;
-      break;
-    }
+  // if desired, take over the
+  // inversion of a 4x4 tensor
+  // from the FullMatrix
+  AssertThrow (false, ExcNotImplemented());
 
-    case 3:
-    {
-      const Number t4 = t[0][0]*t[1][1],
-                   t6 = t[0][0]*t[1][2],
-                   t8 = t[0][1]*t[1][0],
-                   t00 = t[0][2]*t[1][0],
-                   t01 = t[0][1]*t[2][0],
-                   t04 = t[0][2]*t[2][0],
-                   det = (t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
-                          t00*t[2][1]+t01*t[1][2]-t04*t[1][1]),
-                         t07 = 1.0/det;
-      return_tensor[0][0] = (t[1][1]*t[2][2]-t[1][2]*t[2][1])*t07;
-      return_tensor[0][1] = (t[0][2]*t[2][1]-t[0][1]*t[2][2])*t07;
-      return_tensor[0][2] = (t[0][1]*t[1][2]-t[0][2]*t[1][1])*t07;
-      return_tensor[1][0] = (t[1][2]*t[2][0]-t[1][0]*t[2][2])*t07;
-      return_tensor[1][1] = (t[0][0]*t[2][2]-t04)*t07;
-      return_tensor[1][2] = (t00-t6)*t07;
-      return_tensor[2][0] = (t[1][0]*t[2][1]-t[1][1]*t[2][0])*t07;
-      return_tensor[2][1] = (t01-t[0][0]*t[2][1])*t07;
-      return_tensor[2][2] = (t4-t8)*t07;
-
-      break;
-    }
-
-    // if desired, take over the
-    // inversion of a 4x4 tensor
-    // from the FullMatrix
-    default:
-      AssertThrow (false, ExcNotImplemented());
-    }
   return Tensor<2,dim,Number>(return_tensor);
 }
 
 
+#ifndef DOXYGEN
+
+template <typename Number>
+inline
+Tensor<2,1,Number>
+invert (const Tensor<2,1,Number> &t)
+{
+  Number return_tensor [1][1];
+
+  return_tensor[0][0] = 1.0/t[0][0];
+
+  return Tensor<2,1,Number>(return_tensor);
+}
+
+
+template <typename Number>
+inline
+Tensor<2,2,Number>
+invert (const Tensor<2,2,Number> &t)
+{
+  Tensor<2,2,Number> return_tensor;
+
+  // this is Maple output,
+  // thus a bit unstructured
+  const Number inv_det_t = 1.0/(t[0][0]*t[1][1]-t[1][0]*t[0][1]);
+  return_tensor[0][0] = t[1][1];
+  return_tensor[0][1] = -t[0][1];
+  return_tensor[1][0] = -t[1][0];
+  return_tensor[1][1] = t[0][0];
+  return_tensor *= inv_det_t;
+
+  return return_tensor;
+}
+
+
+template <typename Number>
+inline
+Tensor<2,3,Number>
+invert (const Tensor<2,3,Number> &t)
+{
+  Tensor<2,3,Number> return_tensor;
+
+  const Number t4 = t[0][0]*t[1][1],
+               t6 = t[0][0]*t[1][2],
+               t8 = t[0][1]*t[1][0],
+               t00 = t[0][2]*t[1][0],
+               t01 = t[0][1]*t[2][0],
+               t04 = t[0][2]*t[2][0],
+               inv_det_t = 1.0/(t4*t[2][2]-t6*t[2][1]-t8*t[2][2]+
+                                t00*t[2][1]+t01*t[1][2]-t04*t[1][1]);
+  return_tensor[0][0] = t[1][1]*t[2][2]-t[1][2]*t[2][1];
+  return_tensor[0][1] = t[0][2]*t[2][1]-t[0][1]*t[2][2];
+  return_tensor[0][2] = t[0][1]*t[1][2]-t[0][2]*t[1][1];
+  return_tensor[1][0] = t[1][2]*t[2][0]-t[1][0]*t[2][2];
+  return_tensor[1][1] = t[0][0]*t[2][2]-t04;
+  return_tensor[1][2] = t00-t6;
+  return_tensor[2][0] = t[1][0]*t[2][1]-t[1][1]*t[2][0];
+  return_tensor[2][1] = t01-t[0][0]*t[2][1];
+  return_tensor[2][2] = t4-t8;
+  return_tensor *= inv_det_t;
+
+  return return_tensor;
+}
+
+#endif /* DOXYGEN */
+
+
 /**
  * Return the transpose of the given tensor.
  *
@@ -1981,4 +2005,3 @@ DEAL_II_NAMESPACE_CLOSE
 #include <deal.II/base/tensor_deprecated.h>
 
 #endif
-
diff --git a/tests/base/tensor_27.cc b/tests/base/tensor_27.cc
new file mode 100644 (file)
index 0000000..86cfe69
--- /dev/null
@@ -0,0 +1,81 @@
+// ---------------------------------------------------------------------
+//
+// 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 ("Tensor dim 1\n")
+t1 = [2]
+inv(t1)
+
+printf ("Tensor dim 2\n")
+t2 = [2 1; -1 1.5]
+inv(t2)
+
+printf ("Tensor dim 3\n")
+t3 = [2 1 0.5; -1 1.5 0.25; 1.5 -0.75 1.25]
+inv(t3)
+*/
+
+#include "../tests.h"
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/symmetric_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 << "Tensor dim 1" << std::endl;
+  Tensor<2,1> t1;
+  t1[0][0] = 2.0;
+  deallog << invert(t1) << std::endl;
+  Assert((invert(t1)*t1 - unit_symmetric_tensor<1>()).norm() < 1e-12,
+         ExcMessage("Dim 1 inverse tensor definition is incorrect"));
+
+  deallog << "Tensor dim 2" << std::endl;
+  Tensor<2,2> t2;
+  t2[0][0] = 2.0;
+  t2[0][1] = 1.0;
+  t2[1][0] = -1.0;
+  t2[1][1] = 1.5;
+  deallog << invert(t2) << std::endl;
+  Assert((invert(t2)*t2 - unit_symmetric_tensor<2>()).norm() < 1e-12,
+         ExcMessage("Dim 2 inverse tensor definition is incorrect"));
+
+  deallog << "Tensor dim 3" << std::endl;
+  Tensor<2,3> t3;
+  t3[0][0] = 2.0;
+  t3[0][1] = 1.0;
+  t3[0][2] = 0.5;
+  t3[1][0] = -1.0;
+  t3[1][1] = 1.5;
+  t3[1][2] = 0.25;
+  t3[2][0] = 1.5;
+  t3[2][1] = -0.75;
+  t3[2][2] = 1.25;
+  deallog << invert(t3) << std::endl;
+  Assert((invert(t3)*t3 - unit_symmetric_tensor<3>()).norm() < 1e-12,
+         ExcMessage("Dim 3 inverse tensor definition is incorrect"));
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/base/tensor_27.output b/tests/base/tensor_27.output
new file mode 100644 (file)
index 0000000..9a2c62c
--- /dev/null
@@ -0,0 +1,8 @@
+
+DEAL::Tensor dim 1
+DEAL::0.50000
+DEAL::Tensor dim 2
+DEAL::0.37500 -0.25000 0.25000 0.50000
+DEAL::Tensor dim 3
+DEAL::0.41250 -0.32500 -0.10000 0.32500 0.35000 -0.20000 -0.30000 0.60000 0.80000
+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.