// ---------------------------------------------------------------------
//
-// 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.
//
#include <deal.II/base/config.h>
+#include <complex>
+
+
DEAL_II_NAMESPACE_OPEN
template <bool, typename> struct constraint_and_return_value;
+/**
+ * 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 <code>product</code> 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
+ * <code>product</code> would have.
+ *
+ *
+ * <h3>Where is this useful</h3>
+ *
+ * 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<float>, 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
+ * <code>std::complex@<double@></code> 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 <code>Tensor@<1,dim,std::complex@<double@>@></code> because
+ * that's what you get when you multiply a complex number by a
+ * <code>Tensor@<1,dim@></code> (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
+ * <code>Tensor@<1,dim@></code> (because the shape functions have type
+ * <code>Tensor@<1,dim@></code>). Finally, if you store the $U_j$ as
+ * objects of type <code>std::complex@<double@></code> 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
+ * <code>Tensor@<2,dim,std::complex@<double@> @></code>.
+ *
+ * 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 <typename T, typename U>
+struct ProductType
+{
+#ifdef DEAL_II_WITH_CXX11
+ typedef decltype(T() * U()) type;
+#endif
+};
+
+#ifndef DEAL_II_WITH_CXX11
+
+template <typename T>
+struct ProductType<T,double>
+{
+ typedef T type;
+};
+
+template <typename T>
+struct ProductType<T,bool>
+{
+ typedef T type;
+};
+
+template <typename T>
+struct ProductType<double,T>
+{
+ typedef T type;
+};
+
+template <typename T>
+struct ProductType<bool, T>
+{
+ typedef T type;
+};
+
+template <>
+struct ProductType<bool,double>
+{
+ typedef double type;
+};
+
+template <>
+struct ProductType<double,bool>
+{
+ typedef double type;
+};
+
+#endif
+
+
+// Annoyingly, there is no std::complex<T>::operator*(U) for scalars U
+// other than T. Consequently, even with C++11, we need the following
+// specializations:
+template <typename T, typename U>
+struct ProductType<std::complex<T>,std::complex<U> >
+{
+ typedef std::complex<typename ProductType<T,U>::type> type;
+};
+
+template <typename U>
+struct ProductType<double,std::complex<U> >
+{
+ typedef std::complex<typename ProductType<double,U>::type> type;
+};
+
+template <typename T>
+struct ProductType<std::complex<T>,double>
+{
+ typedef std::complex<typename ProductType<T,double>::type> type;
+};
+
+
+template <typename U>
+struct ProductType<float,std::complex<U> >
+{
+ typedef std::complex<typename ProductType<float,U>::type> type;
+};
+
+template <typename T>
+struct ProductType<std::complex<T>,float>
+{
+ typedef std::complex<typename ProductType<T,float>::type> type;
+};
+
+
+
+
// --------------- inline functions -----------------
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <iomanip>
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+#include <typeinfo>
+#include <complex>
+
+#include <deal.II/base/template_constraints.h>
+#include <deal.II/base/tensor.h>
+
+
+template <typename T, typename U, typename CompareType>
+void check()
+{
+ Assert (typeid(typename ProductType<T,U>::type) == typeid(CompareType),
+ ExcInternalError());
+ Assert (typeid(typename ProductType<T,U>::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<double,double,double>();
+ check<float,double,double>();
+ check<double,float,double>();
+
+ check<int,int,int>();
+ check<int,double,double>();
+ check<double,int,double>();
+
+ // check product with Tensor<1,dim>
+ check<Tensor<1,2,double>,double,Tensor<1,2,double> >();
+ check<Tensor<1,2,float>,double,Tensor<1,2,double> >();
+ check<double,Tensor<1,2,float>,Tensor<1,2,double> >();
+
+ // check product with Tensor<2,dim> (which is a different class)
+ check<Tensor<2,2,double>,double,Tensor<2,2,double> >();
+ check<Tensor<2,2,float>,double,Tensor<2,2,double> >();
+ check<double,Tensor<2,2,float>,Tensor<2,2,double> >();
+
+ // check product with std::complex. rather annoyingly, there is no
+ // product between std::complex<double> and float, or the other way
+ // around, so stay within the same type system
+ check<std::complex<double>,double,std::complex<double> >();
+ check<std::complex<float>,float,std::complex<float> >();
+ check<Tensor<1,2>,std::complex<double>,Tensor<1,2,std::complex<double> > >();
+ check<std::complex<double>,Tensor<1,2>,Tensor<1,2,std::complex<double> > >();
+
+ deallog << "OK" << std::endl;
+}