]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix computations in SymmetricTensor::norm() for complex-valued tensors.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 8 Jul 2016 21:04:03 +0000 (16:04 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 8 Jul 2016 21:04:03 +0000 (16:04 -0500)
include/deal.II/base/symmetric_tensor.h

index c4fa23fe41a54f594601c39f7bd3185c9df75a2a..501c0ec94063cae37a2b1c72c8995885ef8852d0 100644 (file)
@@ -18,6 +18,7 @@
 
 
 #include <deal.II/base/tensor.h>
+#include <deal.II/base/numbers.h>
 #include <deal.II/base/table_indices.h>
 #include <deal.II/base/template_constraints.h>
 
@@ -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<Number>::real_type
+  norm () const;
 
   /**
    * Tensors can be unrolled by simply pasting all elements into one long
@@ -1749,70 +1751,77 @@ namespace internal
 {
   template <int dim, typename Number>
   inline
-  Number
+  typename numbers::NumberTraits<Number>::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<Number>::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<Number>::abs_square(data[0]) +
+                         numbers::NumberTraits<Number>::abs_square(data[1]) +
+                         2. * numbers::NumberTraits<Number>::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<Number>::abs_square(data[0]) +
+                         numbers::NumberTraits<Number>::abs_square(data[1]) +
+                         numbers::NumberTraits<Number>::abs_square(data[2]) +
+                         2. * numbers::NumberTraits<Number>::abs_square(data[3]) +
+                         2. * numbers::NumberTraits<Number>::abs_square(data[4]) +
+                         2. * numbers::NumberTraits<Number>::abs_square(data[5]));
+
       default:
-        return_value = Number();
+      {
+        typename numbers::NumberTraits<Number>::real_type return_value
+          = typename numbers::NumberTraits<Number>::real_type();
+
         for (unsigned int d=0; d<dim; ++d)
-          return_value += data[d] * data[d];
+          return_value += numbers::NumberTraits<Number>::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<Number>::abs_square(data[d]);
+
+        return std::sqrt(return_value);
+      }
       }
-    return return_value;
   }
 
 
 
   template <int dim, typename Number>
   inline
-  Number
+  typename numbers::NumberTraits<Number>::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<Number>::abs (data[0][0]);
+
       default:
-        return_value = Number();
+      {
+        typename numbers::NumberTraits<Number>::real_type return_value
+          = typename numbers::NumberTraits<Number>::real_type();
+
+        const unsigned int n_independent_components = data.dimension;
+
         for (unsigned int i=0; i<dim; ++i)
           for (unsigned int j=0; j<dim; ++j)
-            return_value += data[i][j] * data[i][j];
+            return_value += numbers::NumberTraits<Number>::abs_square(data[i][j]);
         for (unsigned int i=0; i<dim; ++i)
           for (unsigned int j=dim; j<n_independent_components; ++j)
-            return_value += Number(2.) * data[i][j] * data[i][j];
+            return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
         for (unsigned int i=dim; i<n_independent_components; ++i)
           for (unsigned int j=0; j<dim; ++j)
-            return_value += Number(2.) * data[i][j] * data[i][j];
+            return_value += 2. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
         for (unsigned int i=dim; i<n_independent_components; ++i)
           for (unsigned int j=dim; j<n_independent_components; ++j)
-            return_value += 4. * data[i][j] * data[i][j];
-        return_value = std::sqrt(return_value);
-      }
+            return_value += 4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
 
-    return return_value;
+        return std::sqrt(return_value);
+      }
+      }
   }
 
 } // end of namespace internal
@@ -1821,7 +1830,7 @@ namespace internal
 
 template <int rank, int dim, typename Number>
 inline
-Number
+typename numbers::NumberTraits<Number>::real_type
 SymmetricTensor<rank,dim,Number>::norm () const
 {
   return internal::compute_norm<dim,Number> (data);

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.