/**
* 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
/**
* 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
* of weights.
*/
const vector<double> & get_weights () const;
-
- /**
- * Exception
- */
- DeclException0 (ExcInternalError);
-
+
protected:
/**
* List of quadrature points. To be filled
* for a description of the orientation of the different faces.
*/
template <int dim>
-class QProjector {
+class QProjector
+{
public:
/**
* Compute the quadrature points on the
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
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);
};
//---------------------------- quadrature.cc ---------------------------
+#include <grid/geometry_info.h>
#include <base/quadrature.h>
#include <cmath>
// 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);
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;
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;
};
template <>
-const vector<double> & Quadrature<0>::get_weights () const {
+const vector<double> & Quadrature<0>::get_weights () const
+{
Assert (false, ExcInternalError());
return weights;
};
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));
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());
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:
{
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());
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());
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());
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());
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:
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:
};
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)
{
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:
};
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)
{
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:
};
+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)
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;
};
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;
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;
// 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);
// 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;
template class QIterated<1>;
template class QIterated<2>;
template class QIterated<3>;
+template class QProjector<2>;
+template class QProjector<3>;