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-2022 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
47namespace 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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...