← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/quadrature/RadauQuadrature.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/RadauQuadrature.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 <cmath>
8
 
#include <dolfin/dolfin_log.h>
9
 
#include <dolfin/Legendre.h>
10
 
#include <dolfin/RadauQuadrature.h>
11
 
 
12
 
using namespace dolfin;
13
 
 
14
 
//-----------------------------------------------------------------------------
15
 
RadauQuadrature::RadauQuadrature(unsigned int n) : GaussianQuadrature(n)
16
 
{
17
 
  init();
18
 
 
19
 
  if ( !check(2*n-2) )
20
 
    error("Radau quadrature not ok, check failed.");
21
 
 
22
 
  //message("Radau quadrature computed for n = %d, check passed.", n);
23
 
}
24
 
//-----------------------------------------------------------------------------
25
 
void RadauQuadrature::disp() const
26
 
{
27
 
  cout << "Radau quadrature points and weights on [-1,1] for n = " 
28
 
       << n << ":" << endl;
29
 
 
30
 
  cout << " i    points                   weights" << endl;
31
 
  cout << "-----------------------------------------------------" << endl;
32
 
 
33
 
  for (unsigned int i = 0; i < n; i++)
34
 
    message("%2d   %.16e   %.16e", i, points[i], weights[i]);
35
 
}
36
 
//-----------------------------------------------------------------------------
37
 
void RadauQuadrature::computePoints()
38
 
{
39
 
  // Compute the Radau quadrature points in [-1,1] as -1 and the zeros
40
 
  // of ( Pn-1(x) + Pn(x) ) / (1+x) where Pn is the n:th Legendre
41
 
  // polynomial. Computation is a little different than for Gauss and
42
 
  // Lobatto quadrature, since we don't know of any good initial
43
 
  // approximation for the Newton iterations.
44
 
  
45
 
  // Special case n = 1
46
 
  if ( n == 1 )
47
 
  {
48
 
    points[0] = -1.0;
49
 
    return;
50
 
  }
51
 
 
52
 
  Legendre p1(n-1), p2(n);
53
 
  real x, dx, step, sign;
54
 
  
55
 
  // Set size of stepping for seeking starting points
56
 
  step = 2.0 / ( real(n-1) * 10.0 );
57
 
  
58
 
  // Set the first nodal point which is -1
59
 
  points[0] = -1.0;
60
 
  
61
 
  // Start at -1 + step
62
 
  x = -1.0 + step;
63
 
  
64
 
  // Set the sign at -1 + epsilon
65
 
  sign = ( (p1(x) + p2(x)) > 0 ? 1.0 : -1.0 );
66
 
  
67
 
  // Compute the rest of the nodes by Newton's method
68
 
  for (unsigned int i = 1; i < n; i++) {
69
 
    
70
 
    // Step to a sign change
71
 
    while ( (p1(x) + p2(x))*sign > 0.0 )
72
 
      x += step;
73
 
    
74
 
    // Newton's method
75
 
    do
76
 
    {
77
 
      dx = - (p1(x) + p2(x)) / (p1.ddx(x) + p2.ddx(x));
78
 
      x  = x + dx;
79
 
    } while ( fabs(dx) > DOLFIN_EPS );
80
 
    
81
 
    // Set the node value
82
 
    points[i] = x;
83
 
    
84
 
    // Fix step so that it's not too large
85
 
    if ( step > (points[i] - points[i-1])/10.0 )
86
 
      step = (points[i] - points[i-1]) / 10.0;
87
 
    
88
 
    // Step forward
89
 
    sign = - sign;
90
 
    x += step;
91
 
    
92
 
  }
93
 
  
94
 
}
95
 
//-----------------------------------------------------------------------------