From 9692686965595340e5f5383bb6791edf40f6f333 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 27 Jul 2017 23:44:06 +0200 Subject: [PATCH] Allow tensor products with different formulas in each direction --- doc/news/changes/minor/20170726DanielArndt-1 | 2 +- include/deal.II/base/quadrature.h | 27 +-- source/base/quadrature.cc | 96 +++++--- source/base/quadrature_lib.cc | 15 +- tests/base/quadrature_check_tensor_product.cc | 218 ++++++++++++++++++ .../quadrature_check_tensor_product.output | 44 ++++ tests/base/quadrature_is_tensor_product.cc | 2 +- .../base/quadrature_is_tensor_product.output | 8 +- 8 files changed, 354 insertions(+), 58 deletions(-) create mode 100644 tests/base/quadrature_check_tensor_product.cc create mode 100644 tests/base/quadrature_check_tensor_product.output 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().
(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> &get_tensor_basis() const; protected: /** @@ -255,10 +248,18 @@ protected: 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 + * 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 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 > 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 + void + set_tensor_basis_1d(const Quadrature &, + std::vector> &) + {} + + template <> + void + set_tensor_basis_1d(const Quadrature<1> &quadrature, + std::vector> &tensor_basis) + { + tensor_basis = std::vector >(1, quadrature); + } +} + + + template <> Quadrature<0>::Quadrature (const unsigned int n_q) : @@ -48,7 +67,9 @@ Quadrature::Quadrature (const unsigned int n_q) quadrature_points (n_q, Point()), weights (n_q, 0), is_tensor_product_flag (dim==1) -{} +{ + set_tensor_basis_1d(*this, this->tensor_basis); +} @@ -74,6 +95,7 @@ Quadrature::Quadrature (const std::vector > &points, { Assert (weights.size() == points.size(), ExcDimensionMismatch(weights.size(), points.size())); + set_tensor_basis_1d(*this, this->tensor_basis); } @@ -87,6 +109,7 @@ Quadrature::Quadrature (const std::vector > &points) { Assert(weights.size() == points.size(), ExcDimensionMismatch(weights.size(), points.size())); + set_tensor_basis_1d(*this, this->tensor_basis); } @@ -97,7 +120,9 @@ Quadrature::Quadrature (const Point &point) quadrature_points(std::vector > (1, point)), weights(std::vector (1, 1.)), is_tensor_product_flag (dim==1) -{} +{ + set_tensor_basis_1d(*this, this->tensor_basis); +} template <> @@ -115,7 +140,7 @@ Quadrature::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::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 >(1,q2)) { unsigned int present_index = 0; for (unsigned int i2=0; i2 >(dim, q); } @@ -253,7 +287,8 @@ Quadrature::Quadrature (const Quadrature &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::operator= (const Quadrature &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::memory_consumption () const template -Quadrature<1> +const std::vector > & 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!")); + 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(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); + return tensor_basis; } @@ -347,6 +359,9 @@ QAnisotropic::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::QAnisotropic(const Quadrature<1> &qx) template QAnisotropic::QAnisotropic(const Quadrature<1> &qx, const Quadrature<1> &qy) - : Quadrature(qx.size() - *qy.size()) + : Quadrature(qx.size()*qy.size()) { Assert (dim==2, ExcImpossibleInDim(dim)); unsigned int k=0; @@ -367,6 +381,10 @@ QAnisotropic::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 QAnisotropic::QAnisotropic(const Quadrature<1> &qx, const Quadrature<1> &qy, const Quadrature<1> &qz) - : Quadrature(qx.size() - *qy.size() - *qz.size()) + : Quadrature(qx.size()*qy.size()*qz.size()) { Assert (dim==3, ExcImpossibleInDim(dim)); unsigned int k=0; @@ -391,6 +407,11 @@ QAnisotropic::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 >(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 >(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(points[i])); this->weights[i] = 0.5*w[i]; } + this->tensor_basis = std::vector >(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 >(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 >(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 >(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 >(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 >(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 >(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 >(1, *this); } @@ -1120,6 +1129,7 @@ QTelles<1>::QTelles ( weights[q] = J * weights[q]; } + this->tensor_basis = std::vector >(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 >(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 >(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 >(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 +#include +#include + +template +void +check_tensor_product(const std::vector > &quadratures, + const std::vector &quadrature_names) +{ + Assert(false, ExcNotImplemented()); +} + +template <> +void +check_tensor_product(const std::vector > &quadratures, + const std::vector &quadrature_names) +{ + for (unsigned int i=0; i &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 +void +check_tensor_product(const std::vector > &quadratures, + const std::vector &quadrature_names) +{ + for (unsigned int i=0; i &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 +void +check_tensor_product(const std::vector > &quadratures, + const std::vector &quadrature_names) +{ + for (unsigned int i=0; i &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 +void +fill_quadrature_vector(std::vector > &quadratures, + std::vector &quadrature_names) +{ + quadratures.push_back(Quadrature()); + quadrature_names.push_back("Quadrature"); + + quadratures.push_back(QIterated (QGauss<1>(2), 2)); + quadrature_names.push_back("QIterated"); + + quadratures.push_back(QGauss (2)); + quadrature_names.push_back("QGauss"); + + quadratures.push_back(QGaussLobatto (2)); + quadrature_names.push_back("QGaussLobatto"); + + quadratures.push_back(QMidpoint ()); + quadrature_names.push_back("QMidPoint"); + + quadratures.push_back(QSimpson ()); + quadrature_names.push_back("QSimpson"); + + quadratures.push_back(QTrapez ()); + quadrature_names.push_back("QTrapez"); + + quadratures.push_back(QMilne ()); + quadrature_names.push_back("QMilne"); + + quadratures.push_back(QWeddle ()); + quadrature_names.push_back("QWeddle"); + + quadratures.push_back(QGaussChebyshev (3)); + quadrature_names.push_back("QGaussChebyshev"); + + quadratures.push_back(QGaussRadauChebyshev (2)); + quadrature_names.push_back("QGaussRadauChebyshev"); + + quadratures.push_back(QGaussLobattoChebyshev (2)); + quadrature_names.push_back("QGaussLobattoChebyshev"); + + quadratures.push_back(QSorted (Quadrature())); + quadrature_names.push_back("QSorted"); + + quadratures.push_back(QTelles (1, Point())); + quadrature_names.push_back("QTelles"); +} + +int main() +{ + initlog(); + deallog << std::boolalpha; + + Quadrature<1> q; + + std::vector > quadratures_1d; + std::vector 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 > quadratures_2d; + std::vector 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 > quadratures_3d; + std::vector 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 -- 2.39.5