From: Guido Kanschat Date: Wed, 31 May 2000 04:42:11 +0000 (+0000) Subject: new functions in QProjector; weeded out doc for Quadrature X-Git-Tag: v8.0.0~20431 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bfa8a32c347fc86d8bd87bf7b099c30192dbf765;p=dealii.git new functions in QProjector; weeded out doc for Quadrature git-svn-id: https://svn.dealii.org/trunk@2980 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/quadrature.h b/deal.II/base/include/base/quadrature.h index d75b92c208..56a2903739 100644 --- a/deal.II/base/include/base/quadrature.h +++ b/deal.II/base/include/base/quadrature.h @@ -21,55 +21,28 @@ /** * 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# 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 class Quadrature @@ -108,13 +81,13 @@ class Quadrature /** * Return the #i#th quadrature point. */ - const Point & quad_point (const unsigned int i) const; + const Point & point (const unsigned int i) const; /** * Return a reference to the whole array of * quadrature points. */ - const vector > & get_quad_points () const; + const vector > & get_points () const; /** * Return the weight of the #i#th @@ -127,12 +100,7 @@ class Quadrature * of weights. */ const vector & get_weights () const; - - /** - * Exception - */ - DeclException0 (ExcInternalError); - + protected: /** * List of quadrature points. To be filled @@ -237,7 +205,8 @@ class QIterated : public Quadrature * for a description of the orientation of the different faces. */ template -class QProjector { +class QProjector +{ public: /** * Compute the quadrature points on the @@ -250,6 +219,14 @@ class QProjector { const unsigned int face_no, vector > &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 &quadrature, + vector > &q_points); + /** * Compute the quadrature points on the * cell if the given quadrature formula @@ -264,9 +241,14 @@ class QProjector { vector > &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 &quadrature, + vector > &q_points); }; diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc index 6c817780a5..54ba4c2ab7 100644 --- a/deal.II/base/source/quadrature.cc +++ b/deal.II/base/source/quadrature.cc @@ -12,6 +12,7 @@ //---------------------------- quadrature.cc --------------------------- +#include #include #include @@ -71,9 +72,9 @@ Quadrature::Quadrature (const Quadrature &q1, // product in the last component for (unsigned int d=0; d::~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 -const Point & Quadrature::quad_point (const unsigned int i) const { +const Point & Quadrature::point (const unsigned int i) const { Assert (i -const vector > & Quadrature<0>::get_quad_points () const { +const vector > & Quadrature<0>::get_points () const { Assert (false, ExcInternalError()); return quadrature_points; }; template -const vector > & Quadrature::get_quad_points () const { +const vector > & Quadrature::get_points () const { return quadrature_points; }; @@ -145,7 +146,8 @@ const vector & Quadrature::get_weights () const { template <> -const vector & Quadrature<0>::get_weights () const { +const vector & Quadrature<0>::get_weights () const +{ Assert (false, ExcInternalError()); return weights; }; @@ -154,7 +156,8 @@ const vector & Quadrature<0>::get_weights () const { template <> void QProjector<2>::project_to_face (const Quadrature<1> &quadrature, const unsigned int face_no, - vector > &q_points) { + vector > &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(quadrature.quad_point(p)(0),0); + q_points[p] = Point(quadrature.point(p)(0),0); break; case 1: - q_points[p] = Point(1,quadrature.quad_point(p)(0)); + q_points[p] = Point(1,quadrature.point(p)(0)); break; case 2: - q_points[p] = Point(quadrature.quad_point(p)(0),1); + q_points[p] = Point(quadrature.point(p)(0),1); break; case 3: - q_points[p] = Point(0,quadrature.quad_point(p)(0)); + q_points[p] = Point(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(quadrature.quad_point(p)(0), + q_points[p] = Point(quadrature.point(p)(0), 0, - quadrature.quad_point(p)(1)); + quadrature.point(p)(1)); break; case 1: - q_points[p] = Point(quadrature.quad_point(p)(0), + q_points[p] = Point(quadrature.point(p)(0), 1, - quadrature.quad_point(p)(1)); + quadrature.point(p)(1)); break; case 2: - q_points[p] = Point(quadrature.quad_point(p)(0), - quadrature.quad_point(p)(1), + q_points[p] = Point(quadrature.point(p)(0), + quadrature.point(p)(1), 0); break; case 3: q_points[p] = Point(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(quadrature.quad_point(p)(0), - quadrature.quad_point(p)(1), + q_points[p] = Point(quadrature.point(p)(0), + quadrature.point(p)(1), 1); break; case 5: q_points[p] = Point(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(quadrature.quad_point(p)(0)/2,0); + = Point(quadrature.point(p)(0)/2,0); break; case 1: q_points[p] - = Point(quadrature.quad_point(p)(0)/2+0.5,0); + = Point(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(1,quadrature.quad_point(p)(0)/2); + q_points[p] = Point(1,quadrature.point(p)(0)/2); break; case 1: - q_points[p] = Point(1,quadrature.quad_point(p)(0)/2+0.5); + q_points[p] = Point(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(quadrature.quad_point(p)(0)/2,1); + q_points[p] = Point(quadrature.point(p)(0)/2,1); break; case 1: - q_points[p] = Point(quadrature.quad_point(p)(0)/2+0.5,1); + q_points[p] = Point(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(0,quadrature.quad_point(p)(0)/2); + q_points[p] = Point(0,quadrature.point(p)(0)/2); break; case 1: - q_points[p] = Point(0,quadrature.quad_point(p)(0)/2+0.5); + q_points[p] = Point(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(quadrature.quad_point(p)(0)/2, + q_points[p] = Point(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(quadrature.quad_point(p)(0)/2, + q_points[p] = Point(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(quadrature.quad_point(p)(0)/2, - quadrature.quad_point(p)(1)/2, + q_points[p] = Point(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(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(quadrature.quad_point(p)(0)/2, - quadrature.quad_point(p)(1)/2, + q_points[p] = Point(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(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 +void +QProjector::project_to_faces (const Quadrature &quadrature, + vector > &q_points) +{ + unsigned int npt = quadrature.n_quadrature_points; + unsigned int nf = GeometryInfo::faces_per_cell; + + q_points.resize (npt*nf); + vector > help(npt); + + unsigned k=0; + for (unsigned int i=0;i +void +QProjector::project_to_subfaces (const Quadrature &quadrature, + vector > &q_points) +{ + unsigned int npt = quadrature.n_quadrature_points; + unsigned int nf = GeometryInfo::faces_per_cell; + unsigned int nc = GeometryInfo::subfaces_per_face; + + q_points.resize (npt*nf*nc); + vector > help(npt); + + unsigned k=0; + for (unsigned int i=0;i 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(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(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(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 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(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>;