From 5508df0ee1292d6e953bb319cea7b8598aaad6b5 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Sun, 30 Jul 2017 15:58:54 +0200 Subject: [PATCH] Convert to std::array --- include/deal.II/base/quadrature.h | 14 +- source/base/quadrature.cc | 164 ++++++++++++++------- source/base/quadrature_lib.cc | 17 +-- tests/base/quadrature_is_tensor_product.cc | 8 +- 4 files changed, 123 insertions(+), 80 deletions(-) diff --git a/include/deal.II/base/quadrature.h b/include/deal.II/base/quadrature.h index 2fb29302b2..4aebe91b02 100644 --- a/include/deal.II/base/quadrature.h +++ b/include/deal.II/base/quadrature.h @@ -21,6 +21,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -129,7 +130,7 @@ public: /** * Copy constructor. */ - Quadrature (const Quadrature &q); + Quadrature (const Quadrature &q) = default; /** * Move constructor. Construct a new quadrature object by transferring the @@ -163,13 +164,13 @@ public: /** * Virtual destructor. */ - virtual ~Quadrature (); + virtual ~Quadrature () = default; /** * Assignment operator. Copies contents of #weights and #quadrature_points * as well as size. */ - Quadrature &operator = (const Quadrature &); + Quadrature &operator = (const Quadrature &) = default; /** * Test for equality of two quadratures. @@ -232,7 +233,8 @@ public: * returns the one-dimensional basis objects. * Otherwise, calling this function is not allowed. */ - const std::vector> &get_tensor_basis() const; + typename std::conditional, dim>,const std::array,dim>&>::type + get_tensor_basis() const; protected: /** @@ -259,7 +261,7 @@ protected: * Stores the one-dimensional tensor basis objects in case this object * can be represented by a tensor product. */ - std::vector > tensor_basis; + typename std::conditional,dim>>::type tensor_basis; }; @@ -442,8 +444,6 @@ Quadrature<0>::Quadrature (const Quadrature<-1> &, const Quadrature<1> &); template <> Quadrature<0>::Quadrature (const Quadrature<1> &); -template <> -Quadrature<0>::~Quadrature (); template <> Quadrature<1>::Quadrature (const Quadrature<0> &, diff --git a/source/base/quadrature.cc b/source/base/quadrature.cc index 4a107dc22a..7e6b090997 100644 --- a/source/base/quadrature.cc +++ b/source/base/quadrature.cc @@ -26,25 +26,6 @@ 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) : @@ -55,12 +36,6 @@ Quadrature<0>::Quadrature (const unsigned int n_q) -template <> -Quadrature<0>::~Quadrature () -{} - - - template Quadrature::Quadrature (const unsigned int n_q) : @@ -68,7 +43,6 @@ Quadrature::Quadrature (const unsigned int n_q) weights (n_q, 0), is_tensor_product_flag (dim==1) { - set_tensor_basis_1d(*this, this->tensor_basis); } @@ -95,7 +69,6 @@ Quadrature::Quadrature (const std::vector > &points, { Assert (weights.size() == points.size(), ExcDimensionMismatch(weights.size(), points.size())); - set_tensor_basis_1d(*this, this->tensor_basis); } @@ -109,7 +82,6 @@ Quadrature::Quadrature (const std::vector > &points) { Assert(weights.size() == points.size(), ExcDimensionMismatch(weights.size(), points.size())); - set_tensor_basis_1d(*this, this->tensor_basis); } @@ -119,9 +91,23 @@ Quadrature::Quadrature (const Point &point) : quadrature_points(std::vector > (1, point)), weights(std::vector (1, 1.)), - is_tensor_product_flag (dim==1) + is_tensor_product_flag (true) +{ + for (unsigned int i=0; i > quad_vec_1d (1, Point<1>(point[i])); + tensor_basis[i] = Quadrature<1>(quad_vec_1d, weights); + } +} + + +template <> +Quadrature<1>::Quadrature (const Point<1> &point) + : + quadrature_points(std::vector > (1, point)), + weights(std::vector (1, 1.)), + is_tensor_product_flag (true) { - set_tensor_basis_1d(*this, this->tensor_basis); } @@ -175,8 +161,9 @@ Quadrature::Quadrature (const SubQuadrature &q1, if (is_tensor_product_flag) { - tensor_basis = q1.get_tensor_basis(); - tensor_basis.push_back(q2); + for (unsigned int i=0; i::Quadrature (const SubQuadrature &, : quadrature_points (q2.size()), weights (q2.size()), - is_tensor_product_flag (true), - tensor_basis(std::vector >(1,q2)) + is_tensor_product_flag (true) { unsigned int present_index = 0; for (unsigned int i2=0; i2 >(dim, q); + tensor_basis[i] = q; } -template +/*template Quadrature::Quadrature (const Quadrature &q) : Subscriptor(), quadrature_points (q.quadrature_points), weights (q.weights), - is_tensor_product_flag (q.is_tensor_product_flag), - tensor_basis (q.tensor_basis) -{} + is_tensor_product_flag (q.is_tensor_product_flag) +{ + if (is_tensor_product_flag) + tensor_basis = q.tensor_basis; +} -template + +template <> +Quadrature<1>::Quadrature (const Quadrature<1> &q) + : + Subscriptor(), + quadrature_points (q.quadrature_points), + weights (q.weights), + is_tensor_product_flag (q.is_tensor_product_flag) +{}*/ + + + +/*template Quadrature & Quadrature::operator= (const Quadrature &q) { @@ -303,6 +303,38 @@ Quadrature::operator= (const Quadrature &q) return *this; } +template <> +Quadrature<1> & +Quadrature<1>::operator= (const Quadrature<1> &q) +{ + std::cout << tensor_basis << std::endl; + + Assert (q.size() == q.weights.size(), ExcInternalError()); + Assert (q.size() == q.quadrature_points.size(), ExcInternalError()); + for (unsigned int i=0; i @@ -316,12 +348,6 @@ Quadrature::operator == (const Quadrature &q) const -template -Quadrature::~Quadrature () -{} - - - template std::size_t Quadrature::memory_consumption () const @@ -333,19 +359,31 @@ Quadrature::memory_consumption () const template -const std::vector > & +typename std::conditional, dim>,const std::array,dim>&>::type 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()); return tensor_basis; } +template <> +std::array, 1> +Quadrature<1>::get_tensor_basis () const +{ + Assert (this->is_tensor_product_flag == true, + ExcMessage("This function only makes sense if " + "this object represents a tensor product!")); + + return std::array, 1> {*this}; +} + + + //--------------------------------------------------------------------------- template QAnisotropic::QAnisotropic(const Quadrature<1> &qx) @@ -360,18 +398,25 @@ QAnisotropic::QAnisotropic(const Quadrature<1> &qx) } Assert (k==this->size(), ExcInternalError()); this->is_tensor_product_flag = true; - this->tensor_basis.resize(1); - this->tensor_basis[0] = qx; + //this->tensor_basis[0] = qx; } template -QAnisotropic::QAnisotropic(const Quadrature<1> &qx, - const Quadrature<1> &qy) - : Quadrature(qx.size()*qy.size()) +QAnisotropic::QAnisotropic(const Quadrature<1> &, + const Quadrature<1> &) { Assert (dim==2, ExcImpossibleInDim(dim)); +} + + + +template <> +QAnisotropic<2>::QAnisotropic(const Quadrature<1> &qx, + const Quadrature<1> &qy) + : Quadrature<2>(qx.size()*qy.size()) +{ unsigned int k=0; for (unsigned int k2=0; k2::QAnisotropic(const Quadrature<1> &qx, } 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; } @@ -396,6 +440,16 @@ QAnisotropic::QAnisotropic(const Quadrature<1> &qx, : Quadrature(qx.size()*qy.size()*qz.size()) { Assert (dim==3, ExcImpossibleInDim(dim)); +} + + + +template <> +QAnisotropic<3>::QAnisotropic(const Quadrature<1> &qx, + const Quadrature<1> &qy, + const Quadrature<1> &qz) + : Quadrature<3>(qx.size()*qy.size()*qz.size()) +{ unsigned int k=0; for (unsigned int k3=0; k3::QAnisotropic(const Quadrature<1> &qx, } 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; @@ -1812,7 +1865,6 @@ 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 8e19826722..23fc5b2ba6 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -159,7 +159,6 @@ QGauss<1>::QGauss (const unsigned int n) this->weights[i-1] = w; this->weights[n-i] = w; } - this->tensor_basis = std::vector >(1, *this); } @@ -180,7 +179,6 @@ 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); } @@ -338,7 +336,6 @@ QMidpoint<1>::QMidpoint () { this->quadrature_points[0] = Point<1>(0.5); this->weights[0] = 1.0; - this->tensor_basis = std::vector >(1, *this); } @@ -356,7 +353,6 @@ QTrapez<1>::QTrapez () this->quadrature_points[i] = Point<1>(xpts[i]); this->weights[i] = wts[i]; }; - this->tensor_basis = std::vector >(1, *this); } @@ -374,7 +370,6 @@ QSimpson<1>::QSimpson () this->quadrature_points[i] = Point<1>(xpts[i]); this->weights[i] = wts[i]; }; - this->tensor_basis = std::vector >(1, *this); } @@ -392,7 +387,6 @@ QMilne<1>::QMilne () this->quadrature_points[i] = Point<1>(xpts[i]); this->weights[i] = wts[i]; }; - this->tensor_basis = std::vector >(1, *this); } @@ -412,7 +406,6 @@ QWeddle<1>::QWeddle () this->quadrature_points[i] = Point<1>(xpts[i]); this->weights[i] = wts[i]; }; - this->tensor_basis = std::vector >(1, *this); } @@ -433,7 +426,6 @@ 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 <> @@ -770,7 +762,6 @@ 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); } @@ -948,6 +939,10 @@ QSorted::QSorted(const Quadrature &quad) : std::placeholders::_1, std::placeholders::_2)); + this->quadrature_points = quad.get_points(); + Assert (quad.size() == quad.get_points().size(), ExcInternalError()); + Assert (quad.size() == quad.get_weights().size(), ExcInternalError()); + Assert (this->quadrature_points.size()==quad.size(), ExcInternalError()); for (unsigned int i=0; iweights[i] = quad.weight(permutation[i]); @@ -1129,7 +1124,6 @@ QTelles<1>::QTelles ( weights[q] = J * weights[q]; } - this->tensor_basis = std::vector >(1, *this); } @@ -1184,7 +1178,6 @@ 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); } @@ -1274,7 +1267,6 @@ 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); } @@ -1339,7 +1331,6 @@ 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_is_tensor_product.cc b/tests/base/quadrature_is_tensor_product.cc index d551d15f6f..0f076b59f8 100644 --- a/tests/base/quadrature_is_tensor_product.cc +++ b/tests/base/quadrature_is_tensor_product.cc @@ -28,7 +28,7 @@ template void print_is_tensor_product() { - Quadrature q_0; + Quadrature q_0 (1); deallog << "Quadrature<" << dim << ">: " << q_0.is_tensor_product() << std::endl; @@ -76,7 +76,7 @@ print_is_tensor_product() deallog << "QGaussLobattoChebyshev<" << dim << ">: " << q_11.is_tensor_product() << std::endl; - QSorted q_16(q_0); + QSorted q_16(QGauss(2)); deallog << "QSorted<" << dim << ">: " << q_16.is_tensor_product() << std::endl; @@ -90,7 +90,7 @@ int main() initlog(); deallog << std::boolalpha; - Quadrature<1> q; + QGauss<1> q(1); print_is_tensor_product<1>(); QAnisotropic<1> q_1(q); @@ -103,7 +103,7 @@ int main() 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>()); + QGaussOneOverR<2> q_15(2, Point<2>()); deallog << "QGaussOneOverR<2>: " << q_15.is_tensor_product() << std::endl; print_is_tensor_product<3>(); -- 2.39.5