]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added QGaussLogR class
authorheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 11 Mar 2009 18:06:11 +0000 (18:06 +0000)
committerheltai <heltai@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 11 Mar 2009 18:06:11 +0000 (18:06 +0000)
git-svn-id: https://svn.dealii.org/trunk@18481 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/quadrature_lib.h
deal.II/base/source/quadrature_lib.cc
deal.II/doc/news/changes.h

index 0606d258ea5973b9f50f96cb2e4101303a8ba50a..fc4429dcb88b85ae40604c4b985263e6b357b031 100644 (file)
@@ -309,24 +309,18 @@ class QWeddle : public Quadrature<dim>
 
 
 /**
- * Gauss Quadrature Formula to integrate <tt>ln[(x-a)/(b-a)]*f(x)</tt> on the
- * interval <tt>[a,b]</tt>, where f is a smooth function without
- * singularities. See <tt>Numerical Recipes</tt> for more information.  The
- * quadrature formula weights are chosen in order to integrate
- * <tt>f(x)*ln[(x-a)/(b-a)]</tt> (that is, the argument of the logarithm goes
- * linearly from 0 to 1 on the interval of integration). So, the only
- * thing one can specify is the function <tt>f(x)</tt>. Using logarithm
- * properties, one can add the integral of <tt>f(x)*ln(b-a)</tt> (calculated
- * with QGauss class formulas) to the integral of <tt>f(x)*ln[(x-a)/(b-a)]</tt>
- * to finally obtain the integral of <tt>f(x)*ln(x-a)</tt>.  It is also
- * possible to integrate with weighting function <tt>ln[(b-x)/(b-a)]</tt> if
- * revert is set to true.
- * 
- * It works up to <tt>n=12</tt>.
+ * Gauss Quadrature Formula with logarithmic weighting function. This
+ * formula is used to to integrate <tt>ln|x|*f(x)</tt> on the interval
+ * <tt>[0,1]</tt>, where f is a smooth function without
+ * singularities. The collection of quadrature points and weights has
+ * been obtained using <tt>Numerical Recipes</tt>.
  *
-//TODO: implement the calculation of points and weights as it is
-//described in Numerical Recipes, instead of passing them to set_*
-//functions, so that n can also be grater than 12.
+ * Notice that only the function <tt>f(x)</tt> should be provided,
+ * i.e., $\int_0^1 f(x) ln|x| dx = \sum_{i=0}^N w_i f(q_i)$. Setting
+ * the @p revert flag to true at construction time swithces the weight
+ * from <tt>ln|x|</tt> to <tt>ln|1-x|</tt>.
+ *
+ * The weights and functions have been tabulated up to order 12.
  *
  */
 template <int dim>
@@ -355,21 +349,54 @@ class QGaussLog : public Quadrature<dim>
   std::vector<double>
   set_quadrature_weights(const unsigned int n) const;
 
-                                    /**  
-                                    * If the singularity is at the rightmost point 
-                                    * of the interval just revert the vectors
-                                    * of points and weights and everything 
-                                    * works.
-                                    */
-  void
-  revert_points_and_weights(std::vector<double> &p,
-                           std::vector<double> &w) const;
-
 };
 
 
 
 
+/**
+ * Gauss Quadrature Formula with arbitrary logarithmic weighting
+ * function. This formula is used to to integrate
+ * <tt>ln(|x-x0|/alpha)*f(x)</tt> on the interval <tt>[0,1]</tt>,
+ * where f is a smooth function without singularities, and x0 and
+ * alpha are given at construction time, and are the location of the
+ * singularity <tt>x0</tt> and an arbitrary scaling factor in the
+ * singularity.
+ *
+ * This quadrature formula is rather expensive, since it uses
+ * internally two Gauss quadrature formulas of order n to integrate
+ * the nonsingular part of the factor, and two GaussLog quadrature
+ * formulas to integrate on the separate segments [0,x0] and
+ * [x0,1]. If the singularity is one of the extremes and the factor
+ * alpha is 1, then this quadrature is the same as QGaussLog.
+ *
+ * Notice again that only the function <tt>f(x)</tt> should be
+ * provided, i.e., $\int_0^1 f(x) ln(|x-x0|/alpha) dx = \sum_{i=0}^N
+ * w_i f(q_i)$. 
+ *
+ * The weights and functions have been tabulated up to order 12.
+ *
+ */
+template<int dim>
+class QGaussLogR : public Quadrature<dim> {
+public:
+    /** The constructor takes three arguments arguments: the order of
+     * the gauss formula on each of the segments [0,x0] and [x0,1],
+     * the actual location of the singularity and the scale factor for
+     * the logarithmic function. */
+    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
+     * the two extremes have been selected. */
+    const double fraction;
+};
+
+
 
 
 /*@}*/
@@ -391,7 +418,6 @@ template <>
 unsigned int 
 QGaussLobatto<1>::gamma(const unsigned int n) const;
 
-template <> void QGaussLog<1>::revert_points_and_weights(std::vector<double> &, std::vector<double> &) const;
 template <> std::vector<double> QGaussLog<1>::set_quadrature_points(const unsigned int) const;
 template <> std::vector<double> QGaussLog<1>::set_quadrature_weights(const unsigned int) const;
 
