From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Wed, 26 Jul 2017 17:26:21 +0000 (+0200)
Subject: Introduce Quadrature::is_tensor_product_flag
X-Git-Tag: v9.0.0-rc1~1367^2~5
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1c63f0ff263f944bbd377243d99fb6a5652cc385;p=dealii.git

Introduce Quadrature::is_tensor_product_flag
---

diff --git a/doc/news/changes/minor/20170726DanielArndt-1 b/doc/news/changes/minor/20170726DanielArndt-1
new file mode 100644
index 0000000000..d5332c597a
--- /dev/null
+++ b/doc/news/changes/minor/20170726DanielArndt-1
@@ -0,0 +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
+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 e0f2e67c8f..ad933339a0 100644
--- a/include/deal.II/base/quadrature.h
+++ b/include/deal.II/base/quadrature.h
@@ -221,6 +221,26 @@ public:
   template <class Archive>
   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.
+   */
+  bool is_tensor_product() const;
+
+  /**
+   * In case the quadrature formula is a tensor product, this function
+   * returns the one-dimensional basis object.
+   * 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;
+
 protected:
   /**
    * List of quadrature points. To be filled by the constructors of derived
@@ -233,6 +253,12 @@ protected:
    * constructors of derived classes.
    */
   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
+   */
+  bool is_tensor_product_flag;
 };
 
 
@@ -381,6 +407,16 @@ Quadrature<dim>::get_weights () const
 
 
 
+template <int dim>
+inline
+bool
+Quadrature<dim>::is_tensor_product () const
+{
+  return is_tensor_product_flag;
+}
+
+
+
 template <int dim>
 template <class Archive>
 inline
diff --git a/source/base/quadrature.cc b/source/base/quadrature.cc
index 159899cf65..ad2252a4b9 100644
--- a/source/base/quadrature.cc
+++ b/source/base/quadrature.cc
@@ -30,7 +30,8 @@ template <>
 Quadrature<0>::Quadrature (const unsigned int n_q)
   :
   quadrature_points (n_q),
-  weights (n_q, 0)
+  weights (n_q, 0),
+  is_tensor_product_flag (false)
 {}
 
 
@@ -45,7 +46,8 @@ template <int dim>
 Quadrature<dim>::Quadrature (const unsigned int n_q)
   :
   quadrature_points (n_q, Point<dim>()),
-  weights (n_q, 0)
+  weights (n_q, 0),
+  is_tensor_product_flag (dim==1)
 {}
 
 
@@ -67,7 +69,8 @@ Quadrature<dim>::Quadrature (const std::vector<Point<dim> > &points,
                              const std::vector<double>      &weights)
   :
   quadrature_points(points),
