orthpoly.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Header functions for functions used to generate orthogonal polynomials
31 //Include guards to prevent multiple inclusion of the header
32 
33 #ifndef OOMPH_ORTHPOLY_HEADER
34 #define OOMPH_ORTHPOLY_HEADER
35 
36 // Config header generated by autoconfig
37 #ifdef HAVE_CONFIG_H
38  #include <oomph-lib-config.h>
39 #endif
40 
41 #ifdef OOMPH_HAS_MPI
42 #include "mpi.h"
43 #endif
44 
45 //Oomph-lib include
46 #include "Vector.h"
47 #include "oomph_utilities.h"
48 
49 #include <cmath>
50 
51 namespace oomph
52 {
53 
54 //Let's put these things in a namespace
55 namespace Orthpoly
56 {
57  const double eps = 1.0e-15;
58 
59  ///\short Calculates Legendre polynomial of degree p at x
60  /// using the three term recurrence relation
61  /// \f$ (n+1) P_{n+1} = (2n+1)xP_{n} - nP_{n-1} \f$
62  inline double legendre(const unsigned &p, const double &x)
63  {
64  //Return the constant value
65  if(p==0) return 1.0;
66  //Return the linear polynomial
67  else if(p==1) return x;
68  //Otherwise use the recurrence relation
69  else{
70  //Initialise the terms in the recurrence relation, we're going
71  //to shift down before using the relation.
72  double L0 = 1.0, L1 = 1.0, L2 = x;
73  //Loop over the remaining polynomials
74  for(unsigned n=1;n<p;n++)
75  {
76  //Shift down the values
77  L0 = L1; L1 = L2;
78  //Use the three term recurrence
79  L2 = ((2*n + 1)*x*L1 - n*L0)/(n+1);
80  }
81  //Once we've finished return the final value
82  return L2;
83  }
84  }
85 
86 
87  /// \short Calculates Legendre polynomial of degree p at x
88  /// using three term recursive formula. Returns all polynomials up to
89  /// order p in the vector
90  inline void legendre_vector(const unsigned &p, const double &x,
91  Vector<double> &polys)
92  {
93  //Set the constant term
94  polys[0] = 1.0;
95  //If we're only asked for the constant term return
96  if(p==0) {return;}
97  //Set the linear polynomial
98  polys[1] = x;
99  //Initialise terms for the recurrence relation
100  double L0 = 1.0, L1 = 1.0, L2 = x;
101  //Loop over the remaining terms
102  for(unsigned n=1;n<p;n++){
103  //Shift down the values
104  L0=L1;
105  L1=L2;
106  //Use the recurrence relation
107  L2 = ((2*n+1)*x*L1 - n*L0)/(n+1);
108  //Set the newly calculated polynomial
109  polys[n+1] = L2;
110  }
111  }
112 
113 
114  /// \short Calculates first derivative of Legendre
115  /// polynomial of degree p at x
116  /// using three term recursive formula.
117  /// \f$ nP_{n+1}^{'} = (2n+1)xP_{n}^{'} - (n+1)P_{n-1}^{'} \f$
118  inline double dlegendre(const unsigned &p, const double &x)
119  {
120  double dL1 = 1.0, dL2 = 3*x, dL3=0.0;
121  if(p==0) return 0.0;
122  else if(p==1) return dL1;
123  else if(p==2) return dL2;
124  else{
125  for(unsigned n=2;n<p;n++){
126  dL3 = 1.0/n*((2.0*n+1.0)*x*dL2-(n+1.0)*dL1);
127  dL1=dL2;
128  dL2=dL3;
129  }
130  return dL3;
131  }
132  }
133 
134  /// \short Calculates second derivative of Legendre
135  /// polynomial of degree p at x
136  /// using three term recursive formula.
137  inline double ddlegendre(const unsigned &p, const double &x)
138  {
139  double ddL2 = 3.0, ddL3 = 15*x, ddL4=0.0;
140  if(p==0) return 0.0;
141  else if(p==1) return 0.0;
142  else if(p==2) return ddL2;
143  else if(p==3) return ddL3;
144  else{
145  for(unsigned i=3;i<p;i++){
146  ddL4 = 1.0/(i-1.0)*((2.0*i+1.0)*x*ddL3-(i+2.0)*ddL2);
147  ddL2=ddL3;
148  ddL3=ddL4;
149  }
150  return ddL4;
151  }
152  }
153 
154  /// \short Calculate the Jacobi polnomials
155  inline double jacobi(const int &alpha, const int &beta, const unsigned &p,
156  const double &x)
157  {
158  double P0 = 1.0;
159  double P1 = 0.5*(alpha-beta+(alpha+beta+2.0)*x);
160  double P2;
161  if(p==0) return P0;
162  else if(p==1) return P1;
163  else{
164  for(unsigned n=1;n<p;n++){
165  double a1n = 2*(n+1)*(n+alpha+beta+1)*(2*n+alpha+beta);
166  double a2n = (2*n+alpha+beta+1)*(alpha*alpha-beta*beta);
167  double a3n = (2*n+alpha+beta)*(2*n+alpha+beta+1)*(2*n+alpha+beta+2);
168  double a4n = 2*(n+alpha)*(n+beta)*(2*n+alpha+beta+2);
169  P2 = ( (a2n+a3n*x)*P1 - a4n*P0 ) / a1n;
170  P0 = P1;
171  P1 = P2;
172  }
173  return P2;
174  }
175  }
176 
177  /// \short Calculate the Jacobi polnomials all in one goe
178  inline void jacobi(const int &alpha, const int &beta, const unsigned &p,
179  const double &x,
180  Vector<double> &polys)
181  {
182  //Set the constant term
183  polys[0] = 1.0;
184  //If we've only been asked for the constant term, bin out
185  if(p==0) {return;}
186  //Set the linear polynomial
187  polys[1] = 0.5*(alpha-beta+(alpha+beta+2.0)*x);
188  //Initialise the terms for the recurrence relation
189  double P0 = 1.0;
190  double P1 = 1.0;
191  double P2 = 0.5*(alpha-beta+(alpha+beta+2.0)*x);
192  //Loop over the remaining terms
193  for(unsigned n=1;n<p;n++)
194  {
195  //Shift down the terms
196  P0 = P1;
197  P1 = P2;
198  //Set the constants
199  double a1n = 2*(n+1)*(n+alpha+beta+1)*(2*n+alpha+beta);
200  double a2n = (2*n+alpha+beta+1)*(alpha*alpha-beta*beta);
201  double a3n = (2*n+alpha+beta)*(2*n+alpha+beta+1)*(2*n+alpha+beta+2);
202  double a4n = 2*(n+alpha)*(n+beta)*(2*n+alpha+beta+2);
203  //Set the latest term
204  P2 = ( (a2n+a3n*x)*P1 - a4n*P0 ) / a1n;
205  //Set the newly calculate polynomial
206  polys[n+1] = P2;
207  }
208  }
209 
210 
211 /// Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1
212  void gll_nodes(const unsigned &Nnode, Vector<double> &x);
213 
214 // This version of gll_nodes calculates the abscissas AND weights
215  void gll_nodes(const unsigned &Nnode, Vector<double> &x, Vector<double> &w);
216 
217 // Calculates the Gauss Legendre abscissas of degree p=Nnode-1
218  void gl_nodes(const unsigned &Nnode, Vector<double> &x);
219 
220 // This version of gl_nodes calculates the abscissas AND weights
221  void gl_nodes(const unsigned &Nnode, Vector<double> &x, Vector<double> &w);
222 
223 } //End of the namespace
224 
225 }
226 
227 #endif
double dlegendre(const unsigned &p, const double &x)
Calculates first derivative of Legendre polynomial of degree p at x using three term recursive formul...
Definition: orthpoly.h:118
cstr elem_len * i
Definition: cfortran.h:607
double legendre(const unsigned &p, const double &x)
Calculates Legendre polynomial of degree p at x using the three term recurrence relation ...
Definition: orthpoly.h:62
void gll_nodes(const unsigned &Nnode, Vector< double > &x)
Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1.
Definition: orthpoly.cc:38
double ddlegendre(const unsigned &p, const double &x)
Calculates second derivative of Legendre polynomial of degree p at x using three term recursive formu...
Definition: orthpoly.h:137
void gl_nodes(const unsigned &Nnode, Vector< double > &x)
Definition: orthpoly.cc:96
void legendre_vector(const unsigned &p, const double &x, Vector< double > &polys)
Calculates Legendre polynomial of degree p at x using three term recursive formula. Returns all polynomials up to order p in the vector.
Definition: orthpoly.h:90
double jacobi(const int &alpha, const int &beta, const unsigned &p, const double &x)
Calculate the Jacobi polnomials.
Definition: orthpoly.h:155
const double eps
Definition: orthpoly.h:57