// $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
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
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);
};
}
+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,