← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/math/Legendre.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/math/Legendre.cpp
1
 
// Copyright (C) 2003-2005 Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// First added:  2003-06-03
5
 
// Last changed: 2005
6
 
 
7
 
#include <cmath>
8
 
#include <dolfin/dolfin_log.h>
9
 
#include <dolfin/Legendre.h>
10
 
 
11
 
using namespace dolfin;
12
 
 
13
 
//-----------------------------------------------------------------------------
14
 
Legendre::Legendre(int n)
15
 
{
16
 
  if ( n < 0 )
17
 
    error("Degree for Legendre polynomial must be non-negative.");
18
 
 
19
 
  this->n = n;
20
 
}
21
 
//-----------------------------------------------------------------------------
22
 
real Legendre::operator() (real x)
23
 
{
24
 
  return eval(n, x);
25
 
}
26
 
//-----------------------------------------------------------------------------
27
 
real Legendre::ddx(real x)
28
 
{
29
 
  return ddx(n, x);
30
 
}
31
 
//-----------------------------------------------------------------------------
32
 
real Legendre::d2dx(real x)
33
 
{
34
 
  return d2dx(n, x);
35
 
}
36
 
//-----------------------------------------------------------------------------
37
 
real Legendre::eval(int n, real x)
38
 
{
39
 
  // Special case n = 0
40
 
  if ( n == 0 )
41
 
    return 1.0;
42
 
 
43
 
  // Special case n = 1
44
 
  if ( n == 1 )
45
 
    return x;
46
 
  
47
 
  // Recurrence, BETA page 254
48
 
  real nn = real(n);
49
 
  return ( (2.0*nn-1.0)*x*eval(n-1, x) - (nn-1.0)*eval(n-2, x) ) / nn;
50
 
}
51
 
//-----------------------------------------------------------------------------
52
 
real Legendre::ddx(int n, real x)
53
 
{
54
 
  // Special case n = 0
55
 
  if ( n == 0 )
56
 
    return 0.0;
57
 
  
58
 
  // Special case n = 1
59
 
  if ( n == 1 )
60
 
    return 1.0;
61
 
  
62
 
  // Avoid division by zero
63
 
  if ( fabs(x - 1.0) < DOLFIN_EPS )
64
 
    x -= 2.0*DOLFIN_EPS;
65
 
  if ( fabs(x + 1.0) < DOLFIN_EPS )
66
 
    x += 2.0*DOLFIN_EPS;
67
 
  
68
 
  // Formula, BETA page 254
69
 
  real nn = real(n);
70
 
  return nn * (x*eval(n, x) - eval(n-1, x)) / (x*x - 1.0);
71
 
}
72
 
//-----------------------------------------------------------------------------
73
 
real Legendre::d2dx(int, real x)
74
 
{
75
 
  // Special case n = 0
76
 
  if ( n == 0 )
77
 
    return 0.0;
78
 
 
79
 
  // Special case n = 1
80
 
  if ( n == 1 )
81
 
    return 0.0;
82
 
 
83
 
  // Avoid division by zero
84
 
  if ( fabs(x - 1.0) < DOLFIN_EPS )
85
 
    x -= 2.0*DOLFIN_EPS;
86
 
  if ( fabs(x + 1.0) < DOLFIN_EPS )
87
 
    x += 2.0*DOLFIN_EPS;
88
 
 
89
 
  // Formula, BETA page 254
90
 
  real nn = real(n);
91
 
  return ( 2.0*x*ddx(n, x) - nn*(nn+1)*eval(n, x) ) / (1.0-x*x);
92
 
}
93
 
//-----------------------------------------------------------------------------