From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Thu, 27 Jul 2017 21:44:06 +0000 (+0200)
Subject: Allow tensor products with different formulas in each direction
X-Git-Tag: v9.0.0-rc1~1367^2~4
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9692686965595340e5f5383bb6791edf40f6f333;p=dealii.git

Allow tensor products with different formulas in each direction
---

diff --git a/doc/news/changes/minor/20170726DanielArndt-1 b/doc/news/changes/minor/20170726DanielArndt-1
index d5332c597a..843b8ec9a5 100644
--- a/doc/news/changes/minor/20170726DanielArndt-1
+++ b/doc/news/changes/minor/20170726DanielArndt-1
@@ -1,6 +1,6 @@
 New: The Quadrature class has a new boolean is_tensor_product_flag
 to indicate if the represented quadrature formula is in fact a
-tensor product of (identical) 1D formulas. The corresponding 1D object
+tensor product of 1D formulas. The corresponding 1D objects
 can be queried using Quadrature::get_tensor_basis().
 <br>
 (Daniel Arndt, 2017/07/26)
diff --git a/include/deal.II/base/quadrature.h b/include/deal.II/base/quadrature.h
index ad933339a0..2fb29302b2 100644
--- a/include/deal.II/base/quadrature.h
+++ b/include/deal.II/base/quadrature.h
@@ -222,24 +222,17 @@ public:
   void serialize (Archive &ar, const unsigned int version);
 
   /**
-   * This quadrature object is a tensor product in case the tensor product
-   * is formed by the first $size^(1/dim)$ quadrature points.
-   * Return if this is the case.
+   * This function returns true if the quadrature object is a tensor product
+   * of one-dimensional formulas.
    */
   bool is_tensor_product() const;
 
   /**
    * In case the quadrature formula is a tensor product, this function
-   * returns the one-dimensional basis object.
+   * returns the one-dimensional basis objects.
    * Otherwise, calling this function is not allowed.
-   *
-   * Each time this function is called the constructor for the
-   * one-dimensional Quadrature is called and a copy returned.
-   *
-   * The weights for the one-dimensional formula are rescaled
-   * but they are not modified in case the sum is zero or infinite.
    */
-  Quadrature<1> get_tensor_basis() const;
+  const std::vector<Quadrature<1>> &get_tensor_basis() const;
 
 protected:
   /**
@@ -255,10 +248,18 @@ protected:
   std::vector<double>      weights;
 
   /**
-   * This quadrature object is a tensor product in case the tensor product
-   * is formed by the first $size^(1/dim)$ quadrature points
+   * Indicates if this object represents quadrature formula that is a tensor
+   * product of one-dimensional formulas.
+   * This flag is set if dim==1 or the constructors taking a Quadrature<1>
+   * (and possibly a Quadrature<dim-1> object) is called.
    */
   bool is_tensor_product_flag;
+
+  /**
+   * Stores the one-dimensional tensor basis objects in case this object
+   * can be represented by a tensor product.
+   */
+  std::vector<Quadrature<1> > tensor_basis;
 };
 
 
diff --git a/source/base/quadrature.cc b/source/base/quadrature.cc
index ad2252a4b9..fc7117ab3d 100644
--- a/source/base/quadrature.cc
+++ b/source/base/quadrature.cc
@@ -26,6 +26,25 @@
 DEAL_II_NAMESPACE_OPEN
 
 
+namespace
+{
+  template <int dim>
+  void
+  set_tensor_basis_1d(const Quadrature<dim> &,
+                      std::vector<Quadrature<1>> &)
+  {}
+
+  template <>
+  void
+  set_tensor_basis_1d(const Quadrature<1> &quadrature,
+                      std::vector<Quadrature<1>> &tensor_basis)
+  {
+    tensor_basis = std::vector<Quadrature<1> >(1, quadrature);
+  }
+}
+
+
+
 template <>
 Quadrature<0>::Quadrature (const unsigned int n_q)
   :
