From 6e4a9eb759e1d5881f61c657c1c6ba26c414785d Mon Sep 17 00:00:00 2001 From: Reza Rastak Date: Mon, 4 May 2020 19:18:19 -0700 Subject: [PATCH] The tolerance parameter is removed from project_tensor --- include/deal.II/base/tensor.h | 43 ++++++++----------- tests/tensors/project_orthogonal.cc | 4 +- ...project_orthogonal.with_lapack=true.output | 8 ++++ 3 files changed, 28 insertions(+), 27 deletions(-) diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index 407f1a78a3..a8d6dd3157 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -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 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 matrix(dim); - LAPACKFullMatrix lapack_matrix(dim); - LAPACKFullMatrix result(dim); + Tensor<2, dim, Number> output_tensor; + FullMatrix matrix(dim); + LAPACKFullMatrix lapack_matrix(dim); + LAPACKFullMatrix 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; } diff --git a/tests/tensors/project_orthogonal.cc b/tests/tensors/project_orthogonal.cc index 463b836220..f62bcf5805 100644 --- a/tests/tensors/project_orthogonal.cc +++ b/tests/tensors/project_orthogonal.cc @@ -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 diff --git a/tests/tensors/project_orthogonal.with_lapack=true.output b/tests/tensors/project_orthogonal.with_lapack=true.output index bf9984b243..cec2818537 100644 --- a/tests/tensors/project_orthogonal.with_lapack=true.output +++ b/tests/tensors/project_orthogonal.with_lapack=true.output @@ -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 -- 2.39.5