From ae38b46b77f7d2dc2603001e080795b3c0da3c7d Mon Sep 17 00:00:00 2001 From: Menno Fraters Date: Fri, 28 Feb 2020 14:42:40 -0800 Subject: [PATCH] adress Wolfgangs comments. --- include/deal.II/base/tensor.h | 33 ++++++++++----------- tests/ad_common_tests/tensor_functions_01.h | 19 ++++++------ tests/tensors/constexpr_tensor.cc | 4 +-- 3 files changed, 27 insertions(+), 29 deletions(-) diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index 81a06ca789..dd4727d2f2 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -2618,31 +2618,29 @@ cofactor(const Tensor<2, dim, Number> &t) /** - * Return the nearest orthogonal matrix using a SVD if the - * deteriminant is more than a tolerance away from one. - * The orthogonalization is done by combining the products - * of the SVD decomposition: $((U*I)*V^T)^T$, where I is the - * idententy matrix and $U$ and $V$ are computed from the SVD - * decomposition: $\mathbf U \cdot \mathbf S \cdot \mathbf V^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 + * combining the products of the SVD decomposition: $U V^T$, where + * $U$ and $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 constexpr Tensor<2, dim, Number> -orthogonalize(const Tensor<2, dim, Number> &tensor, const double tolerance) +project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &tensor, + const double tolerance) { if (std::abs(determinant(tensor) - 1.0) > tolerance) { - LAPACKFullMatrix identity_matrix(dim); - for (size_t i = 0; i < dim; i++) - { - identity_matrix.set(i, i, 1.); - } - Tensor<2, dim, Number> output_tensor; FullMatrix matrix(dim); LAPACKFullMatrix lapack_matrix(dim); LAPACKFullMatrix result(dim); - LAPACKFullMatrix result2(dim); // todo: find or add dealii functionallity to copy in one step. matrix.copy_from(tensor); @@ -2651,12 +2649,11 @@ orthogonalize(const Tensor<2, dim, Number> &tensor, const double tolerance) // now compute the svd of the matrices lapack_matrix.compute_svd(); - // Use the SVD results to orthogonalize: ((U*I)*V^T)^T - lapack_matrix.get_svd_u().mmult(result, identity_matrix); - result.mmult(result2, (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 = result2; + matrix = result; matrix.copy_to(output_tensor); return output_tensor; } diff --git a/tests/ad_common_tests/tensor_functions_01.h b/tests/ad_common_tests/tensor_functions_01.h index 85f4fc853a..d0813fa81a 100644 --- a/tests/ad_common_tests/tensor_functions_01.h +++ b/tests/ad_common_tests/tensor_functions_01.h @@ -55,15 +55,16 @@ test_tensor() const Tensor<2, dim, ADNumberType> C5 = A * a; const Tensor<2, dim, ADNumberType> C6 = A / a; - const ADNumberType det_A = determinant(A); - const ADNumberType tr_A = trace(A); - const Tensor<2, dim, ADNumberType> A_inv = invert(A); - const Tensor<2, dim, ADNumberType> A_T = transpose(A); - const Tensor<2, dim, ADNumberType> A_adj = adjugate(A); - const Tensor<2, dim, ADNumberType> A_cof = cofactor(A); - const Tensor<2, dim, ADNumberType> A_orth = orthogonalize(A); - const ADNumberType A_l1_norm = l1_norm(A, 1e-8); - const ADNumberType A_linf_norm = linfty_norm(A); + const ADNumberType det_A = determinant(A); + const ADNumberType tr_A = trace(A); + const Tensor<2, dim, ADNumberType> A_inv = invert(A); + const Tensor<2, dim, ADNumberType> A_T = transpose(A); + const Tensor<2, dim, ADNumberType> A_adj = adjugate(A); + const Tensor<2, dim, ADNumberType> A_cof = cofactor(A); + const Tensor<2, dim, ADNumberType> A_orth = + project_onto_orthogonal_tensors(A); + const ADNumberType A_l1_norm = l1_norm(A, 1e-8); + const ADNumberType A_linf_norm = linfty_norm(A); const ADNumberType A_ddot_B = double_contract<0, 0, 1, 1>(A, B); const Tensor<2, dim, ADNumberType> A_dot_B = contract<1, 0>(A, B); diff --git a/tests/tensors/constexpr_tensor.cc b/tests/tensors/constexpr_tensor.cc index d1f2f0905e..11b026535f 100644 --- a/tests/tensors/constexpr_tensor.cc +++ b/tests/tensors/constexpr_tensor.cc @@ -181,7 +181,7 @@ main() DEAL_II_CONSTEXPR const auto dummy_8 = cofactor(a); deallog << "Deteriminant before orthogonalization: " << determinant(a) << std::endl; - const auto dummy_9 = orthogonalize(a, 1e-8); + const auto dummy_9 = project_onto_orthogonal_tensors(a, 1e-8); deallog << "Deteriminant after orthogonalization: " << determinant(dummy_9) << std::endl; Assert(determinant(dummy_9) - 1. < 1e-8, ExcInternalError()); @@ -194,7 +194,7 @@ main() deallog << "Deteriminant before orthogonalization: " << determinant(non_orthogonal) << std::endl; - const auto dummy_10 = orthogonalize(non_orthogonal, 1e-8); + const auto dummy_10 = project_onto_orthogonal_tensors(non_orthogonal, 1e-8); deallog << "Deteriminant after orthogonalization: " << determinant(dummy_10) << std::endl; Assert(determinant(dummy_10) - 1. < 1e-8, ExcInternalError()); -- 2.39.5