/home/runner/work/kynema/kynema/kynema/src/math/gll_quadrature.hpp Source File

Kynema API: /home/runner/work/kynema/kynema/kynema/src/math/gll_quadrature.hpp Source File
Kynema API
A flexible multibody structural dynamics code for wind turbines
Loading...
Searching...
No Matches
gll_quadrature.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <cmath>
4#include <stdexcept>
5#include <vector>
6
7namespace kynema::math {
8
9inline std::vector<double> GetGllLocations(size_t order) {
10 switch (order) {
11 case 1UL:
12 return std::vector{-1., 1.};
13 case 2UL:
14 return std::vector{-1., 0., 1.};
15 case 3UL:
16 return std::vector{-1., -std::sqrt(1. / 5.), std::sqrt(1. / 5.), 1.};
17 case 4UL:
18 return std::vector{-1., -std::sqrt(3. / 7.), 0., std::sqrt(3. / 7.), 1.};
19 case 5UL:
20 return std::vector{
21 -1.,
22 -std::sqrt((1. / 3.) + (2. * std::sqrt(7.) / 21.)),
23 -std::sqrt((1. / 3.) - (2. * std::sqrt(7.) / 21.)),
24 std::sqrt((1. / 3.) - (2. * std::sqrt(7.) / 21.)),
25 std::sqrt((1. / 3.) + (2. * std::sqrt(7.) / 21.)),
26 1.
27 };
28 case 6UL:
29 return std::vector{
30 -1.,
31 -std::sqrt((5. / 11.) + ((2. / 11.) * std::sqrt(5. / 3.))),
32 -std::sqrt((5. / 11.) - ((2. / 11.) * std::sqrt(5. / 3.))),
33 0.,
34 std::sqrt((5. / 11.) - ((2. / 11.) * std::sqrt(5. / 3.))),
35 std::sqrt((5. / 11.) + ((2. / 11.) * std::sqrt(5. / 3.))),
36 1.
37 };
38 case 7UL:
39 return std::vector{
40 -1.,
41 -0.87174014850960657,
42 -0.59170018143314229,
43 -0.20929921790247885,
44 0.20929921790247885,
45 0.59170018143314229,
46 0.87174014850960657,
47 1.
48 };
49 case 8UL:
50 return std::vector{-1., -0.89975799541146018, -0.67718627951073773, -0.36311746382617816,
51 0., 0.36311746382617816, 0.67718627951073773, 0.89975799541146018,
52 1.};
53 case 9UL:
54 return std::vector{
55 -1.,
56 -0.91953390816645886,
57 -0.73877386510550502,
58 -0.47792494981044453,
59 -0.16527895766638701,
60 0.16527895766638701,
61 0.47792494981044448,
62 0.73877386510550502,
63 0.91953390816645886,
64 1.
65 };
66 case 10UL:
67 return std::vector{-1.,
68 -0.93400143040805916,
69 -0.78448347366314441,
70 -0.56523532699620493,
71 -0.29575813558693936,
72 0.,
73 0.29575813558693936,
74 0.56523532699620493,
75 0.78448347366314441,
76 0.93400143040805916,
77 1.};
78 case 11UL:
79 return std::vector{
80 -1.,
81 -0.94489927222288217,
82 -0.81927932164400663,
83 -0.63287615303186073,
84 -0.39953094096534891,
85 -0.13655293285492756,
86 0.13655293285492756,
87 0.39953094096534891,
88 0.63287615303186073,
89 0.81927932164400663,
90 0.94489927222288217,
91 1.
92 };
93 case 12UL:
94 return std::vector{
95 -1.,
96 -0.95330984664216400,
97 -0.84634756465187233,
98 -0.68618846908175746,
99 -0.48290982109133623,
100 -0.24928693010623998,
101 0.,
102 0.24928693010623998,
103 0.48290982109133623,
104 0.68618846908175746,
105 0.84634756465187233,
106 0.95330984664216400,
107 1.
108 };
109 case 13UL:
110 return std::vector{
111 -1.,
112 -0.95993504526726092,
113 -0.86780105383034722,
114 -0.72886859909132617,
115 -0.55063940292864699,
116 -0.34272401334271285,
117 -0.11633186888370387,
118 0.11633186888370387,
119 0.34272401334271285,
120 0.55063940292864710,
121 0.72886859909132617,
122 0.86780105383034722,
123 0.95993504526726092,
124 1.
125 };
126 case 14UL:
127 return std::vector{
128 -1.,
129 -0.96524592650383856,
130 -0.88508204422297632,
131 -0.76351968995181518,
132 -0.60625320546984574,
133 -0.42063805471367249,
134 -0.21535395536379423,
135 0.,
136 0.21535395536379423,
137 0.42063805471367249,
138 0.60625320546984574,
139 0.76351968995181518,
140 0.88508204422297632,
141 0.96524592650383856,
142 1.
143 };
144 case 15UL:
145 return std::vector{
146 -1.,
147 -0.96956804627021798,
148 -0.89920053309347214,
149 -0.79200829186181509,
150 -0.65238870288249307,
151 -0.48605942188713763,
152 -0.29983046890076320,
153 -0.10132627352194944,
154 0.10132627352194944,
155 0.29983046890076320,
156 0.48605942188713763,
157 0.65238870288249307,
158 0.79200829186181509,
159 0.89920053309347214,
160 0.96956804627021798,
161 1.
162 };
163 case 16UL:
164 return std::vector{
165 -1.0,
166 -0.9731321766314183,
167 -0.9108799959155736,
168 -0.8156962512217704,
169 -0.6910289806276847,
170 -0.5413853993301015,
171 -0.372174433565477,
172 -0.1895119735183174,
173 0.0,
174 0.1895119735183174,
175 0.372174433565477,
176 0.5413853993301015,
177 0.6910289806276847,
178 0.8156962512217704,
179 0.9108799959155736,
180 0.9731321766314183,
181 1.0
182 };
183 case 17UL:
184 return std::vector{
185 -1.0,
186 -0.9761055574121985,
187 -0.9206491853475339,
188 -0.8355935352180902,
189 -0.7236793292832426,
190 -0.5885048343186617,
191 -0.43441503691212396, // NOLINT
192 -0.266362652878281,
193 -0.08974909348465213,
194 0.08974909348465213,
195 0.266362652878281,
196 0.43441503691212396, // NOLINT
197 0.5885048343186617,
198 0.7236793292832426,
199 0.8355935352180902,
200 0.9206491853475339,
201 0.9761055574121985,
202 1.0
203 };
204 case 18UL:
205 return std::vector{
206 -1.0,
207 -0.9786117662220801,
208 -0.9289015281525863,
209 -0.8524605777966461,
210 -0.751494202552613,
211 -0.6289081372652205,
212 -0.4882292856807135,
213 -0.33350484782449863,
214 -0.16918602340928157,
215 0.0,
216 0.16918602340928157,
217 0.33350484782449863,
218 0.4882292856807135,
219 0.6289081372652205,
220 0.751494202552613,
221 0.8524605777966461,
222 0.9289015281525863,
223 0.9786117662220801,
224 1.0
225 };
226 case 19UL:
227 return std::vector{
228 -1.0,
229 -0.9807437048939142,
230 -0.9359344988126654,
231 -0.8668779780899502,
232 -0.7753682609520559,
233 -0.6637764022903113,
234 -0.5349928640318863,
235 -0.39235318371390926,
236 -0.2395517059229865,
237 -0.08054593723882184,
238 0.08054593723882184,
239 0.2395517059229865,
240 0.39235318371390926,
241 0.5349928640318863,
242 0.6637764022903113,
243 0.7753682609520559,
244 0.8668779780899502,
245 0.9359344988126654,
246 0.9807437048939142,
247 1.0
248 };
249 case 20UL:
250 return std::vector{
251 -1.0,
252 -0.982572296604548,
253 -0.9419762969597455,
254 -0.8792947553235905,
255 -0.7960019260777124,
256 -0.6940510260622232, // NOLINT
257 -0.5758319602618307,
258 -0.4441157832790021,
259 -0.30198985650876486,
260 -0.15278551580218547,
261 0.0,
262 0.15278551580218547,
263 0.30198985650876486,
264 0.4441157832790021,
265 0.5758319602618307,
266 0.6940510260622232, // NOLINT
267 0.7960019260777124,
268 0.8792947553235905,
269 0.9419762969597455,
270 0.982572296604548,
271 1.0
272 };
273 case 100UL:
274 return std::vector{
275 -1.0,
276 -0.9992732578112129,
277 -0.997564432044587,
278 -0.9948806385217545,
279 -0.9912247073833753,
280 -0.9866002474942588,
281 -0.9810117878154722,
282 -0.9744647926317048,
283 -0.9669656606203273,
284 -0.9585217199147787,
285 -0.9491412214185269,
286 -0.9388333309390072,
287 -0.9276081203214385,
288 -0.9154765576522249,
289 -0.9024504965657505,
290 -0.8885426646755374,
291 -0.873766651145992,
292 -0.8581368934193241,
293 -0.8416686631118906,
294 -0.8243780510944396,
295 -0.8062819517712353,
296 -0.7873980465736545,
297 -0.7677447866845291,
298 -0.7473413750101983,
299 -0.7262077474179458,
300 -0.7043645532571894,
301 -0.6818331351834797,
302 -0.6586355083050388,
303 -0.6347943386722272,
304 -0.6103329211309613,
305 -0.5852751565617291,
306 -0.5596455285264476,
307 -0.5334690793459793,
308 -0.5067713856316896,
309 -0.4795785332949497,
310 -0.4519170920590072,
311 -0.4238140894981276,
312 -0.39529698462937374,
313 -0.3663936410828277,
314 -0.33713229987646925,
315 -0.30754155182231335,
316 -0.2776503095907678,
317 -0.24748777946050743,
318 -0.21708343278146724,
319 -0.18646697717883803,
320 -0.15566832752620036,
321 -0.12471757671615918,
322 -0.0936449662570384,
323 -0.06248085672436486,
324 -0.03125569809601377,
325 0.0,
326 0.03125569809601377,
327 0.06248085672436486,
328 0.0936449662570384,
329 0.12471757671615918,
330 0.15566832752620036,
331 0.18646697717883803,
332 0.21708343278146724,
333 0.24748777946050743,
334 0.2776503095907678,
335 0.30754155182231335,
336 0.33713229987646925,
337 0.3663936410828277,
338 0.39529698462937374,
339 0.4238140894981276,
340 0.4519170920590072,
341 0.4795785332949497,
342 0.5067713856316896,
343 0.5334690793459793,
344 0.5596455285264476,
345 0.5852751565617291,
346 0.6103329211309613,
347 0.6347943386722272,
348 0.6586355083050388,
349 0.6818331351834797,
350 0.7043645532571894,
351 0.7262077474179458,
352 0.7473413750101983,
353 0.7677447866845291,
354 0.7873980465736545,
355 0.8062819517712353,
356 0.8243780510944396,
357 0.8416686631118906,
358 0.8581368934193241,
359 0.873766651145992,
360 0.8885426646755374,
361 0.9024504965657505,
362 0.9154765576522249,
363 0.9276081203214385,
364 0.9388333309390072,
365 0.9491412214185269,
366 0.9585217199147787,
367 0.9669656606203273,
368 0.9744647926317048,
369 0.9810117878154722,
370 0.9866002474942588,
371 0.9912247073833753,
372 0.9948806385217545,
373 0.997564432044587,
374 0.9992732578112129,
375 1.0,
376 };
377 default:
378 throw std::runtime_error("Supported orders are 1 - 20");
379 }
380}
381
382inline std::vector<double> GetGllWeights(size_t order) {
383 switch (order) {
384 case 1UL:
385 return std::vector{1., 1.};
386 case 2UL:
387 return std::vector{1. / 3., 4. / 3., 1. / 3.};
388 case 3UL:
389 return std::vector{1. / 6., 5. / 6., 5. / 6., 1. / 6.};
390 case 4UL:
391 return std::vector{1. / 10., 49. / 90., 32. / 45., 49. / 90., 1. / 10.};
392 case 5UL:
393 return std::vector{
394 1. / 15.,
395 (14. - std::sqrt(7.)) / 30.,
396 (14. + std::sqrt(7.)) / 30.,
397 (14. + std::sqrt(7.)) / 30.,
398 (14. - std::sqrt(7.)) / 30.,
399 1. / 15.
400 };
401 case 6UL:
402 return std::vector{1. / 21.,
403 (124. - 7. * std::sqrt(15.)) / 350.,
404 (124. + 7. * std::sqrt(15.)) / 350.,
405 256. / 525.,
406 (124. + 7. * std::sqrt(15.)) / 350.,
407 (124. - 7. * std::sqrt(15.)) / 350,
408 1. / 21.};
409 case 7UL:
410 return std::vector{3.5714285714285712E-002, 0.21070422714350615, 0.34112269248350441,
411 0.41245879465870372, 0.41245879465870372, 0.34112269248350441,
412 0.21070422714350615, 3.5714285714285712E-002};
413 case 8UL:
414 return std::vector{2.7777777777777776E-002, 0.16549536156080552,
415 0.27453871250016160, 0.34642851097304617,
416 0.37151927437641724, 0.34642851097304617,
417 0.27453871250016160, 0.16549536156080552,
418 2.7777777777777776E-002};
419 case 9UL:
420 return std::vector{2.2222222222222223E-002, 0.13330599085107006, 0.22488934206312636,
421 0.29204268367968378, 0.32753976118389744, 0.32753976118389744,
422 0.29204268367968378, 0.22488934206312636, 0.13330599085107006,
423 2.2222222222222223E-002};
424 case 10UL:
425 return std::vector{1.8181818181818181E-002, 0.10961227326699498, 0.18716988178030541,
426 0.24804810426402829, 0.28687912477900801, 0.30021759545569071,
427 0.28687912477900823, 0.24804810426402829, 0.18716988178030541,
428 0.10961227326699498, 1.8181818181818181E-002};
429 case 11UL:
430 return std::vector{1.5151515151515152E-002, 9.1684517413196109E-002,
431 0.15797470556437004, 0.21250841776102122,
432 0.25127560319920106, 0.27140524091069618,
433 0.27140524091069618, 0.25127560319920106,
434 0.21250841776102122, 0.15797470556437004,
435 9.1684517413196109E-002, 1.5151515151515152E-002};
436 case 12UL:
437 return std::vector{1.2820512820512820E-002, 7.7801686746818866E-002,
438 0.13498192668960840, 0.18364686520355006,
439 0.22076779356611012, 0.24401579030667625,
440 0.25193084933344673, 0.24401579030667625,
441 0.22076779356611012, 0.18364686520355006,
442 0.13498192668960840, 7.7801686746818866E-002,
443 1.2820512820512820E-002};
444 case 13UL:
445 return std::vector{1.0989010989010990E-002, 6.6837284497681185E-002, 0.11658665589871173,
446 0.16002185176295217, 0.19482614937341614, 0.21912625300977098,
447 0.23161279446845698, 0.23161279446845698, 0.21912625300977098,
448 0.19482614937341600, 0.16002185176295217, 0.11658665589871173,
449 6.6837284497681185E-002, 1.0989010989010990E-002};
450 case 14UL:
451 return std::vector{
452 9.5238095238095247E-003, 5.8029893028601086E-002, 0.10166007032571801,
453 0.14051169980242798, 0.17278964725360088, 0.19698723596461334,
454 0.21197358592682095, 0.21704811634881566, 0.21197358592682075,
455 0.19698723596461334, 0.17278964725360088, 0.14051169980242798,
456 0.10166007032571801, 5.8029893028601086E-002, 9.5238095238095247E-003
457 };
458 case 15UL:
459 return std::vector{8.3333333333333332E-003, 5.0850361005920039E-002,
460 8.9393697325930832E-002, 0.12425538213251400,
461 0.15402698080716443, 0.17749191339170411,
462 0.19369002382520362, 0.20195830817822993,
463 0.20195830817822993, 0.19369002382520362,
464 0.17749191339170411, 0.15402698080716443,
465 0.12425538213251400, 8.9393697325930832E-002,
466 5.0850361005920039E-002, 8.3333333333333332E-003};
467 case 16UL:
468 return std::vector{0.007352941176470588, 0.04492194054325405, 0.0791982705036871,
469 0.11059290900702815, 0.13798774620192658, 0.16039466199762148,
470 0.1770042535156577, 0.18721633967761928, 0.19066187475346943,
471 0.18721633967761928, 0.1770042535156577, 0.16039466199762148,
472 0.13798774620192658, 0.11059290900702815, 0.0791982705036871,
473 0.04492194054325405, 0.007352941176470588};
474 case 17UL:
475 return std::vector{0.006535947712418301, 0.03997062881091395, 0.07063716688563365,
476 0.09901627171750278, 0.12421053313296708, 0.14541196157380232,
477 0.16193951723760242, 0.17326210948945636, 0.17901586343970305,
478 0.17901586343970305, 0.17326210948945636, 0.16193951723760242,
479 0.14541196157380232, 0.12421053313296708, 0.09901627171750278,
480 0.07063716688563365, 0.03997062881091395, 0.006535947712418301};
481 case 18UL:
482 return std::vector{0.005847953216374269, 0.035793365186176644, 0.06338189176262979,
483 0.08913175709920705, 0.11231534147730492, 0.1322672804487507,
484 0.14841394259593885, 0.1602909240440612, 0.16755658452714278,
485 0.17000191928482725, 0.16755658452714278, 0.1602909240440612,
486 0.14841394259593885, 0.1322672804487507, 0.11231534147730492,
487 0.08913175709920705, 0.06338189176262979, 0.035793365186176644,
488 0.005847953216374269};
489 case 19UL:
490 return std::vector{0.005263157894736842, 0.03223712318848893, 0.0571818021275668,
491 0.08063176399611967, 0.10199149969945068, 0.12070922762867466,
492 0.13630048235872422, 0.14836155407091683, 0.15658010264747543,
493 0.16074328638784569, 0.16074328638784569, 0.15658010264747543,
494 0.14836155407091683, 0.13630048235872422, 0.12070922762867466,
495 0.10199149969945068, 0.08063176399611967, 0.0571818021275668,
496 0.03223712318848893, 0.005263157894736842};
497 case 20UL:
498 return std::vector{0.004761904761904762, 0.029184840098505565, 0.05184316900084964,
499 0.07327391818507417, 0.092985467957886, 0.11051708321912326,
500 0.1254581211908689, 0.1374584628600413, 0.14623686244797737,
501 0.15158757511168136, 0.15338519033217496, 0.15158757511168136,
502 0.14623686244797737, 0.1374584628600413, 0.1254581211908689,
503 0.11051708321912326, 0.092985467957886, 0.07327391818507417,
504 0.05184316900084964, 0.029184840098505565, 0.004761904761904762};
505 case 100UL:
506 return std::vector{0.00019801980198019803, 0.0012204276494756096, 0.0021967379918678047,
507 0.0031703893238251937, 0.004140872623634033, 0.005107292474504979,
508 0.00606871613662751, 0.007024207698695935, 0.00797283485601666,
509 0.008913671273050698, 0.009845797936967065, 0.01076830421785277,
510 0.011680288824616627, 0.012580860715228265, 0.013469139981628431,
511 0.014344258716811612, 0.015205361866815714, 0.016051608068429056,
512 0.01688217047259503, 0.017696237553117427, 0.018493013900093766,
513 0.01927172099742115, 0.02003159798368662, 0.02077190239573547,
514 0.021491910894219514, 0.022190919970428696, 0.022868246633728155,
515 0.02352322907893218, 0.024155227332968956, 0.02476362388020823,
516 0.025347824265839494, 0.025907257676715513, 0.0264413774990937,
517 0.02694966185272911, 0.02743161410080139, 0.027886763335174196,
518 0.028314664836515926, 0.028714900508829515, 0.02908707928797069,
519 0.029430837523751276, 0.02974583933525808, 0.03003177693903739,
520 0.03028837094982707, 0.030515370653539507, 0.030712554252231405,
521 0.030879729080819606, 0.031016731795331418, 0.031123428532505784,
522 0.031199715040589396, 0.03124551678119947, 0.03126078900215414,
523 0.03124551678119947, 0.031199715040589396, 0.031123428532505784,
524 0.031016731795331418, 0.030879729080819606, 0.030712554252231405,
525 0.030515370653539507, 0.03028837094982707, 0.03003177693903739,
526 0.02974583933525808, 0.029430837523751276, 0.02908707928797069,
527 0.028714900508829515, 0.028314664836515926, 0.027886763335174196,
528 0.02743161410080139, 0.02694966185272911, 0.0264413774990937,
529 0.025907257676715513, 0.025347824265839494, 0.02476362388020823,
530 0.024155227332968956, 0.02352322907893218, 0.022868246633728155,
531 0.022190919970428696, 0.021491910894219514, 0.02077190239573547,
532 0.02003159798368662, 0.01927172099742115, 0.018493013900093766,
533 0.017696237553117427, 0.01688217047259503, 0.016051608068429056,
534 0.015205361866815714, 0.014344258716811612, 0.013469139981628431,
535 0.012580860715228265, 0.011680288824616627, 0.01076830421785277,
536 0.009845797936967065, 0.008913671273050698, 0.00797283485601666,
537 0.007024207698695935, 0.00606871613662751, 0.005107292474504979,
538 0.004140872623634033, 0.0031703893238251937, 0.0021967379918678047,
539 0.0012204276494756096, 0.00019801980198019803};
540 default:
541 throw std::runtime_error("Supported orders are 1 - 20");
542 }
543}
544
545} // namespace kynema::math
Definition gl_quadrature.hpp:8
std::vector< double > GetGllLocations(size_t order)
Definition gll_quadrature.hpp:9
std::vector< double > GetGllWeights(size_t order)
Definition gll_quadrature.hpp:382