From: Mathias Mentler Date: Sun, 12 Feb 2017 10:00:15 +0000 (+0100) Subject: Removed use of constructor of Number in SymmetricTensor as this could X-Git-Tag: v8.5.0-rc1~111^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b8ce6632e2d56ea976cb5cf9e5cd23334da7c988;p=dealii.git Removed use of constructor of Number in SymmetricTensor as this could also call the non-existant constructor of VectorizedArray. Replaced by call to static member of template struct defined in vectorization.h and numbers.h --- diff --git a/doc/news/changes/minor/20160214MathiasMentler b/doc/news/changes/minor/20160214MathiasMentler new file mode 100644 index 0000000000..ade94a5a2f --- /dev/null +++ b/doc/news/changes/minor/20160214MathiasMentler @@ -0,0 +1,7 @@ +Fixed: Double contraction of two SymmetricTensor<..,VectorizedArray> now +works. Introduced internal::NumberType with static member +internal::Numbertype::value to be called instead of the constructor Number() +in symmetric_tensor.h. +
+(Mathias Mentler, 2016/02/14) + diff --git a/include/deal.II/base/numbers.h b/include/deal.II/base/numbers.h index 40fd69ccce..83d65d5c40 100644 --- a/include/deal.II/base/numbers.h +++ b/include/deal.II/base/numbers.h @@ -353,7 +353,36 @@ namespace numbers } +namespace internal +{ + /** + * The structs below are needed since VectorizedArray is a POD-type without a constructor and + * can be a template argument for SymmetricTensor<...,T2> where T2 would equal VectorizedArray. + * 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 + * breaks the POD-idioms needed elsewhere. Calls to constructors of T2 subsequently got replaced by a + * call to internal::NumberType 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 + struct NumberType + { + static T value (const T &t) + { + return t; + } + }; + template + struct NumberType > + { + static std::complex value (const T &t) + { + return std::complex(t); + } + }; +} DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index c1767a936a..90285b84f8 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -21,6 +21,7 @@ #include #include #include +#include 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::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::value(2.0); for (unsigned int d=0; d &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::value(2.0) * t.data[3]*t.data[4]*t.data[5] ); default: Assert (false, ExcNotImplemented()); return 0; diff --git a/include/deal.II/base/vectorization.h b/include/deal.II/base/vectorization.h index 7b1f874184..80004246a8 100644 --- a/include/deal.II/base/vectorization.h +++ b/include/deal.II/base/vectorization.h @@ -69,6 +69,31 @@ namespace std DEAL_II_NAMESPACE_OPEN +namespace internal +{ + /** + * The structs below are needed since VectorizedArray is a POD-type without a constructor and + * can be a template argument for SymmetricTensor<...,T2> where T2 would equal VectorizedArray. + * 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 + * breaks the POD-idioms needed elsewhere. Calls to constructors of T2 subsequently got replaced by a + * call to internal::NumberType 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 + struct NumberType > + { + static VectorizedArray value (const T &t) + { + VectorizedArray tmp; + tmp=t; + return tmp; + } + }; +} + + // Enable the EnableIfScalar type trait for VectorizedArray such // that it can be used as a Number type in Tensor, etc. diff --git a/tests/base/symmetric_tensor_39.cc b/tests/base/symmetric_tensor_39.cc new file mode 100644 index 0000000000..8de62a6045 --- /dev/null +++ b/tests/base/symmetric_tensor_39.cc @@ -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 +#include + +int main() +{ + std::string logname = "output"; + std::ofstream logfile(logname); + + const unsigned int dim = 3; + + SymmetricTensor<4,dim,VectorizedArray > I; + SymmetricTensor<2,dim,VectorizedArray > A; + SymmetricTensor<2,dim,VectorizedArray > 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::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 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::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 index 0000000000..d86bac9de5 --- /dev/null +++ b/tests/base/symmetric_tensor_39.output @@ -0,0 +1 @@ +OK diff --git a/tests/base/vectorization_07.cc b/tests/base/vectorization_07.cc new file mode 100644 index 0000000000..d75b88c01d --- /dev/null +++ b/tests/base/vectorization_07.cc @@ -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 + +#include + +int main() +{ + std::string logname = "output"; + std::ofstream logfile(logname); + + if (!std::is_pod >::value) + logfile << "Not OK because VectorizedArray is not POD!" << std::endl; + + if (!std::is_trivial >::value) + logfile << "Not OK because VectorizedArray 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 index 0000000000..d86bac9de5 --- /dev/null +++ b/tests/base/vectorization_07.with_cxx11=on.output @@ -0,0 +1 @@ +OK