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 940a783ecd89b34e958814339fa89ec4a7279008..bbc088fe4a86d095c714be50b07458635f887278 100644 (file)
--- 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 <tt>1/(R)*f(x)</tt> on the reference
+ * can be used to to integrate <tt>1/(R)*f(x)</tt> on the reference
  * element <tt>[0,1]^2</tt>, 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<int dim>
 class QGaussOneOverR : public Quadrature<dim> {
 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:
+     *
+     <code> 
+     QGaussOneOverR singular_quad(order, vertex_id, false);
+     // This will produce the integral of f(x)/R
+     for(unsigned int i=0; i<singular_quad.size(); ++i)
+       integral += f(singular_quad.point(i))*singular_quad.weight(i);
+     
+     // And the same here
+     QGaussOneOverR singular_quad_noR(order, vertex_id, true);
+
+     // This also will produce the integral of f(x)/R, but 1/R has to
+     // be specified.
+     for(unsigned int i=0; i<singular_quad.size(); ++i) {
+        double R = (singular_quad_noR.point(i)-cell->vertex(vertex_id)).norm();
+       integral += f(singular_quad_noR.point(i))*singular_quad_noR.weight(i)/R;
+     }
+     </code>
+     */
     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 16e2e945d25c07cf9338799ea50a8b238c7fa9b3..973ba0d4577d5748adff69f1f8fe9b154c7a6882 100644 (file)
--- a/deal.II/base/source/quadrature_lib.cc
+++ b/deal.II/base/source/quadrature_lib.cc
@@ -13,6 +13,7 @@
 
 
 #include <base/quadrature_lib.h>
+#include <base/geometry_info.h>
 #include <cmath>
 
 #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<Point<2> >     &ps = this->quadrature_points;
     std::vector<double>                &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];