G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsLobattoRule.hpp
Go to the documentation of this file.
1 
14 #pragma once
15 
16 #include <gsCore/gsBasis.h>
17 #include <gsIO/gsOptionList.h>
18 
19 namespace gismo
20 {
21 
22 template<class T> void
23 gsLobattoRule<T>::setNodes( gsVector<index_t> const & numNodes,
24  unsigned digits)
25 {
26  const short_t d = static_cast<short_t>(numNodes.rows());
27  static const T epsilon = std::pow(10.0, -REAL_DIG * 0.85);
28  // Get base rule nodes and weights
29  std::vector<gsVector<T> > nodes(d);
30  std::vector<gsVector<T> > weights(d);
31 
32  if (digits == 0)
33  {
34  for (short_t i = 0; i < d; ++i)
35  {
36  if (!lookupReference(numNodes[i], nodes[i], weights[i]))
37  computeReference(numNodes[i], nodes[i], weights[i],
38  0==REAL_DIG?2:REAL_DIG);
39  if (1!=numNodes[i])
40  nodes[i].last() -= epsilon; //interval may be half-open
41  }
42  }
43  else
44  {
45  for (short_t i = 0; i < d; ++i)
46  {
47  computeReference(numNodes[i], nodes[i], weights[i], digits);
48  if (1!=numNodes[i])
49  nodes[i].last() -= epsilon; //interval may be half-open
50  }
51  }
52 
53  //std::copy(nodes.begin(), nodes.end(), std::ostream_iterator<gsVector<T> >(gsInfo, "\n"));
54 
55  this->computeTensorProductRule(nodes, weights);
56 }
57 
58 template<class T> void
59 gsLobattoRule<T>::computeReference(index_t n, // Number of points
60  gsVector<T> & x, // Quadrature points
61  gsVector<T> & w, // Quadrature weights
62  unsigned digits) // Number of exact decimal digits
63 {
64  GISMO_ASSERT(digits!=0, "Digits cannot be 0");
65  // Allocate space for points and weights.
66  x.resize(n);
67  w.resize(n);
68 
69  index_t i, j;
70  T test;
71  const T tolerance = math::pow(T(0.1), static_cast<int>(digits));
72 
73  if ( n == 1 )
74  {
75  x[0] = -1.0;
76  w[0] = 2.0;
77  return;
78  }
79 
80  // Initial estimate ( Chebyshev-Gauss-Lobatto nodes)
81  for ( i = 0; i < n; i++ )
82  x[i] = math::cos ( (T)(EIGEN_PI) * (T)( i ) / (T)( n - 1 ) );
83 
84  gsVector<T> xold(n);
85  gsVector<T> p(n*n);
86 
87  do
88  {
89  for ( i = 0; i < n; i++ )
90  {
91  xold[i] = p[i+1*n] = x[i];
92  p[i+0*n] = 1.0;
93  }
94 
95  for ( j = 2; j <= n-1; j++ )
96  {
97  for ( i = 0; i < n; i++)
98  {
99  p[i+j*n] = ( (T) ( 2 * j - 1 ) * x[i] * p[i+(j-1)*n]
100  + (T) ( - j + 1 ) * p[i+(j-2)*n] )
101  / (T) ( j );
102  }
103  }
104 
105  for ( i = 0; i < n; i++ )
106  x[i] = xold[i] - ( x[i] * p[i+(n-1)*n] - p[i+(n-2)*n] )
107  / ( (T) ( n ) * p[i+(n-1)*n] );
108 
109  test = 0.0;
110  for ( i = 0; i < n; i++ )
111  test = math::max( test, math::abs ( x[i] - xold[i] ) );
112 
113  } while ( tolerance <= test );
114 
115  x.reverseInPlace();
116 
117  for ( i = 0; i < n; i++ )
118  w[i] = 2.0 / ( (T) ( ( n - 1 ) * n ) * math::pow ( p[i+(n-1)*n], 2 ) );
119 }
120 
121 
122 template<class T> bool
123 gsLobattoRule<T>::lookupReference(index_t n, // Number of points
124  gsVector<T> & x, // Quadrature points
125  gsVector<T> & w) // Quadrature weights
126 {
127  if ( REAL_DIG >= 28 )
128  {
129  // The generated points and weights are only accurate up to ~30 decimal digits (leaving some wiggle room inside the conditional).
130  // If the precision of the number format aliased by real_t is higher than that the lookup is refused.
131  // More precise weights must be computed on-the-fly.
132  return false;
133  }
134 
135  x.resize(n);
136  w.resize(n);
137 
138  switch (n)
139  {
140  case 1 :
141  {
142  x << -1.0;
143  w << 2.0;
144  return true;
145  }
146  case 2 :
147  {
148  x << - 1.0E+00, 1.0E+00;
149  w << 1.0E+00, 1.0E+00;
150  return true;
151  }
152  case 3 :
153  {
154  x << - 1.0E+00, 0.0E+00, 1.0E+00;
155  w << 1.0 / 3.0E+00, 4.0 / 3.0E+00, 1.0 / 3.0E+00;
156  return true;
157  }
158  case 4 :
159  {
160  x << - 1.0E+00,- 0.447213595499957939281834733746E+00, 0.447213595499957939281834733746E+00, 1.0E+00;
161  w << 1.0E+00 / 6.0E+00, 5.0E+00 / 6.0E+00, 5.0E+00 / 6.0E+00, 1.0E+00 / 6.0E+00;
162  return true;
163  }
164  case 5 :
165  {
166  x << - 1.0E+00, - 0.654653670707977143798292456247E+00, 0.0E+00, 0.654653670707977143798292456247E+00, 1.0E+00;
167  w << 9.0E+00 / 90.0E+00, 49.0E+00 / 90.0E+00, 64.0E+00 / 90.0E+0, 49.0E+00 / 90.0E+00, 9.0E+00 / 90.0E+00;
168  return true;
169  }
170  case 6 :
171  {
172  x << - 1.0E+00,- 0.765055323929464692851002973959E+00,- 0.285231516480645096314150994041E+00, 0.285231516480645096314150994041E+00, 0.765055323929464692851002973959E+00, 1.0E+00;
173  w << 0.066666666666666666666666666667E+00, 0.378474956297846980316612808212E+00, 0.554858377035486353016720525121E+00, 0.554858377035486353016720525121E+00, 0.378474956297846980316612808212E+00, 0.066666666666666666666666666667E+00;
174  return true;
175  }
176  case 7 :
177  {
178  x << - 1.0E+00, - 0.830223896278566929872032213967E+00, - 0.468848793470714213803771881909E+00, 0.0E+00, 0.468848793470714213803771881909E+00, 0.830223896278566929872032213967E+00, 1.0E+00;
179  w << 0.476190476190476190476190476190E-01, 0.276826047361565948010700406290E+00, 0.431745381209862623417871022281E+00, 0.487619047619047619047619047619E+00, 0.431745381209862623417871022281E+00, 0.276826047361565948010700406290E+00, 0.476190476190476190476190476190E-01;
180  return true;
181  }
182  case 8 :
183  {
184  x << - 1.0E+00, - 0.871740148509606615337445761221E+00, - 0.591700181433142302144510731398E+00, - 0.209299217902478868768657260345E+00, 0.209299217902478868768657260345E+00, 0.591700181433142302144510731398E+00, 0.871740148509606615337445761221E+00, 1.0E+00;
185  w << 0.357142857142857142857142857143E-01, 0.210704227143506039382991065776E+00, 0.341122692483504364764240677108E+00, 0.412458794658703881567052971402E+00, 0.412458794658703881567052971402E+00, 0.341122692483504364764240677108E+00, 0.210704227143506039382991065776E+00, 0.357142857142857142857142857143E-01;
186  return true;
187  }
188 
189  }
190 
191  switch (n)
192  {
193  case 9 :
194  {
195  x << - 1.0E+00, - 0.899757995411460157312345244418E+00, - 0.677186279510737753445885427091E+00, - 0.363117463826178158710752068709E+00, 0.0E+00, 0.363117463826178158710752068709E+00, 0.677186279510737753445885427091E+00, 0.899757995411460157312345244418E+00, 1.0E+00;
196 
197  w << 0.277777777777777777777777777778E-01, 0.165495361560805525046339720029E+00, 0.274538712500161735280705618579E+00, 0.346428510973046345115131532140E+00, 0.371519274376417233560090702948E+00, 0.346428510973046345115131532140E+00, 0.274538712500161735280705618579E+00, 0.165495361560805525046339720029E+00, 0.277777777777777777777777777778E-01;
198  return true;
199  }
200  case 10 :
201  {
202  x << - 1.0E+00, - 0.919533908166458813828932660822E+00, - 0.738773865105505075003106174860E+00, - 0.477924949810444495661175092731E+00, - 0.165278957666387024626219765958E+00, 0.165278957666387024626219765958E+00, 0.477924949810444495661175092731E+00, 0.738773865105505075003106174860E+00, 0.919533908166458813828932660822E+00, 1.0E+00;
203 
204  w << 0.222222222222222222222222222222E-01, 0.133305990851070111126227170755E+00, 0.224889342063126452119457821731E+00, 0.292042683679683757875582257374E+00, 0.327539761183897456656510527917E+00, 0.327539761183897456656510527917E+00, 0.292042683679683757875582257374E+00, 0.224889342063126452119457821731E+00, 0.133305990851070111126227170755E+00, 0.222222222222222222222222222222E-01;
205  return true;
206  }
207  case 11 :
208  {
209  x << - 1.0E+00, - 0.934001430408059134332274136099E+00, - 0.784483473663144418622417816108E+00, - 0.565235326996205006470963969478E+00, - 0.295758135586939391431911515559E+00, 0.0E+00, 0.295758135586939391431911515559E+00, 0.565235326996205006470963969478E+00, 0.784483473663144418622417816108E+00, 0.934001430408059134332274136099E+00, 1.0E+00;
210 
211  w << 0.181818181818181818181818181818E-01, 0.109612273266994864461403449580E+00, 0.187169881780305204108141521899E+00, 0.248048104264028314040084866422E+00, 0.286879124779008088679222403332E+00, 0.300217595455690693785931881170E+00, 0.286879124779008088679222403332E+00, 0.248048104264028314040084866422E+00, 0.187169881780305204108141521899E+00, 0.109612273266994864461403449580E+00, 0.181818181818181818181818181818E-01;
212 
213  return true;
214  }
215  case 12 :
216  {
217  x << - 1.0E+00, - 0.944899272222882223407580138303E+00, - 0.819279321644006678348641581717E+00, - 0.632876153031869677662404854444E+00, - 0.399530940965348932264349791567E+00, - 0.136552932854927554864061855740E+00, 0.136552932854927554864061855740E+00, 0.399530940965348932264349791567E+00, 0.632876153031869677662404854444E+00, 0.819279321644006678348641581717E+00, 0.944899272222882223407580138303E+00, 1.0E+00;
218 
219  w << 0.151515151515151515151515151515E-01, 0.916845174131961306683425941341E-01, 0.157974705564370115164671062700E+00, 0.212508417761021145358302077367E+00, 0.251275603199201280293244412148E+00, 0.271405240910696177000288338500E+00, 0.271405240910696177000288338500E+00, 0.251275603199201280293244412148E+00, 0.212508417761021145358302077367E+00, 0.157974705564370115164671062700E+00, 0.916845174131961306683425941341E-01, 0.151515151515151515151515151515E-01;
220  return true;
221  }
222  case 13 :
223  {
224  x << - 1.0E+00, - 0.953309846642163911896905464755E+00, - 0.846347564651872316865925607099E+00, - 0.686188469081757426072759039566E+00, - 0.482909821091336201746937233637E+00, - 0.249286930106239992568673700374E+00, 0.0E+00, 0.249286930106239992568673700374E+00, 0.482909821091336201746937233637E+00, 0.686188469081757426072759039566E+00, 0.846347564651872316865925607099E+00, 0.953309846642163911896905464755E+00, 1.0E+00;
225 
226  w << 0.128205128205128205128205128205E-01, 0.778016867468189277935889883331E-01, 0.134981926689608349119914762589E+00, 0.183646865203550092007494258747E+00, 0.220767793566110086085534008379E+00, 0.244015790306676356458578148360E+00, 0.251930849333446736044138641541E+00, 0.244015790306676356458578148360E+00, 0.220767793566110086085534008379E+00, 0.183646865203550092007494258747E+00, 0.134981926689608349119914762589E+00, 0.778016867468189277935889883331E-01, 0.128205128205128205128205128205E-01;
227  return true;
228  }
229  case 14 :
230  {
231  x << - 1.0E+00, - 0.959935045267260901355100162015E+00, - 0.867801053830347251000220202908E+00, - 0.728868599091326140584672400521E+00, - 0.550639402928647055316622705859E+00, - 0.342724013342712845043903403642E+00, - 0.116331868883703867658776709736E+00, 0.116331868883703867658776709736E+00, 0.342724013342712845043903403642E+00, 0.550639402928647055316622705859E+00, 0.728868599091326140584672400521E+00, 0.867801053830347251000220202908E+00, 0.959935045267260901355100162015E+00, 1.0E+00;
232 
233  w << 0.109890109890109890109890109890E-01, 0.668372844976812846340706607461E-01, 0.116586655898711651540996670655E+00, 0.160021851762952142412820997988E+00, 0.194826149373416118640331778376E+00, 0.219126253009770754871162523954E+00, 0.231612794468457058889628357293E+00, 0.231612794468457058889628357293E+00, 0.219126253009770754871162523954E+00, 0.194826149373416118640331778376E+00, 0.160021851762952142412820997988E+00, 0.116586655898711651540996670655E+00, 0.668372844976812846340706607461E-01, 0.109890109890109890109890109890E-01;
234  return true;
235  }
236  case 15 :
237  {
238  x << - 1.0E+00, - 0.965245926503838572795851392070E+00, - 0.885082044222976298825401631482E+00, - 0.763519689951815200704118475976E+00, - 0.606253205469845711123529938637E+00, - 0.420638054713672480921896938739E+00, - 0.215353955363794238225679446273E+00, 0.0E+00, 0.215353955363794238225679446273E+00, 0.420638054713672480921896938739E+00, 0.606253205469845711123529938637E+00, 0.763519689951815200704118475976E+00, 0.885082044222976298825401631482E+00, 0.965245926503838572795851392070E+00, 1.0E+00;
239 
240  w << 0.952380952380952380952380952381E-02, 0.580298930286012490968805840253E-01, 0.101660070325718067603666170789E+00, 0.140511699802428109460446805644E+00, 0.172789647253600949052077099408E+00, 0.196987235964613356092500346507E+00, 0.211973585926820920127430076977E+00, 0.217048116348815649514950214251E+00, 0.211973585926820920127430076977E+00, 0.196987235964613356092500346507E+00, 0.172789647253600949052077099408E+00, 0.140511699802428109460446805644E+00, 0.101660070325718067603666170789E+00, 0.580298930286012490968805840253E-01, 0.952380952380952380952380952381E-02;
241  return true;
242  }
243  case 16 :
244  {
245  x << - 1.0E+00, - 0.969568046270217932952242738367E+00, - 0.899200533093472092994628261520E+00, - 0.792008291861815063931088270963E+00, - 0.652388702882493089467883219641E+00, - 0.486059421887137611781890785847E+00, - 0.299830468900763208098353454722E+00, - 0.101326273521949447843033005046E+00, 0.101326273521949447843033005046E+00, 0.299830468900763208098353454722E+00, 0.486059421887137611781890785847E+00, 0.652388702882493089467883219641E+00, 0.792008291861815063931088270963E+00, 0.899200533093472092994628261520E+00, 0.969568046270217932952242738367E+00, 1.0E+00;
246 
247  w << 0.833333333333333333333333333333E-02, 0.508503610059199054032449195655E-01, 0.893936973259308009910520801661E-01, 0.124255382132514098349536332657E+00, 0.154026980807164280815644940485E+00, 0.177491913391704125301075669528E+00, 0.193690023825203584316913598854E+00, 0.201958308178229871489199125411E+00, 0.201958308178229871489199125411E+00, 0.193690023825203584316913598854E+00, 0.177491913391704125301075669528E+00, 0.154026980807164280815644940485E+00, 0.124255382132514098349536332657E+00, 0.893936973259308009910520801661E-01, 0.508503610059199054032449195655E-01, 0.833333333333333333333333333333E-02;
248  return true;
249  }
250  case 17 :
251  {
252  x << - 1.0E+00, - 0.973132176631418314156979501874E+00, - 0.910879995915573595623802506398E+00, - 0.815696251221770307106750553238E+00, - 0.691028980627684705394919357372E+00, - 0.541385399330101539123733407504E+00, - 0.372174433565477041907234680735E+00, - 0.189511973518317388304263014753E+00, 0.0E+00, 0.189511973518317388304263014753E+00, 0.372174433565477041907234680735E+00, 0.541385399330101539123733407504E+00, 0.691028980627684705394919357372E+00, 0.815696251221770307106750553238E+00, 0.910879995915573595623802506398E+00, 0.973132176631418314156979501874E+00, 1.0E+00;
253 
254  w << 0.735294117647058823529411764706E-02, 0.449219405432542096474009546232E-01, 0.791982705036871191902644299528E-01, 0.110592909007028161375772705220E+00, 0.137987746201926559056201574954E+00, 0.160394661997621539516328365865E+00, 0.177004253515657870436945745363E+00, 0.187216339677619235892088482861E+00, 0.190661874753469433299407247028E+00, 0.187216339677619235892088482861E+00, 0.177004253515657870436945745363E+00, 0.160394661997621539516328365865E+00, 0.137987746201926559056201574954E+00, 0.110592909007028161375772705220E+00, 0.791982705036871191902644299528E-01, 0.449219405432542096474009546232E-01, 0.735294117647058823529411764706E-02;
255  return true;
256  }
257  case 18 :
258  {
259  x << - 1.0E+00, - 0.976105557412198542864518924342E+00, - 0.920649185347533873837854625431E+00, - 0.835593535218090213713646362328E+00, - 0.723679329283242681306210365302E+00, - 0.588504834318661761173535893194E+00, - 0.434415036912123975342287136741E+00, - 0.266362652878280984167665332026E+00, - 0.897490934846521110226450100886E-01, 0.897490934846521110226450100886E-01, 0.266362652878280984167665332026E+00, 0.434415036912123975342287136741E+00,0.588504834318661761173535893194E+00,0.723679329283242681306210365302E+00,0.835593535218090213713646362328E+00,0.920649185347533873837854625431E+00,0.976105557412198542864518924342E+00,1.0E+00;
260 
261  w << 0.653594771241830065359477124183E-02, 0.399706288109140661375991764101E-01, 0.706371668856336649992229601678E-01, 0.990162717175028023944236053187E-01, 0.124210533132967100263396358897E+00, 0.145411961573802267983003210494E+00, 0.161939517237602489264326706700E+00, 0.173262109489456226010614403827E+00, 0.179015863439703082293818806944E+00, 0.179015863439703082293818806944E+00, 0.173262109489456226010614403827E+00, 0.161939517237602489264326706700E+00, 0.145411961573802267983003210494E+00, 0.124210533132967100263396358897E+00, 0.990162717175028023944236053187E-01, 0.706371668856336649992229601678E-01, 0.399706288109140661375991764101E-01, 0.653594771241830065359477124183E-02;
262  return true;
263  }
264  case 19 :
265  {
266  x << - 1.0E+00, - 0.978611766222080095152634063110E+00, - 0.928901528152586243717940258797E+00, - 0.852460577796646093085955970041E+00, - 0.751494202552613014163637489634E+00, - 0.628908137265220497766832306229E+00, - 0.488229285680713502777909637625E+00, - 0.333504847824498610298500103845E+00, - 0.169186023409281571375154153445E+00, 0.0E+00, 0.169186023409281571375154153445E+00, 0.333504847824498610298500103845E+00,0.488229285680713502777909637625E+00,0.628908137265220497766832306229E+00,0.751494202552613014163637489634E+00,0.852460577796646093085955970041E+00,0.928901528152586243717940258797E+00,0.978611766222080095152634063110E+00,1.0E+00;
267 
268  w << 0.584795321637426900584795321637E-02, 0.357933651861764771154255690351E-01, 0.633818917626297368516956904183E-01, 0.891317570992070844480087905562E-01, 0.112315341477305044070910015464E+00, 0.132267280448750776926046733910E+00, 0.148413942595938885009680643668E+00, 0.160290924044061241979910968184E+00, 0.167556584527142867270137277740E+00, 0.170001919284827234644672715617E+00, 0.167556584527142867270137277740E+00, 0.160290924044061241979910968184E+00, 0.148413942595938885009680643668E+00, 0.132267280448750776926046733910E+00, 0.112315341477305044070910015464E+00, 0.891317570992070844480087905562E-01, 0.633818917626297368516956904183E-01, 0.357933651861764771154255690351E-01, 0.584795321637426900584795321637E-02;
269  return true;
270  }
271  case 20 :
272  {
273  x << - 1.0E+00, - 0.980743704893914171925446438584E+00, - 0.935934498812665435716181584931E+00, - 0.866877978089950141309847214616E+00, - 0.775368260952055870414317527595E+00, - 0.663776402290311289846403322971E+00, - 0.534992864031886261648135961829E+00, - 0.392353183713909299386474703816E+00, - 0.239551705922986495182401356927E+00,- 0.805459372388218379759445181596E-01, 0.805459372388218379759445181596E-01, 0.239551705922986495182401356927E+00,0.392353183713909299386474703816E+00,0.534992864031886261648135961829E+00,0.663776402290311289846403322971E+00,0.775368260952055870414317527595E+00,0.866877978089950141309847214616E+00,0.935934498812665435716181584931E+00,0.980743704893914171925446438584E+00, 1.0E+00;
274 
275  w << 0.526315789473684210526315789474E-02, 0.322371231884889414916050281173E-01, 0.571818021275668260047536271732E-01, 0.806317639961196031447768461137E-01, 0.101991499699450815683781205733E+00, 0.120709227628674725099429705002E+00, 0.136300482358724184489780792989E+00, 0.148361554070916825814713013734E+00, 0.156580102647475487158169896794E+00, 0.160743286387845749007726726449E+00, 0.160743286387845749007726726449E+00, 0.156580102647475487158169896794E+00, 0.148361554070916825814713013734E+00, 0.136300482358724184489780792989E+00, 0.120709227628674725099429705002E+00, 0.101991499699450815683781205733E+00, 0.806317639961196031447768461137E-01, 0.571818021275668260047536271732E-01, 0.322371231884889414916050281173E-01, 0.526315789473684210526315789474E-02;
276  return true;
277  }
278  default:
279  {
280  //gsWarn << " Illegal value of N = " << n << "\n";
281  gsWarn << "Precomputed Lobatto rule (1,..,20) not found for N="<<n<<".\n";
282  return false;
283  }
284 
285  }//switch
286 
287 }// lookupReference
288 
289 } // namespace gismo
short_t
#define short_t
Definition: gsConfig.h:35
gsBasis.h
Provides declaration of Basis abstract interface.
index_t
#define index_t
Definition: gsConfig.h:32
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
gsOptionList.h
Provides a list of labeled parameters/options that can be set and accessed easily.
gsWarn
#define gsWarn
Definition: gsDebug.h:50
gismo::gsLobattoRule::lookupReference
static bool lookupReference(index_t n, gsVector< T > &x, gsVector< T > &w)
Look up function for the Lobatto quadrature rule in the interval [-1,1].
Definition: gsLobattoRule.hpp:123
gismo::expr::abs
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4486
gismo::gsLobattoRule::setNodes
void setNodes(gsVector< index_t > const &numNodes, unsigned digits=0)
Initialize quadrature rule with numNodes number of quadrature points per integration variable...
Definition: gsLobattoRule.hpp:23
gismo::gsLobattoRule::computeReference
static void computeReference(index_t n, gsVector< T > &x, gsVector< T > &w, unsigned digits=0)
Computes the Lobatto quadrature rule with n nodes in the interval [-1,1].
Definition: gsLobattoRule.hpp:59