@@ -407,6 +433,7 @@ template <> QSimpson<1>::QSimpson ();
 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);
 
 
 
index d05d70bc67a7f3497f69d5f7609a812e5c5ea864..5723a507c1ee5298006136e7936181dd002fa914 100644 (file)
@@ -645,16 +645,16 @@ QGaussLog<1>::QGaussLog(const unsigned int n,
 
   std::vector<double> p=set_quadrature_points(n);
   std::vector<double> w=set_quadrature_weights(n);
-  
-  if (revert == true)
-      revert_points_and_weights(p,w);
 
   for (unsigned int i=0; i<this->size(); ++i)
     {
-      this->quadrature_points[i] = Point<1>(p[i]);
-      this->weights[i]           = w[i];
+       // Using the change of variables x=1-t, it's possible to show
+       // that int f(x)ln|1-x| = int f(1-t) ln|t|, which implies that
+       // we can use this quadrature formula also with weight ln|1-x|.
+       this->quadrature_points[i] = revert ? Point<1>(1-p[n-1-i]) : Point<1>(p[i]);
+       this->weights[i]           = revert ? w[n-1-i] : w[i];
     }
-
+  
 }
 
 template <>
@@ -925,29 +925,63 @@ QGaussLog<1>::set_quadrature_weights(const unsigned int n) const
 }
 
 
-template <> 
-void
-QGaussLog<1>::revert_points_and_weights(std::vector<double> &p,
-                                       std::vector<double> &w) const
+template<>
+QGaussLogR<1>::QGaussLogR(const unsigned int n, 
+                         const Point<1> origin,
+                         const double alpha) :
+    Quadrature<1>( ( (origin[0] == 0) || (origin[0] == 1) ) ? 
+                  (alpha == 1 ? n : 2*n ) : 4*n ),
+    fraction( ( (origin[0] == 0) || (origin[0] == 1.) ) ? 1. : origin[0] )
 {
+    // The three quadrature formulas that make this one up. There are
+    // at most two when the origin is one of the extremes, and there is
+    // only one if the origin is one of the extremes and alpha is
+    // equal to one.
+    //
+    // If alpha is different from one, then we need a correction which
+    // is performed with a standard Gauss quadrature rule on each
+    // segment. This is not needed in the standard case where alpha is
+    // equal to one and the origin is on one of the extremes. We
+    // integrate with weight ln|(x-o)/alpha|. In the easy cases, we
+    // only need n quadrature points. In the most difficult one, we
+    // need 2*n points for the first segment, and 2*n points for the
+    // second segment.
+    QGaussLog<1> quad1(n, origin[0] != 0);
+    QGaussLog<1> quad2(n);
+    QGauss<1> quad(n);
+    
+    // Check that the origin is inside 0,1
+    Assert( (fraction >= 0) && (fraction <= 1),
+           ExcMessage("Origin is outside [0,1]."));
 
-  double temp;
-  const unsigned int sz = this->size();
+    // Non singular offset. This is the start of non singular quad
+    // points.
+    unsigned int ns_offset = (fraction == 1) ? n : 2*n;
+    
+    for(unsigned int i=0, j=ns_offset; i<n; ++i, ++j) {
+       // The first i quadrature points are the same as quad1, and
+       // are by default singular.
+       this->quadrature_points[i] = quad1.point(i)*fraction;
+       this->weights[i] = quad1.weight(i)*fraction;
+       
+       // We need to scale with -log|fraction*alpha|
+       if( (alpha != 1) || (fraction != 1) ) {
+           this->quadrature_points[j] = quad.point(i)*fraction;
+           this->weights[j] = -log(alpha/fraction)*quad.weight(i)*fraction;
+       }
 
-  for (unsigned int i=0; i<sz/2; i++)
-    {
-      temp = p[sz-1-i];
-      p[sz-1-i] = p[i];
-      p[i]=temp;
-      
-      temp = w[sz-1-i];
-      w[sz-1-i]=w[i];
-      w[i]=temp;      
+       // In case we need the second quadrature as well, do it now.
+       if(fraction != 1) {
+           this->quadrature_points[i+n] = quad2.point(i)*(1-fraction)+Point<1>(fraction);
+           this->weights[i+n] = quad2.weight(i)*(1-fraction);
+           
+           // We need to scale with -log|fraction*alpha|
+           this->quadrature_points[j+n] = quad.point(i)*(1-fraction)+Point<1>(fraction);
+           this->weights[j+n] = -log(alpha/(1-fraction))*quad.weight(i)*(1-fraction);
+       }
     }
-
 }
-
-
+       
 
 
 // construct the quadrature formulae in higher dimensions by
index 50dd32393a256ab0116f9cd5464fc20740369a54..edd544096bc065dcb219cb882e50fb94e4a519c9 100644 (file)
@@ -400,6 +400,14 @@ inconvenience this causes.
 <h3>base</h3>
 
 <ol>
+  <li>
+  <p>
+  New: There is now a new QGaussLogR class, that generalizes the QGaussLog class to 
+  allow for arbitrary location of singularities, and singularity factors.
+  <br> 
+  (Luca Heltai 2009/03/11)
+  </p>
+
   <li>
   <p>
   New: The FunctionParser class now supports the fparser library's interface to use

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.