G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
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
19namespace gismo
20{
21
22template<class T> void
23gsLobattoRule<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
58template<class T> void
59gsLobattoRule<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
122template<class T> bool
123gsLobattoRule<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
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
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::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gsBasis.h
Provides declaration of Basis abstract interface.
short_t
#define short_t
Definition gsConfig.h:35
index_t
#define index_t
Definition gsConfig.h:32
gsWarn
#define gsWarn
Definition gsDebug.h:50
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.
gismo
The G+Smo namespace, containing all definitions for the library.