From ecffc9afe0dc4a5c6096795d7673baa82d04e3b4 Mon Sep 17 00:00:00 2001 From: heltai Date: Sat, 21 Mar 2009 15:14:23 +0000 Subject: [PATCH] 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 --- deal.II/base/include/base/quadrature_lib.h | 30 ++++++++- deal.II/base/source/quadrature_lib.cc | 77 +++++++++++++++++++++- deal.II/doc/news/changes.h | 11 ++++ 3 files changed, 115 insertions(+), 3 deletions(-) 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 x0 = Point(), const double alpha = 1); - - double weight_correction(const Point &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 1/(R)*f(x) on the reference + * element [0,1]^2, 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 +class QGaussOneOverR : public Quadrature { +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 > &ps = this->quadrature_points; + std::vector &ws = this->weights; + double pi4 = numbers::PI/4; + + for(unsigned int q=0; q &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; qbase
    +
  1. +

    + 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. +
    + (Luca Heltai 2009/03/11) +

    +
  2. New: There is now a new QGaussLogR class, that generalizes the QGaussLog class to -- 2.39.5