Functions | Variables
oomph::Orthpoly Namespace Reference

Functions

void gll_nodes (const unsigned &Nnode, Vector< double > &x)
 Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1. More...
 
void gll_nodes (const unsigned &Nnode, Vector< double > &x, Vector< double > &w)
 
void gl_nodes (const unsigned &Nnode, Vector< double > &x)
 
void gl_nodes (const unsigned &Nnode, Vector< double > &x, Vector< double > &w)
 
double legendre (const unsigned &p, const double &x)
 
void legendre_vector (const unsigned &p, const double &x, Vector< double > &polys)
 
double dlegendre (const unsigned &p, const double &x)
 
double ddlegendre (const unsigned &p, const double &x)
 
double jacobi (const int &alpha, const int &beta, const unsigned &p, const double &x)
 Calculate the Jacobi polnomials. More...
 
void jacobi (const int &alpha, const int &beta, const unsigned &p, const double &x, Vector< double > &polys)
 Calculate the Jacobi polnomials all in one goe. More...
 

Variables

const double eps = 1.0e-15
 

Function Documentation

◆ ddlegendre()

double oomph::Orthpoly::ddlegendre ( const unsigned &  p,
const double &  x 
)
inline

Calculates second derivative of Legendre polynomial of degree p at x using three term recursive formula.

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  }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
p
float * p
Definition: Tutorial_Map_using.cpp:9
plotDoE.x
list x
Definition: plotDoE.py:28

References i, p, and plotDoE::x.

Referenced by gll_nodes(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::oc_hang_helper(), oomph::OneDimensionalLegendreDShape< NNODE_1D >::OneDimensionalLegendreDShape(), oomph::OneDLegendreDShapeParam::OneDLegendreDShapeParam(), and oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::quad_hang_helper().

◆ dlegendre()

double oomph::Orthpoly::dlegendre ( const unsigned &  p,
const double &  x 
)
inline

Calculates first derivative of Legendre polynomial of degree p at x using three term recursive formula. \( nP_{n+1}^{'} = (2n+1)xP_{n}^{'} - (n+1)P_{n-1}^{'} \)

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  }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11

References n, p, and plotDoE::x.

Referenced by gl_nodes(), gll_nodes(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::oc_hang_helper(), oomph::OneDimensionalLegendreDShape< NNODE_1D >::OneDimensionalLegendreDShape(), oomph::OneDimensionalLegendreShape< NNODE_1D >::OneDimensionalLegendreShape(), oomph::OneDimensionalModalDShape::OneDimensionalModalDShape(), oomph::OneDLegendreDShapeParam::OneDLegendreDShapeParam(), oomph::OneDLegendreShapeParam::OneDLegendreShapeParam(), and oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::quad_hang_helper().

◆ gl_nodes() [1/2]

void oomph::Orthpoly::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  }
cos
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
BiharmonicTestFunctions2::Pi
double Pi
Definition: two_d_biharmonic.cc:235
CRBond_Bessel::eps
double eps
Definition: crbond_bessel.cc:24
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
oomph::Orthpoly::dlegendre
double dlegendre(const unsigned &p, const double &x)
Definition: orthpoly.h:121
oomph::Orthpoly::legendre
double legendre(const unsigned &p, const double &x)
Definition: orthpoly.h:57
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References cos(), dlegendre(), eps, boost::multiprecision::fabs(), j, legendre(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, p, oomph::MathematicalConstants::Pi, sqrt(), and plotDoE::x.

Referenced by oomph::GaussLegendre< 1, NPTS_1D >::GaussLegendre(), oomph::GaussLegendre< 2, NPTS_1D >::GaussLegendre(), oomph::GaussLegendre< 3, NPTS_1D >::GaussLegendre(), and gl_nodes().

◆ gl_nodes() [2/2]

void oomph::Orthpoly::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  }
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
oomph::Orthpoly::gl_nodes
void gl_nodes(const unsigned &Nnode, Vector< double > &x, Vector< double > &w)
Definition: orthpoly.cc:148

References dlegendre(), gl_nodes(), i, p, w, and plotDoE::x.

◆ gll_nodes() [1/2]

void oomph::Orthpoly::gll_nodes ( const unsigned &  Nnode,
Vector< double > &  x 
)

Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1.

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  }
oomph::Orthpoly::ddlegendre
double ddlegendre(const unsigned &p, const double &x)
Definition: orthpoly.h:144

References cos(), ddlegendre(), dlegendre(), eps, boost::multiprecision::fabs(), j, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, p, oomph::MathematicalConstants::Pi, sqrt(), and plotDoE::x.

Referenced by oomph::OneDimensionalLegendreShape< NNODE_1D >::calculate_nodal_positions(), oomph::OneDLegendreShapeParam::calculate_nodal_positions(), oomph::GaussLobattoLegendre< 1, NPTS_1D >::GaussLobattoLegendre(), oomph::GaussLobattoLegendre< 2, NPTS_1D >::GaussLobattoLegendre(), oomph::GaussLobattoLegendre< 3, NPTS_1D >::GaussLobattoLegendre(), oomph::PRefineableQElement< 1, INITIAL_NNODE_1D >::get_node_at_local_coordinate(), oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::get_node_at_local_coordinate(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::get_node_at_local_coordinate(), gll_nodes(), oomph::PRefineableQElement< 3, INITIAL_NNODE_1D >::oc_hang_helper(), and oomph::PRefineableQElement< 2, INITIAL_NNODE_1D >::quad_hang_helper().

◆ gll_nodes() [2/2]

void oomph::Orthpoly::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  }
oomph::Orthpoly::gll_nodes
void gll_nodes(const unsigned &Nnode, Vector< double > &x, Vector< double > &w)
Definition: orthpoly.cc:90

References gll_nodes(), i, legendre(), p, w, and plotDoE::x.

◆ jacobi() [1/2]

double oomph::Orthpoly::jacobi ( const int &  alpha,
const int &  beta,
const unsigned &  p,
const double &  x 
)
inline

Calculate the Jacobi polnomials.

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  }
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:151
beta
Scalar beta
Definition: level2_cplx_impl.h:36
Problem_Parameter::P0
double P0
Definition: two_dim.cc:101

