← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/quadrature/GaussQuadrature.cpp

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/kernel/quadrature/GaussQuadrature.cpp
1
 
// Copyright (C) 2003-2006 Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// First added:  2003-06-03
5
 
// Last changed: 2006-10-23
6
 
 
7
 
#include <stdio.h>
8
 
#include <cmath>
9
 
#include <dolfin/dolfin_log.h>
10
 
#include <dolfin/Legendre.h>
11
 
#include <dolfin/GaussQuadrature.h>
12
 
 
13
 
using namespace dolfin;
14
 
 
15
 
//-----------------------------------------------------------------------------
16
 
GaussQuadrature::GaussQuadrature(unsigned int n) : GaussianQuadrature(n)
17
 
{
18
 
  init();
19
 
 
20
 
  if ( !check(2*n-1) )
21
 
    error("Gauss quadrature not ok, check failed.");
22
 
  
23
 
  //message("Gauss quadrature computed for n = %d, check passed.", n);
24
 
}
25
 
//-----------------------------------------------------------------------------
26
 
void GaussQuadrature::disp() const
27
 
{
28
 
  cout << "Gauss quadrature points and weights on [-1,1] for n = " 
29
 
       << n << ":" << endl;
30
 
 
31
 
  cout << " i    points                   weights" << endl;
32
 
  cout << "-----------------------------------------------------" << endl;
33
 
  
34
 
  for (unsigned int i = 0; i < n; i++)
35
 
    message("%2d   %.16e   %.16e", i, points[i], weights[i]);
36
 
}
37
 
//-----------------------------------------------------------------------------
38
 
void GaussQuadrature::computePoints()
39
 
{
40
 
  // Compute Gauss quadrature points on [-1,1] as the
41
 
  // as the zeroes of the Legendre polynomials using Newton's method
42
 
  
43
 
  // Special case n = 1
44
 
  if ( n == 1 )
45
 
  {
46
 
    points[0] = 0.0;
47
 
    return;
48
 
  }
49
 
 
50
 
  Legendre p(n);
51
 
  real x, dx;
52
 
  
53
 
  // Compute the points by Newton's method
54
 
  for (unsigned int i = 0; i <= ((n-1)/2); i++) {
55
 
    
56
 
    // Initial guess
57
 
    x = cos(DOLFIN_PI*(real(i+1)-0.25)/(real(n)+0.5));
58
 
    
59
 
    // Newton's method
60
 
    do {
61
 
      dx = - p(x) / p.ddx(x);
62
 
      x  = x + dx;
63
 
    } while ( fabs(dx) > DOLFIN_EPS );
64
 
    
65
 
    // Save the value using the symmetry of the points
66
 
    points[i] = - x;
67
 
    points[n-1-i] = x;
68
 
    
69
 
  }
70
 
  
71
 
  // Set middle node
72
 
  if ( (n % 2) != 0 )
73
 
    points[n/2] = 0.0;
74
 
}
75
 
//-----------------------------------------------------------------------------