]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Removed use of constructor of Number in SymmetricTensor as this could 3967/head
authorMathias Mentler <mathiasmentler@web.de>
Sun, 12 Feb 2017 10:00:15 +0000 (11:00 +0100)
committerMathias Mentler <mathiasmentler@web.de>
Tue, 14 Feb 2017 20:58:25 +0000 (21:58 +0100)
also call the non-existant constructor of VectorizedArray. Replaced
by call to static member of template struct defined in
vectorization.h and numbers.h

doc/news/changes/minor/20160214MathiasMentler [new file with mode: 0644]
include/deal.II/base/numbers.h
include/deal.II/base/symmetric_tensor.h
include/deal.II/base/vectorization.h
tests/base/symmetric_tensor_39.cc [new file with mode: 0644]
tests/base/symmetric_tensor_39.output [new file with mode: 0644]
tests/base/vectorization_07.cc [new file with mode: 0644]
tests/base/vectorization_07.with_cxx11=on.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20160214MathiasMentler b/doc/news/changes/minor/20160214MathiasMentler
new file mode 100644 (file)
index 0000000..ade94a5
--- /dev/null
@@ -0,0 +1,7 @@
+Fixed: Double contraction of two SymmetricTensor<..,VectorizedArray<T>> now
+works. Introduced internal::NumberType<T> with static member 
+internal::Numbertype::value to be called instead of the constructor Number()
+in symmetric_tensor.h.
+<br>
+(Mathias Mentler, 2016/02/14)
+
index 40fd69ccce1399807b43c90783556616782f6cb2..83d65d5c40a124a1ea8037648e7ec06e7d7078c1 100644 (file)
@@ -353,7 +353,36 @@ namespace numbers
 
 }
 
+namespace internal
+{
+  /**
+  * The structs below are needed since VectorizedArray<T1> is a POD-type without a constructor and
+  * can be a template argument for SymmetricTensor<...,T2> where T2 would equal VectorizedArray<T1>.
+  * Internally, in previous versions of deal.II, SymmetricTensor<...,T2> would make use of the constructor
+  * of T2 leading to a compile-time error. However simply adding a constructor for VectorizedArray<T1>
+  * breaks the POD-idioms needed elsewhere. Calls to constructors of T2 subsequently got replaced by a
+  * call to internal::NumberType<T2> which then determines the right function to use by template deduction.
+  * A detailled discussion can be found at https://github.com/dealii/dealii/pull/3967 . Also see
+  * numbers.h for another specialization.
+  */
+  template <typename T>
+  struct NumberType
+  {
+    static T value (const T &t)
+    {
+      return t;
+    }
+  };
 
+  template <typename T>
+  struct NumberType<std::complex<T> >
+  {
+    static std::complex<T> value (const T &t)
+    {
+      return std::complex<T>(t);
+    }
+  };
+}
 
 DEAL_II_NAMESPACE_CLOSE
 
index c1767a936a0528d394b70430dfe7b97fc4075086..90285b84f806f4b7d762716ae463041f96db67e1 100644 (file)
@@ -21,6 +21,7 @@
 #include <deal.II/base/numbers.h>
 #include <deal.II/base/table_indices.h>
 #include <deal.II/base/template_constraints.h>
+#include <deal.II/base/vectorization.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -1217,14 +1218,14 @@ namespace internal
       case 2:
         return (data[0] * sdata[0] +
                 data[1] * sdata[1] +
-                Number(2.) * data[2] * sdata[2]);
+                NumberType<Number>::value(2.0) * data[2] * sdata[2]);
       default:
         // Start with the non-diagonal part to avoid some multiplications by
         // 2.
         Number sum = data[dim] * sdata[dim];
         for (unsigned int d=dim+1; d<(dim*(dim+1)/2); ++d)
           sum += data[d] * sdata[d];
-        sum *= 2.;
+        sum *= NumberType<Number>::value(2.0);
         for (unsigned int d=0; d<dim; ++d)
           sum += data[d] * sdata[d];
         return sum;
@@ -2165,7 +2166,7 @@ Number determinant (const SymmetricTensor<2,dim,Number> &t)
                -t.data[0]*t.data[5]*t.data[5]
                -t.data[1]*t.data[4]*t.data[4]
                -t.data[2]*t.data[3]*t.data[3]
-               +Number(2.) * t.data[3]*t.data[4]*t.data[5] );
+               +internal::NumberType<Number>::value(2.0) * t.data[3]*t.data[4]*t.data[5] );
     default:
       Assert (false, ExcNotImplemented());
       return 0;
index 7b1f874184be2a14807b34c59578471b64fc953b..80004246a85c1d0dbf6c39343636c79365c313dd 100644 (file)
@@ -69,6 +69,31 @@ namespace std
 DEAL_II_NAMESPACE_OPEN
 
 
