]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
new functions in QProjector; weeded out doc for Quadrature
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 31 May 2000 04:42:11 +0000 (04:42 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 31 May 2000 04:42:11 +0000 (04:42 +0000)
git-svn-id: https://svn.dealii.org/trunk@2980 0785d39b-7218-0410-832d-ea1e28bc413d

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

index d75b92c208c6853c23e56976557c1c62b5a0d89b..56a2903739e6b7d5b9273aa3ad6dcdd38bf403ef 100644 (file)
 /**
  * Base class for quadrature formulae in arbitrary dimensions. This class
  * stores quadrature points and weights on the unit line [0,1], unit
- * square [0,1]x[0,1], etc. This information is used together with
- * objects of the \Ref{FiniteElement} class to compute the values stored
- * in the \Ref{FEValues} objects.
+ * square [0,1]x[0,1], etc.
  *
- * There are a number of derived classes, denoting concrete integration
- * formulae. These are named by a prefixed #Q#, the name of the formula
- * (e.g. #Midpoint# or #Gauss#) and finally (for Gauss integration 
- * formulae) the number of quadrature points.
+ * There are a number of derived classes, denoting concrete
+ * integration formulae. Their names names prefixed by #Q#. By now,
+ * there are several Newton-Cotes formulae, @ref{QMidpoint},
+ * @ref{QTrapez} and @ref{QSimpson}, as well as N-point Gauss formulae
+ * @p{QGaussN}. The names refer to the one-dimensional formulae. The
+ * schemes for higher dimensions are tensor products of
+ * these. Therefore, a three-dimensional @ref{QGauss5} formula has 125
+ * quadrature points.
  *
- * For each quadrature formula
- * there exists a number #m#, that denotes the maximal degree of polynomials
- * the formula of integration is exact for. This number is given 
- * in the documentation of each formula. The order of integration is then
- * given by #m+1#, that means that the error representation of the quadrature 
- * formula includes the $(m+1).$ derivative of the function to be integrated.
- * As the $(m+1).$ derivate of polynomials of degree #m# is 0, the order of
- * integration is always one larger than the degree of polynomials the
- * quadrature formula is exact for. For example, the #Midpoint# quadrature
- * formula is exact for polynomials of degree 1 (linear polynomials) and its order
- * of integration is 2, i.e. `The midpoint-formula is of order 2'.
- * 
- * Note the special case of Gauss integration formulae: The number #n# in the
- * n-Point-Gauss Quadrature formula #QGaussn# denotes the number of quadrature
- * points of the formula in one dimension. This formula is exact for polynomials
- * of degree #2n-1# and its order of integration is #2n#. For example,
- * #QGauss2<1># denotes the 2-Point-Gauss quadrature formula in 1 dimension.
- * It is exact for polynomials of degree 3 and its order of integration is 4.
+ * For each quadrature formula we denote by #m#, the maximal degree of
+ * polynomials integrated exactly. This number is given in the
+ * documentation of each formula. The order of the integration error
+ * is #m+1#, that is, the error is the size of the cell two the #m+1#
+ * by the Bramble-Hilbert Lemma. The number #m# is to be found in the
+ * documentation of each concrete formula. For the optimal formulae
+ * @p{QGaussN} we have $m = 2N-1$. The tensor product formulae are
+ * exact on tensor product polynomials of degree #m# in each space
+ * direction, but they are still only of #m+1#st order.
  *
- * 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.
- * For example, the #QGauss2<dim># formulae includes $2^dim$ quadrature points
- * in #dim# dimensions but is still exact for polynomials of degree 3 and its
- * order of integration is 4.
- *
- * For some programs it is necessary to have a quadrature object for faces.
- * These programs fail to link if compiled for only one space dimension,
- * since there quadrature rules for faces just don't make no sense. In
- * order to allow these programs to be linked anyway, for class #Quadrature<0>#
- * all functions are provided in the #quadrature.cc# file, but they will
- * throw exceptions if actually called. The only function which is allowed
- * to be called is the constructor taking one integer, which in this case
- * ignores its parameter, and of course the destructor. Besides this, it is
- * necessary to provide a class #Point<0># to make the compiler happy. This
- * class also does nothing.
- *
- * @author Wolfgang Bangerth, 1998, documentation: Ralf Hartmann, 1999
+ * @author Wolfgang Bangerth, 1998, 1999, 2000
  */
 template <int dim>
 class Quadrature
@@ -108,13 +81,13 @@ class Quadrature
                                     /**
                                      * Return the #i#th quadrature point.
                                      */
-    const Point<dim> & quad_point (const unsigned int i) const;
+    const Point<dim> & point (const unsigned int i) const;
 
                                     /**
                                      * Return a reference to the whole array of
                                      * quadrature points.
                                      */
-    const vector<Point<dim> > & get_quad_points () const;
+    const vector<Point<dim> > & get_points () const;
     
                                     /**
                                      * Return the weight of the #i#th
@@ -127,12 +100,7 @@ class Quadrature
                                      * of weights.
                                      */
     const vector<double> & get_weights () const;
-    
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInternalError);
-    
+
   protected:
                                     /**
                                      * List of quadrature points. To be filled
@@ -237,7 +205,8 @@ class QIterated : public Quadrature<dim>
  *  for a description of the orientation of the different faces.
  */
 template <int dim>
