--- /dev/null
+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)
+
}
+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
#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
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;
-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;
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.
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+}
projects / dealii.git / commitdiff