@@ -48,7 +67,9 @@ Quadrature<dim>::Quadrature (const unsigned int n_q)
   quadrature_points (n_q, Point<dim>()),
   weights (n_q, 0),
   is_tensor_product_flag (dim==1)
-{}
+{
+  set_tensor_basis_1d(*this, this->tensor_basis);
+}
 
 
 
@@ -74,6 +95,7 @@ Quadrature<dim>::Quadrature (const std::vector<Point<dim> > &points,
 {
   Assert (weights.size() == points.size(),
           ExcDimensionMismatch(weights.size(), points.size()));
+  set_tensor_basis_1d(*this, this->tensor_basis);
 }
 
 
@@ -87,6 +109,7 @@ Quadrature<dim>::Quadrature (const std::vector<Point<dim> > &points)
 {
   Assert(weights.size() == points.size(),
          ExcDimensionMismatch(weights.size(), points.size()));
+  set_tensor_basis_1d(*this, this->tensor_basis);
 }
 
 
@@ -97,7 +120,9 @@ Quadrature<dim>::Quadrature (const Point<dim> &point)
   quadrature_points(std::vector<Point<dim> > (1, point)),
   weights(std::vector<double> (1, 1.)),
   is_tensor_product_flag (dim==1)
-{}
+{
+  set_tensor_basis_1d(*this, this->tensor_basis);
+}
 
 
 template <>
@@ -115,7 +140,7 @@ Quadrature<dim>::Quadrature (const SubQuadrature &q1,
   :
   quadrature_points (q1.size() * q2.size()),
   weights (q1.size() * q2.size()),
-  is_tensor_product_flag (true)
+  is_tensor_product_flag (q1.is_tensor_product())
 {
   unsigned int present_index = 0;
   for (unsigned int i2=0; i2<q2.size(); ++i2)
@@ -147,6 +172,12 @@ Quadrature<dim>::Quadrature (const SubQuadrature &q1,
       Assert ((sum>0.999999) && (sum<1.000001), ExcInternalError());
     }
 #endif
+
+  if (is_tensor_product_flag)
+    {
+      tensor_basis = q1.get_tensor_basis();
+      tensor_basis.push_back(q2);
+    }
 }
 
 
@@ -157,7 +188,8 @@ Quadrature<1>::Quadrature (const SubQuadrature &,
   :
   quadrature_points (q2.size()),
   weights (q2.size()),
-  is_tensor_product_flag (true)
+  is_tensor_product_flag (true),
+  tensor_basis(std::vector<Quadrature<1> >(1,q2))
 {
   unsigned int present_index = 0;
   for (unsigned int i2=0; i2<q2.size(); ++i2)
@@ -243,6 +275,8 @@ Subscriptor(),
             weights[k] *= q.weight(i2);
           ++k;
         }
+  for (unsigned int i=0; i<dim; ++i)
+    tensor_basis = std::vector<Quadrature<1> >(dim, q);
 }
 
 
@@ -253,7 +287,8 @@ Quadrature<dim>::Quadrature (const Quadrature<dim> &q)
   Subscriptor(),
   quadrature_points (q.quadrature_points),
   weights (q.weights),
-  is_tensor_product_flag (q.is_tensor_product_flag)
+  is_tensor_product_flag (q.is_tensor_product_flag),
+  tensor_basis (q.tensor_basis)
 {}
 
 
@@ -264,6 +299,7 @@ Quadrature<dim>::operator= (const Quadrature<dim> &q)
   weights = q.weights;
   quadrature_points = q.quadrature_points;
   is_tensor_product_flag = q.is_tensor_product_flag;
+  tensor_basis = q.tensor_basis;
   return *this;
 }
 
@@ -297,39 +333,15 @@ Quadrature<dim>::memory_consumption () const
 
 
 template <int dim>
-Quadrature<1>
+const std::vector<Quadrature<1> > &
 Quadrature<dim>::get_tensor_basis () const
 {
   Assert (this->is_tensor_product_flag == true,
           ExcMessage("This function only makes sense if "
                      "this object represents a tensor product!"));
+  Assert (tensor_basis.size()==dim, ExcInternalError());
 
-  // Just take the first components of the first quadrature points
-  // and rescale the weights accordingly.
-  const unsigned int n_q_points = this->size();
-  const unsigned int n_q_points_1d
-    = static_cast<unsigned int>(std::round(std::pow(n_q_points, 1./dim)));
-  Assert(Utilities::fixed_power<dim>(n_q_points_1d) == n_q_points,
-         ExcInternalError());
-
-  std::vector<Point<1> > q_points_1d (n_q_points_1d);
-  std::vector<double> weights_1d (n_q_points_1d);
-
-  double sum = 0.;
-  for (unsigned int i=0; i<n_q_points_1d; ++i)
-    {
-      sum += this->weight(i);
-      q_points_1d[i](0) = this->point(i)(0);
-    }
-
-  if (!std::isinf(sum) || sum>1.e-10)
-    for (unsigned int i=0; i<n_q_points_1d; ++i)
-      weights_1d[i] = this->weight(i)/sum;
-  else
-    for (unsigned int i=0; i<n_q_points_1d; ++i)
-      weights_1d[i] = this->weight(i);
-
-  return Quadrature<1> (q_points_1d, weights_1d);
+  return tensor_basis;
 }
 
 
@@ -347,6 +359,9 @@ QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx)
       this->weights[k++] = qx.weight(k1);
     }
   Assert (k==this->size(), ExcInternalError());
