From af61f22919e3ba4112e9ad1df07d227daac578ab Mon Sep 17 00:00:00 2001 From: Reza Rastak Date: Fri, 1 May 2020 13:36:17 -0700 Subject: [PATCH] Test for project_onto_orthogonal_tensors improved --- tests/tensors/constexpr_tensor.cc | 34 --------- tests/tensors/constexpr_tensor.output | 4 -- tests/tensors/project_orthogonal.cc | 69 +++++++++++++++++++ ...project_orthogonal.with_lapack=true.output | 50 ++++++++++++++ 4 files changed, 119 insertions(+), 38 deletions(-) create mode 100644 tests/tensors/project_orthogonal.cc create mode 100644 tests/tensors/project_orthogonal.with_lapack=true.output 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 -#include - #include "../tests.h" template @@ -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 unrolled_a(9); - a.unroll(unrolled_a); - Vector unrolled_dummy_9(9); - dummy_9.unroll(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 + +#include + +#include "../tests.h" + +template +void +print_type(const Tensor<2, dim, number> &) +{ + if (std::is_same::value) + deallog << " Tensor<2, " << dim << ", float>" << std::endl; + else if (std::is_same::value) + deallog << " Tensor<2, " << dim << ", double>" << std::endl; + else + Assert(false, ExcNotImplemented()); +} + +template +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 -- 2.39.5