+namespace internal
+{
+  /**
+  * The structs below are needed since VectorizedArray<T1> is a POD-type without a constructor and
+  * can be a template argument for SymmetricTensor<...,T2> where T2 would equal VectorizedArray<T1>.
+  * Internally, in previous versions of deal.II, SymmetricTensor<...,T2> would make use of the constructor
+  * of T2 leading to a compile-time error. However simply adding a constructor for VectorizedArray<T1>
+  * breaks the POD-idioms needed elsewhere. Calls to constructors of T2 subsequently got replaced by a
+  * call to internal::NumberType<T2> which then determines the right function to use by template deduction.
+  * A detailled discussion can be found at https://github.com/dealii/dealii/pull/3967 . Also see numbers.h
+  * for other specializations.
+  */
+  template <typename T>
+  struct NumberType<VectorizedArray<T> >
+  {
+    static VectorizedArray<T> value (const T &t)
+    {
+      VectorizedArray<T> tmp;
+      tmp=t;
+      return tmp;
+    }
+  };
+}
+
+
 // Enable the EnableIfScalar type trait for VectorizedArray<Number> such
 // that it can be used as a Number type in Tensor<rank,dim,Number>, etc.
 
diff --git a/tests/base/symmetric_tensor_39.cc b/tests/base/symmetric_tensor_39.cc
new file mode 100644 (file)
index 0000000..8de62a6
--- /dev/null
@@ -0,0 +1,61 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2015 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include "../tests.h"
+
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/vectorization.h>
+
+int main()
+{
+  std::string logname = "output";
+  std::ofstream logfile(logname);
+
+  const unsigned int dim = 3;
+
+  SymmetricTensor<4,dim,VectorizedArray<double> > I;
+  SymmetricTensor<2,dim,VectorizedArray<double> > A;
+  SymmetricTensor<2,dim,VectorizedArray<double> > B;
+
+  // I^sym = 0.5(d_ik*d_jl + d_il*d_jk) -> I^sym : A = A^sym
+  for (unsigned int i = 0; i < dim; i++)
+    for (unsigned int j = 0; j < dim; j++)
+      for (unsigned int k = 0; k < dim; k++)
+        for (unsigned int l = 0; l < dim; l++)
+          I[i][j][k][l] = ( (i == k && j== l && i == l && j == k) ? make_vectorized_array(1.0) : ( (i == k && j== l) || (i == l && j == k) ? make_vectorized_array(0.5) : make_vectorized_array(0.0)));
+
+  double counter = 0.0;
+  for (unsigned int i = 0; i < dim; ++i)
+    for (unsigned int j = 0; j < dim; ++j)
+      for (unsigned int v = 0; v < VectorizedArray<double>::n_array_elements; v++)
+        {
+          A[i][j][v] = counter;
+          counter += 1.0;
+        }
+
+  B = I*A;
+  B -= A;
+
+  // Note that you cannot use B.norm() here even with something
+  // like VectorizedArray<double> frob_norm = B.norm() -> Maybe a TODO?
+  for (unsigned int i = 0; i < dim; ++i)
+    for (unsigned int j = 0; j < dim; ++j)
+      for (unsigned int v = 0; v < VectorizedArray<double>::n_array_elements; ++v)
+        if (B[i][j][v] != 0.0)
+          logfile << "Not OK" << std::endl;
+
+  logfile << "OK" << std::endl;
+  logfile.close();
+}
diff --git a/tests/base/symmetric_tensor_39.output b/tests/base/symmetric_tensor_39.output
new file mode 100644 (file)
index 0000000..d86bac9
--- /dev/null
@@ -0,0 +1 @@
+OK
diff --git a/tests/base/vectorization_07.cc b/tests/base/vectorization_07.cc
new file mode 100644 (file)
index 0000000..d75b88c
--- /dev/null
@@ -0,0 +1,34 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2015 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include "../tests.h"
+#include <type_traits>
+
+#include <deal.II/base/vectorization.h>
+
+int main()
+{
+  std::string logname = "output";
+  std::ofstream logfile(logname);
+
+  if (!std::is_pod<VectorizedArray<double> >::value)
+    logfile << "Not OK because VectorizedArray<double> is not POD!" << std::endl;
+
+  if (!std::is_trivial<VectorizedArray<double> >::value)
+    logfile << "Not OK because VectorizedArray<double> is not trivial!" << std::endl;
+
+  logfile << "OK" << std::endl;
+  logfile.close();
+}
diff --git a/tests/base/vectorization_07.with_cxx11=on.output b/tests/base/vectorization_07.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..d86bac9
--- /dev/null
@@ -0,0 +1 @@
+OK

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.