← Back to branch summary

~ubuntu-branches/ubuntu/vivid/ffc/vivid

« back to all changes in this revision

Viewing changes to ffc/ext/Legendre.cpp

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2014-01-10 13:56:45 UTC
  • mfrom: (1.1.14)
  • Revision ID: package-import@ubuntu.com-20140110135645-4ozcd71y1oggj44z
Tags: 1.3.0-1
* New upstream release.
* debian/watch: Update URL for move to Bitbucket.
* debian/docs: README -> README.rst and remove TODO.
* debian/control:
  - Add python-numpy to Build-Depends.
  - Replace python-all with python-all-dev in Build-Depends.
  - Add ${shlibs:Depends} to Depends.
  - Change to Architecture: any.
  - Bump Standards-Version to 3.9.5 (no changes needed).
* debian/rules: Call dh_numpy in override_dh_python2.

Show diffs side-by-side

added added

removed removed

Lines of Context:
ffc/ext/Legendre.cpp
 
1
// Copyright (C) 2003-2008 Anders Logg
 
2
//
 
3
// This file is part of FFC.
 
4
//
 
5
// FFC is free software: you can redistribute it and/or modify
 
6
// it under the terms of the GNU Lesser General Public License as published by
 
7
// the Free Software Foundation, either version 3 of the License, or
 
8
// (at your option) any later version.
 
9
//
 
10
// FFC is distributed in the hope that it will be useful,
 
11
// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
// GNU Lesser General Public License for more details.
 
14
//
 
15
// You should have received a copy of the GNU Lesser General Public License
 
16
// along with FFC.  If not, see <http://www.gnu.org/licenses/>.
 
17
//
 
18
// Modified by Benjamin Kehlet 2011-2012
 
19
//
 
20
// First added:  2012-08-20
 
21
// Last changed: 2012-09-05
 
22
 
 
23
 
 
24
#include "Legendre.h"
 
25
#include <cmath>
 
26
 
 
27
//-----------------------------------------------------------------------------
 
28
Legendre::Legendre(unsigned int n) : n(n), cache_x(0.0), cache(n + 1)
 
29
{
 
30
  cache[0] = 1.0; //constant value
 
31
 
 
32
  // eval to initialize cache
 
33
  eval(n, -1.0);
 
34
}
 
35
//-----------------------------------------------------------------------------
 
36
double Legendre::operator() (double x)
 
37
{
 
38
  return eval(n, x);
 
39
}
 
40
//-----------------------------------------------------------------------------
 
41
double Legendre::ddx(double x)
 
42
{
 
43
  return ddx(n, x);
 
44
}
 
45
//-----------------------------------------------------------------------------
 
46
double Legendre::d2dx(double x)
 
47
{
 
48
  return d2dx(n, x);
 
49
}
 
50
//-----------------------------------------------------------------------------
 
51
double Legendre::eval(unsigned int nn, double x)
 
52
{
 
53
  //recursive formula, BETA page 254
 
54
  //return ( (2.0*nn-1.0)*x*eval(nn-1, x) - (nn-1.0)*eval(nn-2, x) ) / nn;
 
55
 
 
56
  //The special cases
 
57
  if (n == 0)
 
58
    return 1.0;
 
59
  else if (n == 1)
 
60
    return x;
 
61
 
 
62
  //check cache
 
63
  if (x != cache_x)
 
64
  {
 
65
    cache[1] = x;
 
66
    for (unsigned int i = 2; i <= n; ++i)
 
67
    {
 
68
      double ii = i;
 
69
      cache[i] = ( (2.0*ii-1.0)*x*cache[i-1] - (ii-1.0)*cache[i-2] ) / ii;
 
70
    }
 
71
    cache_x = x;
 
72
  }
 
73
 
 
74
  return cache[nn];
 
75
}
 
76
//-----------------------------------------------------------------------------
 
77
double Legendre::ddx(unsigned int n, double x)
 
78
{
 
79
  // Special cases
 
80
  if (n == 0)
 
81
    return 0.0;
 
82
  else if (n == 1)
 
83
    return 1.0;
 
84
 
 
85
  // Avoid division by zero
 
86
  if (std::abs(x - 1.0) < EPSILON)
 
87
    x -= 2.0*EPSILON;
 
88
 
 
89
  if (std::abs(x + 1.0) < EPSILON)
 
90
    x += 2.0*EPSILON;
 
91
 
 
92
  // Formula, BETA page 254
 
93
  const double nn = n;
 
94
  return nn * (x*eval(n, x) - eval(n-1, x)) / (x*x - 1.0);
 
95
}
 
96
//-----------------------------------------------------------------------------
 
97
double Legendre::d2dx(unsigned int, double x)
 
98
{
 
99
  // Special case n = 0
 
100
  if (n == 0)
 
101
    return 0.0;
 
102
 
 
103
  // Special case n = 1
 
104
  if (n == 1)
 
105
    return 0.0;
 
106
 
 
107
  // Avoid division by zero
 
108
  if (std::abs(x - 1.0) < EPSILON)
 
109
    x -= 2.0*EPSILON;
 
110
  if (std::abs(x + 1.0) < EPSILON)
 
111
    x += 2.0*EPSILON;
 
112
 
 
113
  // Formula, BETA page 254
 
114
  const double nn = double(n);
 
115
  return (2.0*x*ddx(n, x) - nn*(nn+1)*eval(n, x)) / (1.0-x*x);
 
116
}
 
117
//-----------------------------------------------------------------------------