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