--- /dev/null
+//---------------------------- symmetric_tensor_25.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2006, 2008, 2009, 2011 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- symmetric_tensor_25.cc ---------------------------
+
+// Test the invariants of tensors using the Cayley-Hamilton theorem,
+// see http://en.wikipedia.org/wiki/Invariants_of_tensors
+
+#include "../tests.h"
+#include <base/symmetric_tensor.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <fstream>
+#include <iomanip>
+
+
+template <int dim>
+void check ()
+{
+ for (unsigned int round=0; round<10; ++round)
+ {
+ SymmetricTensor<2,dim> S;
+ for (unsigned int i=0; i<S.n_independent_components; ++i)
+ S[S.unrolled_to_component_indices (i)] = std::rand () % 10;
+
+ deallog << "S = " << S << std::endl;
+ deallog << "first invariant = " << first_invariant(S) << std::endl;
+ deallog << "second invariant = " << second_invariant(S) << std::endl;
+ deallog << "third invariant = " << third_invariant(S) << std::endl;
+
+ Tensor<2,dim> S_cubed;
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ for (unsigned int f=0; f<dim; ++f)
+ for (unsigned int g=0; g<dim; ++g)
+ S_cubed[d][e] += S[d][f] * S[f][g] * S[g][e];
+
+ Tensor<2,dim> S_squared;
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ for (unsigned int f=0; f<dim; ++f)
+ S_squared[d][e] += S[d][f] * S[f][e];
+
+ Tensor<2,dim> R = S_cubed
+ - first_invariant(S) * S_squared
+ + second_invariant(S) * S
+ - third_invariant(S) * unit_symmetric_tensor<dim> ();
+ deallog << R << std::endl;
+
+ Assert (R.norm() < 1e-10, ExcInternalError());
+ }
+}
+
+
+int main ()
+{
+ std::ofstream logfile("symmetric_tensor_25/output");
+ deallog << std::setprecision(3);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ check<3> ();
+
+ deallog << "OK" << std::endl;
+}
--- /dev/null
+
+DEAL::S = 3.00 5.00 3.00 5.00 6.00 5.00 3.00 5.00 7.00
+DEAL::first invariant = 16.0
+DEAL::second invariant = 22.0
+DEAL::third invariant = -28.0
+DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
+DEAL::S = 6.00 1.00 2.00 1.00 2.00 7.00 2.00 7.00 9.00
+DEAL::first invariant = 17.0
+DEAL::second invariant = 30.0
+DEAL::third invariant = -175.
+DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
+DEAL::S = 0.00 6.00 0.00 6.00 9.00 6.00 0.00 6.00 3.00
+DEAL::first invariant = 12.0
+DEAL::second invariant = -45.0
+DEAL::third invariant = -108.
+DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
+DEAL::S = 2.00 8.00 7.00 8.00 6.00 9.00 7.00 9.00 1.00
+DEAL::first invariant = 9.00
+DEAL::second invariant = -174.
+DEAL::third invariant = 500.
+DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
+DEAL::S = 2.00 3.00 7.00 3.00 0.00 5.00 7.00 5.00 2.00
+DEAL::first invariant = 4.00
+DEAL::second invariant = -79.0
+DEAL::third invariant = 142.
+DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
+DEAL::S = 9.00 8.00 9.00 8.00 2.00 7.00 9.00 7.00 2.00
+DEAL::first invariant = 13.0
+DEAL::second invariant = -154.
+DEAL::third invariant = 313.
+DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
+DEAL::S = 3.00 2.00 9.00 2.00 6.00 3.00 9.00 3.00 1.00
+DEAL::first invariant = 10.0
+DEAL::second invariant = -67.0
+DEAL::third invariant = -391.
+DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
+DEAL::S = 1.00 7.00 8.00 7.00 9.00 4.00 8.00 4.00 4.00
+DEAL::first invariant = 14.0
+DEAL::second invariant = -80.0
+DEAL::third invariant = -304.
+DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
+DEAL::S = 5.00 6.00 1.00 6.00 0.00 0.00 1.00 0.00 3.00
+DEAL::first invariant = 8.00
+DEAL::second invariant = -22.0
+DEAL::third invariant = -108.
+DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
+DEAL::S = 6.00 0.00 6.00 0.00 3.00 1.00 6.00 1.00 2.00
+DEAL::first invariant = 11.0
+DEAL::second invariant = -1.00
+DEAL::third invariant = -78.0
+DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
+DEAL::OK