orthpoly.cc
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
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#include "orthpoly.h"
27
28namespace oomph
29{
30 namespace Orthpoly
31 {
32 // Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1
33 void gll_nodes(const unsigned& Nnode, Vector<double>& x)
34 {
35 double z, zold, del;
36 unsigned p = Nnode - 1;
37 x.resize(Nnode);
38 if (Nnode < 2)
39 {
40 throw OomphLibError("Invalid number of nodes",
41 OOMPH_CURRENT_FUNCTION,
42 OOMPH_EXCEPTION_LOCATION);
43 }
44 else if (Nnode == 2)
45 {
46 x[0] = -1.0;
47 x[1] = 1.0;
48 }
49 else if (Nnode == 3)
50 {
51 x[0] = -1;
52 x[1] = 0.0;
53 x[2] = 1.0;
54 }
55 else if (Nnode == 4)
56 {
57 x[0] = -1.0;
58 x[1] = -std::sqrt(5.0) / 5.0;
59 x[2] = -x[1];
60 x[3] = 1.0;
61 }
62 else
63 {
64 unsigned mid;
65 if (p % 2 == 0)
66 {
67 mid = p / 2;
68 x[mid] = 0.0;
69 }
70 else
71 mid = p / 2 + 1;
72 for (unsigned j = 1; j < mid; j++)
73 {
74 z = -std::cos(j * MathematicalConstants::Pi / double(p));
75 do
76 {
77 del = dlegendre(p, z) / ddlegendre(p, z);
78 zold = z;
79 z = zold - del;
80 } while (std::fabs(z - zold) > eps);
81 x[j] = z;
82 x[p - j] = -z;
83 }
84 x[0] = -1.0;
85 x[p] = 1.0;
86 }
87 }
88
89 // This version of gll_nodes calculates the abscissas AND weights
90 void gll_nodes(const unsigned& Nnode, Vector<double>& x, Vector<double>& w)
91 {
92 gll_nodes(Nnode, x);
93 // Now calculate the corresponding weights
94 double l_z;
95 unsigned p = Nnode - 1;
96 for (unsigned i = 0; i < p + 1; i++)
97 {
98 l_z = legendre(p, x[i]);
99 w[i] = 2.0 / (p * (p + 1) * l_z * l_z);
100 }
101 }
102
103
104 // Calculates the Gauss Legendre abscissas of degree p=Nnode-1
105 void gl_nodes(const unsigned& Nnode, Vector<double>& x)
106 {
107 double z, zold, del;
108 unsigned p = Nnode - 1;
109 x.resize(Nnode);
110 if (Nnode < 2)
111 {
112 throw OomphLibError("Invalid number of nodes",
113 OOMPH_CURRENT_FUNCTION,
114 OOMPH_EXCEPTION_LOCATION);
115 }
116 else if (Nnode == 2)
117 {
118 x[0] = -1.0 / 3.0 * std::sqrt(3.0);
119 x[1] = -x[0];
120 }
121 else
122 {
123 unsigned mid;
124 if (p % 2 == 0)
125 {
126 mid = p / 2;
127 x[mid] = 0.0;
128 }
129 else
130 mid = p / 2 + 1;
131 for (unsigned j = 0; j < mid; j++)
132 {
133 z = -std::cos((2.0 * j + 1.0) * MathematicalConstants::Pi /
134 (2 * double(p) + 2.0));
135 do
136 {
137 del = legendre(p + 1, z) / dlegendre(p + 1, z);
138 zold = z;
139 z = zold - del;
140 } while (std::fabs(z - zold) > eps);
141 x[j] = z;
142 x[p - j] = -z;
143 }
144 }
145 }
146
147 // This version of gl_nodes calculates the abscissas AND weights
148 void gl_nodes(const unsigned& Nnode, Vector<double>& x, Vector<double>& w)
149 {
150 gl_nodes(Nnode, x);
151 // Now calculate the corresponding weights
152 double dl_z;
153 unsigned p = Nnode - 1;
154 for (unsigned i = 0; i < p + 1; i++)
155 {
156 dl_z = dlegendre(p + 1, x[i]);
157 w[i] = 2.0 / (1 - x[i] * x[i]) / (dl_z * dl_z);
158 }
159 }
160 } // namespace Orthpoly
161
162} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
An OomphLibError object which should be thrown when an run-time error is encountered....
const double Pi
50 digits from maple
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 gll_nodes(const unsigned &Nnode, Vector< double > &x)
Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1.
Definition: orthpoly.cc:33
//////////////////////////////////////////////////////////////////// ////////////////////////////////...