From ebb96b1ac04b944ea03949a9149ef43c8e5fcce6 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Thu, 24 Aug 2017 17:21:51 +0200 Subject: [PATCH] Reimplement ProductType do better deal with qualified number types. With reference to #4950. --- .../20170824Jean-PaulPelteretWolfgangBangerth | 6 + include/deal.II/base/template_constraints.h | 103 ++++++++++++------ .../differentiation/ad/sacado_product_types.h | 81 +++++++------- include/deal.II/fe/fe_values.h | 34 +++--- source/fe/fe_values.cc | 40 ++++--- 5 files changed, 155 insertions(+), 109 deletions(-) create mode 100644 doc/news/changes/incompatibilities/20170824Jean-PaulPelteretWolfgangBangerth diff --git a/doc/news/changes/incompatibilities/20170824Jean-PaulPelteretWolfgangBangerth b/doc/news/changes/incompatibilities/20170824Jean-PaulPelteretWolfgangBangerth new file mode 100644 index 0000000000..850f9c0a3b --- /dev/null +++ b/doc/news/changes/incompatibilities/20170824Jean-PaulPelteretWolfgangBangerth @@ -0,0 +1,6 @@ +Changed: Specialization of the ProductType class are now implemented through +specialization of the class internal::ProductTypeImpl . This was done in order +to ensure that product operations performed with qualified number types do not +result in the intended specializations being overlooked by the compiler. +
+(Jean-Paul Pelteret, Wolfgang Bangerth, 2017/08/24) diff --git a/include/deal.II/base/template_constraints.h b/include/deal.II/base/template_constraints.h index 6171fd5ab3..5f2a01149d 100644 --- a/include/deal.II/base/template_constraints.h +++ b/include/deal.II/base/template_constraints.h @@ -326,6 +326,31 @@ struct types_are_equal : std::is_same +namespace internal +{ + + /** + * A struct that implements the default product type resulting from the + * multiplication of two types. + * + * @note Care should be taken when @p T or @p U have qualifiers (@p const or + * @p volatile) or are @p lvalue or @p rvalue references! It is recommended + * that specialization of this class is only made for unqualified (fully + * stripped) types and that the ProductType class be used to determine the + * result of operating with (potentially) qualified types. + * + * @author Wolfgang Bangerth, Jean-Paul Pelteret, 2017 + */ + template + struct ProductTypeImpl + { + typedef decltype(std::declval() * std::declval()) type; + }; + +} + + + /** * A class with a local typedef that represents the type that results from the * product of two variables of type @p T and @p U. In other words, we would @@ -372,54 +397,60 @@ struct types_are_equal : std::is_same * used for the result of computing the product of unknowns and the values, * gradients, or other properties of shape functions. * - * @author Wolfgang Bangerth, 2015 + * @author Wolfgang Bangerth, 2015, 2017 */ template struct ProductType { - typedef decltype(std::declval() * std::declval()) type; + typedef typename internal::ProductTypeImpl< + typename std::decay::type, typename std::decay::type>::type type; }; -// Annoyingly, there is no std::complex::operator*(U) for scalars U -// other than T (not even in C++11, or C++14). We provide our own overloads -// in base/complex_overloads.h, but in order for them to work, we have to -// manually specify all products we want to allow: - -template -struct ProductType,std::complex > +namespace internal { - typedef std::complex type; -}; -template -struct ProductType,std::complex > -{ - typedef std::complex::type> type; -}; + // Annoyingly, there is no std::complex::operator*(U) for scalars U + // other than T (not even in C++11, or C++14). We provide our own overloads + // in base/complex_overloads.h, but in order for them to work, we have to + // manually specify all products we want to allow: -template -struct ProductType > -{ - typedef std::complex::type> type; -}; + template + struct ProductTypeImpl,std::complex > + { + typedef std::complex type; + }; -template -struct ProductType,double> -{ - typedef std::complex::type> type; -}; + template + struct ProductTypeImpl,std::complex > + { + typedef std::complex::type> type; + }; -template -struct ProductType > -{ - typedef std::complex::type> type; -}; + template + struct ProductTypeImpl > + { + typedef std::complex::type> type; + }; -template -struct ProductType,float> -{ - typedef std::complex::type> type; -}; + template + struct ProductTypeImpl,double> + { + typedef std::complex::type> type; + }; + + template + struct ProductTypeImpl > + { + typedef std::complex::type> type; + }; + + template + struct ProductTypeImpl,float> + { + typedef std::complex::type> type; + }; + +} diff --git a/include/deal.II/differentiation/ad/sacado_product_types.h b/include/deal.II/differentiation/ad/sacado_product_types.h index 54de3762b1..0ff6ce32fc 100644 --- a/include/deal.II/differentiation/ad/sacado_product_types.h +++ b/include/deal.II/differentiation/ad/sacado_product_types.h @@ -29,47 +29,52 @@ DEAL_II_ENABLE_EXTRA_DIAGNOSTICS DEAL_II_NAMESPACE_OPEN -template -struct ProductType, float> -{ - typedef Sacado::Fad::DFad type; -}; - -template -struct ProductType > -{ - typedef Sacado::Fad::DFad type; -}; - -template -struct ProductType, double> -{ - typedef Sacado::Fad::DFad type; -}; - -template -struct ProductType > +namespace internal { - typedef Sacado::Fad::DFad type; -}; -template -struct ProductType, int> -{ - typedef Sacado::Fad::DFad type; -}; + template + struct ProductTypeImpl, float> + { + typedef Sacado::Fad::DFad type; + }; + + template + struct ProductTypeImpl > + { + typedef Sacado::Fad::DFad type; + }; + + template + struct ProductTypeImpl, double> + { + typedef Sacado::Fad::DFad type; + }; + + template + struct ProductTypeImpl > + { + typedef Sacado::Fad::DFad type; + }; + + template + struct ProductTypeImpl, int> + { + typedef Sacado::Fad::DFad type; + }; + + template + struct ProductTypeImpl > + { + typedef Sacado::Fad::DFad type; + }; + + template + struct ProductTypeImpl, Sacado::Fad::DFad > + { + typedef Sacado::Fad::DFad::type > type; + }; -template -struct ProductType > -{ - typedef Sacado::Fad::DFad type; -}; - -template -struct ProductType, Sacado::Fad::DFad > -{ - typedef Sacado::Fad::DFad::type > type; -}; +} template struct EnableIfScalar > diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h index 3839ea98d0..6faeb04ade 100644 --- a/include/deal.II/fe/fe_values.h +++ b/include/deal.II/fe/fe_values.h @@ -178,31 +178,31 @@ namespace FEValuesViews * A typedef for the data type of the product of a @p Number and the * values of the view the Scalar class. */ - typedef typename ProductType::type, typename Scalar::value_type>::type value_type; + typedef typename ProductType::value_type>::type value_type; /** * A typedef for the data type of the product of a @p Number and the * gradients of the view the Scalar class. */ - typedef typename ProductType::type, typename Scalar::gradient_type>::type gradient_type; + typedef typename ProductType::gradient_type>::type gradient_type; /** * A typedef for the data type of the product of a @p Number and the * laplacians of the view the Scalar class. */ - typedef typename ProductType::type, typename Scalar::value_type>::type laplacian_type; + typedef typename ProductType::value_type>::type laplacian_type; /** * A typedef for the data type of the product of a @p Number and the * hessians of the view the Scalar class. */ - typedef typename ProductType::type, typename Scalar::hessian_type>::type hessian_type; + typedef typename ProductType::hessian_type>::type hessian_type; /** * A typedef for the data type of the product of a @p Number and the * third derivatives of the view the Scalar class. */ - typedef typename ProductType::type, typename Scalar::third_derivative_type>::type third_derivative_type; + typedef typename ProductType::third_derivative_type>::type third_derivative_type; }; /** @@ -595,49 +595,49 @@ namespace FEValuesViews * A typedef for the data type of the product of a @p Number and the * values of the view the Vector class. */ - typedef typename ProductType::type, typename Vector::value_type>::type value_type; + typedef typename ProductType::value_type>::type value_type; /** * A typedef for the data type of the product of a @p Number and the * gradients of the view the Vector class. */ - typedef typename ProductType::type, typename Vector::gradient_type>::type gradient_type; + typedef typename ProductType::gradient_type>::type gradient_type; /** * A typedef for the data type of the product of a @p Number and the * symmetric gradients of the view the Vector class. */ - typedef typename ProductType::type, typename Vector::symmetric_gradient_type>::type symmetric_gradient_type; + typedef typename ProductType::symmetric_gradient_type>::type symmetric_gradient_type; /** * A typedef for the data type of the product of a @p Number and the * divergences of the view the Vector class. */ - typedef typename ProductType::type, typename Vector::divergence_type>::type divergence_type; + typedef typename ProductType::divergence_type>::type divergence_type; /** * A typedef for the data type of the product of a @p Number and the * laplacians of the view the Vector class. */ - typedef typename ProductType::type, typename Vector::value_type>::type laplacian_type; + typedef typename ProductType::value_type>::type laplacian_type; /** * A typedef for the data type of the product of a @p Number and the * curls of the view the Vector class. */ - typedef typename ProductType::type, typename Vector::curl_type>::type curl_type; + typedef typename ProductType::curl_type>::type curl_type; /** * A typedef for the data type of the product of a @p Number and the * hessians of the view the Vector class. */ - typedef typename ProductType::type, typename Vector::hessian_type>::type hessian_type; + typedef typename ProductType::hessian_type>::type hessian_type; /** * A typedef for the data type of the product of a @p Number and the * third derivatives of the view the Vector class. */ - typedef typename ProductType::type, typename Vector::third_derivative_type>::type third_derivative_type; + typedef typename ProductType::third_derivative_type>::type third_derivative_type; }; /** @@ -1157,13 +1157,13 @@ namespace FEValuesViews * A typedef for the data type of the product of a @p Number and the * values of the view the SymmetricTensor class. */ - typedef typename ProductType::type, typename SymmetricTensor<2,dim,spacedim>::value_type>::type value_type; + typedef typename ProductType::value_type>::type value_type; /** * A typedef for the data type of the product of a @p Number and the * divergences of the view the SymmetricTensor class. */ - typedef typename ProductType::type, typename SymmetricTensor<2,dim,spacedim>::divergence_type>::type divergence_type; + typedef typename ProductType::divergence_type>::type divergence_type; }; /** @@ -1413,13 +1413,13 @@ namespace FEValuesViews * A typedef for the data type of the product of a @p Number and the * values of the view the Tensor class. */ - typedef typename ProductType::type, typename Tensor<2,dim,spacedim>::value_type>::type value_type; + typedef typename ProductType::value_type>::type value_type; /** * A typedef for the data type of the product of a @p Number and the * divergences of the view the Tensor class. */ - typedef typename ProductType::type, typename Tensor<2,dim,spacedim>::divergence_type>::type divergence_type; + typedef typename ProductType::divergence_type>::type divergence_type; }; /** diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index 0e2ac8ae79..f1ae74645b 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -462,14 +462,15 @@ namespace FEValuesViews do_function_values (const ArrayView &dof_values, const Table<2,double> &shape_values, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector::type,double>::type> &values) + std::vector::type> &values) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? shape_values.n_cols() : values.size(); AssertDimension (values.size(), n_quadrature_points); - std::fill (values.begin(), values.end(), dealii::internal::NumberType::type>::value(0.0)); + std::fill (values.begin(), values.end(), + dealii::internal::NumberType::type>::value(0.0)); for (unsigned int shape_function=0; shape_function &dof_values, const Table<2,dealii::Tensor > &shape_derivatives, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector::type,dealii::Tensor >::type> &derivatives) + std::vector >::type> &derivatives) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? @@ -503,7 +504,7 @@ namespace FEValuesViews AssertDimension (derivatives.size(), n_quadrature_points); std::fill (derivatives.begin(), derivatives.end(), - typename ProductType::type,dealii::Tensor >::type()); + typename ProductType >::type()); for (unsigned int shape_function=0; shape_function &dof_values, const Table<2,double> &shape_values, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector::type,dealii::Tensor<1,spacedim> >::type> &values) + std::vector >::type> &values) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? shape_values.n_cols() : values.size(); AssertDimension (values.size(), n_quadrature_points); - std::fill (values.begin(), values.end(), typename ProductType::type,dealii::Tensor<1,spacedim> >::type()); + std::fill (values.begin(), values.end(), + typename ProductType >::type()); for (unsigned int shape_function=0; shape_function &dof_values, const Table<2,dealii::Tensor > &shape_derivatives, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector::type,dealii::Tensor >::type> &derivatives) + std::vector >::type> &derivatives) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? @@ -618,7 +620,7 @@ namespace FEValuesViews AssertDimension (derivatives.size(), n_quadrature_points); std::fill (derivatives.begin(), derivatives.end(), - typename ProductType::type,dealii::Tensor >::type()); + typename ProductType >::type()); for (unsigned int shape_function=0; shape_function &dof_values, const Table<2,dealii::Tensor<1,spacedim> > &shape_gradients, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector::type,dealii::SymmetricTensor<2,spacedim> >::type> &symmetric_gradients) + std::vector >::type> &symmetric_gradients) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? @@ -672,7 +674,7 @@ namespace FEValuesViews AssertDimension (symmetric_gradients.size(), n_quadrature_points); std::fill (symmetric_gradients.begin(), symmetric_gradients.end(), - typename ProductType::type,dealii::SymmetricTensor<2,spacedim> >::type()); + typename ProductType >::type()); for (unsigned int shape_function=0; shape_function &dof_values, const Table<2,dealii::Tensor<1,spacedim> > &shape_gradients, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector::type,typename dealii::internal::CurlType::type>::type> &curls) + std::vector::type>::type> &curls) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? @@ -776,7 +778,7 @@ namespace FEValuesViews AssertDimension (curls.size(), n_quadrature_points); std::fill (curls.begin(), curls.end(), - typename ProductType::type,typename dealii::internal::CurlType::type>::type()); + typename ProductType::type>::type()); switch (spacedim) { @@ -1015,7 +1017,7 @@ namespace FEValuesViews do_function_values (const ArrayView &dof_values, const dealii::Table<2,double> &shape_values, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector::type,dealii::SymmetricTensor<2,spacedim> >::type> &values) + std::vector >::type> &values) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? @@ -1023,7 +1025,7 @@ namespace FEValuesViews AssertDimension (values.size(), n_quadrature_points); std::fill (values.begin(), values.end(), - typename ProductType::type,dealii::SymmetricTensor<2,spacedim> >::type()); + typename ProductType >::type()); for (unsigned int shape_function=0; shape_function &dof_values, const dealii::Table<2,double> &shape_values, const std::vector::ShapeFunctionData> &shape_function_data, - std::vector::type,dealii::Tensor<2,spacedim> >::type> &values) + std::vector >::type> &values) { const unsigned int dofs_per_cell = dof_values.size(); const unsigned int n_quadrature_points = dofs_per_cell > 0 ? @@ -1165,7 +1167,7 @@ namespace FEValuesViews AssertDimension (values.size(), n_quadrature_points); std::fill (values.begin(), values.end(), - typename ProductType::type,dealii::Tensor<2,spacedim> >::type()); + typename ProductType >::type()); for (unsigned int shape_function=0; shape_function::type>::value(0.0)); + std::fill_n (values.begin(), n_quadrature_points, + dealii::internal::NumberType::type>::value(0.0)); // add up contributions of trial functions. note that here we deal with // scalar finite elements, so no need to check for non-primitivity of @@ -2936,7 +2939,8 @@ namespace internal AssertDimension(laplacians.size(), n_quadrature_points); // initialize with zero - std::fill_n (laplacians.begin(), n_quadrature_points, dealii::internal::NumberType::type>::value(0.0)); + std::fill_n (laplacians.begin(), n_quadrature_points, + dealii::internal::NumberType::type>::value(0.0)); // add up contributions of trial functions. note that here we deal with // scalar finite elements and also note that the Laplacian is -- 2.39.5