+  this->is_tensor_product_flag = true;
+  this->tensor_basis.resize(1);
+  this->tensor_basis[0] = qx;
 }
 
 
@@ -354,8 +369,7 @@ QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx)
 template <int dim>
 QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx,
                                 const Quadrature<1> &qy)
-  : Quadrature<dim>(qx.size()
-                    *qy.size())
+  :  Quadrature<dim>(qx.size()*qy.size())
 {
   Assert (dim==2, ExcImpossibleInDim(dim));
   unsigned int k=0;
@@ -367,6 +381,10 @@ QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx,
         this->weights[k++] = qx.weight(k1) * qy.weight(k2);
       }
   Assert (k==this->size(), ExcInternalError());
+  this->is_tensor_product_flag = true;
+  this->tensor_basis.resize(2);
+  this->tensor_basis[0] = qx;
+  this->tensor_basis[1] = qy;
 }
 
 
@@ -375,9 +393,7 @@ template <int dim>
 QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx,
                                 const Quadrature<1> &qy,
                                 const Quadrature<1> &qz)
-  : Quadrature<dim>(qx.size()
-                    *qy.size()
-                    *qz.size())
+  : Quadrature<dim>(qx.size()*qy.size()*qz.size())
 {
   Assert (dim==3, ExcImpossibleInDim(dim));
   unsigned int k=0;
@@ -391,6 +407,11 @@ QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx,
           this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
         }
   Assert (k==this->size(), ExcInternalError());
+  this->is_tensor_product_flag = true;
+  this->tensor_basis.resize(3);
+  this->tensor_basis[0] = qx;
+  this->tensor_basis[1] = qy;
+  this->tensor_basis[2] = qz;
 }
 
 
