From: Reza Rastak <rastak@stanford.edu>
Date: Fri, 1 May 2020 20:36:17 +0000 (-0700)
Subject: Test for project_onto_orthogonal_tensors improved
X-Git-Tag: v9.2.0-rc1~141^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=af61f22919e3ba4112e9ad1df07d227daac578ab;p=dealii.git

Test for project_onto_orthogonal_tensors improved
---

diff --git a/tests/tensors/constexpr_tensor.cc b/tests/tensors/constexpr_tensor.cc
index 80844d985e..49748c3ae7 100644
--- a/tests/tensors/constexpr_tensor.cc
+++ b/tests/tensors/constexpr_tensor.cc
@@ -17,8 +17,6 @@
 
 #include <deal.II/base/tensor.h>
 
-#include <deal.II/lac/full_matrix.h>
-
 #include "../tests.h"
 
 template <int rank, int dim, typename Number>
@@ -179,38 +177,6 @@ main()
     DEAL_II_CONSTEXPR const auto dummy_6 = contract3(a, middle, a);
     DEAL_II_CONSTEXPR const auto dummy_7 = adjugate(a);
     DEAL_II_CONSTEXPR const auto dummy_8 = cofactor(a);
-    deallog << "Determinant before orthogonalization: " << determinant(a)
-            << std::endl;
-    const auto dummy_9 = project_onto_orthogonal_tensors(a, 0);
-    deallog << "Determinant after  orthogonalization: " << determinant(dummy_9)
-            << std::endl;
-    Assert(determinant(dummy_9) - 1. < 1e-8, ExcInternalError());
-
-    // check wheter the output is the same as the input. For example, it could
-    // have been transposed.
-    Vector<double> unrolled_a(9);
-    a.unroll<double>(unrolled_a);
-    Vector<double> unrolled_dummy_9(9);
-    dummy_9.unroll<double>(unrolled_dummy_9);
-    for (size_t i = 0; i < unrolled_a.size(); i++)
-      {
-        Assert(std::fabs(unrolled_a[i] - unrolled_dummy_9[i]) < 1e-8,
-               ExcInternalError());
-      }
-
-
-
-    constexpr double       non_orthogonal_init[3][3] = {{1., 2., 3.},
-                                                  {1., 1.5, 2.},
-                                                  {2., 0., 1}};
-    constexpr Tensor<2, 3> non_orthogonal{non_orthogonal_init};
-
-    deallog << "Determinant before orthogonalization: "
-            << determinant(non_orthogonal) << std::endl;
-    const auto dummy_10 = project_onto_orthogonal_tensors(non_orthogonal, 1e-8);
-    deallog << "Determinant after  orthogonalization: " << determinant(dummy_10)
-            << std::endl;
-    Assert(determinant(dummy_10) - 1. < 1e-8, ExcInternalError());
   }
 
   {
diff --git a/tests/tensors/constexpr_tensor.output b/tests/tensors/constexpr_tensor.output
index 31dd851cc2..5b07862e37 100644
--- a/tests/tensors/constexpr_tensor.output
+++ b/tests/tensors/constexpr_tensor.output
@@ -166,8 +166,4 @@ DEAL::Using Tensor within constexpr functions
 DEAL::15.2000
 DEAL::116.000
 DEAL::Used memory: 72
-DEAL::Determinant before orthogonalization: 1.00000
-DEAL::Determinant after  orthogonalization: 1.00000
-DEAL::Determinant before orthogonalization: -1.50000
-DEAL::Determinant after  orthogonalization: -1.00000
 DEAL::OK
diff --git a/tests/tensors/project_orthogonal.cc b/tests/tensors/project_orthogonal.cc
new file mode 100644
index 0000000000..463b836220
--- /dev/null
+++ b/tests/tensors/project_orthogonal.cc
@@ -0,0 +1,69 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test the function projection_onto_orthogonal_tensors
+
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/lac/full_matrix.h>
+
+#include "../tests.h"
+
+template <int dim, typename number>
+void
+print_type(const Tensor<2, dim, number> &)
+{
+  if (std::is_same<number, float>::value)
+    deallog << " Tensor<2, " << dim << ", float>" << std::endl;
+  else if (std::is_same<number, double>::value)
+    deallog << " Tensor<2, " << dim << ", double>" << std::endl;
+  else
+    Assert(false, ExcNotImplemented());
+}
+
+template <int dim, typename number>
+void
+test(const Tensor<2, dim, number> &A)
+{
+  print_type(A);
+  deallog << "original tensor A = " << A << std::endl;
+  deallog << "det A = " << determinant(A) << std::endl;
+  deallog << "A A^T = " << (A * transpose(A)) << std::endl;
+  const Tensor<2, dim, number> Q = project_onto_orthogonal_tensors(A, 1.e-10);
+  deallog << "projected to Q = " << Q << std::endl;
+  deallog << "det Q = " << determinant(Q) << std::endl;
+  deallog << "Q Q^T = " << (Q * transpose(Q)) << std::endl;
+  deallog << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+  // not orthogonal
+  test(Tensor<2, 3, double>{{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}});
+  // already orthogonal
+  test(Tensor<2, 3, double>{{{1, 0, 0}, {0, 0, 1}, {0, 1, 0}}});
+  // 2D not orthogonal
+  test(Tensor<2, 2, double>{{{1, 2}, {3, 4}}});
+  // 2D already orthogonal
+  test(Tensor<2, 2, double>{{{0, 1}, {1, 0}}});
+  // 1D not orthogonal
+  test(Tensor<2, 1, double>{{{2.5}}});
+  // 1D already orthogonal
+  test(Tensor<2, 1, double>{{{1.}}});
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/tensors/project_orthogonal.with_lapack=true.output b/tests/tensors/project_orthogonal.with_lapack=true.output
new file mode 100644
index 0000000000..bf9984b243
--- /dev/null
+++ b/tests/tensors/project_orthogonal.with_lapack=true.output
@@ -0,0 +1,50 @@
+
+DEAL:: Tensor<2, 3, double>
+DEAL::original tensor A = 1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 9.00000
+DEAL::det A = 0.00000
+DEAL::A A^T = 14.0000 32.0000 50.0000 32.0000 77.0000 122.000 50.0000 122.000 194.000
+DEAL::projected to Q = -0.752720 0.389148 0.531015 0.389148 -0.387594 0.835664 0.531015 0.835664 0.140313
+DEAL::det Q = 1.00000
+DEAL::Q Q^T = 1.00000 5.55112e-17 2.77556e-16 5.55112e-17 1.00000 5.55112e-17 2.77556e-16 5.55112e-17 1.00000
+DEAL::
+DEAL:: Tensor<2, 3, double>
+DEAL::original tensor A = 1.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 1.00000 0.00000
+DEAL::det A = -1.00000
+DEAL::A A^T = 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000
+DEAL::projected to Q = 1.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 1.00000 0.00000
+DEAL::det Q = -1.00000
+DEAL::Q Q^T = 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000
+DEAL::
+DEAL:: Tensor<2, 2, double>
+DEAL::original tensor A = 1.00000 2.00000 3.00000 4.00000
+DEAL::det A = -2.00000
+DEAL::A A^T = 5.00000 11.0000 11.0000 25.0000
+DEAL::projected to Q = -0.514496 0.857493 0.857493 0.514496
+DEAL::det Q = -1.00000
+DEAL::Q Q^T = 1.00000 2.77556e-16 2.77556e-16 1.00000
+DEAL::
+DEAL:: Tensor<2, 2, double>
+DEAL::original tensor A = 0.00000 1.00000 1.00000 0.00000
+DEAL::det A = -1.00000
+DEAL::A A^T = 1.00000 0.00000 0.00000 1.00000
+DEAL::projected to Q = 0.00000 1.00000 1.00000 0.00000
+DEAL::det Q = -1.00000
+DEAL::Q Q^T = 1.00000 0.00000 0.00000 1.00000
+DEAL::
+DEAL:: Tensor<2, 1, double>
+DEAL::original tensor A = 2.50000
+DEAL::det A = 2.50000
+DEAL::A A^T = 6.25000
+DEAL::projected to Q = 1.00000
+DEAL::det Q = 1.00000
+DEAL::Q Q^T = 1.00000
+DEAL::
+DEAL:: Tensor<2, 1, double>
+DEAL::original tensor A = 1.00000
+DEAL::det A = 1.00000
+DEAL::A A^T = 1.00000
+DEAL::projected to Q = 1.00000
+DEAL::det Q = 1.00000
+DEAL::Q Q^T = 1.00000
+DEAL::
+DEAL::OK