From d9581d2eea4ec7b0f3b93aad4c9908c9fad53b78 Mon Sep 17 00:00:00 2001 From: Reza Rastak Date: Thu, 18 Jul 2019 11:48:17 -0700 Subject: [PATCH] constexpr_symmetric_tensor test improved --- doc/news/changes/minor/20190713RezaRastak | 2 +- include/deal.II/base/symmetric_tensor.h | 5 +- tests/tensors/constexpr_symmetric_tensor.cc | 80 ++++++++++++++----- .../tensors/constexpr_symmetric_tensor.output | 13 +-- 4 files changed, 72 insertions(+), 28 deletions(-) diff --git a/doc/news/changes/minor/20190713RezaRastak b/doc/news/changes/minor/20190713RezaRastak index 0ecbf23d08..96207401e4 100644 --- a/doc/news/changes/minor/20190713RezaRastak +++ b/doc/news/changes/minor/20190713RezaRastak @@ -1,3 +1,3 @@ -Improved: constexpr evaluation is enabled in SymmetricTensor and TableIndices +Improved: constexpr evaluation is enabled in SymmetricTensor and TableIndices.
(Reza Rastak, 2019/07/13) \ No newline at end of file diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index 5be47f7cab..e8bd996113 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -1800,9 +1800,8 @@ namespace internal default: // to do the rest, sort our indices before comparing { - TableIndices<2> sorted_indices(indices); - sorted_indices.sort(); - + TableIndices<2> sorted_indices(std::min(indices[0], indices[1]), + std::max(indices[0], indices[1])); for (unsigned int d = 0, c = 0; d < dim; ++d) for (unsigned int e = d + 1; e < dim; ++e, ++c) if ((sorted_indices[0] == d) && (sorted_indices[1] == e)) diff --git a/tests/tensors/constexpr_symmetric_tensor.cc b/tests/tensors/constexpr_symmetric_tensor.cc index e34d0d130d..6d5fca5ac7 100644 --- a/tests/tensors/constexpr_symmetric_tensor.cc +++ b/tests/tensors/constexpr_symmetric_tensor.cc @@ -16,9 +16,54 @@ // create and manipulate constexpr SymmetricTensor objects #include +#include #include "../tests.h" +template +void +test_symmetric_tensor() +{ + deallog << "*** Test constexpr SymmetricTensor functions, " + << "dim = " << Utilities::to_string(dim) << std::endl; + + constexpr Number a = 1.0; + constexpr Tensor<1, dim, Number> v{}; + DEAL_II_CONSTEXPR const SymmetricTensor<2, dim, Number> A( + unit_symmetric_tensor()); + DEAL_II_CONSTEXPR const SymmetricTensor<2, dim, Number> B( + unit_symmetric_tensor()); + DEAL_II_CONSTEXPR const Tensor<2, dim, Number> A_ns( + unit_symmetric_tensor()); + DEAL_II_CONSTEXPR const SymmetricTensor<4, dim, Number> HH( + identity_tensor()); + + DEAL_II_CONSTEXPR const SymmetricTensor<2, dim, Number> C1 = A + B; + DEAL_II_CONSTEXPR const SymmetricTensor<2, dim, Number> C2 = A - B; + DEAL_II_CONSTEXPR const SymmetricTensor<2, dim, Number> C4 = a * A; + DEAL_II_CONSTEXPR const SymmetricTensor<2, dim, Number> C5 = A * a; + DEAL_II_CONSTEXPR const SymmetricTensor<2, dim, Number> C6 = A / a; + + DEAL_II_CONSTEXPR const Number det_A = determinant(A); + DEAL_II_CONSTEXPR const Number tr_A = trace(A); + DEAL_II_CONSTEXPR const Number I1_A = first_invariant(A); + DEAL_II_CONSTEXPR const Number I2_A = second_invariant(A); + DEAL_II_CONSTEXPR const Number I3_A = third_invariant(A); + + DEAL_II_CONSTEXPR const SymmetricTensor<2, dim, Number> A_inv = invert(A); + DEAL_II_CONSTEXPR const SymmetricTensor<2, dim, Number> A_T = transpose(A); + DEAL_II_CONSTEXPR const SymmetricTensor<2, dim, Number> A_dev = deviator(A); + + DEAL_II_CONSTEXPR const Number A_ddot_B_1 = A * B; + DEAL_II_CONSTEXPR const Number sp_A_B = scalar_product(A, B); + DEAL_II_CONSTEXPR const Tensor<4, dim, Number> op_A_B = outer_product(A, B); + + DEAL_II_CONSTEXPR const Tensor<1, dim, Number> v3 = A * v; + DEAL_II_CONSTEXPR const Tensor<1, dim, Number> v4 = v * A; + DEAL_II_CONSTEXPR const Tensor<2, dim, Number> C7 = A * A_ns; + DEAL_II_CONSTEXPR const Tensor<2, dim, Number> C8 = A_ns * A; +} + DEAL_II_CONSTEXPR SymmetricTensor<2, 2> get_tensor_2() { @@ -46,29 +91,26 @@ main() { initlog(); - DEAL_II_CONSTEXPR const auto A = get_tensor_2(); { - LogStream::Prefix p("SymmetricTensor<2,2>"); - DEAL_II_CONSTEXPR const double calculation = A[0][0] * A[1][0] - A[1][1]; - DEAL_II_CONSTEXPR const double invariants[] = {first_invariant(A), - second_invariant(A), - third_invariant(A)}; - DEAL_II_CONSTEXPR const double norm_square = scalar_product(A, A); - deallog << "calculation result = " << calculation << std::endl; - deallog << "invariants = " << invariants[0] << ", " << invariants[1] << ", " - << invariants[2] << std::endl; - deallog << "norm square = " << norm_square << std::endl; + LogStream::Prefix p("float"); + test_symmetric_tensor<1, float>(); + test_symmetric_tensor<2, float>(); + test_symmetric_tensor<3, float>(); } - DEAL_II_CONSTEXPR const auto B = get_tensor_4(); { - LogStream::Prefix p("SymmetricTensor<4,2>"); - DEAL_II_CONSTEXPR const double calculation = - B[0][0][0][0] * B[1][0][0][1] - B[1][1][0][0]; - DEAL_II_CONSTEXPR const auto B_times_A = B * A; - deallog << "calculation result = " << calculation << std::endl; - deallog << "B times A (0,0) = " << B_times_A[0][0] << std::endl; + LogStream::Prefix p("double"); + test_symmetric_tensor<1, double>(); + test_symmetric_tensor<2, double>(); + test_symmetric_tensor<3, double>(); } + + DEAL_II_CONSTEXPR const auto A = get_tensor_2(); + deallog << "SymmetricTensor<2,2> = " << A << std::endl; + + DEAL_II_CONSTEXPR const auto B = get_tensor_4(); + deallog << "SymmetricTensor<4,2> = " << B << std::endl; + deallog << "OK" << std::endl; return 0; -} \ No newline at end of file +} diff --git a/tests/tensors/constexpr_symmetric_tensor.output b/tests/tensors/constexpr_symmetric_tensor.output index 6f69b6bb62..6920777316 100644 --- a/tests/tensors/constexpr_symmetric_tensor.output +++ b/tests/tensors/constexpr_symmetric_tensor.output @@ -1,7 +1,10 @@ -DEAL:SymmetricTensor<2,2>::calculation result = -8.00000 -DEAL:SymmetricTensor<2,2>::invariants = 4.00000, -22.0000, -22.0000 -DEAL:SymmetricTensor<2,2>::norm square = 60.0000 -DEAL:SymmetricTensor<4,2>::calculation result = 0.200000 -DEAL:SymmetricTensor<4,2>::B times A (0,0) = 1.00000 +DEAL:float::*** Test constexpr SymmetricTensor functions, dim = 1 +DEAL:float::*** Test constexpr SymmetricTensor functions, dim = 2 +DEAL:float::*** Test constexpr SymmetricTensor functions, dim = 3 +DEAL:double::*** Test constexpr SymmetricTensor functions, dim = 1 +DEAL:double::*** Test constexpr SymmetricTensor functions, dim = 2 +DEAL:double::*** Test constexpr SymmetricTensor functions, dim = 3 +DEAL::SymmetricTensor<2,2> = 1.00000 -5.00000 -5.00000 3.00000 +DEAL::SymmetricTensor<4,2> = 1.00000 0.00000 0.00000 0.00000 0.00000 0.200000 0.200000 0.00000 0.00000 0.200000 0.200000 0.00000 0.00000 0.00000 0.00000 2.50000 DEAL::OK -- 2.39.5