From 85ca7bf15714112659a08ede45410ac811a78c55 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Wed, 7 Dec 2016 10:05:24 +0100 Subject: [PATCH] Refactor definition of rank-2 Tensor inverse function; add test case This silences warnings produced by `-Warray-bounds`. --- include/deal.II/base/tensor.h | 119 ++++++++++++++++++++-------------- tests/base/tensor_27.cc | 81 +++++++++++++++++++++++ tests/base/tensor_27.output | 8 +++ 3 files changed, 160 insertions(+), 48 deletions(-) create mode 100644 tests/base/tensor_27.cc create mode 100644 tests/base/tensor_27.output diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index d5f600d615..78f857f80d 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -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 +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 +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 +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 #endif - diff --git a/tests/base/tensor_27.cc b/tests/base/tensor_27.cc new file mode 100644 index 0000000000..86cfe69b7c --- /dev/null +++ b/tests/base/tensor_27.cc @@ -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 +#include +#include +#include +#include + +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 index 0000000000..9a2c62c96f --- /dev/null +++ b/tests/base/tensor_27.output @@ -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 -- 2.39.5