@@ -1791,6 +1812,7 @@ QIterated<1>::QIterated (const Quadrature<1> &base_quadrature,
   Assert (std::fabs(sum_of_weights-1) < 1e-13,
           ExcInternalError());
 #endif
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc
index 6e7d9d1960..8e19826722 100644
--- a/source/base/quadrature_lib.cc
+++ b/source/base/quadrature_lib.cc
@@ -159,6 +159,7 @@ QGauss<1>::QGauss (const unsigned int n)
       this->weights[i-1] = w;
       this->weights[n-i] = w;
     }
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
@@ -179,6 +180,7 @@ QGaussLobatto<1>::QGaussLobatto (const unsigned int n)
       this->quadrature_points[i] = Point<1>(0.5 + 0.5*static_cast<double>(points[i]));
       this->weights[i]           = 0.5*w[i];
     }
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
@@ -336,6 +338,7 @@ QMidpoint<1>::QMidpoint ()
 {
   this->quadrature_points[0] = Point<1>(0.5);
   this->weights[0] = 1.0;
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
@@ -353,6 +356,7 @@ QTrapez<1>::QTrapez ()
       this->quadrature_points[i] = Point<1>(xpts[i]);
       this->weights[i] = wts[i];
     };
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
@@ -370,6 +374,7 @@ QSimpson<1>::QSimpson ()
       this->quadrature_points[i] = Point<1>(xpts[i]);
       this->weights[i] = wts[i];
     };
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
@@ -387,6 +392,7 @@ QMilne<1>::QMilne ()
       this->quadrature_points[i] = Point<1>(xpts[i]);
       this->weights[i] = wts[i];
     };
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
@@ -406,6 +412,7 @@ QWeddle<1>::QWeddle ()
       this->quadrature_points[i] = Point<1>(xpts[i]);
       this->weights[i] = wts[i];
     };
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
@@ -426,6 +433,7 @@ QGaussLog<1>::QGaussLog(const unsigned int n,
       this->quadrature_points[i] = revert ? Point<1>(1-p[n-1-i]) : Point<1>(p[i]);
       this->weights[i]           = revert ? w[n-1-i] : w[i];
     }
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 template <>
@@ -762,6 +770,7 @@ QGaussLogR<1>::QGaussLogR(const unsigned int n,
                            "factor out the singularity, which is zero at one point."));
         this->weights[i] /= denominator;
       }
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
@@ -1120,6 +1129,7 @@ QTelles<1>::QTelles (
       weights[q] = J * weights[q];
 
     }
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
@@ -1174,7 +1184,7 @@ QGaussChebyshev<1>::QGaussChebyshev(const unsigned int n)
       this->quadrature_points[i] = Point<1>(p[i]);
       this->weights[i]           = w[i];
     }
-
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
@@ -1264,6 +1274,7 @@ QGaussRadauChebyshev<1>::QGaussRadauChebyshev(const unsigned int n,
       this->quadrature_points[i] = Point<1>(p[i]);
       this->weights[i]           = w[i];
     }
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
@@ -1328,7 +1339,7 @@ QGaussLobattoChebyshev<1>::QGaussLobattoChebyshev(const unsigned int n)
       this->quadrature_points[i] = Point<1>(p[i]);
       this->weights[i]           = w[i];
     }
-
+  this->tensor_basis = std::vector<Quadrature<1> >(1, *this);
 }
 
 
