From: Wolfgang Bangerth Date: Fri, 8 Jul 2016 21:04:03 +0000 (-0500) Subject: Fix computations in SymmetricTensor::norm() for complex-valued tensors. X-Git-Tag: v8.5.0-rc1~896^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f1bbe281185026973fa03ed69f3996f2ad1daa21;p=dealii.git Fix computations in SymmetricTensor::norm() for complex-valued tensors. --- diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index c4fa23fe41..501c0ec940 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -18,6 +18,7 @@ #include +#include #include #include @@ -727,7 +728,8 @@ public: * upper right as well as lower left entries, not just one of them, although * they are equal for symmetric tensors). */ - Number norm () const; + typename numbers::NumberTraits::real_type + norm () const; /** * Tensors can be unrolled by simply pasting all elements into one long @@ -1749,70 +1751,77 @@ namespace internal { template inline - Number + typename numbers::NumberTraits::real_type compute_norm (const typename SymmetricTensorAccessors::StorageType<2,dim,Number>::base_tensor_type &data) { - Number return_value; switch (dim) { case 1: - return_value = std::abs(data[0]); - break; + return numbers::NumberTraits::abs(data[0]); + case 2: - return_value = std::sqrt(data[0]*data[0] + data[1]*data[1] + - Number(2.) * data[2]*data[2]); - break; + return std::sqrt(numbers::NumberTraits::abs_square(data[0]) + + numbers::NumberTraits::abs_square(data[1]) + + 2. * numbers::NumberTraits::abs_square(data[2])); + case 3: - return_value = std::sqrt(data[0]*data[0] + data[1]*data[1] + - data[2]*data[2] + - Number(2.) * data[3]*data[3] + - Number(2.) * data[4]*data[4] + - Number(2.) * data[5]*data[5]); - break; + return std::sqrt(numbers::NumberTraits::abs_square(data[0]) + + numbers::NumberTraits::abs_square(data[1]) + + numbers::NumberTraits::abs_square(data[2]) + + 2. * numbers::NumberTraits::abs_square(data[3]) + + 2. * numbers::NumberTraits::abs_square(data[4]) + + 2. * numbers::NumberTraits::abs_square(data[5])); + default: - return_value = Number(); + { + typename numbers::NumberTraits::real_type return_value + = typename numbers::NumberTraits::real_type(); + for (unsigned int d=0; d::abs_square(data[d]); for (unsigned int d=dim; d<(dim*dim+dim)/2; ++d) - return_value += Number(2.) * data[d] * data[d]; - return_value = std::sqrt(return_value); + return_value += 2. * numbers::NumberTraits::abs_square(data[d]); + + return std::sqrt(return_value); + } } - return return_value; } template inline - Number + typename numbers::NumberTraits::real_type compute_norm (const typename SymmetricTensorAccessors::StorageType<4,dim,Number>::base_tensor_type &data) { - Number return_value; - const unsigned int n_independent_components = data.dimension; - switch (dim) { case 1: - return_value = std::abs (data[0][0]); - break; + return numbers::NumberTraits::abs (data[0][0]); + default: - return_value = Number(); + { + typename numbers::NumberTraits::real_type return_value + = typename numbers::NumberTraits::real_type(); + + const unsigned int n_independent_components = data.dimension; + for (unsigned int i=0; i::abs_square(data[i][j]); for (unsigned int i=0; i::abs_square(data[i][j]); for (unsigned int i=dim; i::abs_square(data[i][j]); for (unsigned int i=dim; i::abs_square(data[i][j]); - return return_value; + return std::sqrt(return_value); + } + } } } // end of namespace internal @@ -1821,7 +1830,7 @@ namespace internal template inline -Number +typename numbers::NumberTraits::real_type SymmetricTensor::norm () const { return internal::compute_norm (data);