From: heltai Date: Sun, 5 Apr 2009 17:24:59 +0000 (+0000) Subject: Added flag to remove singularity from the weighting functions X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9fd5438f57fb8ea448a9314055a3453a3a11fffa;p=dealii-svn.git Added flag to remove singularity from the weighting functions of QGaussOneOverR. git-svn-id: https://svn.dealii.org/trunk@18549 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/quadrature_lib.h b/deal.II/base/include/base/quadrature_lib.h index 940a783ecd..bbc088fe4a 100644 --- a/deal.II/base/include/base/quadrature_lib.h +++ b/deal.II/base/include/base/quadrature_lib.h @@ -397,10 +397,11 @@ protected: /** * Gauss Quadrature Formula with 1/R weighting function. This formula - * is used to to integrate 1/(R)*f(x) on the reference + * can be used to to integrate 1/(R)*f(x) on the reference * element [0,1]^2, where f is a smooth function without * singularities, and R is the distance from the point x to the vertex - * xi, given at construction time by specifying its index. + * xi, given at construction time by specifying its index. Notice that + * this distance is evaluated in the reference element. * * This quadrature formula is obtained from two QGauss quadrature * formula, upon transforming them into polar coordinate system @@ -410,16 +411,53 @@ protected: * the reference element is transformed into a triangle by collapsing * one of the side adjacent to the singularity. The Jacobian of this * transformation contains R, which is removed before scaling the - * original quadrature, and this process is repeated for the next + * original quadrature, and this process is repeated for the next half + * element. + * + * Upon construction it is possible to specify wether we want the + * singularity removed, or not. In other words, this quadrature can be + * used to integrate f(x) = 1/R*g(x), or simply g(x), with the 1/R + * factor already included in the quadrature weights. */ template class QGaussOneOverR : public Quadrature { public: - /** The constructor takes two arguments arguments: the order of - * the gauss formula, and the index of the vertex where the - * singularity is located. */ + /** The constructor takes three arguments: the order of the gauss + * formula, the index of the vertex where the singularity is + * located, and whether we include the weighting singular function + * inside the quadrature, or we leave it in the user function to + * be integrated. + * + * Traditionally, quadrature formulas include their weighting + * function, and the last argument is set to false by + * default. There are cases, however, where this is undesirable + * (for example when you only know that your singularity has the + * same order of 1/R, but cannot be written exactly in this + * way). + * + * In other words, you can use this function in either of + * the following way, obtaining the same result: + * + + QGaussOneOverR singular_quad(order, vertex_id, false); + // This will produce the integral of f(x)/R + for(unsigned int i=0; ivertex(vertex_id)).norm(); + integral += f(singular_quad_noR.point(i))*singular_quad_noR.weight(i)/R; + } + + */ QGaussOneOverR(const unsigned int n, - const unsigned int vertex_index); + const unsigned int vertex_index, + const bool factor_out_singular_weight=false); }; @@ -459,7 +497,7 @@ template <> QMilne<1>::QMilne (); template <> QWeddle<1>::QWeddle (); template <> QGaussLog<1>::QGaussLog (const unsigned int n, const bool revert); template <> QGaussLogR<1>::QGaussLogR (const unsigned int n, const Point<1> x0, const double alpha); -template <> QGaussOneOverR<2>::QGaussOneOverR (const unsigned int n, const unsigned int index); +template <> QGaussOneOverR<2>::QGaussOneOverR (const unsigned int n, const unsigned int index, const bool flag); diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index 16e2e945d2..973ba0d457 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -13,6 +13,7 @@ #include +#include #include #ifdef HAVE_STD_NUMERIC_LIMITS @@ -986,7 +987,8 @@ QGaussLogR<1>::QGaussLogR(const unsigned int n, template<> QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n, - const unsigned int vertex_index) : + const unsigned int vertex_index, + const bool factor_out_singularity) : Quadrature<2>(2*n*n) { // Start with the gauss quadrature formula on the (u,v) reference @@ -1011,7 +1013,8 @@ QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n, // // J = pi/4 R / cos(pi/4 v) // - // And we get rid of R to take into account the singularity. + // And we get rid of R to take into account the singularity, + // unless specified differently in the constructor. std::vector > &ps = this->quadrature_points; std::vector &ws = this->weights; double pi4 = numbers::PI/4; @@ -1021,6 +1024,8 @@ QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n, ps[q][0] = gp[0]; ps[q][1] = gp[0] * std::tan(pi4*gp[1]); ws[q] = gauss.weight(q)*pi4/std::cos(pi4 *gp[1]); + if(factor_out_singularity) + ws[q] *= (ps[q]-GeometryInfo<2>::unit_cell_vertex(0)).norm(); // The other half of the quadrilateral is symmetric with // respect to xy plane. ws[gauss.size()+q] = ws[q];