From 1c63f0ff263f944bbd377243d99fb6a5652cc385 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Wed, 26 Jul 2017 19:26:21 +0200 Subject: [PATCH] Introduce Quadrature::is_tensor_product_flag --- doc/news/changes/minor/20170726DanielArndt-1 | 6 + include/deal.II/base/quadrature.h | 36 ++++++ source/base/quadrature.cc | 70 +++++++++-- tests/base/quadrature_is_tensor_product.cc | 113 ++++++++++++++++++ .../base/quadrature_is_tensor_product.output | 49 ++++++++ 5 files changed, 264 insertions(+), 10 deletions(-) create mode 100644 doc/news/changes/minor/20170726DanielArndt-1 create mode 100644 tests/base/quadrature_is_tensor_product.cc create mode 100644 tests/base/quadrature_is_tensor_product.output 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(). +
+(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 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 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::get_weights () const +template +inline +bool +Quadrature::is_tensor_product () const +{ + return is_tensor_product_flag; +} + + + template template 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 Quadrature::Quadrature (const unsigned int n_q) : quadrature_points (n_q, Point()), - weights (n_q, 0) + weights (n_q, 0), + is_tensor_product_flag (dim==1) {} @@ -67,7 +69,8 @@ Quadrature::Quadrature (const std::vector > &points, const std::vector &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 Quadrature::Quadrature (const std::vector > &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 Quadrature::Quadrature (const Point &point) : quadrature_points(std::vector > (1, point)), - weights(std::vector (1, 1.)) + weights(std::vector (1, 1.)), + is_tensor_product_flag (dim==1) {} @@ -109,7 +114,8 @@ Quadrature::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::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::Quadrature (const Quadrature<1> &) : Subscriptor(), // quadrature_points(1), - weights(1,1.) + weights(1,1.), + is_tensor_product_flag (false) {} @@ -209,7 +217,8 @@ Quadrature::Quadrature (const Quadrature &q) : Subscriptor(), quadrature_points (Utilities::fixed_power(q.size())), - weights (Utilities::fixed_power(q.size())) + weights (Utilities::fixed_power(q.size())), + is_tensor_product_flag(true) { Assert (dim <= 3, ExcNotImplemented()); @@ -243,7 +252,8 @@ Quadrature::Quadrature (const Quadrature &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::operator= (const Quadrature &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::memory_consumption () const } + +template +Quadrature<1> +Quadrature::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(std::round(std::pow(n_q_points, 1./dim))); + Assert(Utilities::fixed_power(n_q_points_1d) == n_q_points, + ExcInternalError()); + + std::vector > q_points_1d (n_q_points_1d); + std::vector weights_1d (n_q_points_1d); + + double sum = 0.; + for (unsigned int i=0; iweight(i); + q_points_1d[i](0) = this->point(i)(0); + } + + if (!std::isinf(sum) || sum>1.e-10) + for (unsigned int i=0; iweight(i)/sum; + else + for (unsigned int i=0; iweight(i); + + return Quadrature<1> (q_points_1d, weights_1d); +} + + + //--------------------------------------------------------------------------- template QAnisotropic::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 +#include +#include + +template +void +print_is_tensor_product() +{ + Quadrature q_0; + deallog << "Quadrature<" << dim << ">: " + << q_0.is_tensor_product() << std::endl; + + QIterated q_1 (QGauss<1>(1), 1); + deallog << "QIterated<" << dim << ">: " + << q_1.is_tensor_product() << std::endl; + + QGauss q_2 (1); + deallog << "QGauss<" << dim << ">: " + << q_2.is_tensor_product() << std::endl; + + QGaussLobatto q_3 (2); + deallog << "QGaussLobatto<" << dim << ">: " + << q_3.is_tensor_product() << std::endl; + + QMidpoint q_4; + deallog << "QMidpoint<" << dim << ">: " + << q_4.is_tensor_product() << std::endl; + + QSimpson q_5; + deallog << "QSimpson<" << dim << ">: " + << q_5.is_tensor_product() << std::endl; + + QTrapez q_6; + deallog << "QTrapez<" << dim << ">: " + << q_6.is_tensor_product() << std::endl; + + QMilne q_7; + deallog << "QMilne<" << dim << ">: " + << q_7.is_tensor_product() << std::endl; + + QWeddle q_8; + deallog << "QWeddle<" << dim << ">: " + << q_8.is_tensor_product() << std::endl; + + QGaussChebyshev q_9(3); + deallog << "QGaussChebyshev<" << dim << ">: " + << q_9.is_tensor_product() << std::endl; + + QGaussRadauChebyshev q_10(1); + deallog << "QGaussRadauChebyshev<" << dim << ">: " + << q_10.is_tensor_product() << std::endl; + + QGaussLobattoChebyshev q_11(2); + deallog << "QGaussLobattoChebyshev<" << dim << ">: " + << q_11.is_tensor_product() << std::endl; + + QSorted q_16(q_0); + deallog << "QSorted<" << dim << ">: " + << q_16.is_tensor_product() << std::endl; + + QTelles q_17(1, Point()); + 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 -- 2.39.5