]> https://gitweb.dealii.org/ - dealii.git/commitdiff
The tolerance parameter is removed from project_tensor 10033/head
authorReza Rastak <rastak@stanford.edu>
Tue, 5 May 2020 02:18:19 +0000 (19:18 -0700)
committerReza Rastak <rastak@stanford.edu>
Tue, 5 May 2020 02:18:19 +0000 (19:18 -0700)
include/deal.II/base/tensor.h
tests/tensors/project_orthogonal.cc
tests/tensors/project_orthogonal.with_lapack=true.output

index 407f1a78a3572bc95f6c5e535f6de965806c8bf0..a8d6dd315766fd70b43a305adb9c1823ccc42c82 100644 (file)
@@ -2632,47 +2632,38 @@ cofactor(const Tensor<2, dim, Number> &t)
 
 
 /**
- * Return the nearest orthogonal matrix using a SVD if the determinant is
- * more than a tolerance away from one. The orthogonalization is done by
+ * Return the nearest orthogonal matrix by
  * combining the products of the SVD decomposition: $\mathbf U \mathbf{V}^T$,
  * where $\mathbf U$ and $\mathbf V$ are computed from the SVD decomposition:
  * $\mathbf U  \mathbf S \mathbf V^T$,
  * effectively replacing $\mathbf S$ with the identity matrix.
  * @param tensor The tensor which to find the closest orthogonal
  * tensor to.
- * @param tolerance If the $\text{determinant} - 1$ is smaller than
- * this value, it will just return the current tensor.
- * Otherwise it will return the nearest orthogonal tensor.
  * @relatesalso Tensor
  */
 template <int dim, typename Number>
 Tensor<2, dim, Number>
-project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &tensor,
-                                const double                  tolerance)
+project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &tensor)
 {
-  if (std::abs(determinant(tensor) - 1.0) > tolerance)
-    {
-      Tensor<2, dim, Number>   output_tensor;
-      FullMatrix<Number>       matrix(dim);
-      LAPACKFullMatrix<Number> lapack_matrix(dim);
-      LAPACKFullMatrix<Number> result(dim);
+  Tensor<2, dim, Number>   output_tensor;
+  FullMatrix<Number>       matrix(dim);
+  LAPACKFullMatrix<Number> lapack_matrix(dim);
+  LAPACKFullMatrix<Number> result(dim);
 
-      // todo: find or add dealii functionallity to copy in one step.
-      matrix.copy_from(tensor);
-      lapack_matrix.copy_from(matrix);
+  // todo: find or add dealii functionallity to copy in one step.
+  matrix.copy_from(tensor);
+  lapack_matrix.copy_from(matrix);
 
-      // now compute the svd of the matrices
-      lapack_matrix.compute_svd();
+  // now compute the svd of the matrices
+  lapack_matrix.compute_svd();
 
-      // Use the SVD results to orthogonalize: $U V^T$
-      lapack_matrix.get_svd_u().mmult(result, lapack_matrix.get_svd_vt());
+  // Use the SVD results to orthogonalize: $U V^T$
+  lapack_matrix.get_svd_u().mmult(result, lapack_matrix.get_svd_vt());
 
-      // todo: find or add dealii functionallity to copy in one step.
-      matrix = result;
-      matrix.copy_to(output_tensor);
-      return output_tensor;
-    }
-  return tensor;
+  // todo: find or add dealii functionallity to copy in one step.
+  matrix = result;
+  matrix.copy_to(output_tensor);
+  return output_tensor;
 }
 
 
index 463b836220704adac8bc05baa0bc0713b9bb6a3b..f62bcf5805f742a5799403410c75b01aa0845d29 100644 (file)
@@ -41,7 +41,7 @@ test(const Tensor<2, dim, number> &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);
+  const Tensor<2, dim, number> Q = project_onto_orthogonal_tensors(A);
   deallog << "projected to Q = " << Q << std::endl;
   deallog << "det Q = " << determinant(Q) << std::endl;
   deallog << "Q Q^T = " << (Q * transpose(Q)) << std::endl;
@@ -56,6 +56,8 @@ main()
   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}}});
+  // not orthogonal, but det = 1
+  test(Tensor<2, 3, double>{{{1, 2, 0}, {0, 1, 0}, {0, 0, 1}}});
   // 2D not orthogonal
   test(Tensor<2, 2, double>{{{1, 2}, {3, 4}}});
   // 2D already orthogonal
index bf9984b2439149be8d1f15f544c5df3a23fdda2d..cec2818537fae9fe69a3b301fd18044b322f1470 100644 (file)
@@ -15,6 +15,14 @@ DEAL::projected to Q = 1.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 1
 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, 3, double>
+DEAL::original tensor A = 1.00000 2.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000
+DEAL::det A = 1.00000
+DEAL::A A^T = 5.00000 2.00000 0.00000 2.00000 1.00000 0.00000 0.00000 0.00000 1.00000
+DEAL::projected to Q = 0.707107 0.707107 0.00000 -0.707107 0.707107 0.00000 0.00000 0.00000 1.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

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.