From e81865815564cecac10c5d57cb82a432c3cc9ae7 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 4 Feb 2015 08:45:58 -0600 Subject: [PATCH] Introduce a new type ProductType to evaluate the result of a product. --- doc/news/changes.h | 6 + include/deal.II/base/template_constraints.h | 141 +++++++++++++++++++- tests/base/product_type_01.cc | 76 +++++++++++ tests/base/product_type_01.output | 2 + 4 files changed, 224 insertions(+), 1 deletion(-) create mode 100644 tests/base/product_type_01.cc create mode 100644 tests/base/product_type_01.output diff --git a/doc/news/changes.h b/doc/news/changes.h index c0b4a4cb06..b193e76b9f 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -330,6 +330,12 @@ inconvenience this causes.

Specific improvements

    +
  1. New: There is now a new class ProductType that can be used + to infer the type of the product of two objects. +
    + (Wolfgang Bangerth, 2015/02/04) +
    2. +
  2. New: The Tensor classes now have copy constructors and copy operators that allow assignment from other tensors with different underlying scalar types. diff --git a/include/deal.II/base/template_constraints.h b/include/deal.II/base/template_constraints.h index 0411ea6164..acd529d168 100644 --- a/include/deal.II/base/template_constraints.h +++ b/include/deal.II/base/template_constraints.h @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2003 - 2014 by the deal.II authors +// Copyright (C) 2003 - 2015 by the deal.II authors // // This file is part of the deal.II library. // @@ -19,6 +19,9 @@ #include +#include + + DEAL_II_NAMESPACE_OPEN template struct constraint_and_return_value; @@ -297,6 +300,142 @@ struct types_are_equal +/** + * 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 like to infer the type of the product variable + * in code like this: + * @code + * T t; + * U u; + * auto product = t*u; + * @endcode + * The local typedef of this structure represents the type the variable + * product would have. + * + * + *

    Where is this useful

    + * + * The purpose of this class is principally to represent the type one needs + * to use to represent the values or gradients of finite element fields at + * quadrature points. For example, assume you are storing the values $U_j$ of + * unknowns in a Vector, then evaluating + * $u_h(x_q) = \sum_j U_j \varphi_j(x_q)$ at quadrature points results + * in values $u_h(x_q)$ that need to be stored as @p double variables + * because the $U_j$ are @p float values and the $\varphi_j(x_q)$ are + * computed as @p double values, and the product are then @p double + * values. On the other hand, if you store your unknowns $U_j$ as + * std::complex@ values and you try to evaluate + * $\nabla u_h(x_q) = \sum_j U_j \nabla\varphi_j(x_q)$ at quadrature points, + * then the gradients $\nabla u_h(x_q)$ need to be stored as objects of + * type Tensor@<1,dim,std::complex@@> because + * that's what you get when you multiply a complex number by a + * Tensor@<1,dim@> (the type used to represent the gradient + * of shape functions of scalar finite elements). + * + * Likewise, if you are using a vector valued element (with dim components) + * and the $U_j$ are stored as @p double variables, then + * $u_h(x_q) = \sum_j U_j \varphi_j(x_q)$ needs to have type + * Tensor@<1,dim@> (because the shape functions have type + * Tensor@<1,dim@>). Finally, if you store the $U_j$ as + * objects of type std::complex@ and you have a + * vector valued element, then the gradients + * $\nabla u_h(x_q) = \sum_j U_j \nabla\varphi_j(x_q)$ will result in + * objects of type + * Tensor@<2,dim,std::complex@ @>. + * + * In all of these cases, this type is used to identify which type needs + * to be used for the result of computing the product of unknowns + * and the values, gradients, or other properties of shape functions. + * + * @author Wolfgang Bangerth, 2015 + */ +template +struct ProductType +{ +#ifdef DEAL_II_WITH_CXX11 + typedef decltype(T() * U()) type; +#endif +}; + +#ifndef DEAL_II_WITH_CXX11 + +template +struct ProductType +{ + typedef T type; +}; + +template +struct ProductType +{ + typedef T type; +}; + +template +struct ProductType +{ + typedef T type; +}; + +template +struct ProductType +{ + typedef T type; +}; + +template <> +struct ProductType +{ + typedef double type; +}; + +template <> +struct ProductType +{ + typedef double type; +}; + +#endif + + +// Annoyingly, there is no std::complex::operator*(U) for scalars U +// other than T. Consequently, even with C++11, we need the following +// specializations: +template +struct ProductType,std::complex > +{ + typedef std::complex::type> type; +}; + +template +struct ProductType > +{ + typedef std::complex::type> type; +}; + +template +struct ProductType,double> +{ + typedef std::complex::type> type; +}; + + +template +struct ProductType > +{ + typedef std::complex::type> type; +}; + +template +struct ProductType,float> +{ + typedef std::complex::type> type; +}; + + + + // --------------- inline functions ----------------- diff --git a/tests/base/product_type_01.cc b/tests/base/product_type_01.cc new file mode 100644 index 0000000000..67bf7dd2ad --- /dev/null +++ b/tests/base/product_type_01.cc @@ -0,0 +1,76 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + + +// test ProductType + +#include "../tests.h" +#include +#include +#include +#include +#include +#include + +#include +#include + + +template +void check() +{ + Assert (typeid(typename ProductType::type) == typeid(CompareType), + ExcInternalError()); + Assert (typeid(typename ProductType::type) == typeid(T() * U()), + ExcInternalError()); +} + + +int main() +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + // check product of scalars + check(); + check(); + check(); + + check(); + check(); + check(); + + // check product with Tensor<1,dim> + check,double,Tensor<1,2,double> >(); + check,double,Tensor<1,2,double> >(); + check,Tensor<1,2,double> >(); + + // check product with Tensor<2,dim> (which is a different class) + check,double,Tensor<2,2,double> >(); + check,double,Tensor<2,2,double> >(); + check,Tensor<2,2,double> >(); + + // check product with std::complex. rather annoyingly, there is no + // product between std::complex and float, or the other way + // around, so stay within the same type system + check,double,std::complex >(); + check,float,std::complex >(); + check,std::complex,Tensor<1,2,std::complex > >(); + check,Tensor<1,2>,Tensor<1,2,std::complex > >(); + + deallog << "OK" << std::endl; +} diff --git a/tests/base/product_type_01.output b/tests/base/product_type_01.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/base/product_type_01.output @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5