diff --git a/tests/base/quadrature_check_tensor_product.cc b/tests/base/quadrature_check_tensor_product.cc
new file mode 100644
index 0000000000..44c39010bc
--- /dev/null
+++ b/tests/base/quadrature_check_tensor_product.cc
@@ -0,0 +1,218 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the boolean is_tensor_product for all the quadrature classes
+
+
+#include "../tests.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/quadrature.h>
+
+template <int dim>
+void
+check_tensor_product(const std::vector<Quadrature<dim> > &quadratures,
+                     const std::vector<std::string> &quadrature_names)
+{
+  Assert(false, ExcNotImplemented());
+}
+
+template <>
+void
+check_tensor_product(const std::vector<Quadrature<1> > &quadratures,
+                     const std::vector<std::string> &quadrature_names)
+{
+  for (unsigned int i=0; i<quadratures.size(); ++i)
+    {
+      const Quadrature<1> &quadrature = quadratures[i];
+      if (quadrature.is_tensor_product())
+        {
+          deallog << "1D " << quadrature_names[i];
+          const auto &q_basis = quadrature.get_tensor_basis();
+          AssertThrow(q_basis.size()==1, ExcInternalError());
+          const auto &q_points = quadrature.get_points();
+          const auto &q_weights = quadrature.get_weights();
+          AssertThrow(q_basis[0].size() == q_points.size(),
+                      ExcInternalError());
+          for (unsigned int q=0; q<quadrature.size(); ++q)
+            {
+              std::cout << q_points[q] << " " << q_basis[0].get_points()[q] << std::endl;
+              AssertThrow(std::abs((q_points[q]-q_basis[0].get_points()[q]).norm())<1.e-10,
+                          ExcInternalError());
+              AssertThrow(std::abs(q_weights[q]-q_basis[0].get_weights()[q])<1.e-10,
+                          ExcInternalError());
+            }
+          deallog << " OK" << std::endl;
+        }
+    }
+}
+
+template <>
+void
+check_tensor_product(const std::vector<Quadrature<2> > &quadratures,
+                     const std::vector<std::string> &quadrature_names)
+{
+  for (unsigned int i=0; i<quadratures.size(); ++i)
+    {
+      const Quadrature<2> &quadrature = quadratures[i];
+      if (quadrature.is_tensor_product())
+        {
+          deallog << "2D " << quadrature_names[i];
+          const auto &q_basis = quadrature.get_tensor_basis();
+          AssertThrow(q_basis.size()==2, ExcInternalError());
+          const auto &q_points = quadrature.get_points();
+          const auto &q_weights = quadrature.get_weights();
+          AssertThrow(q_basis[0].size()*q_basis[1].size() == q_points.size(),
+                      ExcInternalError());
+          unsigned int q=0;
+          for (unsigned int q2=0; q2<q_basis[0].size(); ++q2)
+            for (unsigned int q1=0; q1<q_basis[1].size(); ++q1)
+              {
+                AssertThrow(std::abs(q_points[q][0]-q_basis[0].get_points()[q1][0])<1.e-10,
+                            ExcInternalError());
+                AssertThrow(std::abs(q_points[q][1]-q_basis[1].get_points()[q2][0])<1.e-10,
+                            ExcInternalError());
+                AssertThrow(std::abs((q_weights[q]-q_basis[0].get_weights()[q1]*q_basis[1].get_weights()[q2]))<1.e-10,
+                            ExcInternalError());
+                ++q;
+              }
+          deallog << " OK" << std::endl;
+        }
+    }
+}
+
+template <>
+void
+check_tensor_product(const std::vector<Quadrature<3> > &quadratures,
+                     const std::vector<std::string> &quadrature_names)
+{
+  for (unsigned int i=0; i<quadratures.size(); ++i)
+    {
+      const Quadrature<3> &quadrature = quadratures[i];
+      if (quadrature.is_tensor_product())
+        {
+          deallog << "3D " << quadrature_names[i];
+          const auto &q_basis = quadrature.get_tensor_basis();
+          AssertThrow(q_basis.size()==3, ExcInternalError());
+          const auto &q_points = quadrature.get_points();
+          const auto &q_weights = quadrature.get_weights();
+          AssertThrow(q_basis[0].size()*q_basis[1].size()*q_basis[2].size()
+                      == q_points.size(), ExcInternalError());
+          unsigned int q=0;
+          for (unsigned int q3=0; q3<q_basis[2].size(); ++q3)
+            for (unsigned int q2=0; q2<q_basis[1].size(); ++q2)
+              for (unsigned int q1=0; q1<q_basis[0].size(); ++q1)
+                {
+                  AssertThrow(std::abs(q_points[q][0]-q_basis[0].get_points()[q1][0])<1.e-10,
+                              ExcInternalError());
+                  AssertThrow(std::abs(q_points[q][1]-q_basis[1].get_points()[q2][0])<1.e-10,
+                              ExcInternalError());
+                  AssertThrow(std::abs(q_points[q][2]-q_basis[2].get_points()[q3][0])<1.e-10,
+                              ExcInternalError());
+                  AssertThrow(std::abs(q_weights[q]-q_basis[0].get_weights()[q1]*q_basis[1].get_weights()[q2]*q_basis[2].get_weights()[q3])<1.e-10,
+                              ExcInternalError());
+                  ++q;
+                }
+          deallog << " OK" << std::endl;
+        }
+    }
+}
+
+template <int dim>
+void
+fill_quadrature_vector(std::vector<Quadrature<dim> > &quadratures,
+                       std::vector<std::string> &quadrature_names)
+{
+  quadratures.push_back(Quadrature<dim>());
+  quadrature_names.push_back("Quadrature");
+
+  quadratures.push_back(QIterated<dim> (QGauss<1>(2), 2));
+  quadrature_names.push_back("QIterated");
+
+  quadratures.push_back(QGauss<dim> (2));
+  quadrature_names.push_back("QGauss");
+
+  quadratures.push_back(QGaussLobatto<dim> (2));
+  quadrature_names.push_back("QGaussLobatto");
+
+  quadratures.push_back(QMidpoint<dim> ());
+  quadrature_names.push_back("QMidPoint");
+
+  quadratures.push_back(QSimpson<dim> ());
+  quadrature_names.push_back("QSimpson");
+
+  quadratures.push_back(QTrapez<dim> ());
+  quadrature_names.push_back("QTrapez");
+
+  quadratures.push_back(QMilne<dim> ());
+  quadrature_names.push_back("QMilne");
+
+  quadratures.push_back(QWeddle<dim> ());
+  quadrature_names.push_back("QWeddle");
+
+  quadratures.push_back(QGaussChebyshev<dim> (3));
+  quadrature_names.push_back("QGaussChebyshev");
+
+  quadratures.push_back(QGaussRadauChebyshev<dim> (2));
+  quadrature_names.push_back("QGaussRadauChebyshev");
+
+  quadratures.push_back(QGaussLobattoChebyshev<dim> (2));
+  quadrature_names.push_back("QGaussLobattoChebyshev");
+
+  quadratures.push_back(QSorted<dim> (Quadrature<dim>()));
+  quadrature_names.push_back("QSorted");
+
+  quadratures.push_back(QTelles<dim> (1, Point<dim>()));
+  quadrature_names.push_back("QTelles");
+}
+
+int main()
+{
+  initlog();
+  deallog << std::boolalpha;
+
+  Quadrature<1> q;
+
+  std::vector<Quadrature<1> > quadratures_1d;
+  std::vector<std::string> quadrature_names_1d;
+  fill_quadrature_vector(quadratures_1d, quadrature_names_1d);
+  quadratures_1d.push_back(QAnisotropic<1>(q));
+  quadrature_names_1d.push_back("QAnisotropic");
+  quadratures_1d.push_back(QGaussLog<1> (1));
+  quadrature_names_1d.push_back("QGaussLog");
+  quadratures_1d.push_back(QGaussLogR<1> (1));
+  quadrature_names_1d.push_back("QGaussLogR");
+  check_tensor_product(quadratures_1d, quadrature_names_1d);
+
+  std::vector<Quadrature<2> > quadratures_2d;
+  std::vector<std::string> quadrature_names_2d;
+  fill_quadrature_vector(quadratures_2d, quadrature_names_2d);
+  quadratures_2d.push_back(QAnisotropic<2>(q,q));
+  quadrature_names_2d.push_back("QAnisotropic");
+  quadratures_2d.push_back(QGaussOneOverR<2> (1, Point<2>()));
+  quadrature_names_2d.push_back("QGaussOneOverR");
+  check_tensor_product(quadratures_2d, quadrature_names_2d);
+
+  std::vector<Quadrature<3> > quadratures_3d;
+  std::vector<std::string> quadrature_names_3d;
+  fill_quadrature_vector(quadratures_3d, quadrature_names_3d);
+  quadratures_3d.push_back(QAnisotropic<3>(q,q,q));
+  quadrature_names_3d.push_back("QAnisotropic");
+  check_tensor_product(quadratures_3d, quadrature_names_3d);
+}
+
diff --git a/tests/base/quadrature_check_tensor_product.output b/tests/base/quadrature_check_tensor_product.output
new file mode 100644
index 0000000000..47df9efd25
--- /dev/null
+++ b/tests/base/quadrature_check_tensor_product.output
@@ -0,0 +1,44 @@
+
+DEAL::1D Quadrature OK
+DEAL::1D QIterated OK
+DEAL::1D QGauss OK
+DEAL::1D QGaussLobatto OK
+DEAL::1D QMidPoint OK
+DEAL::1D QSimpson OK
+DEAL::1D QTrapez OK
+DEAL::1D QMilne OK
+DEAL::1D QWeddle OK
+DEAL::1D QGaussChebyshev OK
+DEAL::1D QGaussRadauChebyshev OK
+DEAL::1D QGaussLobattoChebyshev OK
+DEAL::1D QSorted OK
+DEAL::1D QTelles OK
+DEAL::1D QAnisotropic OK
+DEAL::1D QGaussLog OK
+DEAL::1D QGaussLogR OK
+DEAL::2D QIterated OK
+DEAL::2D QGauss OK
+DEAL::2D QGaussLobatto OK
+DEAL::2D QMidPoint OK
+DEAL::2D QSimpson OK
+DEAL::2D QTrapez OK
+DEAL::2D QMilne OK
+DEAL::2D QWeddle OK
+DEAL::2D QGaussChebyshev OK
+DEAL::2D QGaussRadauChebyshev OK
+DEAL::2D QGaussLobattoChebyshev OK
+DEAL::2D QTelles OK
+DEAL::2D QAnisotropic OK
+DEAL::3D QIterated OK
+DEAL::3D QGauss OK
+DEAL::3D QGaussLobatto OK
+DEAL::3D QMidPoint OK
+DEAL::3D QSimpson OK
+DEAL::3D QTrapez OK
+DEAL::3D QMilne OK
+DEAL::3D QWeddle OK
+DEAL::3D QGaussChebyshev OK
+DEAL::3D QGaussRadauChebyshev OK
+DEAL::3D QGaussLobattoChebyshev OK
+DEAL::3D QTelles OK
+DEAL::3D QAnisotropic OK
diff --git a/tests/base/quadrature_is_tensor_product.cc b/tests/base/quadrature_is_tensor_product.cc
index a9fbe1c06f..d551d15f6f 100644
--- a/tests/base/quadrature_is_tensor_product.cc
+++ b/tests/base/quadrature_is_tensor_product.cc
@@ -15,7 +15,7 @@
 
 
 
