]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduce a new type ProductType to evaluate the result of a product.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 4 Feb 2015 14:45:58 +0000 (08:45 -0600)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 12 Feb 2015 13:31:42 +0000 (07:31 -0600)
doc/news/changes.h
include/deal.II/base/template_constraints.h
tests/base/product_type_01.cc [new file with mode: 0644]
tests/base/product_type_01.output [new file with mode: 0644]

index c0b4a4cb060591ff57830d0e734ad9a03772ba27..b193e76b9ff8be7389364cb545d79a1d94a455ed 100644 (file)
@@ -330,6 +330,12 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> New: There is now a new class ProductType that can be used
+  to infer the type of the product of two objects.
+  <br>
+  (Wolfgang Bangerth, 2015/02/04)
+  </li>
+
   <li> New: The Tensor classes now have copy constructors and copy
   operators that allow assignment from other tensors with different
   underlying scalar types.
index 0411ea61648e7a112ecef9c47531f1c969bb1521..acd529d168ad5a357d35fae63efd722a4e06e464 100644 (file)
@@ -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 <deal.II/base/config.h>
 
+#include <complex>
+
+
 DEAL_II_NAMESPACE_OPEN
 
 template <bool, typename> struct constraint_and_return_value;
@@ -297,6 +300,142 @@ struct types_are_equal<T,T>
 
 
 
+/**
+ * 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 -----------------
 
 
diff --git a/tests/base/product_type_01.cc b/tests/base/product_type_01.cc
new file mode 100644 (file)
index 0000000..67bf7dd
--- /dev/null
@@ -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 <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;
+}
diff --git a/tests/base/product_type_01.output b/tests/base/product_type_01.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.