References alpha, beta, n, p, Problem_Parameter::P0, and plotDoE::x.

◆ jacobi() [2/2]

void oomph::Orthpoly::jacobi ( const int &  alpha,
const int &  beta,
const unsigned &  p,
const double &  x,
Vector< double > &  polys 
)
inline

Calculate the Jacobi polnomials all in one goe.

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  }

References alpha, beta, n, p, Problem_Parameter::P0, and plotDoE::x.

◆ legendre()

double oomph::Orthpoly::legendre ( const unsigned &  p,
const double &  x 
)
inline

Calculates Legendre polynomial of degree p at x using the three term recurrence relation \( (n+1) P_{n+1} = (2n+1)xP_{n} - nP_{n-1} \)

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  }

References n, p, and plotDoE::x.

Referenced by gl_nodes(), gll_nodes(), oomph::OneDimensionalLegendreDShape< NNODE_1D >::OneDimensionalLegendreDShape(), oomph::OneDimensionalLegendreShape< NNODE_1D >::OneDimensionalLegendreShape(), oomph::OneDimensionalModalDShape::OneDimensionalModalDShape(), oomph::OneDimensionalModalShape::OneDimensionalModalShape(), oomph::OneDLegendreDShapeParam::OneDLegendreDShapeParam(), and oomph::OneDLegendreShapeParam::OneDLegendreShapeParam().

◆ legendre_vector()

void oomph::Orthpoly::legendre_vector ( const unsigned &  p,
const double &  x,
Vector< double > &  polys 
)
inline

Calculates Legendre polynomial of degree p at x using three term recursive formula. Returns all polynomials up to order p in the vector

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  }

References n, p, and plotDoE::x.

Variable Documentation

◆ eps

const double oomph::Orthpoly::eps = 1.0e-15

Referenced by gl_nodes(), gll_nodes(), oomph::OneDimensionalLegendreDShape< NNODE_1D >::OneDimensionalLegendreDShape(), oomph::OneDimensionalLegendreShape< NNODE_1D >::OneDimensionalLegendreShape(), oomph::OneDLegendreDShapeParam::OneDLegendreDShapeParam(), and oomph::OneDLegendreShapeParam::OneDLegendreShapeParam().