-  weights(weights)
+  weights(weights),
+  is_tensor_product_flag (dim==1)
 {
   Assert (weights.size() == points.size(),
           ExcDimensionMismatch(weights.size(), points.size()));
@@ -79,7 +82,8 @@ template <int dim>
 Quadrature<dim>::Quadrature (const std::vector<Point<dim> > &points)
   :
   quadrature_points(points),
-  weights(points.size(), std::atof("Inf"))
+  weights(points.size(), std::atof("Inf")),
+  is_tensor_product_flag (dim==1)
 {
   Assert(weights.size() == points.size(),
          ExcDimensionMismatch(weights.size(), points.size()));
@@ -91,7 +95,8 @@ template <int dim>
 Quadrature<dim>::Quadrature (const Point<dim> &point)
   :
   quadrature_points(std::vector<Point<dim> > (1, point)),
-  weights(std::vector<double> (1, 1.))
+  weights(std::vector<double> (1, 1.)),
+  is_tensor_product_flag (dim==1)
 {}
 
 
@@ -109,7 +114,8 @@ Quadrature<dim>::Quadrature (const SubQuadrature &q1,
                              const Quadrature<1> &q2)
   :
   quadrature_points (q1.size() * q2.size()),
-  weights (q1.size() * q2.size())
+  weights (q1.size() * q2.size()),
+  is_tensor_product_flag (true)
 {
   unsigned int present_index = 0;
   for (unsigned int i2=0; i2<q2.size(); ++i2)
@@ -150,7 +156,8 @@ Quadrature<1>::Quadrature (const SubQuadrature &,
                            const Quadrature<1> &q2)
   :
   quadrature_points (q2.size()),
-  weights (q2.size())
+  weights (q2.size()),
+  is_tensor_product_flag (true)
 {
   unsigned int present_index = 0;
   for (unsigned int i2=0; i2<q2.size(); ++i2)
@@ -187,7 +194,8 @@ Quadrature<0>::Quadrature (const Quadrature<1> &)
   :
   Subscriptor(),
 //              quadrature_points(1),
-  weights(1,1.)
+  weights(1,1.),
+  is_tensor_product_flag (false)
 {}
 
 
@@ -209,7 +217,8 @@ Quadrature<dim>::Quadrature (const Quadrature<dim != 1 ? 1 : 0> &q)
 :
 Subscriptor(),
             quadrature_points (Utilities::fixed_power<dim>(q.size())),
-            weights (Utilities::fixed_power<dim>(q.size()))
+            weights (Utilities::fixed_power<dim>(q.size())),
+            is_tensor_product_flag(true)
 {
   Assert (dim <= 3, ExcNotImplemented());
 
@@ -243,7 +252,8 @@ Quadrature<dim>::Quadrature (const Quadrature<dim> &q)
   :
   Subscriptor(),
   quadrature_points (q.quadrature_points),
-  weights (q.weights)
+  weights (q.weights),
+  is_tensor_product_flag (q.is_tensor_product_flag)
 {}
 
 
@@ -253,6 +263,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;
   return *this;
 }
 
@@ -284,6 +295,45 @@ Quadrature<dim>::memory_consumption () const
 }
 
 
+
+template <int dim>
+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!"));
+
+  // 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);
+}
+
+
+
 //---------------------------------------------------------------------------
 template <int dim>
 QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx)
diff --git a/tests/base/quadrature_is_tensor_product.cc b/tests/base/quadrature_is_tensor_product.cc
new file mode 100644
index 0000000000..a9fbe1c06f
--- /dev/null
+++ b/tests/base/quadrature_is_tensor_product.cc
@@ -0,0 +1,113 @@
+// ---------------------------------------------------------------------
+//
+// 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 type_trait 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
+print_is_tensor_product()
+{
+  Quadrature<dim> q_0;
+  deallog << "Quadrature<" << dim << ">: "
+          << q_0.is_tensor_product() << std::endl;
+
+  QIterated<dim> q_1 (QGauss<1>(1), 1);
+  deallog << "QIterated<" << dim << ">: "
+          << q_1.is_tensor_product() << std::endl;
+
+  QGauss<dim> q_2 (1);
+  deallog << "QGauss<" << dim << ">: "
+          << q_2.is_tensor_product() << std::endl;
+
+  QGaussLobatto<dim> q_3 (2);
+  deallog << "QGaussLobatto<" << dim << ">: "
+          << q_3.is_tensor_product() << std::endl;
+
+  QMidpoint<dim> q_4;
+  deallog << "QMidpoint<" << dim << ">: "
+          << q_4.is_tensor_product() << std::endl;
+
+  QSimpson<dim> q_5;
+  deallog << "QSimpson<" << dim << ">: "
+          << q_5.is_tensor_product() << std::endl;
+
+  QTrapez<dim> q_6;
+  deallog << "QTrapez<" << dim << ">: "
+          << q_6.is_tensor_product() << std::endl;
+
+  QMilne<dim> q_7;
+  deallog << "QMilne<" << dim << ">: "
+          << q_7.is_tensor_product() << std::endl;
+
+  QWeddle<dim> q_8;
+  deallog << "QWeddle<"                << dim << ">: "
+          << q_8.is_tensor_product() << std::endl;
+
+  QGaussChebyshev<dim> q_9(3);
+  deallog << "QGaussChebyshev<"        << dim << ">: "
+          << q_9.is_tensor_product() << std::endl;
+
+  QGaussRadauChebyshev<dim> q_10(1);
+  deallog << "QGaussRadauChebyshev<" << dim << ">: "
+          << q_10.is_tensor_product() << std::endl;
+
+  QGaussLobattoChebyshev<dim> q_11(2);
+  deallog << "QGaussLobattoChebyshev<" << dim << ">: "
+          << q_11.is_tensor_product() << std::endl;
+
+  QSorted<dim> q_16(q_0);
+  deallog << "QSorted<" << dim << ">: "
+          << q_16.is_tensor_product() << std::endl;
+
+  QTelles<dim> q_17(1, Point<dim>());
+  deallog << "QTelles<" << dim << ">: "
+          << q_17.is_tensor_product() << std::endl;
+}
+
+int main()
+{
+  initlog();
+  deallog << std::boolalpha;
+
+  Quadrature<1> q;
+
+  print_is_tensor_product<1>();
+  QAnisotropic<1> q_1(q);
+  deallog << "QAnisotropic<1>: " << q_1.is_tensor_product() << std::endl;
+  QGaussLog<1> q_13(1);
+  deallog << "QGaussLog<1>: " << q_13.is_tensor_product() << std::endl;
+  QGaussLogR<1> q_14(1);
+  deallog << "QGaussLogR<1>: " << q_14.is_tensor_product() << std::endl;
+
+  print_is_tensor_product<2>();
+  QAnisotropic<2> q_2(q, q);
+  deallog << "QAnisotropic<2>: " << q_2.is_tensor_product() << std::endl;
+  QGaussOneOverR<2> q_15(1, Point<2>());
+  deallog << "QGaussOneOverR<2>: " << q_15.is_tensor_product() << std::endl;
+
+  print_is_tensor_product<3>();
+  QAnisotropic<3> q_3(q, q, q);
+  deallog << "QAnisotropic<3>: " << q_3.is_tensor_product() << std::endl;
+}
+
diff --git a/tests/base/quadrature_is_tensor_product.output b/tests/base/quadrature_is_tensor_product.output
new file mode 100644
index 0000000000..c9c4e21637
--- /dev/null
+++ b/tests/base/quadrature_is_tensor_product.output
@@ -0,0 +1,49 @@
+
+DEAL::Quadrature<1>: true
+DEAL::QIterated<1>: true
+DEAL::QGauss<1>: true
+DEAL::QGaussLobatto<1>: true
+DEAL::QMidpoint<1>: true
+DEAL::QSimpson<1>: true
+DEAL::QTrapez<1>: true
+DEAL::QMilne<1>: true
+DEAL::QWeddle<1>: true
+DEAL::QGaussChebyshev<1>: true
+DEAL::QGaussRadauChebyshev<1>: true
+DEAL::QGaussLobattoChebyshev<1>: true
+DEAL::QSorted<1>: true
+DEAL::QTelles<1>: true
+DEAL::QAnisotropic<1>: true
+DEAL::QGaussLog<1>: true
+DEAL::QGaussLogR<1>: true
+DEAL::Quadrature<2>: false
+DEAL::QIterated<2>: true
+DEAL::QGauss<2>: true
+DEAL::QGaussLobatto<2>: true
+DEAL::QMidpoint<2>: true
+DEAL::QSimpson<2>: true
+DEAL::QTrapez<2>: true
+DEAL::QMilne<2>: true
+DEAL::QWeddle<2>: true
+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::QGaussOneOverR<2>: false
+DEAL::Quadrature<3>: false
+DEAL::QIterated<3>: true
+DEAL::QGauss<3>: true
+DEAL::QGaussLobatto<3>: true
+DEAL::QMidpoint<3>: true
+DEAL::QSimpson<3>: true
+DEAL::QTrapez<3>: true
+DEAL::QMilne<3>: true
+DEAL::QWeddle<3>: true
+DEAL::QGaussChebyshev<3>: true
+DEAL::QGaussRadauChebyshev<3>: true
+DEAL::QGaussLobattoChebyshev<3>: true
+DEAL::QSorted<3>: false
+DEAL::QTelles<3>: false
+DEAL::QAnisotropic<3>: false