]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add implementation of the determinant() function for rank-2 tensors of all sizes.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 28 Jun 2009 18:03:47 +0000 (18:03 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 28 Jun 2009 18:03:47 +0000 (18:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@18978 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/tensor.h
deal.II/doc/news/changes.h
tests/base/full_tensor_10.cc [new file with mode: 0644]
tests/base/full_tensor_10/cmp/generic [new file with mode: 0644]

index 9ad9e4630860943d58688c2d03256785dce3856c..63703cc944adf57e3a299b1ccfc36dcefaffe219 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -1401,6 +1401,45 @@ double determinant (const Tensor<2,3> &t)
 
 
 
+/**
+ * Compute the determinant of a tensor or rank 2, here for all dimensions for
+ * which no explicit specialization is available above.
+ *
+ * @relates Tensor
+ * @author Wolfgang Bangerth, 2009
+ */
+template <int dim>
+inline
+double determinant (const Tensor<2,dim> &t)
+{
+                                  // compute the determinant using the
+                                  // Laplace expansion of the
+                                  // determinant. this may not be the most
+                                  // efficient algorithm, but it does for
+                                  // small n.
+                                  //
+                                  // for some algorithmic simplicity, we use
+                                  // the expansion along the last row
+  double det = 0;
+  
+  for (unsigned int k=0; k<dim; ++k)
+    {
+      Tensor<2,dim-1> minor;
+      for (unsigned int i=0; i<dim-1; ++i)
+       for (unsigned int j=0; j<dim-1; ++j)
+         minor[i][j] = t[i][j<k ? j : j+1];
+
+      const double cofactor = std::pow (-1, k+1) *
+                             determinant (minor);
+
+      det += t[dim-1][k] * cofactor;
+    }
+  
+  return std::pow (-1, dim) * det;
+}
+
+
+
 /**
  * Compute and return the trace of a tensor of rank 2, i.e. the sum of
  * its diagonal entries.
index 71e59cc700a59ddc5b704170bc2197f8bd126e51..2ec5653fe0504a5c006656fbfcbd4546fd21944e 100644 (file)
@@ -93,6 +93,16 @@ inconvenience this causes.
 <h3>base</h3>
 
 <ol>
+  <li>
+  <p>
+  New: The determinant() function is now implemented for rank-2 Tensor
+  arguments of all sizes. The implementation is not efficient for very large
+  matrix sizes, however.
+  <br>
+  (WB 2009/06/28)
+  </p>
+  </li>
+
   <li>
   <p>
   Improved: The QGaussLobatto::gamma function now returns a long double
diff --git a/tests/base/full_tensor_10.cc b/tests/base/full_tensor_10.cc
new file mode 100644 (file)
index 0000000..61f6e42
--- /dev/null
@@ -0,0 +1,60 @@
+//----------------------------  full_tensor_10.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005, 2009 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.
+//
+//----------------------------  full_tensor_10.cc  ---------------------------
+
+// test the determinant code for n>3
+
+#include "../tests.h"
+#include <base/tensor.h>
+#include <base/logstream.h>
+#include <fstream>
+#include <iomanip>
+
+
+template <int dim>
+void test ()
+{
+  Tensor<2,dim> t;
+
+                                  // choose the same symmetric tensor
+                                  // as in symmetric_tensor_10
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=0; j<dim; ++j)
+      t[i][j] = 1. * rand() / RAND_MAX;
+
+  for (unsigned int i=0; i<dim; ++i)
+    for (unsigned int j=0; j<dim; ++j)
+      deallog << "A[" << i+1 << ',' << j+1 << "] := " << t[i][j] << ';'
+             << std::endl;
+  
+  deallog << determinant(t) << std::endl;
+}
+
+  
+
+
+int main ()
+{
+  std::ofstream logfile("full_tensor_10/output");
+  deallog << std::setprecision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test<4> ();
+  test<5> ();
+  test<6> ();
+  test<7> ();
+  test<8> ();
+  
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/base/full_tensor_10/cmp/generic b/tests/base/full_tensor_10/cmp/generic
new file mode 100644 (file)
index 0000000..5d5b86c
--- /dev/null
@@ -0,0 +1,197 @@
+
+DEAL::A[1,1] := 0.840;
+DEAL::A[1,2] := 0.394;
+DEAL::A[1,3] := 0.783;
+DEAL::A[1,4] := 0.798;
+DEAL::A[2,1] := 0.912;
+DEAL::A[2,2] := 0.198;
+DEAL::A[2,3] := 0.335;
+DEAL::A[2,4] := 0.768;
+DEAL::A[3,1] := 0.278;
+DEAL::A[3,2] := 0.554;
+DEAL::A[3,3] := 0.477;
+DEAL::A[3,4] := 0.629;
+DEAL::A[4,1] := 0.365;
+DEAL::A[4,2] := 0.513;
+DEAL::A[4,3] := 0.952;
+DEAL::A[4,4] := 0.916;
+DEAL::0.0542
+DEAL::A[1,1] := 0.636;
+DEAL::A[1,2] := 0.717;
+DEAL::A[1,3] := 0.142;
+DEAL::A[1,4] := 0.607;
+DEAL::A[1,5] := 0.0163;
+DEAL::A[2,1] := 0.243;
+DEAL::A[2,2] := 0.137;
+DEAL::A[2,3] := 0.804;
+DEAL::A[2,4] := 0.157;
+DEAL::A[2,5] := 0.401;
+DEAL::A[3,1] := 0.130;
+DEAL::A[3,2] := 0.109;
+DEAL::A[3,3] := 0.999;
+DEAL::A[3,4] := 0.218;
+DEAL::A[3,5] := 0.513;
+DEAL::A[4,1] := 0.839;
+DEAL::A[4,2] := 0.613;
+DEAL::A[4,3] := 0.296;
+DEAL::A[4,4] := 0.638;
+DEAL::A[4,5] := 0.524;
+DEAL::A[5,1] := 0.494;
+DEAL::A[5,2] := 0.973;
+DEAL::A[5,3] := 0.293;
+DEAL::A[5,4] := 0.771;
+DEAL::A[5,5] := 0.527;
+DEAL::-0.0166
+DEAL::A[1,1] := 0.770;
+DEAL::A[1,2] := 0.400;
+DEAL::A[1,3] := 0.892;
+DEAL::A[1,4] := 0.283;
+DEAL::A[1,5] := 0.352;
+DEAL::A[1,6] := 0.808;
+DEAL::A[2,1] := 0.919;
+DEAL::A[2,2] := 0.0698;
+DEAL::A[2,3] := 0.949;
+DEAL::A[2,4] := 0.526;
+DEAL::A[2,5] := 0.0861;
+DEAL::A[2,6] := 0.192;
+DEAL::A[3,1] := 0.663;
+DEAL::A[3,2] := 0.890;
+DEAL::A[3,3] := 0.349;
+DEAL::A[3,4] := 0.0642;
+DEAL::A[3,5] := 0.0200;
+DEAL::A[3,6] := 0.458;
+DEAL::A[4,1] := 0.0631;
+DEAL::A[4,2] := 0.238;
+DEAL::A[4,3] := 0.971;
+DEAL::A[4,4] := 0.902;
+DEAL::A[4,5] := 0.851;
+DEAL::A[4,6] := 0.267;
+DEAL::A[5,1] := 0.540;
+DEAL::A[5,2] := 0.375;
+DEAL::A[5,3] := 0.760;
+DEAL::A[5,4] := 0.513;
+DEAL::A[5,5] := 0.668;
+DEAL::A[5,6] := 0.532;
+DEAL::A[6,1] := 0.0393;
+DEAL::A[6,2] := 0.438;
+DEAL::A[6,3] := 0.932;
+DEAL::A[6,4] := 0.931;
+DEAL::A[6,5] := 0.721;
+DEAL::A[6,6] := 0.284;
+DEAL::0.0101
+DEAL::A[1,1] := 0.739;
+DEAL::A[1,2] := 0.640;
+DEAL::A[1,3] := 0.354;
+DEAL::A[1,4] := 0.688;
+DEAL::A[1,5] := 0.166;
+DEAL::A[1,6] := 0.440;
+DEAL::A[1,7] := 0.880;
+DEAL::A[2,1] := 0.829;
+DEAL::A[2,2] := 0.330;
+DEAL::A[2,3] := 0.229;
+DEAL::A[2,4] := 0.893;
+DEAL::A[2,5] := 0.350;
+DEAL::A[2,6] := 0.687;
+DEAL::A[2,7] := 0.956;
+DEAL::A[3,1] := 0.589;
+DEAL::A[3,2] := 0.657;
+DEAL::A[3,3] := 0.859;
+DEAL::A[3,4] := 0.440;
+DEAL::A[3,5] := 0.924;
+DEAL::A[3,6] := 0.398;
+DEAL::A[3,7] := 0.815;
+DEAL::A[4,1] := 0.684;
+DEAL::A[4,2] := 0.911;
+DEAL::A[4,3] := 0.482;
+DEAL::A[4,4] := 0.216;
+DEAL::A[4,5] := 0.950;
+DEAL::A[4,6] := 0.920;
+DEAL::A[4,7] := 0.148;
+DEAL::A[5,1] := 0.881;
+DEAL::A[5,2] := 0.641;
+DEAL::A[5,3] := 0.432;
+DEAL::A[5,4] := 0.620;
+DEAL::A[5,5] := 0.281;
+DEAL::A[5,6] := 0.786;
+DEAL::A[5,7] := 0.307;
+DEAL::A[6,1] := 0.447;
+DEAL::A[6,2] := 0.226;
+DEAL::A[6,3] := 0.188;
+DEAL::A[6,4] := 0.276;
+DEAL::A[6,5] := 0.556;
+DEAL::A[6,6] := 0.417;
+DEAL::A[6,7] := 0.170;
+DEAL::A[7,1] := 0.907;
+DEAL::A[7,2] := 0.103;
+DEAL::A[7,3] := 0.126;
+DEAL::A[7,4] := 0.495;
+DEAL::A[7,5] := 0.760;
+DEAL::A[7,6] := 0.985;
+DEAL::A[7,7] := 0.935;
+DEAL::0.0109
+DEAL::A[1,1] := 0.684;
+DEAL::A[1,2] := 0.383;
+DEAL::A[1,3] := 0.750;
+DEAL::A[1,4] := 0.369;
+DEAL::A[1,5] := 0.294;
+DEAL::A[1,6] := 0.232;
+DEAL::A[1,7] := 0.584;
+DEAL::A[1,8] := 0.244;
+DEAL::A[2,1] := 0.152;
+DEAL::A[2,2] := 0.732;
+DEAL::A[2,3] := 0.125;
+DEAL::A[2,4] := 0.793;
+DEAL::A[2,5] := 0.164;
+DEAL::A[2,6] := 0.745;
+DEAL::A[2,7] := 0.0745;
+DEAL::A[2,8] := 0.950;
+DEAL::A[3,1] := 0.0525;
+DEAL::A[3,2] := 0.522;
+DEAL::A[3,3] := 0.176;
+DEAL::A[3,4] := 0.240;
+DEAL::A[3,5] := 0.798;
+DEAL::A[3,6] := 0.733;
+DEAL::A[3,7] := 0.657;
+DEAL::A[3,8] := 0.967;
+DEAL::A[4,1] := 0.639;
+DEAL::A[4,2] := 0.760;
+DEAL::A[4,3] := 0.0935;
+DEAL::A[4,4] := 0.135;
+DEAL::A[4,5] := 0.520;
+DEAL::A[4,6] := 0.0782;
+DEAL::A[4,7] := 0.0699;
+DEAL::A[4,8] := 0.205;
+DEAL::A[5,1] := 0.461;
+DEAL::A[5,2] := 0.820;
+DEAL::A[5,3] := 0.573;
+DEAL::A[5,4] := 0.756;
+DEAL::A[5,5] := 0.0519;
+DEAL::A[5,6] := 0.158;
+DEAL::A[5,7] := 1.00;
+DEAL::A[5,8] := 0.204;
+DEAL::A[6,1] := 0.890;
+DEAL::A[6,2] := 0.125;
+DEAL::A[6,3] := 0.998;
+DEAL::A[6,4] := 0.0541;
+DEAL::A[6,5] := 0.871;
+DEAL::A[6,6] := 0.0723;
+DEAL::A[6,7] := 0.00416;
+DEAL::A[6,8] := 0.923;
+DEAL::A[7,1] := 0.594;
+DEAL::A[7,2] := 0.180;
+DEAL::A[7,3] := 0.163;
+DEAL::A[7,4] := 0.392;
+DEAL::A[7,5] := 0.913;
+DEAL::A[7,6] := 0.820;
+DEAL::A[7,7] := 0.359;
+DEAL::A[7,8] := 0.552;
+DEAL::A[8,1] := 0.579;
+DEAL::A[8,2] := 0.453;
+DEAL::A[8,3] := 0.687;
+DEAL::A[8,4] := 0.0996;
+DEAL::A[8,5] := 0.531;
+DEAL::A[8,6] := 0.757;
+DEAL::A[8,7] := 0.304;
+DEAL::A[8,8] := 0.992;
+DEAL::-0.0318
+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.