From d88188eba7036dbe6997c94133d320047437cfb3 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Wed, 7 Dec 2016 11:27:00 +0100 Subject: [PATCH] Implemented the inverse of rank-2 SymmetricTensors. Includes a test case, gratis. --- doc/news/changes.h | 7 ++ include/deal.II/base/symmetric_tensor.h | 116 ++++++++++++++++++++++++ tests/base/symmetric_tensor_35.cc | 106 ++++++++++++++++++++++ tests/base/symmetric_tensor_35.output | 8 ++ 4 files changed, 237 insertions(+) create mode 100644 tests/base/symmetric_tensor_35.cc create mode 100644 tests/base/symmetric_tensor_35.output diff --git a/doc/news/changes.h b/doc/news/changes.h index 4aeeeed31f..001c32791d 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -223,6 +223,13 @@ inconvenience this causes.
    +
  1. 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. +
    + (Jean-Paul Pelteret, 2016/12/07) +
    2. +
  2. 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 dd754ec942..f5f9b69505 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -32,6 +32,8 @@ template SymmetricTensor<4,dim,Number> deviator_tensor (); template SymmetricTensor<4,dim,Number> identity_tensor (); +template SymmetricTensor<2,dim,Number> +invert (const SymmetricTensor<2,dim,Number> &); template SymmetricTensor<4,dim,Number> invert (const SymmetricTensor<4,dim,Number> &); template Number @@ -823,6 +825,9 @@ private: template friend SymmetricTensor<4,dim2,Number2> identity_tensor (); + template + friend SymmetricTensor<2,dim2,Number2> invert (const SymmetricTensor<2,dim2,Number2> &); + template 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 +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 +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 +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 +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 index 0000000000..21ec1f24b2 --- /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 +#include +#include +#include +#include + +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 >(invert(t1))*static_cast >(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 >(invert(t2))*static_cast >(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 >(invert(t3))*static_cast >(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 index 0000000000..dcd51a7435 --- /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 -- 2.39.5