-// check the type_trait is_tensor_product for all the quadrature classes
+// check the boolean is_tensor_product for all the quadrature classes
 
 
 #include "../tests.h"
diff --git a/tests/base/quadrature_is_tensor_product.output b/tests/base/quadrature_is_tensor_product.output
index c9c4e21637..ab53fc52fa 100644
--- a/tests/base/quadrature_is_tensor_product.output
+++ b/tests/base/quadrature_is_tensor_product.output
@@ -29,8 +29,8 @@ DEAL::QGaussChebyshev<2>: true
 DEAL::QGaussRadauChebyshev<2>: true
 DEAL::QGaussLobattoChebyshev<2>: true
 DEAL::QSorted<2>: false
-DEAL::QTelles<2>: false
-DEAL::QAnisotropic<2>: false
+DEAL::QTelles<2>: true
+DEAL::QAnisotropic<2>: true
 DEAL::QGaussOneOverR<2>: false
 DEAL::Quadrature<3>: false
 DEAL::QIterated<3>: true
@@ -45,5 +45,5 @@ DEAL::QGaussChebyshev<3>: true
 DEAL::QGaussRadauChebyshev<3>: true
 DEAL::QGaussLobattoChebyshev<3>: true
 DEAL::QSorted<3>: false
-DEAL::QTelles<3>: false
-DEAL::QAnisotropic<3>: false
+DEAL::QTelles<3>: true
+DEAL::QAnisotropic<3>: true