From: Luca Heltai Date: Wed, 21 Jul 2010 14:09:29 +0000 (+0000) Subject: Added new feature to QGaussOneOverR X-Git-Tag: v8.0.0~5782 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3e61d73b7ff33d9b62a63639065165d5c00193de;p=dealii.git Added new feature to QGaussOneOverR git-svn-id: https://svn.dealii.org/trunk@21551 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 ffc6ee4ee1..78efbce072 100644 --- a/deal.II/base/include/base/quadrature_lib.h +++ b/deal.II/base/include/base/quadrature_lib.h @@ -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 class QGaussOneOverR : public Quadrature { 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; ivertex(vertex_id)).norm(); + * integral += f(singular_quad_noR.point(i))*singular_quad_noR.weight(i)/R; + * } + * @endcode + */ + QGaussOneOverR(const unsigned int n, + const Point 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 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 singularity, + const unsigned int n); }; diff --git a/deal.II/base/source/quadrature_lib.cc b/deal.II/base/source/quadrature_lib.cc index 3d7b4cdbb9..3331fa173a 100644 --- a/deal.II/base/source/quadrature_lib.cc +++ b/deal.II/base/source/quadrature_lib.cc @@ -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 > quads; + std::vector > 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 &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, diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index cf94deca73..9716af98f1 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -145,8 +145,13 @@ inconvenience this causes.

base

+
    -
  1. None. +
  2. New: QGaussOneOverR now has a new constructor for arbitrary quadrature + points and not only the vertices of the reference cell. +
    + (Luca Heltai 2010/07/21) +