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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #include "orthpoly.h"
31 
32 namespace oomph
33 {
34 
35 namespace Orthpoly
36 {
37 // Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1
38 void gll_nodes(const unsigned &Nnode, Vector<double> &x)
39 {
40  double z,zold,del;
41  unsigned p=Nnode-1;
42  x.resize(Nnode);
43  if(Nnode<2)
44  {
45  throw OomphLibError("Invalid number of nodes",
46  OOMPH_CURRENT_FUNCTION,
47  OOMPH_EXCEPTION_LOCATION);
48  }
49  else if (Nnode==2){
50  x[0]=-1.0;x[1]=1.0;
51  }
52  else if (Nnode==3){
53  x[0]=-1;x[1]=0.0;x[2]=1.0;
54  }
55  else if (Nnode==4){
56  x[0]=-1.0;x[1]=-std::sqrt(5.0)/5.0;x[2]=-x[1];x[3]=1.0;
57  }
58  else
59  {
60  unsigned mid;
61  if(p%2==0){
62  mid = p/2;
63  x[mid]=0.0;
64  }
65  else mid=p/2+1;
66  for(unsigned j=1;j<mid;j++)
67  {
68  z=-std::cos(j*MathematicalConstants::Pi/double(p));
69  do{
70  del= dlegendre(p,z)/ddlegendre(p,z);
71  zold=z;
72  z=zold-del;
73  } while(std::fabs(z-zold) > eps);
74  x[j]=z;
75  x[p-j]=-z;
76  }
77  x[0]=-1.0;x[p]=1.0;
78  }
79 }
80 
81 // This version of gll_nodes calculates the abscissas AND weights
82 void gll_nodes(const unsigned &Nnode, Vector<double> &x, Vector<double> &w)
83 {
84  gll_nodes(Nnode,x);
85  // Now calculate the corresponding weights
86  double l_z;
87  unsigned p = Nnode-1;
88  for (unsigned i=0;i<p+1;i++){
89  l_z=legendre(p,x[i]);
90  w[i]=2.0/(p*(p+1)*l_z*l_z);
91  }
92 }
93 
94 
95 // Calculates the Gauss Legendre abscissas of degree p=Nnode-1
96 void gl_nodes(const unsigned &Nnode, Vector<double> &x)
97 {
98  double z,zold,del;
99  unsigned p=Nnode-1;
100  x.resize(Nnode);
101  if(Nnode<2)
102  {
103  throw OomphLibError("Invalid number of nodes",
104  OOMPH_CURRENT_FUNCTION,
105  OOMPH_EXCEPTION_LOCATION);
106  }
107  else if (Nnode==2){
108  x[0]=-1.0/3.0*std::sqrt(3.0);x[1]=-x[0];
109  }
110  else
111  {
112  unsigned mid;
113  if(p%2==0){
114  mid = p/2;
115  x[mid]=0.0;
116  }
117  else mid=p/2+1;
118  for(unsigned j=0;j<mid;j++)
119  {
120  z=-std::cos((2.0*j+1.0)*MathematicalConstants::Pi/(2*double(p)+2.0));
121  do{
122  del = legendre(p+1,z)/dlegendre(p+1,z);
123  zold=z;
124  z=zold-del;
125  } while(std::fabs(z-zold)>eps);
126  x[j]=z;
127  x[p-j]=-z;
128  }
129  }
130 }
131 
132 // This version of gl_nodes calculates the abscissas AND weights
133 void gl_nodes(const unsigned &Nnode, Vector<double> &x, Vector<double> &w)
134 {
135  gl_nodes(Nnode,x);
136  // Now calculate the corresponding weights
137  double dl_z;
138  unsigned p = Nnode-1;
139  for (unsigned i=0;i<p+1;i++){
140  dl_z=dlegendre(p+1,x[i]);
141  w[i]=2.0/(1-x[i]*x[i])/(dl_z*dl_z);
142  }
143 }
144 }
145 
146 }
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:118
cstr elem_len * i
Definition: cfortran.h:607
const double Pi
50 digits from maple
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:62
void gll_nodes(const unsigned &Nnode, Vector< double > &x)
Calculates the Gauss Lobatto Legendre abscissas for degree p = NNode-1.
Definition: orthpoly.cc:38
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:137
void gl_nodes(const unsigned &Nnode, Vector< double > &x)
Definition: orthpoly.cc:96
const double eps
Definition: orthpoly.h:57