From: wolf Date: Sat, 5 Jun 1999 14:05:57 +0000 (+0000) Subject: Check in a class which is able to iterate a given quadrature formula n times in each... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=06ebffea7f3b1a7ddb362da5bfb2825a102e6567;p=dealii-svn.git Check in a class which is able to iterate a given quadrature formula n times in each direction. git-svn-id: https://svn.dealii.org/trunk@1369 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/quadrature.h b/deal.II/base/include/base/quadrature.h index 5c78399d41..6a17905eff 100644 --- a/deal.II/base/include/base/quadrature.h +++ b/deal.II/base/include/base/quadrature.h @@ -126,6 +126,59 @@ class Quadrature { +/** + * Quadrature formula constructed by iteration of another quadrature formula in + * each direction. In more than one space dimension, the resulting quadrature + * formula is constructed in the usual way by building the tensor product of + * the respective iterated quadrature formula in one space dimension. + * + * In one space dimension, the given base formula is copied and scaled onto + * a given number of subintervals of length #1/n_copies#. If the quadrature + * formula uses both end points of theunit interval, then in the interior + * of the iterated quadrature formula there would be quadrature points which + * are used twice; we merge them into one with a weight which is the sum + * of the weights of the left- and the rightmost quadrature point. + * + * Since all dimensions higher than one are built up by tensor products of + * one dimensional and #dim-1# dimensional quadrature formulae, the + * argument given to the constructor needs to be a quadrature formula in + * one space dimension, rather than in #dim# dimensions. + * + * @author Wolfgang Bangerth 1999 + */ +template +class QIterated : public Quadrature +{ + public: + /** + * Constructor. Iterate the given + * quadrature formula $n_copies$ times in + * each direction. + */ + QIterated (const Quadrature<1> &base_quadrature, + const unsigned int n_copies); + + /** + * Exception + */ + DeclException0 (ExcSumOfWeightsNotOne); + /** + * Exception + */ + DeclException0 (ExcInvalidQuadratureFormula); + + private: + /** + * check whether the given + * quadrature formula has quadrature + * points at the left and right end points + * of the interval + */ + static bool uses_both_endpoints (const Quadrature<1> &base_quadrature); +}; + + + /** * This class is a helper class to facilitate the usage of quadrature formulae diff --git a/deal.II/base/source/quadrature.cc b/deal.II/base/source/quadrature.cc index fe33ddfb4b..6c70bc76f9 100644 --- a/deal.II/base/source/quadrature.cc +++ b/deal.II/base/source/quadrature.cc @@ -3,6 +3,8 @@ #include +#include + template <> @@ -448,7 +450,134 @@ void QProjector<3>::project_to_subface (const Quadrature<2> &quadrature, +template <> +bool +QIterated<1>::uses_both_endpoints (const Quadrature<1> &base_quadrature) +{ + bool at_left = false, + at_right = false; + for (unsigned int i=0; i(0.0)) + at_left = true; + if (base_quadrature.quad_point(i) == Point<1>(1.0)) + at_right = true; + }; + + return (at_left && at_right); +}; + + + +template <> +QIterated<1>::QIterated (const Quadrature<1> &base_quadrature, + const unsigned int n_copies) : + Quadrature<1> (uses_both_endpoints(base_quadrature) ? + (base_quadrature.n_quadrature_points-1) * n_copies + 1 : + base_quadrature.n_quadrature_points * n_copies) +{ + // we don't have to skip some points + // in order to get a reasonable quadrature + // formula + if (!uses_both_endpoints(base_quadrature)) + { + unsigned int next_point = 0; + for (unsigned int copy=0; copy(base_quadrature.quad_point(q_point)(0) / n_copies + + + (1.0*copy)/n_copies); + weights[next_point] = base_quadrature.weight(q_point) / n_copies; + + ++next_point; + }; + } + else + // skip doubly available points + { + unsigned int next_point = 0; + + // first find out the weights of + // the left and the right boundary + // points. note that these usually + // are but need not necessarily be + // the same + double double_point_weight = 0; + unsigned int n_end_points = 0; + for (unsigned int i=0; i(0.0)) || + (base_quadrature.quad_point(i) == Point<1>(1.0))) + { + double_point_weight += base_quadrature.weight(i); + ++n_end_points; + }; + // scale the weight correctly + double_point_weight /= n_copies; + + // make sure the base quadrature formula + // has only one quadrature point + // per end point + Assert (n_end_points == 2, ExcInvalidQuadratureFormula()); + + + for (unsigned int copy=0; copy 0) && + (base_quadrature.quad_point(q_point) == Point<1>(0.0))) + continue; + + quadrature_points[next_point] = Point<1>(base_quadrature.quad_point(q_point)(0) / n_copies + + + (1.0*copy)/n_copies); + + // if this is the rightmost point of + // 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))) + weights[next_point] = double_point_weight; + else + weights[next_point] = base_quadrature.weight(q_point) / n_copies; + + ++next_point; + }; + }; + +#if DEBUG + double sum_of_weights = 0; + for (unsigned int i=0; i +QIterated::QIterated (const Quadrature<1> &base_quadrature, + const unsigned int N): + Quadrature (QIterated(base_quadrature, N), + QIterated<1>(base_quadrature, N)) +{}; + + + + // explicit instantiations; note: we need them all for all dimensions template class Quadrature<1>; template class Quadrature<2>; template class Quadrature<3>; +template class QIterated<1>; +template class QIterated<2>; +template class QIterated<3>;