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-2023 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 #include "orthpoly.h"
27 
28 namespace 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
//////////////////////////////////////////////////////////////////// ////////////////////////////////...