]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added new feature to QGaussOneOverR
authorLuca Heltai <luca.heltai@sissa.it>
Wed, 21 Jul 2010 14:09:29 +0000 (14:09 +0000)
committerLuca Heltai <luca.heltai@sissa.it>
Wed, 21 Jul 2010 14:09:29 +0000 (14:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@21551 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 ffc6ee4ee118a574d8c650cf50eb7d343c71650b..78efbce0727dd68c48532567670cd67a6d891b7c 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -456,12 +456,51 @@ template<int dim>
 class QGaussOneOverR : public Quadrature<dim>
 {
   public:
+    /**
+     * This constructor takes three arguments: the order of the Gauss
+     * formula, the point of the reference element in which 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, q_point, 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, q_point, 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;
+     * }
+     * @endcode
+     */
+    QGaussOneOverR(const unsigned int n, 
+                  const Point<dim> singularity,
+                  const bool factor_out_singular_weight=false);
     /**
      * 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.
+     * be integrated. Notice that this is a specialized version of the
+     * previous constructor which works only for the vertices of the
+     * quadrilateral.
      *
      * Traditionally, quadrature formulas include their weighting
      * function, and the last argument is set to false by
@@ -493,6 +532,13 @@ class QGaussOneOverR : public Quadrature<dim>
     QGaussOneOverR(const unsigned int n, 
                   const unsigned int vertex_index,
                   const bool factor_out_singular_weight=false);
+  private:
+   /** Given a quadrature point and a degree n, this function returns
+     * the size of the singular quadrature rule, considering whether
+     * the point is inside the cell, on an edge of the cell, or on a
+     * corner of the cell. */
+    static unsigned int quad_size(const Point<dim> singularity,
+                                 const unsigned int n);
 };
 
 
index 3d7b4cdbb9a53954a05253128548fd817adb16ba..3331fa173a9d0cefa0a4c907c1c88f53f76ea38f 100644 (file)
@@ -990,6 +990,75 @@ QGaussLogR<1>::QGaussLogR(const unsigned int n,
 }
 
 
+template<>
+unsigned int QGaussOneOverR<2>::quad_size(const Point<2> singularity,
+                                                const unsigned int n)
+{
+  double eps=1e-8;
+  bool on_edge=false;
+  bool on_vertex=false;
+  for(unsigned int i=0; i<2; ++i)
+    if( ( std::abs(singularity[i]  ) < eps ) ||
+       ( std::abs(singularity[i]-1) < eps ) )
+      on_edge = true;
+  if(on_edge && (std::abs( (singularity-Point<2>(.5, .5)).square()-.5)
+                < eps) )
+    on_vertex = true;
+  if(on_vertex) return (2*n*n);
+  if(on_edge) return (4*n*n);
+  return (8*n*n);
+}
+
+template<>
+QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n,
+                                 const Point<2> singularity,
+                                 const bool factor_out_singularity) :
+               Quadrature<2>(quad_size(singularity, n))
+{
+                                  // We treat all the cases in the
+                                  // same way. Split the element in 4
+                                  // pieces, measure the area, if
+                                  // it's relevant, add the
+                                  // quadrature connected to that
+                                  // singularity.
+  std::vector<QGaussOneOverR<2> > quads;
+  std::vector<Point<2> > origins;
+                                  // Id of the corner with a
+                                  // singularity
+  quads.push_back(QGaussOneOverR(n, 3, factor_out_singularity));
+  quads.push_back(QGaussOneOverR(n, 2, factor_out_singularity));
+  quads.push_back(QGaussOneOverR(n, 1, factor_out_singularity));
+  quads.push_back(QGaussOneOverR(n, 0, factor_out_singularity));
+
+  origins.push_back(Point<2>(0.,0.));
+  origins.push_back(Point<2>(singularity[0],0.));
+  origins.push_back(Point<2>(0.,singularity[1]));
+  origins.push_back(singularity);
+  
+                                  // Lexycographical ordering.
+      
+  double eps = 1e-8;
+  unsigned int q_id = 0; // Current quad point index.
+  double area = 0;
+  Point<2> dist;
+
+  for(unsigned int box=0; box<4; ++box) 
+    {
+      dist = (singularity-GeometryInfo<2>::unit_cell_vertex(box));
+      dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));      
+      area = dist[0]*dist[1];
+      if(area > eps) 
+       for(unsigned int q=0; q<quads[box].size(); ++q, ++q_id) 
+         {
+           const Point<2> &qp = quads[box].point(q);
+           this->quadrature_points[q_id] =
+             origins[box]+
+             Point<2>(dist[0]*qp[0], dist[1]*qp[1]);
+           this->weights[q_id] = quads[box].weight(q)*area;
+         }
+    }
+}
+
 
 template<>
 QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n,
index cf94deca73ff0b354a6606252dbafd5384a70b07..9716af98f1bcb7331e099df715e3dd977de1b604 100644 (file)
@@ -145,8 +145,13 @@ inconvenience this causes.
 <a name="base"></a>
 <h3>base</h3>
 
+
 <ol>
-  <li>None.
+  <li><p> New: QGaussOneOverR now has a new constructor for arbitrary quadrature
+  points and not only the vertices of the reference cell.
+  <br>
+  (Luca Heltai 2010/07/21)
+  </p>
 </ol>
 
 

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.