1
// Copyright (C) 2003-2006 Anders Logg.
2
// Licensed under the GNU LGPL Version 2.1.
4
// First added: 2003-06-03
5
// Last changed: 2006-10-23
8
#include <dolfin/dolfin_log.h>
9
#include <dolfin/Legendre.h>
10
#include <dolfin/RadauQuadrature.h>
12
using namespace dolfin;
14
//-----------------------------------------------------------------------------
15
RadauQuadrature::RadauQuadrature(unsigned int n) : GaussianQuadrature(n)
20
error("Radau quadrature not ok, check failed.");
22
//message("Radau quadrature computed for n = %d, check passed.", n);
24
//-----------------------------------------------------------------------------
25
void RadauQuadrature::disp() const
27
cout << "Radau quadrature points and weights on [-1,1] for n = "
30
cout << " i points weights" << endl;
31
cout << "-----------------------------------------------------" << endl;
33
for (unsigned int i = 0; i < n; i++)
34
message("%2d %.16e %.16e", i, points[i], weights[i]);
36
//-----------------------------------------------------------------------------
37
void RadauQuadrature::computePoints()
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.
52
Legendre p1(n-1), p2(n);
53
real x, dx, step, sign;
55
// Set size of stepping for seeking starting points
56
step = 2.0 / ( real(n-1) * 10.0 );
58
// Set the first nodal point which is -1
64
// Set the sign at -1 + epsilon
65
sign = ( (p1(x) + p2(x)) > 0 ? 1.0 : -1.0 );
67
// Compute the rest of the nodes by Newton's method
68
for (unsigned int i = 1; i < n; i++) {
70
// Step to a sign change
71
while ( (p1(x) + p2(x))*sign > 0.0 )
77
dx = - (p1(x) + p2(x)) / (p1.ddx(x) + p2.ddx(x));
79
} while ( fabs(dx) > DOLFIN_EPS );
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;
95
//-----------------------------------------------------------------------------