From: Wolfgang Bangerth Date: Thu, 2 Jul 1998 08:41:37 +0000 (+0000) Subject: Add functionality to get quadrature formulae in higher dimensions by tensor product... X-Git-Tag: v8.0.0~22826 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2826a02c023462363bfb4aeae7f4e38cc33c2f2b;p=dealii.git Add functionality to get quadrature formulae in higher dimensions by tensor product of formulae in lower dimensions. git-svn-id: https://svn.dealii.org/trunk@428 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/quadrature.h b/deal.II/base/include/base/quadrature.h index c2b6a5b88b..94f53529f4 100644 --- a/deal.II/base/include/base/quadrature.h +++ b/deal.II/base/include/base/quadrature.h @@ -1,5 +1,5 @@ /*---------------------------- quadrature.h ---------------------------*/ -/* */ +/* $Id$ */ #ifndef __quadrature_H #define __quadrature_H /*---------------------------- quadrature.h ---------------------------*/ @@ -25,6 +25,12 @@ * third order exact. In general, a formula of order #n# exactly * integrates polynomials of order #2n-1#. * + * Most integration formulae in more than one space dimension are tensor + * products of quadrature formulae in one space dimension, or more + * generally the tensor product of a formula in #(dim-1)# dimensions and + * one in one dimension. There is a special constructor to generate a + * quadrature formula from two others. + * * @author Wolfgang Bangerth, 1998 */ template @@ -39,6 +45,15 @@ class Quadrature { * Constructor. */ Quadrature (const unsigned int n_quadrature_points); + + /** + * Build this quadrature formula as the + * tensor product of a formula in a + * dimension one less than the present and + * a formula in one dimension. + */ + Quadrature (const Quadrature &, + const Quadrature<1> &); /** * Virtual destructor needed, since someone @@ -78,7 +93,11 @@ class Quadrature { int, int, << "The index " << arg1 << " is out of range, it should be less than " << arg2); - + /** + * Exception + */ + DeclException0 (ExcInternalError); + protected: /** * List of quadrature points. To be filled diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc index 63182b380b..957bda9c6e 100644 --- a/deal.II/base/source/quadrature.cc +++ b/deal.II/base/source/quadrature.cc @@ -10,6 +10,49 @@ Quadrature::Quadrature (const unsigned int n_q) : weights (n_q, 0) {}; +template <> +Quadrature<1>::Quadrature (const Quadrature<0> &, + const Quadrature<1> &) { + Assert (false, ExcInternalError()); +}; + + + +template +Quadrature::Quadrature (const Quadrature &q1, + const Quadrature<1> &q2) : + n_quadrature_points (q1.n_quadrature_points * + q2.n_quadrature_points), + quadrature_points (n_quadrature_points), + weights (n_quadrature_points, 0) +{ + unsigned int present_index = 0; + for (unsigned int i=0; i Quadrature::~Quadrature () {}; diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index 595e72cb5b..b2e89198db 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -319,328 +319,90 @@ QTrapez<1>::QTrapez () : -#if deal_II_dimension >= 2 -template <> -QGauss2<2>::QGauss2 () : - Quadrature<2> (4) -{ - // points on [-1,1] - static const double xpts_normal[] = { -sqrt(1./3.), sqrt(1./3.) }; - // weights on [-1,1] - static const double wts_normal[] = { 1., 1. }; +#if deal_II_dimension >= 2 - // points and weights on [0,1]^2 - static const double xpts[] = { (xpts_normal[0]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[1]+1)/2. }; - static const double ypts[] = { (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2. }; - static const double wts[] = { wts_normal[0]/4., - wts_normal[1]/4., - wts_normal[0]/4., - wts_normal[1]/4. }; +// construct the quadrature formulae in higher dimensions by +// tensor product of lower dimensions - for (unsigned int i=0; i(xpts[i], ypts[i]); - weights[i] = wts[i]; - }; +// note: the tensor product code has not been tested for more than +// two space dimensions, so we leave in these Assert'ions until +// that is done +template +QGauss2::QGauss2 () : + Quadrature (QGauss2(), QGauss2<1>()){ + Assert (dim<3, ExcInternalError()); }; - - -template <> -QGauss3<2>::QGauss3 () : - Quadrature<2> (9) -{ - // points on [-1,1] - static const double xpts_normal[] = { -sqrt(3./5.), - 0., - sqrt(3./5.) }; - // weights on [-1,1] - static const double wts_normal[] = { 5./9., - 8./9., - 5./9. }; - - // points and weights on [0,1]^2 - static const double xpts[] = { (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2. }; - static const double ypts[] = { (xpts_normal[0]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[2]+1)/2. }; - static const double wts[] = { wts_normal[0]/2.*wts_normal[0]/2., - wts_normal[1]/2.*wts_normal[0]/2., - wts_normal[2]/2.*wts_normal[0]/2., - wts_normal[0]/2.*wts_normal[1]/2., - wts_normal[1]/2.*wts_normal[1]/2., - wts_normal[2]/2.*wts_normal[1]/2., - wts_normal[0]/2.*wts_normal[2]/2., - wts_normal[1]/2.*wts_normal[2]/2., - wts_normal[2]/2.*wts_normal[2]/2. }; - - for (unsigned int i=0; i(xpts[i], ypts[i]); - weights[i] = wts[i]; - }; +template +QGauss3::QGauss3 () : + Quadrature (QGauss3(), QGauss3<1>()){ + Assert (dim<3, ExcInternalError()); }; - - -template <> -QGauss4<2>::QGauss4 () : - Quadrature<2> (16) -{ - // points on [-1,1] - static const double xpts_normal[] = { -sqrt(1./7.*(3-4*sqrt(0.3))), - -sqrt(1./7.*(3+4*sqrt(0.3))), - +sqrt(1./7.*(3-4*sqrt(0.3))), - +sqrt(1./7.*(3+4*sqrt(0.3))) }; - // weights on [-1,1] - static const double wts_normal[] = { 1./2. + 1./12.*sqrt(10./3.), - 1./2. - 1./12.*sqrt(10./3.), - 1./2. + 1./12.*sqrt(10./3.), - 1./2. - 1./12.*sqrt(10./3.) }; - - // points and weights on [0,1]^2 - static const double xpts[] = { (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[3]+1)/2. }; - static const double ypts[] = { (xpts_normal[0]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[3]+1)/2.}; - static const double wts[] = { wts_normal[0]/2.*wts_normal[0]/2., - wts_normal[1]/2.*wts_normal[0]/2., - wts_normal[2]/2.*wts_normal[0]/2., - wts_normal[3]/2.*wts_normal[0]/2., - wts_normal[0]/2.*wts_normal[1]/2., - wts_normal[1]/2.*wts_normal[1]/2., - wts_normal[2]/2.*wts_normal[1]/2., - wts_normal[3]/2.*wts_normal[1]/2., - wts_normal[0]/2.*wts_normal[2]/2., - wts_normal[1]/2.*wts_normal[2]/2., - wts_normal[2]/2.*wts_normal[2]/2., - wts_normal[3]/2.*wts_normal[2]/2., - wts_normal[0]/2.*wts_normal[3]/2., - wts_normal[1]/2.*wts_normal[3]/2., - wts_normal[2]/2.*wts_normal[3]/2., - wts_normal[3]/2.*wts_normal[3]/2. }; - - for (unsigned int i=0; i(xpts[i], ypts[i]); - weights[i] = wts[i]; - }; +template +QGauss4::QGauss4 () : + Quadrature (QGauss4(), QGauss4<1>()){ + Assert (dim<3, ExcInternalError()); }; - - -template <> -QGauss5<2>::QGauss5 () : - Quadrature<2> (25) -{ - // points on [-1,1] - static const double xpts_normal[] = { -sqrt(1./9.*(5.-2*sqrt(10./7.))), - -sqrt(1./9.*(5.+2*sqrt(10./7.))), - 0, - +sqrt(1./9.*(5.-2*sqrt(10./7.))), - +sqrt(1./9.*(5.+2*sqrt(10./7.))) }; - // weights on [-1,1] - static const double wts_normal[] = { 0.3*(-0.7+5.*sqrt(0.7))/(-2.+5.*sqrt(0.7)), - 0.3*(+0.7+5.*sqrt(0.7))/(+2.+5.*sqrt(0.7)), - 128./225., - 0.3*(-0.7+5.*sqrt(0.7))/(-2.+5.*sqrt(0.7)), - 0.3*(+0.7+5.*sqrt(0.7))/(+2.+5.*sqrt(0.7)) }; - - // points and weights on [0,1]^2 - static const double xpts[] = { (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[4]+1)/2., - - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[4]+1)/2., - - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[4]+1)/2., - - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[4]+1)/2., - - (xpts_normal[0]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[4]+1)/2. }; - - static const double ypts[] = { (xpts_normal[0]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[0]+1)/2., - (xpts_normal[0]+1)/2., - - (xpts_normal[1]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[1]+1)/2., - (xpts_normal[1]+1)/2., - - (xpts_normal[2]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[2]+1)/2., - (xpts_normal[2]+1)/2., - - (xpts_normal[3]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[3]+1)/2., - (xpts_normal[3]+1)/2., - - (xpts_normal[4]+1)/2., - (xpts_normal[4]+1)/2., - (xpts_normal[4]+1)/2., - (xpts_normal[4]+1)/2., - (xpts_normal[4]+1)/2. }; - - static const double wts[] = { wts_normal[0]/2.*wts_normal[0]/2., - wts_normal[1]/2.*wts_normal[0]/2., - wts_normal[2]/2.*wts_normal[0]/2., - wts_normal[3]/2.*wts_normal[0]/2., - wts_normal[4]/2.*wts_normal[0]/2., - - wts_normal[0]/2.*wts_normal[1]/2., - wts_normal[1]/2.*wts_normal[1]/2., - wts_normal[2]/2.*wts_normal[1]/2., - wts_normal[3]/2.*wts_normal[1]/2., - wts_normal[4]/2.*wts_normal[1]/2., - - wts_normal[0]/2.*wts_normal[2]/2., - wts_normal[1]/2.*wts_normal[2]/2., - wts_normal[2]/2.*wts_normal[2]/2., - wts_normal[3]/2.*wts_normal[2]/2., - wts_normal[4]/2.*wts_normal[2]/2., - - wts_normal[0]/2.*wts_normal[3]/2., - wts_normal[1]/2.*wts_normal[3]/2., - wts_normal[2]/2.*wts_normal[3]/2., - wts_normal[3]/2.*wts_normal[3]/2., - wts_normal[4]/2.*wts_normal[3]/2., - - wts_normal[0]/2.*wts_normal[4]/2., - wts_normal[1]/2.*wts_normal[4]/2., - wts_normal[2]/2.*wts_normal[4]/2., - wts_normal[3]/2.*wts_normal[4]/2., - wts_normal[4]/2.*wts_normal[4]/2. }; - - for (unsigned int i=0; i(xpts[i], ypts[i]); - weights[i] = wts[i]; - }; +template +QGauss5::QGauss5 () : + Quadrature (QGauss5(), QGauss5<1>()){ + Assert (dim<3, ExcInternalError()); }; +template +QGauss6::QGauss6 () : + Quadrature (QGauss6(), QGauss6<1>()){ + Assert (dim<3, ExcInternalError()); +}; +template +QGauss7::QGauss7 () : + Quadrature (QGauss7(), QGauss7<1>()){ + Assert (dim<3, ExcInternalError()); +}; -template <> -QMidpoint<2>::QMidpoint () : - Quadrature<2>(1) -{ - quadrature_points[0] = Point<2> (1./2., 1./2.); - weights[0] = 1.0; +template +QGauss8::QGauss8 () : + Quadrature (QGauss8(), QGauss8<1>()){ + Assert (dim<3, ExcInternalError()); }; +template +QMidpoint::QMidpoint () : + Quadrature (QMidpoint(), QMidpoint<1>()){ + Assert (dim<3, ExcInternalError()); +}; -template <> -QSimpson<2>::QSimpson () : - Quadrature<2> (9) -{ - static const double xpts[] = { 0.0, 0.0, 0.0, - 0.5, 0.5, 0.5, - 1.0, 1.0, 1.0 }; - static const double ypts[] = { 0.0, 0.5, 1.0, - 0.0, 0.5, 1.0, - 0.0, 0.5, 1.0 }; - static const double wts[] = { 1./6.*1./6., 2./3.*1./6., 1./6.*1./6., - 1./6.*2./3., 2./3.*2./3., 1./6.*2./3., - 1./6.*1./6., 2./3.*1./6., 1./6.*1./6. }; +template +QSimpson::QSimpson () : + Quadrature (QSimpson(), QSimpson<1>()){ + Assert (dim<3, ExcInternalError()); +}; - for (unsigned int i=0; i(xpts[i], ypts[i]); - weights[i] = wts[i]; - }; +template +QTrapez::QTrapez () : + Quadrature (QTrapez(), QTrapez<1>()){ + Assert (dim<3, ExcInternalError()); }; -template <> -QTrapez<2>::QTrapez () : - Quadrature<2> (4) -{ - static const double xpts[] = { 0.0, 0.0, 1.0, 1.0 }; - static const double ypts[] = { 0.0, 1.0, 0.0, 1.0 }; - static const double wts[] = { 0.25, 0.25, 0.25, 0.25 }; - for (unsigned int i=0; i(xpts[i], ypts[i]); - weights[i] = wts[i]; - }; -}; +// explicite specialization +template class QGauss2; +template class QGauss3; +template class QGauss4; +template class QGauss5; +template class QGauss6; +template class QGauss7; +template class QGauss8; +template class QMidpoint; +template class QSimpson; +template class QTrapez; + #endif