-class QProjector {
+class QProjector
+{
   public:
                                     /**
                                      * Compute the quadrature points on the
@@ -250,6 +219,14 @@ class QProjector {
                                 const unsigned int       face_no,
                                 vector<Point<dim> >     &q_points);
 
+                                    /**
+                                     * Projection to all faces.
+                                     * Generate a formula that integrates
+                                     * over all faces at the same time.
+                                     */
+    static void project_to_faces (const Quadrature<dim-1> &quadrature,
+                                 vector<Point<dim> >     &q_points);
+    
                                     /**
                                      * Compute the quadrature points on the
                                      * cell if the given quadrature formula
@@ -264,9 +241,14 @@ class QProjector {
                                    vector<Point<dim> >     &q_points);
 
                                     /**
-                                     * Exception
+                                     * Projection to all child faces.
+                                     * Project to the children of all
+                                     * faces at the same time. The
+                                     * ordering is first by face,
+                                     * then by subface
                                      */
-    DeclException0 (ExcInternalError);
+    static void project_to_subfaces (const Quadrature<dim-1> &quadrature,
+                                    vector<Point<dim> >     &q_points);
 };
 
 
index 6c817780a5975597b6aa2a84c1971c4bbda5dd62..54ba4c2ab76cd99724248833f44be30b08821ccc 100644 (file)
@@ -12,6 +12,7 @@
 //----------------------------  quadrature.cc  ---------------------------
 
 
+#include <grid/geometry_info.h>
 #include <base/quadrature.h>
 #include <cmath>
 
@@ -71,9 +72,9 @@ Quadrature<dim>::Quadrature (const Quadrature<dim-1> &q1,
                                         // product in the last component
        for (unsigned int d=0; d<dim-1; ++d)
          quadrature_points[present_index](d)
-           = q1.quad_point(i)(d);
+           = q1.point(i)(d);
        quadrature_points[present_index](dim-1)
-         = q2.quad_point(j)(0);
+         = q2.point(j)(0);
                                               
        weights[present_index] = q1.weight(i) * q2.weight(j);
 
@@ -97,7 +98,7 @@ Quadrature<dim>::~Quadrature () {};
 
 
 template <>
-const Point<0> & Quadrature<0>::quad_point (const unsigned int) const {
+const Point<0> & Quadrature<0>::point (const unsigned int) const {
   Assert (false, ExcInternalError());
   static const Point<0> dummy;
   return dummy;
@@ -105,21 +106,21 @@ const Point<0> & Quadrature<0>::quad_point (const unsigned int) const {
 
 
 template <int dim>
-const Point<dim> & Quadrature<dim>::quad_point (const unsigned int i) const {
+const Point<dim> & Quadrature<dim>::point (const unsigned int i) const {
   Assert (i<n_quadrature_points, ExcIndexRange(i, 0, n_quadrature_points));
   return quadrature_points[i];
 };
 
 
 template <>
-const vector<Point<0> > & Quadrature<0>::get_quad_points () const {
+const vector<Point<0> > & Quadrature<0>::get_points () const {
   Assert (false, ExcInternalError());
   return quadrature_points;
 };
 
 
 template <int dim>
-const vector<Point<dim> > & Quadrature<dim>::get_quad_points () const {
+const vector<Point<dim> > & Quadrature<dim>::get_points () const {
   return quadrature_points;
 };
 
@@ -145,7 +146,8 @@ const vector<double> & Quadrature<dim>::get_weights () const {
 
 
 template <>
-const vector<double> & Quadrature<0>::get_weights () const {
+const vector<double> & Quadrature<0>::get_weights () const
+{
   Assert (false, ExcInternalError());
   return weights;
 };
@@ -154,7 +156,8 @@ const vector<double> & Quadrature<0>::get_weights () const {
 template <>
 void QProjector<2>::project_to_face (const Quadrature<1> &quadrature,
                                     const unsigned int   face_no,
-                                    vector<Point<2> >   &q_points) {
+                                    vector<Point<2> >   &q_points)
+{
   const unsigned int dim=2;
   Assert (face_no<2*dim, ExcIndexRange (face_no, 0, 2*dim));
   
@@ -162,16 +165,16 @@ void QProjector<2>::project_to_face (const Quadrature<1> &quadrature,
     switch (face_no)
       {
        case 0:
-             q_points[p] = Point<dim>(quadrature.quad_point(p)(0),0);
+             q_points[p] = Point<dim>(quadrature.point(p)(0),0);
              break;       
        case 1:
-             q_points[p] = Point<dim>(1,quadrature.quad_point(p)(0));
+             q_points[p] = Point<dim>(1,quadrature.point(p)(0));
              break;       
        case 2:
-             q_points[p] = Point<dim>(quadrature.quad_point(p)(0),1);
+             q_points[p] = Point<dim>(quadrature.point(p)(0),1);
              break;       
        case 3:
-             q_points[p] = Point<dim>(0,quadrature.quad_point(p)(0));
+             q_points[p] = Point<dim>(0,quadrature.point(p)(0));
              break;
        default:
              Assert (false, ExcInternalError());
@@ -190,34 +193,34 @@ void QProjector<3>::project_to_face (const Quadrature<2> &quadrature,
     switch (face_no)
       {
        case 0:
-             q_points[p] = Point<dim>(quadrature.quad_point(p)(0),
+             q_points[p] = Point<dim>(quadrature.point(p)(0),
                                       0,
-                                      quadrature.quad_point(p)(1));
+                                      quadrature.point(p)(1));
              break;       
        case 1:
-             q_points[p] = Point<dim>(quadrature.quad_point(p)(0),
+             q_points[p] = Point<dim>(quadrature.point(p)(0),
                                       1,
-                                      quadrature.quad_point(p)(1));
+                                      quadrature.point(p)(1));
              break;       
        case 2:
-             q_points[p] = Point<dim>(quadrature.quad_point(p)(0),
-                                      quadrature.quad_point(p)(1),
+             q_points[p] = Point<dim>(quadrature.point(p)(0),
+                                      quadrature.point(p)(1),
                                       0);
              break;
        case 3:
              q_points[p] = Point<dim>(1,
-                                      quadrature.quad_point(p)(0),
-                                      quadrature.quad_point(p)(1));
+                                      quadrature.point(p)(0),
+                                      quadrature.point(p)(1));
              break;
        case 4:
-             q_points[p] = Point<dim>(quadrature.quad_point(p)(0),
-                                      quadrature.quad_point(p)(1),
+             q_points[p] = Point<dim>(quadrature.point(p)(0),
+                                      quadrature.point(p)(1),
                                       1);
              break;
        case 5:
              q_points[p] = Point<dim>(0,
-                                      quadrature.quad_point(p)(0),
-                                      quadrature.quad_point(p)(1));
+                                      quadrature.point(p)(0),
+                                      quadrature.point(p)(1));
              break;      
              
        default:
@@ -243,11 +246,11 @@ void QProjector<2>::project_to_subface (const Quadrature<1> &quadrature,
                {
                  case 0:
                        q_points[p]
-                         = Point<dim>(quadrature.quad_point(p)(0)/2,0);
+                         = Point<dim>(quadrature.point(p)(0)/2,0);
                        break;
                  case 1:
                        q_points[p]
-                         = Point<dim>(quadrature.quad_point(p)(0)/2+0.5,0);
+                         = Point<dim>(quadrature.point(p)(0)/2+0.5,0);
                        break;
                  default:
                        Assert (false, ExcInternalError());
@@ -257,10 +260,10 @@ void QProjector<2>::project_to_subface (const Quadrature<1> &quadrature,
              switch (subface_no) 
                {
                  case 0:
-                       q_points[p] = Point<dim>(1,quadrature.quad_point(p)(0)/2);
+                       q_points[p] = Point<dim>(1,quadrature.point(p)(0)/2);
                        break;
                  case 1:
-                       q_points[p] = Point<dim>(1,quadrature.quad_point(p)(0)/2+0.5);
+                       q_points[p] = Point<dim>(1,quadrature.point(p)(0)/2+0.5);
                        break;
                  default:
                        Assert (false, ExcInternalError());
@@ -270,10 +273,10 @@ void QProjector<2>::project_to_subface (const Quadrature<1> &quadrature,
              switch (subface_no) 
                {
                  case 0:
-                       q_points[p] = Point<dim>(quadrature.quad_point(p)(0)/2,1);
+                       q_points[p] = Point<dim>(quadrature.point(p)(0)/2,1);
                        break;
                  case 1:
-                       q_points[p] = Point<dim>(quadrature.quad_point(p)(0)/2+0.5,1);
+                       q_points[p] = Point<dim>(quadrature.point(p)(0)/2+0.5,1);
                        break;
                  default:
                        Assert (false, ExcInternalError());
@@ -283,10 +286,10 @@ void QProjector<2>::project_to_subface (const Quadrature<1> &quadrature,
              switch (subface_no) 
                {
                  case 0:
-                       q_points[p] = Point<dim>(0,quadrature.quad_point(p)(0)/2);
+                       q_points[p] = Point<dim>(0,quadrature.point(p)(0)/2);
                        break;
                  case 1:
-                       q_points[p] = Point<dim>(0,quadrature.quad_point(p)(0)/2+0.5);
+                       q_points[p] = Point<dim>(0,quadrature.point(p)(0)/2+0.5);
                        break;
                  default:
                        Assert (false, ExcInternalError());
@@ -315,9 +318,9 @@ void QProjector<3>::project_to_subface (const Quadrature<2> &quadrature,
     switch (face_no)
       {
        case 0:
-             q_points[p] = Point<dim>(quadrature.quad_point(p)(0)/2,
+             q_points[p] = Point<dim>(quadrature.point(p)(0)/2,
                                       0,
-                                      quadrature.quad_point(p)(1)/2);
+                                      quadrature.point(p)(1)/2);
              switch (subface_no) 
                {
                  case 0:
@@ -338,9 +341,9 @@ void QProjector<3>::project_to_subface (const Quadrature<2> &quadrature,
              
              break;       
        case 1:
-             q_points[p] = Point<dim>(quadrature.quad_point(p)(0)/2,
+             q_points[p] = Point<dim>(quadrature.point(p)(0)/2,
                                       1,
-                                      quadrature.quad_point(p)(1)/2);
+                                      quadrature.point(p)(1)/2);
              switch (subface_no) 
                {
                  case 0:
@@ -360,8 +363,8 @@ void QProjector<3>::project_to_subface (const Quadrature<2> &quadrature,
                };
              break;       
        case 2:
-             q_points[p] = Point<dim>(quadrature.quad_point(p)(0)/2,
-                                      quadrature.quad_point(p)(1)/2,
+             q_points[p] = Point<dim>(quadrature.point(p)(0)/2,
+                                      quadrature.point(p)(1)/2,
                                       0);
              switch (subface_no) 
                {
@@ -383,8 +386,8 @@ void QProjector<3>::project_to_subface (const Quadrature<2> &quadrature,
              break;
        case 3:
              q_points[p] = Point<dim>(1,
-                                      quadrature.quad_point(p)(0)/2,
-                                      quadrature.quad_point(p)(1)/2);
+                                      quadrature.point(p)(0)/2,
+                                      quadrature.point(p)(1)/2);
              switch (subface_no) 
                {
                  case 0:
@@ -404,8 +407,8 @@ void QProjector<3>::project_to_subface (const Quadrature<2> &quadrature,
                };
              break;
        case 4:
-             q_points[p] = Point<dim>(quadrature.quad_point(p)(0)/2,
-                                      quadrature.quad_point(p)(1)/2,
+             q_points[p] = Point<dim>(quadrature.point(p)(0)/2,
+                                      quadrature.point(p)(1)/2,
                                       1);
              switch (subface_no) 
                {
@@ -427,8 +430,8 @@ void QProjector<3>::project_to_subface (const Quadrature<2> &quadrature,
              break;
        case 5:
              q_points[p] = Point<dim>(0,
-                                      quadrature.quad_point(p)(0)/2,
-                                      quadrature.quad_point(p)(1)/2);
+                                      quadrature.point(p)(0)/2,
+                                      quadrature.point(p)(1)/2);
              switch (subface_no) 
                {
                  case 0:
@@ -453,6 +456,55 @@ void QProjector<3>::project_to_subface (const Quadrature<2> &quadrature,
 };
 
 
+template <int dim>
+void
+QProjector<dim>::project_to_faces (const Quadrature<dim-1> &quadrature,
+                                  vector<Point<dim> >     &q_points)
+{
+  unsigned int npt = quadrature.n_quadrature_points;
+  unsigned int nf = GeometryInfo<dim>::faces_per_cell;
+  
+  q_points.resize (npt*nf);
+  vector <Point<dim> > help(npt);
+  
+  unsigned k=0;
+  for (unsigned int i=0;i<nf;++i)
+    {
+      project_to_face(quadrature, i, help);
+      for (unsigned int j=0;j<npt;++j)
+       q_points[k++] = help[j];
+    }
+}
+
+
+
+template <int dim>
+void
+QProjector<dim>::project_to_subfaces (const Quadrature<dim-1> &quadrature,
+                                  vector<Point<dim> >     &q_points)
+{
+  unsigned int npt = quadrature.n_quadrature_points;
+  unsigned int nf = GeometryInfo<dim>::faces_per_cell;
+  unsigned int nc = GeometryInfo<dim>::subfaces_per_face;
+  
+  q_points.resize (npt*nf*nc);
+  vector <Point<dim> > help(npt);
+  
+  unsigned k=0;
+  for (unsigned int i=0;i<nf;++i)
+    for (unsigned int c=0;c<nc;++c)
+      {
+       project_to_subface(quadrature, i, c, help);
+       for (unsigned int j=0;j<npt;++j)
+         q_points[k++] = help[j];
+      }
+}
+
+
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
+
+
 template <>
 bool
 QIterated<1>::uses_both_endpoints (const Quadrature<1> &base_quadrature)
@@ -461,9 +513,9 @@ QIterated<1>::uses_both_endpoints (const Quadrature<1> &base_quadrature)
       at_right = false;
   for (unsigned int i=0; i<base_quadrature.n_quadrature_points; ++i)
     {
-      if (base_quadrature.quad_point(i) == Point<1>(0.0))
+      if (base_quadrature.point(i) == Point<1>(0.0))
        at_left = true;
-      if (base_quadrature.quad_point(i) == Point<1>(1.0))
+      if (base_quadrature.point(i) == Point<1>(1.0))
        at_right = true;
     };
 
@@ -487,7 +539,7 @@ QIterated<1>::QIterated (const Quadrature<1> &base_quadrature,
       for (unsigned int copy=0; copy<n_copies; ++copy)
        for (unsigned int q_point=0; q_point<base_quadrature.n_quadrature_points; ++q_point)
          {
-           quadrature_points[next_point] = Point<1>(base_quadrature.quad_point(q_point)(0) / n_copies
+           quadrature_points[next_point] = Point<1>(base_quadrature.point(q_point)(0) / n_copies
                                                     +
                                                     (1.0*copy)/n_copies);
            weights[next_point]           = base_quadrature.weight(q_point) / n_copies;
@@ -510,8 +562,8 @@ QIterated<1>::QIterated (const Quadrature<1> &base_quadrature,
       for (unsigned int i=0; i<base_quadrature.n_quadrature_points; ++i)
                                         // add up the weight if this
                                         // is an endpoint
-       if ((base_quadrature.quad_point(i) == Point<1>(0.0)) ||
-           (base_quadrature.quad_point(i) == Point<1>(1.0)))
+       if ((base_quadrature.point(i) == Point<1>(0.0)) ||
+           (base_quadrature.point(i) == Point<1>(1.0)))
          {
            double_point_weight += base_quadrature.weight(i);
            ++n_end_points;
@@ -532,10 +584,10 @@ for (unsigned int copy=0; copy<n_copies; ++copy)
                                             // since we have already entered it the
                                             // last time
            if ((copy > 0) &&
-               (base_quadrature.quad_point(q_point) == Point<1>(0.0)))
+               (base_quadrature.point(q_point) == Point<1>(0.0)))
              continue;
            
-           quadrature_points[next_point] = Point<1>(base_quadrature.quad_point(q_point)(0) / n_copies
+           quadrature_points[next_point] = Point<1>(base_quadrature.point(q_point)(0) / n_copies
                                                     +
                                                     (1.0*copy)/n_copies);
 
@@ -543,7 +595,7 @@ for (unsigned int copy=0; copy<n_copies; ++copy)
                                             // one of the non-last copies: give
                                             // it the double weight
            if ((copy != n_copies-1) &&
-               (base_quadrature.quad_point(q_point) == Point<1>(1.0)))
+               (base_quadrature.point(q_point) == Point<1>(1.0)))
              weights[next_point] = double_point_weight;
            else
              weights[next_point] = base_quadrature.weight(q_point) / n_copies;
@@ -579,3 +631,5 @@ template class Quadrature<3>;
 template class QIterated<1>;
 template class QIterated<2>;
 template class QIterated<3>;
+template class QProjector<2>;
+template class QProjector<3>;

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.