From: heltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sat, 21 Mar 2009 15:14:23 +0000 (+0000)
Subject: Added QGaussOneOverR to integrate 1/R singularities
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ecffc9afe0dc4a5c6096795d7673baa82d04e3b4;p=dealii-svn.git

Added QGaussOneOverR to integrate 1/R singularities
on the vertices of reference element.


git-svn-id: https://svn.dealii.org/trunk@18494 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 fc4429dcb8..d00281f8a3 100644
--- a/deal.II/base/include/base/quadrature_lib.h
+++ b/deal.II/base/include/base/quadrature_lib.h
@@ -387,8 +387,6 @@ public:
     QGaussLogR(const unsigned int n, 
 	       const Point<dim> x0 = Point<dim>(), 
 	       const double alpha = 1);
-    
-    double weight_correction(const Point<dim> &q);
 
 protected:
     /** This is the length of interval (0,origin), or 1 if either of
@@ -397,6 +395,33 @@ 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
+ * 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.
+ *
+ * This quadrature formula is obtained from two QGauss quadrature
+ * formula, upon transforming them into polar coordinate system
+ * centered at the singularity, and then again into another reference
+ * element. This allows for the singularity to be cancelled by part of
+ * the Jacobian of the transformation, which contains R. In practice
+ * 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
+ */
+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. */
+    QGaussOneOverR(const unsigned int n, 
+		   const unsigned int vertex_index);
+};
+
 
 
 /*@}*/
@@ -434,6 +459,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);
 
 
 
diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc
index 5723a507c1..16e2e945d2 100644
--- a/deal.II/base/source/quadrature_lib.cc
+++ b/deal.II/base/source/quadrature_lib.cc
@@ -981,7 +981,82 @@ QGaussLogR<1>::QGaussLogR(const unsigned int n,
 	}
     }
 }
-	
+
+
+
+template<>
+QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n, 
+				  const unsigned int vertex_index) :
+    Quadrature<2>(2*n*n)
+{
+    // Start with the gauss quadrature formula on the (u,v) reference
+    // element.
+    QGauss<2> gauss(n);
+    
+    Assert(gauss.size() == n*n, ExcInternalError());
+    
+    // For the moment we only implemented this for the vertices of a
+    // quadrilateral. We are planning to do this also for the support
+    // points of arbitrary FE_Q elements, to allow the use of this
+    // class in boundary element programs with higher order mappings.
+    Assert(vertex_index < 4, ExcIndexRange(vertex_index, 0, 4));
+    
+    // We create only the first one. All other pieces are rotation of
+    // this one.
+    // In this case the transformation is
+    // 
+    // (x,y) = (u, u tan(pi/4 v))
+    //
+    // with Jacobian
+    //
+    // J = pi/4 R / cos(pi/4 v)
+    //
+    // And we get rid of R to take into account the singularity.
+    std::vector<Point<2> >	&ps = this->quadrature_points;
+    std::vector<double>		&ws = this->weights;
+    double pi4 = numbers::PI/4;
+    
+    for(unsigned int q=0; q<gauss.size(); ++q) {
+	const Point<2> &gp = gauss.point(q);
+	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]);
+	// The other half of the quadrilateral is symmetric with
+	// respect to xy plane.
+	ws[gauss.size()+q]    = ws[q];
+	ps[gauss.size()+q][0] = ps[q][1];
+	ps[gauss.size()+q][1] = ps[q][0];
+    }
+    
+    // Now we distribute these vertices in the correct manner
+    double theta = 0;
+    switch(vertex_index) {
+    case 0:
+	theta = 0;
+	break;
+    case 1:
+	// 
+	theta = numbers::PI/2;
+	break;
+    case 2:
+	theta = -numbers::PI/2;
+	break;
+    case 3:
+	theta = numbers::PI;
+	break;
+    }
+    
+    double R00 =  std::cos(theta), R01 = -std::sin(theta);
+    double R10 =  std::sin(theta), R11 =  std::cos(theta);
+    
+    if(vertex_index != 0)
+	for(unsigned int q=0; q<size(); ++q) {
+	    double x = ps[q][0]-.5,  y = ps[q][1]-.5;
+	    
+	    ps[q][0] = R00*x + R01*y + .5;
+	    ps[q][1] = R10*x + R11*y + .5;
+	}
+}
 
 
 // construct the quadrature formulae in higher dimensions by
diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 4d1a9ad516..df758d9108 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -400,6 +400,17 @@ inconvenience this causes.
 <h3>base</h3>
 
 <ol>
+  <li>
+  <p>
+  New: There is now a new QGaussOneOverR class, that allows for integration
+  on the two dimensional reference element of arbitrary polynomial functions
+  with weight 1/R. This class is only instantiated for dim=2, and it is intended
+  for use with collocation type boundary element methods of order 1, where the 
+  singularities are collocated on the vertices of the quadrilaterals. 
+  <br> 
+  (Luca Heltai 2009/03/11)
+  </p>
+
   <li>
   <p>
   New: There is now a new QGaussLogR class, that generalizes the QGaussLog class to