← Back to branch summary

~ubuntu-branches/ubuntu/vivid/deal.ii/vivid

« back to all changes in this revision

Viewing changes to base/source/quadrature_lib.cc

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV, Adam C. Powell, IV, Denis Barbier
  • Date: 2010-07-29 13:47:01 UTC
  • mfrom: (3.1.3 sid)
  • Revision ID: james.westby@ubuntu.com-20100729134701-akb8jb3stwge8tcm
Tags: 6.3.1-1
http://bugs.debian.org/562332
[ Adam C. Powell, IV ]
* Changed to source format 3.0 (quilt).
* Changed maintainer to debian-science with Adam Powell as uploader.
* Added source lintian overrides about Adam Powell's name.
* Added Vcs info on git repository.
* Bumped Standards-Version.
* Changed stamp-patch to patch target and fixed its application criterion.
* Moved make_dependencies and expand_instantiations to a versioned directory
  to avoid shlib package conflicts.

[ Denis Barbier ]
* New upstream release (closes: #562332).
  + Added libtbb support.
  + Forward-ported all patches.
* Updates for new PETSc version, including workaround for different versions
  of petsc and slepc.
* Add debian/watch.
* Update to debhelper 7.
* Added pdebuild patch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
base/source/quadrature_lib.cc
1
1
//---------------------------------------------------------------------------
2
 
//    $Id: quadrature_lib.cc 18647 2009-04-18 11:02:27Z heltai $
 
2
//    $Id: quadrature_lib.cc 21206 2010-06-14 20:34:45Z bangerth $
3
3
//    Version: $Name$
4
4
//
5
 
//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2007, 2008 by the deal.II authors
 
5
//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
6
6
//
7
7
//    This file is subject to QPL and may not be  distributed
8
8
//    without copyright and license information. Please refer
71
71
{
72
72
  if (n == 0)
73
73
    return;
74
 
  
 
74
 
75
75
  const unsigned int m = (n+1)/2;
76
76
 
77
77
                                   // tolerance for the Newton
165
165
      double x = .5*z;
166
166
      this->quadrature_points[i-1] = Point<1>(.5-x);
167
167
      this->quadrature_points[n-i] = Point<1>(.5+x);
168
 
      
 
168
 
169
169
      double w = 1./((1.-z*z)*pp*pp);
170
170
      this->weights[i-1] = w;
171
171
      this->weights[n-i] = w;
174
174
 
175
175
 
176
176
template <>
177
 
QGaussLobatto<1>::QGaussLobatto (unsigned int n)
 
177
QGaussLobatto<1>::QGaussLobatto (const unsigned int n)
178
178
                :
179
179
                Quadrature<1> (n)
180
180
{
181
181
  Assert (n >= 2, ExcNotImplemented());
182
 
  
 
182
 
183
183
  std::vector<long double> points  = compute_quadrature_points(n, 1, 1);
184
184
  std::vector<long double> w       = compute_quadrature_weights(points, 0, 0);
185
185
 
202
202
{
203
203
  const unsigned int m = q-2; // no. of inner points
204
204
  std::vector<long double> x(m);
205
 
  
 
205
 
206
206
                                   // compute quadrature points with
207
207
                                   // a Newton algorithm.
208
208
 
229
229
       :
230
230
       double_eps * 5
231
231
       );
232
 
  
 
232
 
233
233
                                   // we take the zeros of the Chebyshev
234
234
                                   // polynomial (alpha=beta=-0.5) as
235
235
                                   // initial values:
236
236
  for (unsigned int i=0; i<m; ++i)
237
 
    x[i] = - std::cos( (long double) (2*i+1)/(2*m) * M_PI );
238
 
  
 
237
    x[i] = - std::cos( (long double) (2*i+1)/(2*m) * numbers::PI );
 
238
 
239
239
  long double r, s, J_x, f, delta;
240
240
  for (unsigned int k=0; k<m; ++k)
241
241
    {
243
243
      if (k>0)
244
244
        r = (r + x[k-1])/2;
245
245
 
246
 
      do 
 
246
      do
247
247
        {
248
248
          s = 1.;
249
249
          for (unsigned int i=0; i<k; ++i)
255
255
          r += delta;
256
256
        }
257
257
      while (std::fabs(delta) >= tolerance);
258
 
      
 
258
 
259
259
      x[k] = r;
260
260
    } // for
261
261
 
277
277
  const unsigned int q = x.size();
278
278
  std::vector<long double> w(q);
279
279
  long double s = 0.L;
280
 
  
 
280
 
281
281
  const long double factor = std::pow(2., alpha+beta+1) *
282
282
                             gamma(alpha+q) *
283
283
                             gamma(beta+q) /
305
305
                                   // using a recursion formula.
306
306
  std::vector<long double> p(n+1);
307
307
  int v, a1, a2, a3, a4;
308
 
  
 
308
 
309
309
                                   // initial values P_0(x), P_1(x):
310
310
  p[0] = 1.0L;
311
311
  if (n==0) return p[0];
312
312
  p[1] = ((alpha+beta+2)*x + (alpha-beta))/2;
313
313
  if (n==1) return p[1];
314
 
  
 
314
 
315
315
  for (unsigned int i=1; i<=(n-1); ++i)
316
316
    {
317
317
      v  = 2*i + alpha + beta;
328
328
 
329
329
 
330
330
template <>
331
 
unsigned int QGaussLobatto<1>::gamma(const unsigned int n) const
 
331
long double QGaussLobatto<1>::gamma(const unsigned int n) const
332
332
{
333
 
  unsigned int result = n - 1;
 
333
  long double result = n - 1;
334
334
  for (int i=n-2; i>1; --i)
335
335
    result *= i;
336
336
  return result;
354
354
  static const double wts[]  = { wts_normal[0]/2.,
355
355
                                 wts_normal[1]/2.  };
356
356
 
357
 
  for (unsigned int i=0; i<this->size(); ++i) 
 
357
  for (unsigned int i=0; i<this->size(); ++i)
358
358
    {
359
359
      this->quadrature_points[i] = Point<1>(xpts[i]);
360
360
      this->weights[i] = wts[i];
385
385
                                 wts_normal[1]/2.,
386
386
                                 wts_normal[2]/2. };
387
387
 
388
 
  for (unsigned int i=0; i<this->size(); ++i) 
 
388
  for (unsigned int i=0; i<this->size(); ++i)
389
389
    {
390
390
      this->quadrature_points[i] = Point<1>(xpts[i]);
391
391
      this->weights[i] = wts[i];
395
395
 
396
396
 
397
397
template <>
398
 
QGauss4<1>::QGauss4 () 
 
398
QGauss4<1>::QGauss4 ()
399
399
                :
400
400
                Quadrature<1> (4)
401
401
{
420
420
                                 wts_normal[2]/2.,
421
421
                                 wts_normal[3]/2. };
422
422
 
423
 
  for (unsigned int i=0; i<this->size(); ++i) 
 
423
  for (unsigned int i=0; i<this->size(); ++i)
424
424
    {
425
425
      this->quadrature_points[i] = Point<1>(xpts[i]);
426
426
      this->weights[i] = wts[i];
458
458
                                 wts_normal[3]/2.,
459
459
                                 wts_normal[4]/2. };
460
460
 
461
 
  for (unsigned int i=0; i<this->size(); ++i) 
 
461
  for (unsigned int i=0; i<this->size(); ++i)
462
462
    {
463
463
      this->quadrature_points[i] = Point<1>(xpts[i]);
464
464
      this->weights[i] = wts[i];
501
501
                                 wts_normal[4]/2.,
502
502
                                 wts_normal[5]/2. };
503
503
 
504
 
  for (unsigned int i=0; i<this->size(); ++i) 
 
504
  for (unsigned int i=0; i<this->size(); ++i)
505
505
    {
506
506
      this->quadrature_points[i] = Point<1>(xpts[i]);
507
507
      this->weights[i] = wts[i];
548
548
                                 wts_normal[5]/2.,
549
549
                                 wts_normal[6]/2. };
550
550
 
551
 
  for (unsigned int i=0; i<this->size(); ++i) 
 
551
  for (unsigned int i=0; i<this->size(); ++i)
552
552
    {
553
553
      this->quadrature_points[i] = Point<1>(xpts[i]);
554
554
      this->weights[i] = wts[i];
576
576
  static const double xpts[] = { 0.0, 1.0 };
577
577
  static const double wts[]  = { 0.5, 0.5 };
578
578
 
579
 
  for (unsigned int i=0; i<this->size(); ++i) 
 
579
  for (unsigned int i=0; i<this->size(); ++i)
580
580
    {
581
581
      this->quadrature_points[i] = Point<1>(xpts[i]);
582
582
      this->weights[i] = wts[i];
593
593
  static const double xpts[] = { 0.0, 0.5, 1.0 };
594
594
  static const double wts[]  = { 1./6., 2./3., 1./6. };
595
595
 
596
 
  for (unsigned int i=0; i<this->size(); ++i) 
 
596
  for (unsigned int i=0; i<this->size(); ++i)
597
597
    {
598
598
      this->quadrature_points[i] = Point<1>(xpts[i]);
599
599
      this->weights[i] = wts[i];
610
610
  static const double xpts[] = { 0.0, .25, .5, .75, 1.0 };
611
611
  static const double wts[]  = { 7./90., 32./90., 12./90., 32./90., 7./90. };
612
612
 
613
 
  for (unsigned int i=0; i<this->size(); ++i) 
 
613
  for (unsigned int i=0; i<this->size(); ++i)
614
614
    {
615
615
      this->quadrature_points[i] = Point<1>(xpts[i]);
616
616
      this->weights[i] = wts[i];
629
629
                                 27./840., 216./840., 41./840.
630
630
  };
631
631
 
632
 
  for (unsigned int i=0; i<this->size(); ++i) 
 
632
  for (unsigned int i=0; i<this->size(); ++i)
633
633
    {
634
634
      this->quadrature_points[i] = Point<1>(xpts[i]);
635
635
      this->weights[i] = wts[i];
639
639
 
640
640
template <>
641
641
QGaussLog<1>::QGaussLog(const unsigned int n,
642
 
                        const bool revert) 
643
 
                 : 
 
642
                        const bool revert)
 
643
                 :
644
644
                 Quadrature<1> (n)
645
645
{
646
646
 
655
655
        this->quadrature_points[i] = revert ? Point<1>(1-p[n-1-i]) : Point<1>(p[i]);
656
656
        this->weights[i]           = revert ? w[n-1-i] : w[i];
657
657
    }
658
 
  
 
658
 
659
659
}
660
660
 
661
661
template <>
662
662
std::vector<double>
663
663
QGaussLog<1>::set_quadrature_points(const unsigned int n) const
664
664
{
665
 
 
 
665
 
666
666
  std::vector<double> points(n);
667
667
 
668
668
  switch (n)
669
 
    {    
 
669
    {
670
670
    case 1:
671
671
      points[0] = 0.3333333333333333;
672
 
      break; 
 
672
      break;
673
673
 
674
674
    case 2:
675
675
      points[0] = 0.1120088061669761;
676
676
      points[1] = 0.6022769081187381;
677
 
      break;    
 
677
      break;
678
678
 
679
679
    case 3:
680
680
      points[0] = 0.06389079308732544;
694
694
      points[1] = 0.1739772133208974;
695
695
      points[2] =  0.4117025202849029;
696
696
      points[3] = 0.6773141745828183;
697
 
      points[4] = 0.89477136103101;  
 
697
      points[4] = 0.89477136103101;
698
698
      break;
699
699
 
700
700
    case 6:
767
767
      points[9] = 0.914193921640008;
768
768
      points[10] = 0.973860256264123;
769
769
      break;
770
 
      
771
 
    case 12:   
 
770
 
 
771
    case 12:
772
772
      points[0] = 0.006548722279080035;
773
773
      points[1] = 0.03894680956045022;
774
774
      points[2] = 0.0981502631060046;
781
781
      points[9] = 0.850240024421055;
782
782
      points[10] = 0.926749682988251;
783
783
      points[11] = 0.977756129778486;
784
 
      break;   
 
784
      break;
785
785
 
786
786
    default:
787
787
      Assert(false, ExcNotImplemented());
788
788
      break;
789
789
    }
790
 
    
 
790
 
791
791
  return points;
792
792
}
793
793
 
796
796
std::vector<double>
797
797
QGaussLog<1>::set_quadrature_weights(const unsigned int n) const
798
798
{
799
 
 
 
799
 
800
800
  std::vector<double> weights(n);
801
801
 
802
802
  switch (n)
803
 
    {    
 
803
    {
804
804
    case 1:
805
805
      weights[0] = -1.0;
806
 
      break; 
 
806
      break;
807
807
    case 2:
808
808
      weights[0] = -0.7185393190303845;
809
809
      weights[1] = -0.2814606809696154;
898
898
      weights[7] = -0.04054160079499477;
899
899
      weights[8] = -0.01943540249522013;
900
900
      weights[9] = -0.006737429326043388;
901
 
      weights[10] = -0.001152486965101561;      
 
901
      weights[10] = -0.001152486965101561;
902
902
      break;
903
 
      
904
 
    case 12:   
 
903
 
 
904
    case 12:
905
905
      weights[0] = -0.09319269144393;
906
906
      weights[1] = -0.1497518275763289;
907
907
      weights[2] = -0.166557454364573;
913
913
      weights[8] = -0.03007108900074863;
914
914
      weights[9] = -0.01424924540252916;
915
915
      weights[10] = -0.004899924710875609;
916
 
      weights[11] = -0.000834029009809656;     
 
916
      weights[11] = -0.000834029009809656;
917
917
      break;
918
918
 
919
919
    default:
920
920
      Assert(false, ExcNotImplemented());
921
 
      break;  
 
921
      break;
922
922
    }
923
923
 
924
924
  return weights;
927
927
 
928
928
 
929
929
template<>
930
 
QGaussLogR<1>::QGaussLogR(const unsigned int n, 
 
930
QGaussLogR<1>::QGaussLogR(const unsigned int n,
931
931
                          const Point<1> origin,
932
932
                          const double alpha,
933
933
                          const bool factor_out_singularity) :
934
 
    Quadrature<1>( ( (origin[0] == 0) || (origin[0] == 1) ) ? 
 
934
    Quadrature<1>( ( (origin[0] == 0) || (origin[0] == 1) ) ?
935
935
                   (alpha == 1 ? n : 2*n ) : 4*n ),
936
936
    fraction( ( (origin[0] == 0) || (origin[0] == 1.) ) ? 1. : origin[0] )
937
937
{
951
951
    QGaussLog<1> quad1(n, origin[0] != 0);
952
952
    QGaussLog<1> quad2(n);
953
953
    QGauss<1> quad(n);
954
 
    
 
954
 
955
955
    // Check that the origin is inside 0,1
956
956
    Assert( (fraction >= 0) && (fraction <= 1),
957
957
            ExcMessage("Origin is outside [0,1]."));
959
959
    // Non singular offset. This is the start of non singular quad
960
960
    // points.
961
961
    unsigned int ns_offset = (fraction == 1) ? n : 2*n;
962
 
    
 
962
 
963
963
    for(unsigned int i=0, j=ns_offset; i<n; ++i, ++j) {
964
964
        // The first i quadrature points are the same as quad1, and
965
965
        // are by default singular.
966
966
        this->quadrature_points[i] = quad1.point(i)*fraction;
967
967
        this->weights[i] = quad1.weight(i)*fraction;
968
 
        
 
968
 
969
969
        // We need to scale with -log|fraction*alpha|
970
970
        if( (alpha != 1) || (fraction != 1) ) {
971
971
            this->quadrature_points[j] = quad.point(i)*fraction;
972
972
            this->weights[j] = -std::log(alpha/fraction)*quad.weight(i)*fraction;
973
973
        }
974
 
 
975
974
        // In case we need the second quadrature as well, do it now.
976
975
        if(fraction != 1) {
977
976
            this->quadrature_points[i+n] = quad2.point(i)*(1-fraction)+Point<1>(fraction);
978
 
            this->weights[i+n] = quad2.weight(i)*(1-fraction); 
979
 
            
 
977
            this->weights[i+n] = quad2.weight(i)*(1-fraction);
 
978
 
980
979
            // We need to scale with -log|fraction*alpha|
981
980
            this->quadrature_points[j+n] = quad.point(i)*(1-fraction)+Point<1>(fraction);
982
981
            this->weights[j+n] = -std::log(alpha/(1-fraction))*quad.weight(i)*(1-fraction);
983
982
        }
984
983
    }
985
 
    if(factor_out_singularity == true) 
 
984
    if(factor_out_singularity == true)
986
985
        for(unsigned int i=0; i<size(); ++i) {
987
986
            Assert( this->quadrature_points[i] != origin,
988
987
                    ExcMessage("The singularity cannot be on a Gauss point of the same order!") );
993
992
 
994
993
 
995
994
template<>
996
 
QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n, 
 
995
QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n,
997
996
                                  const unsigned int vertex_index,
998
997
                                  const bool factor_out_singularity) :
999
998
    Quadrature<2>(2*n*n)
1001
1000
    // Start with the gauss quadrature formula on the (u,v) reference
1002
1001
    // element.
1003
1002
    QGauss<2> gauss(n);
1004
 
    
 
1003
 
1005
1004
    Assert(gauss.size() == n*n, ExcInternalError());
1006
 
    
 
1005
 
1007
1006
    // For the moment we only implemented this for the vertices of a
1008
1007
    // quadrilateral. We are planning to do this also for the support
1009
1008
    // points of arbitrary FE_Q elements, to allow the use of this
1010
1009
    // class in boundary element programs with higher order mappings.
1011
1010
    Assert(vertex_index < 4, ExcIndexRange(vertex_index, 0, 4));
1012
 
    
 
1011
 
1013
1012
    // We create only the first one. All other pieces are rotation of
1014
1013
    // this one.
1015
1014
    // In this case the transformation is
1016
 
    // 
 
1015
    //
1017
1016
    // (x,y) = (u, u tan(pi/4 v))
1018
1017
    //
1019
1018
    // with Jacobian
1025
1024
    std::vector<Point<2> >      &ps = this->quadrature_points;
1026
1025
    std::vector<double>         &ws = this->weights;
1027
1026
    double pi4 = numbers::PI/4;
1028
 
    
 
1027
 
1029
1028
    for(unsigned int q=0; q<gauss.size(); ++q) {
1030
1029
        const Point<2> &gp = gauss.point(q);
1031
1030
        ps[q][0] = gp[0];
1039
1038
        ps[gauss.size()+q][0] = ps[q][1];
1040
1039
        ps[gauss.size()+q][1] = ps[q][0];
1041
1040
    }
1042
 
    
 
1041
 
1043
1042
    // Now we distribute these vertices in the correct manner
1044
1043
    double theta = 0;
1045
1044
    switch(vertex_index) {
1047
1046
        theta = 0;
1048
1047
        break;
1049
1048
    case 1:
1050
 
        // 
 
1049
        //
1051
1050
        theta = numbers::PI/2;
1052
1051
        break;
1053
1052
    case 2:
1057
1056
        theta = numbers::PI;
1058
1057
        break;
1059
1058
    }
1060
 
    
 
1059
 
1061
1060
    double R00 =  std::cos(theta), R01 = -std::sin(theta);
1062
1061
    double R10 =  std::sin(theta), R11 =  std::cos(theta);
1063
 
    
 
1062
 
1064
1063
    if(vertex_index != 0)
1065
1064
        for(unsigned int q=0; q<size(); ++q) {
1066
1065
            double x = ps[q][0]-.5,  y = ps[q][1]-.5;
1067
 
            
 
1066
 
1068
1067
            ps[q][0] = R00*x + R01*y + .5;
1069
1068
            ps[q][1] = R10*x + R11*y + .5;
1070
1069
        }