]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add functionality to get quadrature formulae in higher dimensions by tensor product...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 2 Jul 1998 08:41:37 +0000 (08:41 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 2 Jul 1998 08:41:37 +0000 (08:41 +0000)
git-svn-id: https://svn.dealii.org/trunk@428 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/quadrature.h
deal.II/base/source/quadrature.cc
deal.II/base/source/quadrature_lib.cc

index c2b6a5b88bc5337568fa0726e4d75343050c3489..94f53529f4da273a29b65c4fd78a620cceebdf72 100644 (file)
@@ -1,5 +1,5 @@
 /*----------------------------   quadrature.h     ---------------------------*/
-/*      <Id:>                 */
+/*      $Id$                 */
 #ifndef __quadrature_H
 #define __quadrature_H
 /*----------------------------   quadrature.h     ---------------------------*/
  * 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 <int dim>
@@ -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<dim-1> &,
+               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
index 63182b380b32197b60c6c5eeb5b099dd9074fcf9..957bda9c6e3690b262adfd6ecd112eec2a4c9d75 100644 (file)
@@ -10,6 +10,49 @@ Quadrature<dim>::Quadrature (const unsigned int n_q) :
                weights (n_q, 0) {};
 
 
+template <>
+Quadrature<1>::Quadrature (const Quadrature<0> &,
+                          const Quadrature<1> &) {
+  Assert (false, ExcInternalError());
+};
+
+
+
+template <int dim>
+Quadrature<dim>::Quadrature (const Quadrature<dim-1> &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<q1.n_quadrature_points; ++i)
+    for (unsigned int j=0; j<q2.n_quadrature_points; ++j)
+      {
+                                        // compose coordinates of
+                                        // new quadrature point by tensor
+                                        // product in the last component
+       for (unsigned int d=0; d<dim-1; ++d)
+         quadrature_points[present_index](d)
+           = q1.quad_point(i)(d);
+       quadrature_points[present_index](dim-1)
+         = q2.quad_point(j)(0);
+                                              
+       weights[present_index] = q1.weight(i) * q2.weight(j);
+
+       ++present_index;
+      };
+
+#ifdef DEBUG
+  double sum = 0;
+  for (unsigned int i=0; i<n_quadrature_points; ++i)
+    sum += weights[i];
+  Assert (sum==1, ExcInternalError());
+#endif
+};
+
+
 
 template <int dim>
 Quadrature<dim>::~Quadrature () {};
index 595e72cb5be4669029504dfb203877bfaa82c7e8..b2e89198db9d2c14bd71ea4e0552ce6d3487bde2 100644 (file)
@@ -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<n_quadrature_points; ++i) 
-    {
-      quadrature_points[i] = Point<2>(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 <int dim>
+QGauss2<dim>::QGauss2 () :
+               Quadrature<dim> (QGauss2<dim-1>(), 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<n_quadrature_points; ++i) 
-    {
-      quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
-      weights[i] = wts[i];
-    };
+template <int dim>
+QGauss3<dim>::QGauss3 () :
+               Quadrature<dim> (QGauss3<dim-1>(), 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<n_quadrature_points; ++i) 
-    {
-      quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
-      weights[i] = wts[i];
-    };
+template <int dim>
+QGauss4<dim>::QGauss4 () :
+               Quadrature<dim> (QGauss4<dim-1>(), 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<n_quadrature_points; ++i) 
-    {
-      quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
-      weights[i] = wts[i];
-    };
+template <int dim>
+QGauss5<dim>::QGauss5 () :
+               Quadrature<dim> (QGauss5<dim-1>(), QGauss5<1>()){
+  Assert (dim<3, ExcInternalError());
 };
 
+template <int dim>
+QGauss6<dim>::QGauss6 () :
+               Quadrature<dim> (QGauss6<dim-1>(), QGauss6<1>()){
+  Assert (dim<3, ExcInternalError());
+};
 
+template <int dim>
+QGauss7<dim>::QGauss7 () :
+               Quadrature<dim> (QGauss7<dim-1>(), 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 <int dim>
+QGauss8<dim>::QGauss8 () :
+               Quadrature<dim> (QGauss8<dim-1>(), QGauss8<1>()){
+  Assert (dim<3, ExcInternalError());
 };
 
 
+template <int dim>
+QMidpoint<dim>::QMidpoint () :
+               Quadrature<dim> (QMidpoint<dim-1>(), 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 <int dim>
+QSimpson<dim>::QSimpson () :
+               Quadrature<dim> (QSimpson<dim-1>(), QSimpson<1>()){
+  Assert (dim<3, ExcInternalError());
+};
 
-  for (unsigned int i=0; i<n_quadrature_points; ++i) 
-    {
-      quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
-      weights[i] = wts[i];
-    };
+template <int dim>
+QTrapez<dim>::QTrapez () :
+               Quadrature<dim> (QTrapez<dim-1>(), 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<n_quadrature_points; ++i) 
-    {
-      quadrature_points[i] = Point<2>(xpts[i], ypts[i]);
-      weights[i] = wts[i];
-    };
-};
+// explicite specialization
+template class QGauss2<deal_II_dimension>;
+template class QGauss3<deal_II_dimension>;
+template class QGauss4<deal_II_dimension>;
+template class QGauss5<deal_II_dimension>;
+template class QGauss6<deal_II_dimension>;
+template class QGauss7<deal_II_dimension>;
+template class QGauss8<deal_II_dimension>;
+template class QMidpoint<deal_II_dimension>;
+template class QSimpson<deal_II_dimension>;
+template class QTrapez<deal_II_dimension>;
+
 
 #endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.