Serac  0.1
Serac is an implicit thermal strucural mechanics simulation code.
polynomials.hpp
Go to the documentation of this file.
1 // Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
2 // other Serac Project Developers. See the top-level LICENSE file for
3 // details.
4 //
5 // SPDX-License-Identifier: (BSD-3-Clause)
6 
13 #pragma once
14 
15 #include "tensor.hpp"
16 #include "geometry.hpp"
17 
18 namespace serac {
19 
20 // clang-format off
27 template <int n, typename T = double >
28 SERAC_HOST_DEVICE constexpr tensor<T, n> GaussLobattoNodes(T a = T(0), T b = T(1)) {
29  if constexpr (n == 2) return {a, b};
30  if constexpr (n == 3) return {a, a + 0.5000000000000000 * (b-a), b};
31  if constexpr (n == 4) return {a, a + 0.2763932022500210 * (b-a), a + 0.7236067977499790 * (b-a), b};
32  return tensor<double, n>{};
33 };
34 
42 template <int n, mfem::Geometry::Type geom>
43 SERAC_HOST_DEVICE constexpr auto GaussLegendreNodes() {
44 
45  if constexpr (geom == mfem::Geometry::SEGMENT) {
46  if constexpr (n == 1) return tensor<double,n>{0.50000000000000000};
47  if constexpr (n == 2) return tensor<double,n>{0.2113248654051871, 0.7886751345948129};
48  if constexpr (n == 3) return tensor<double,n>{0.1127016653792583, 0.500000000000000, 0.887298334620742};
49  if constexpr (n == 4) return tensor<double,n>{0.0694318442029737, 0.330009478207572, 0.669990521792428, 0.930568155797026};
50  if constexpr (n == 5) return tensor<double,n>{0.0469100770306680, 0.230765344947158, 0.500000000000000, 0.769234655052842, 0.953089922969332};
51  if constexpr (n == 6) return tensor<double,n>{0.0337652428984240, 0.169395306766868, 0.380690406958402, 0.619309593041598, 0.830604693233132, 0.966234757101576};
52  if constexpr (n == 7) return tensor<double,n>{0.0254460438286207, 0.129234407200303, 0.297077424311301, 0.500000000000000, 0.702922575688699, 0.87076559279970, 0.97455395617138};
53  if constexpr (n == 8) return tensor<double,n>{0.0198550717512319, 0.101666761293187, 0.237233795041836, 0.408282678752175, 0.591717321247825, 0.76276620495816, 0.89833323870681, 0.98014492824877};
54  }
55 
56 
57  if constexpr (geom == mfem::Geometry::TRIANGLE) {
58  using output_t = tensor< double, n * (n + 1) / 2, 2 >;
59 
60  if constexpr (n == 1) {
61  return output_t{{{0.333333333333333, 0.333333333333333}}};
62  }
63 
64  if constexpr (n == 2) {
65  return output_t{{{0.166666666666667, 0.166666666666667},
66  {0.166666666666667, 0.666666666666667},
67  {0.666666666666667, 0.166666666666667}}};
68  }
69 
70  if constexpr (n == 3) {
71  return output_t{{{0.09157621350978, 0.09157621350978},
72  {0.09157621350978, 0.81684757298044},
73  {0.81684757298044, 0.09157621350978},
74  {0.108103018168071, 0.445948490915964},
75  {0.445948490915964, 0.108103018168071},
76  {0.445948490915964, 0.445948490915964}}};
77  }
78 
79  if constexpr (n == 4) {
80  return output_t{{{0.055564052669793, 0.055564052669793},
81  {0.055564052669793, 0.888871894660413},
82  {0.888871894660413, 0.055564052669793},
83  {0.070255540518384, 0.634210747745723},
84  {0.634210747745723, 0.070255540518384},
85  {0.634210747745723, 0.295533711735893},
86  {0.070255540518384, 0.295533711735893},
87  {0.295533711735893, 0.070255540518384},
88  {0.295533711735893, 0.634210747745723},
89  {0.333333333333333, 0.333333333333333}}};
90  }
91 
92  if constexpr (n == 5) {
93  return output_t{{
94  {0.035870877695734000590,0.035870877695734000590},
95  {0.035870877695734000590,0.92825824460853301190},
96  {0.92825824460853301190,0.035870877695734000590},
97  {0.24172939576796700911,0.24172939576796700911},
98  {0.24172939576796700911,0.51654120846406603729},
99  {0.51654120846406603729,0.24172939576796700911},
100  {0.051382424445842997396,0.47430878777707902172},
101  {0.47430878777707902172,0.051382424445842997396},
102  {0.47430878777707902172,0.47430878777707902172},
103  {0.047312487011716003460,0.75118363110648400660},
104  {0.75118363110648400660,0.047312487011716003460},
105  {0.75118363110648400660,0.20150388188179998994},
106  {0.047312487011716003460,0.20150388188179998994},
107  {0.20150388188179998994,0.047312487011716003460},
108  {0.20150388188179998994,0.75118363110648400660}
109  }};
110  }
111 
112 
113  }
114 
115  if constexpr (geom == mfem::Geometry::TETRAHEDRON) {
116  using output_t = tensor<double, (n * (n + 1) * (n + 2)) / 6, 3>;
117  if constexpr (n == 1) {
118  return output_t{{{0.25, 0.25, 0.25}}};
119  }
120 
121  if constexpr (n == 2) {
122  return output_t{{{0.585410196624967936, 0.138196601125011008, 0.138196601125011008},
123  {0.138196601125011008, 0.585410196624967936, 0.138196601125011008},
124  {0.138196601125011008, 0.138196601125011008, 0.585410196624967936},
125  {0.138196601125011008, 0.138196601125011008, 0.138196601125011008}}};
126  }
127 
128  if constexpr (n == 3) {
129  return output_t{{{0.77849529482132992, 0.0738349017262233984, 0.0738349017262233984},
130  {0.0738349017262233984, 0.77849529482132992, 0.0738349017262233984},
131  {0.0738349017262233984, 0.0738349017262233984, 0.77849529482132992},
132  {0.0738349017262233984, 0.0738349017262233984, 0.0738349017262233984},
133  {0.406244343884051008, 0.406244343884051008, 0.0937556561159491072},
134  {0.406244343884051008, 0.0937556561159491072, 0.406244343884051008},
135  {0.406244343884051008, 0.0937556561159491072, 0.0937556561159491072},
136  {0.0937556561159491072, 0.406244343884051008, 0.406244343884051008},
137  {0.0937556561159491072, 0.406244343884051008, 0.0937556561159491072},
138  {0.0937556561159491072, 0.0937556561159491072, 0.406244343884051008}}};
139  }
140 
141  if constexpr (n == 4) {
142  return output_t{{{0.902942215818267904, 0.0323525947272438976, 0.0323525947272438976},
143  {0.0323525947272438976, 0.902942215818267904, 0.0323525947272438976},
144  {0.0323525947272438976, 0.0323525947272438976, 0.902942215818267904},
145  {0.0323525947272438976, 0.0323525947272438976, 0.0323525947272438976},
146  {0.262682583887778976, 0.616596533061936896, 0.0603604415251420928},
147  {0.616596533061936896, 0.262682583887778976, 0.0603604415251420928},
148  {0.262682583887778976, 0.0603604415251420928, 0.616596533061936896},
149  {0.616596533061936896, 0.0603604415251420928, 0.262682583887778976},
150  {0.262682583887778976, 0.0603604415251420928, 0.0603604415251420928},
151  {0.616596533061936896, 0.0603604415251420928, 0.0603604415251420928},
152  {0.0603604415251420928, 0.262682583887778976, 0.616596533061936896},
153  {0.0603604415251420928, 0.616596533061936896, 0.262682583887778976},
154  {0.0603604415251420928, 0.262682583887778976, 0.0603604415251420928},
155  {0.0603604415251420928, 0.616596533061936896, 0.0603604415251420928},
156  {0.0603604415251420928, 0.0603604415251420928, 0.262682583887778976},
157  {0.0603604415251420928, 0.0603604415251420928, 0.616596533061936896},
158  {0.309769304272861952, 0.309769304272861952, 0.309769304272861952},
159  {0.309769304272861952, 0.309769304272861952, 0.0706920871814129024},
160  {0.309769304272861952, 0.0706920871814129024, 0.309769304272861952},
161  {0.0706920871814129024, 0.309769304272861952, 0.309769304272861952}}};
162  }
163 
164  if constexpr (n == 5) {
165  return output_t{{{0.91978967333688, 0.0267367755543734976, 0.0267367755543734976},
166  {0.0267367755543734976, 0.91978967333688, 0.0267367755543734976},
167  {0.0267367755543734976, 0.0267367755543734976, 0.91978967333688},
168  {0.0267367755543734976, 0.0267367755543734976, 0.0267367755543734976},
169  {0.174035630246894016, 0.747759888481809024, 0.0391022406356488},
170  {0.747759888481809024, 0.174035630246894016, 0.0391022406356488},
171  {0.174035630246894016, 0.0391022406356488, 0.747759888481809024},
172  {0.747759888481809024, 0.0391022406356488, 0.174035630246894016},
173  {0.174035630246894016, 0.0391022406356488, 0.0391022406356488},
174  {0.747759888481809024, 0.0391022406356488, 0.0391022406356488},
175  {0.0391022406356488, 0.174035630246894016, 0.747759888481809024},
176  {0.0391022406356488, 0.747759888481809024, 0.174035630246894016},
177  {0.0391022406356488, 0.174035630246894016, 0.0391022406356488},
178  {0.0391022406356488, 0.747759888481809024, 0.0391022406356488},
179  {0.0391022406356488, 0.0391022406356488, 0.174035630246894016},
180  {0.0391022406356488, 0.0391022406356488, 0.747759888481809024},
181  {0.454754599984483008, 0.454754599984483008, 0.0452454000155171968},
182  {0.454754599984483008, 0.0452454000155171968, 0.454754599984483008},
183  {0.454754599984483008, 0.0452454000155171968, 0.0452454000155171968},
184  {0.0452454000155171968, 0.454754599984483008, 0.454754599984483008},
185  {0.0452454000155171968, 0.454754599984483008, 0.0452454000155171968},
186  {0.0452454000155171968, 0.0452454000155171968, 0.454754599984483008},
187  {0.503118645014598016, 0.223201037962315008, 0.223201037962315008},
188  {0.223201037962315008, 0.503118645014598016, 0.223201037962315008},
189  {0.223201037962315008, 0.223201037962315008, 0.503118645014598016},
190  {0.503118645014598016, 0.223201037962315008, 0.050479279060772},
191  {0.223201037962315008, 0.503118645014598016, 0.050479279060772},
192  {0.223201037962315008, 0.223201037962315008, 0.050479279060772},
193  {0.503118645014598016, 0.050479279060772, 0.223201037962315008},
194  {0.223201037962315008, 0.050479279060772, 0.503118645014598016},
195  {0.223201037962315008, 0.050479279060772, 0.223201037962315008},
196  {0.050479279060772, 0.503118645014598016, 0.223201037962315008},
197  {0.050479279060772, 0.223201037962315008, 0.503118645014598016},
198  {0.050479279060772, 0.223201037962315008, 0.223201037962315008},
199  {0.25, 0.25, 0.25}}};
200  }
201 
202  if constexpr (n == 6) {
203  return output_t{{{0.955143804540822016, 0.0149520651530592, 0.0149520651530592},
204  {0.0149520651530592, 0.955143804540822016, 0.0149520651530592},
205  {0.0149520651530592, 0.0149520651530592, 0.955143804540822016},
206  {0.0149520651530592, 0.0149520651530592, 0.0149520651530592},
207  {0.779976008441539968, 0.151831949165936992, 0.0340960211962614976},
208  {0.151831949165936992, 0.779976008441539968, 0.0340960211962614976},
209  {0.779976008441539968, 0.0340960211962614976, 0.151831949165936992},
210  {0.151831949165936992, 0.0340960211962614976, 0.779976008441539968},
211  {0.779976008441539968, 0.0340960211962614976, 0.0340960211962614976},
212  {0.151831949165936992, 0.0340960211962614976, 0.0340960211962614976},
213  {0.0340960211962614976, 0.779976008441539968, 0.151831949165936992},
214  {0.0340960211962614976, 0.151831949165936992, 0.779976008441539968},
215  {0.0340960211962614976, 0.779976008441539968, 0.0340960211962614976},
216  {0.0340960211962614976, 0.151831949165936992, 0.0340960211962614976},
217  {0.0340960211962614976, 0.0340960211962614976, 0.779976008441539968},
218  {0.0340960211962614976, 0.0340960211962614976, 0.151831949165936992},
219  {0.354934056063979008, 0.552655643106017024, 0.0462051504150017024},
220  {0.552655643106017024, 0.354934056063979008, 0.0462051504150017024},
221  {0.354934056063979008, 0.0462051504150017024, 0.552655643106017024},
222  {0.552655643106017024, 0.0462051504150017024, 0.354934056063979008},
223  {0.354934056063979008, 0.0462051504150017024, 0.0462051504150017024},
224  {0.552655643106017024, 0.0462051504150017024, 0.0462051504150017024},
225  {0.0462051504150017024, 0.354934056063979008, 0.552655643106017024},
226  {0.0462051504150017024, 0.552655643106017024, 0.354934056063979008},
227  {0.0462051504150017024, 0.354934056063979008, 0.0462051504150017024},
228  {0.0462051504150017024, 0.552655643106017024, 0.0462051504150017024},
229  {0.0462051504150017024, 0.0462051504150017024, 0.354934056063979008},
230  {0.0462051504150017024, 0.0462051504150017024, 0.552655643106017024},
231  {0.538104322888002048, 0.228190461068760992, 0.228190461068760992},
232  {0.228190461068760992, 0.538104322888002048, 0.228190461068760992},
233  {0.228190461068760992, 0.228190461068760992, 0.538104322888002048},
234  {0.538104322888002048, 0.228190461068760992, 0.00551475497447749952},
235  {0.228190461068760992, 0.538104322888002048, 0.00551475497447749952},
236  {0.228190461068760992, 0.228190461068760992, 0.00551475497447749952},
237  {0.538104322888002048, 0.00551475497447749952, 0.228190461068760992},
238  {0.228190461068760992, 0.00551475497447749952, 0.538104322888002048},
239  {0.228190461068760992, 0.00551475497447749952, 0.228190461068760992},
240  {0.00551475497447749952, 0.538104322888002048, 0.228190461068760992},
241  {0.00551475497447749952, 0.228190461068760992, 0.538104322888002048},
242  {0.00551475497447749952, 0.228190461068760992, 0.228190461068760992},
243  {0.19618375957456, 0.352305260087993984, 0.352305260087993984},
244  {0.352305260087993984, 0.19618375957456, 0.352305260087993984},
245  {0.352305260087993984, 0.352305260087993984, 0.19618375957456},
246  {0.19618375957456, 0.352305260087993984, 0.0992057202494530048},
247  {0.352305260087993984, 0.19618375957456, 0.0992057202494530048},
248  {0.352305260087993984, 0.352305260087993984, 0.0992057202494530048},
249  {0.19618375957456, 0.0992057202494530048, 0.352305260087993984},
250  {0.352305260087993984, 0.0992057202494530048, 0.19618375957456},
251  {0.352305260087993984, 0.0992057202494530048, 0.352305260087993984},
252  {0.0992057202494530048, 0.19618375957456, 0.352305260087993984},
253  {0.0992057202494530048, 0.352305260087993984, 0.19618375957456},
254  {0.0992057202494530048, 0.352305260087993984, 0.352305260087993984},
255  {0.59656499562101696, 0.134478334792994016, 0.134478334792994016},
256  {0.134478334792994016, 0.59656499562101696, 0.134478334792994016},
257  {0.134478334792994016, 0.134478334792994016, 0.59656499562101696},
258  {0.134478334792994016, 0.134478334792994016, 0.134478334792994016}}};
259  }
260  }
261 };
262 
271 template <int n, mfem::Geometry::Type geom>
272 SERAC_HOST_DEVICE constexpr auto GaussLegendreWeights()
273 {
274  if constexpr (geom == mfem::Geometry::SEGMENT) {
275  if constexpr (n == 1) return tensor<double, n>{1.000000000000000};
276  if constexpr (n == 2) return tensor<double, n>{0.500000000000000, 0.500000000000000};
277  if constexpr (n == 3) return tensor<double, n>{0.277777777777778, 0.444444444444444, 0.277777777777778};
278  if constexpr (n == 4) return tensor<double, n>{0.173927422568727, 0.326072577431273, 0.326072577431273, 0.173927422568727};
279  if constexpr (n == 5) return tensor<double, n>{0.118463442528095, 0.239314335249683, 0.284444444444444, 0.239314335249683, 0.118463442528095};
280  if constexpr (n == 6) return tensor<double, n>{0.085662246189585, 0.180380786524069, 0.233956967286346, 0.233956967286346, 0.180380786524069, 0.085662246189585};
281  if constexpr (n == 7) return tensor<double, n>{0.0647424830844348, 0.139852695744638, 0.190915025252559, 0.208979591836735, 0.190915025252559, 0.139852695744638, 0.0647424830844348};
282  if constexpr (n == 8) return tensor<double, n>{0.0506142681451881, 0.111190517226687, 0.156853322938944, 0.181341891689181, 0.181341891689181, 0.156853322938944, 0.111190517226687, 0.0506142681451881};
283  }
284 
285  if constexpr (geom == mfem::Geometry::TRIANGLE) {
286  using output_t = tensor< double, n * (n + 1) / 2 >;
287  if constexpr (n == 1) return output_t{0.5};
288  if constexpr (n == 2) return output_t{0.166666666666667, 0.166666666666667, 0.166666666666667};
289  if constexpr (n == 3) {
290  return output_t{
291  0.0549758718276665, 0.0549758718276665, 0.0549758718276665,
292  0.111690794839, 0.111690794839, 0.111690794839
293  };
294  }
295  if constexpr (n == 4) {
296  return output_t{
297  0.0209777564983245, 0.0209777564983245, 0.0209777564983245,
298  0.0560492060354435, 0.0560492060354435, 0.0560492060354435, 0.0560492060354435, 0.0560492060354435, 0.0560492060354435,
299  0.100771494292365
300  };
301  }
302 
303  if constexpr (n == 5) {
304  return output_t{
305  0.0089577275061514995830, 0.0089577275061514995830, 0.0089577275061514995830,
306  0.063856097940632502996, 0.063856097940632502996, 0.063856097940632502996,
307  0.038103031192767498891, 0.038103031192767498891, 0.038103031192767498891,
308  0.027874905013557500777, 0.027874905013557500777, 0.027874905013557500777,
309  0.027874905013557500777, 0.027874905013557500777, 0.027874905013557500777
310  };
311  }
312  }
313 
314  if constexpr (geom == mfem::Geometry::TETRAHEDRON) {
315  using output_t = tensor< double, (n * (n + 1) * (n + 2)) / 6 >;
316 
317  if constexpr (n == 1) {
318  return output_t{0.166666666666666656};
319  }
320  if constexpr (n == 2) {
321  return output_t{0.0416666666666666624, 0.0416666666666666624, 0.0416666666666666624, 0.0416666666666666624};
322  }
323  if constexpr (n == 3) {
324  return output_t{0.00793885580720148224, 0.00793885580720148224, 0.00793885580720148224, 0.00793885580720148224,
325  0.022485207239643504, 0.022485207239643504, 0.022485207239643504, 0.022485207239643504,
326  0.022485207239643504, 0.022485207239643504};
327  }
328  if constexpr (n == 4) {
329  return output_t{0.00117784579907825008, 0.00117784579907825008, 0.00117784579907825008, 0.00117784579907825008,
330  0.0078331114953146176, 0.0078331114953146176, 0.0078331114953146176, 0.0078331114953146176,
331  0.0078331114953146176, 0.0078331114953146176, 0.0078331114953146176, 0.0078331114953146176,
332  0.0078331114953146176, 0.0078331114953146176, 0.0078331114953146176, 0.0078331114953146176,
333  0.0169894863816446688, 0.0169894863816446688, 0.0169894863816446688, 0.0169894863816446688};
334  }
335 
336  if constexpr (n == 5) {
337  return output_t{
338  0.000365007732756466624, 0.000365007732756466624, 0.000365007732756466624, 0.000365007732756466624,
339  0.00238992783629441664, 0.00238992783629441664, 0.00238992783629441664, 0.00238992783629441664,
340  0.00238992783629441664, 0.00238992783629441664, 0.00238992783629441664, 0.00238992783629441664,
341  0.00238992783629441664, 0.00238992783629441664, 0.00238992783629441664, 0.00238992783629441664,
342  0.0041717565947791008, 0.0041717565947791008, 0.0041717565947791008, 0.0041717565947791008,
343  0.0041717565947791008, 0.0041717565947791008, 0.00799732221762589952, 0.00799732221762589952,
344  0.00799732221762589952, 0.00799732221762589952, 0.00799732221762589952, 0.00799732221762589952,
345  0.00799732221762589952, 0.00799732221762589952, 0.00799732221762589952, 0.00799732221762589952,
346  0.00799732221762589952, 0.00799732221762589952, 0.0155290955199223328};
347  }
348  if constexpr (n == 6) {
349  return output_t{
350  0.000172885205602333344, 0.000172885205602333344, 0.000172885205602333344, 0.000172885205602333344,
351  0.00160027742332466656, 0.00160027742332466656, 0.00160027742332466656, 0.00160027742332466656,
352  0.00160027742332466656, 0.00160027742332466656, 0.00160027742332466656, 0.00160027742332466656,
353  0.00160027742332466656, 0.00160027742332466656, 0.00160027742332466656, 0.00160027742332466656,
354  0.00274156627997053312, 0.00274156627997053312, 0.00274156627997053312, 0.00274156627997053312,
355  0.00274156627997053312, 0.00274156627997053312, 0.00274156627997053312, 0.00274156627997053312,
356  0.00274156627997053312, 0.00274156627997053312, 0.00274156627997053312, 0.00274156627997053312,
357  0.00256246277522183328, 0.00256246277522183328, 0.00256246277522183328, 0.00256246277522183328,
358  0.00256246277522183328, 0.00256246277522183328, 0.00256246277522183328, 0.00256246277522183328,
359  0.00256246277522183328, 0.00256246277522183328, 0.00256246277522183328, 0.00256246277522183328,
360  0.00489200197292049984, 0.00489200197292049984, 0.00489200197292049984, 0.00489200197292049984,
361  0.00489200197292049984, 0.00489200197292049984, 0.00489200197292049984, 0.00489200197292049984,
362  0.00489200197292049984, 0.00489200197292049984, 0.00489200197292049984, 0.00489200197292049984,
363  0.00610485610675179904, 0.00610485610675179904, 0.00610485610675179904, 0.00610485610675179904};
364  }
365  }
366 
367 };
368 // clang-format on
369 
374 SERAC_HOST_DEVICE constexpr int factorial(int n)
375 {
376  int nfactorial = 1;
377  for (int i = 2; i <= n; i++) {
378  nfactorial *= i;
379  }
380  return nfactorial;
381 }
382 
388 template <int n, typename T>
389 SERAC_HOST_DEVICE constexpr tensor<T, n> powers(T x)
390 {
391  tensor<T, n> values{};
392  values[0] = T(1.0);
393  for (int i = 1; i < n; i++) {
394  values[i] = values[i - 1] * x;
395  }
396  return values;
397 }
398 
406 template <int n, typename S>
407 SERAC_HOST_DEVICE constexpr tensor<S, n> ChebyshevT(S x)
408 {
409  tensor<S, n> T{};
410 
411  if (0 < n) T[0] = 1.0;
412  if (1 < n) T[1] = x;
413  for (int i = 2; i < n; i++) {
414  T[i] = 2 * x * T[i - 1] - T[i - 2];
415  }
416 
417  return T;
418 }
419 
427 template <int n, typename T>
428 SERAC_HOST_DEVICE tensor<T, n> ChebyshevU(T x)
429 {
430  tensor<T, n> U{};
431 
432  if (0 < n) U[0] = 1.0;
433  if (1 < n) U[1] = 2 * x;
434  for (int i = 2; i < n; i++) {
435  U[i] = 2 * x * U[i - 1] - U[i - 2];
436  }
437 
438  return U;
439 }
440 
447 template <int n, typename T>
448 SERAC_HOST_DEVICE tensor<T, n> Legendre(T x)
449 {
450  tensor<T, n> P{};
451 
452  if (0 < n) P[0] = 1.0;
453  if (1 < n) P[1] = x;
454  for (int i = 2; i < n; i++) {
455  P[i] = ((2 * i - 1) * x * P[i - 1] - (i - 1) * P[i - 2]) / T(i);
456  }
457 
458  return P;
459 }
460 
466 template <int n, typename T>
467 SERAC_HOST_DEVICE tensor<T, n> Bernstein(T s)
468 {
469  tensor<T, n> B;
470 
471  T t = 1.0 - s;
472 
473  T f = 1.0;
474  for (int i = 0; i < n; i++) {
475  B[i] = f;
476  f *= s / (i + 1);
477  }
478 
479  f = factorial(n - 1);
480  for (int i = 0; i < n; i++) {
481  B[n - i - 1] *= f;
482  f *= t / (i + 1);
483  }
484 
485  return B;
486 }
487 
493 template <int n, typename T>
494 SERAC_HOST_DEVICE constexpr tensor<T, n> GaussLobattoInterpolation([[maybe_unused]] T x)
495 {
496  static_assert(1 <= n && n <= 4, "error: invalid polynomial order in GaussLobattoInterpolation");
497  if constexpr (n == 1) {
498  return {1.0};
499  }
500  if constexpr (n == 2) {
501  return {1.0 - x, x};
502  }
503  if constexpr (n == 3) {
504  return {(-1.0 + x) * (-1.0 + 2.0 * x), -4.0 * (-1.0 + x) * x, x * (-1.0 + 2.0 * x)};
505  }
506  if constexpr (n == 4) {
507  constexpr double sqrt5 = 2.23606797749978981;
508  return {-(-1.0 + x) * (1.0 + 5.0 * (-1.0 + x) * x), -0.5 * sqrt5 * (5.0 + sqrt5 - 10.0 * x) * (-1.0 + x) * x,
509  -0.5 * sqrt5 * (-1.0 + x) * x * (-5.0 + sqrt5 + 10.0 * x), x * (1.0 + 5.0 * (-1.0 + x) * x)};
510  }
511 
512  return tensor<T, n>{};
513 }
514 
520 template <int n, typename T>
521 SERAC_HOST_DEVICE constexpr tensor<T, n> GaussLobattoInterpolationDerivative([[maybe_unused]] T x)
522 {
523  static_assert(1 <= n && n <= 4, "error: invalid polynomial order in GaussLobattoInterpolationDerivative");
524  if constexpr (n == 1) {
525  return {0.0};
526  }
527  if constexpr (n == 2) {
528  return {-1, 1};
529  }
530  if constexpr (n == 3) {
531  return {-3.0 + 4.0 * x, 4.0 - 8.0 * x, -1.0 + 4.0 * x};
532  }
533  if constexpr (n == 4) {
534  constexpr double sqrt5 = 2.23606797749978981;
535  return {-6.0 + 5.0 * (4.0 - 3.0 * x) * x, 2.5 * (1.0 + sqrt5 + 2.0 * x * (-1.0 - 3.0 * sqrt5 + 3.0 * sqrt5 * x)),
536  -2.5 * (-1.0 + sqrt5 + 2.0 * x * (1.0 - 3.0 * sqrt5 + 3.0 * sqrt5 * x)), 1.0 + 5.0 * x * (-2.0 + 3.0 * x)};
537  }
538 
539  return tensor<T, n>{};
540 }
541 
557 template <int n, typename T>
558 SERAC_HOST_DEVICE constexpr tensor<T, n> GaussLegendreInterpolation([[maybe_unused]] T x)
559 {
560  if constexpr (n == 1) return {1};
561  if constexpr (n == 2)
562  return {1.3660254037844386467637231708 - 1.732050807568877293527446342 * x,
563  -0.3660254037844386467637231708 + 1.732050807568877293527446342 * x};
564  if constexpr (n == 3)
565  return {1.47883055770123614752987757 + x * (-4.6243277820691389617264218 + 3.33333333333333333333333333 * x),
566  -0.666666666666666666666666667 + (6.66666666666666666666666667 - 6.666666666666666666666666667 * x) * x,
567  0.1878361089654305191367891 + x * (-2.04233888459752770494024487 + 3.33333333333333333333333333 * x)};
568  if constexpr (n == 4)
569  return {
570  1.5267881254572667869843283 +
571  x * (-8.546023607872198765973019 + (14.325858354171888152966621 - 7.420540068038946105200642 * x) * x),
572  -0.8136324494869272605618981 +
573  x * (13.807166925689577066158695 + x * (-31.388222363446060212058231 + 18.795449407555060811261716 * x)),
574  0.400761520311650404800281777 + x * (-7.41707042146263907582738061 +
575  (24.9981258592191222217269164 - 18.79544940755506081126171563 * x) * x),
576  -0.113917196281989931222711973 +
577  x * (2.15592710364526077564170438 + x * (-7.9357618499449501626353065 + 7.4205400680389461052006424 * x))};
578  return tensor<T, n>{};
579 }
580 
594 template <int n, typename T>
595 SERAC_HOST_DEVICE constexpr tensor<T, n> GaussLegendreInterpolationDerivative([[maybe_unused]] T x)
596 {
597  if constexpr (n == 1) return {0};
598  if constexpr (n == 2) return {-1.7320508075688772935274463415, 1.7320508075688772935274463415};
599  if constexpr (n == 3)
600  return {-4.6243277820691389617264218 + 6.6666666666666666666666667 * x,
601  6.66666666666666666666666667 - 13.3333333333333333333333333 * x,
602  -2.04233888459752770494024487 + 6.66666666666666666666666667 * x};
603  if constexpr (n == 4)
604  return {-8.5460236078721987659730185 + (28.651716708343776305933241 - 22.261620204116838315601927 * x) * x,
605  13.8071669256895770661586947 + x * (-62.776444726892120424116461 + 56.386348222665182433785147 * x),
606  -7.4170704214626390758273806 + (49.996251718438244443453833 - 56.386348222665182433785147 * x) * x,
607  2.15592710364526077564170438 + x * (-15.871523699889900325270613 + 22.2616202041168383156019272 * x)};
608  return tensor<T, n>{};
609 }
610 
611 } // namespace serac
SERAC_HOST_DEVICE
#define SERAC_HOST_DEVICE
Macro that evaluates to __host__ __device__ when compiling with nvcc and does nothing on a host compi...
Definition: accelerator.hpp:38
serac
Accelerator functionality.
Definition: serac.cpp:38
serac::Bernstein
SERAC_HOST_DEVICE tensor< T, n > Bernstein(T s)
Bernstein Polynomials on the domain [0, 1].
Definition: polynomials.hpp:467
serac::GaussLegendreInterpolationDerivative
constexpr SERAC_HOST_DEVICE tensor< T, n > GaussLegendreInterpolationDerivative([[maybe_unused]] T x)
Derivatives of the Lagrange Interpolating polynomials for nodes at Gauss-Legendre points on the inter...
Definition: polynomials.hpp:595
serac::GaussLobattoNodes
constexpr SERAC_HOST_DEVICE tensor< T, n > GaussLobattoNodes(T a=T(0), T b=T(1))
The positions (in 1D space) of Gauss-Lobatto points.
Definition: polynomials.hpp:28
serac::powers
constexpr SERAC_HOST_DEVICE tensor< T, n > powers(T x)
compute the first n powers of x
Definition: polynomials.hpp:389
serac::GaussLegendreWeights
constexpr SERAC_HOST_DEVICE auto GaussLegendreWeights()
The weights associated with each Gauss-Legendre point.
Definition: polynomials.hpp:272
serac::ChebyshevU
SERAC_HOST_DEVICE tensor< T, n > ChebyshevU(T x)
Chebyshev polynomials of the second kind Satisfying: sin(t) U_n(cos(t)) == sin((n+1)*t)
Definition: polynomials.hpp:428
serac::GaussLobattoInterpolationDerivative
constexpr SERAC_HOST_DEVICE tensor< T, n > GaussLobattoInterpolationDerivative([[maybe_unused]] T x)
Derivatives of the Lagrange Interpolating polynomials for nodes at Gauss-Lobatto points on the interv...
Definition: polynomials.hpp:521
serac::ChebyshevT
constexpr SERAC_HOST_DEVICE tensor< S, n > ChebyshevT(S x)
Chebyshev polynomials of the first kind Satisfying: T_n(cos(t)) == cos(n*t)
Definition: polynomials.hpp:407
serac::Legendre
SERAC_HOST_DEVICE tensor< T, n > Legendre(T x)
Legendre Polynomials, orthogonal on the domain (-1, 1) with unit weight function.
Definition: polynomials.hpp:448
serac::GaussLobattoInterpolation
constexpr SERAC_HOST_DEVICE tensor< T, n > GaussLobattoInterpolation([[maybe_unused]] T x)
Lagrange Interpolating polynomials for nodes at Gauss-Lobatto points on the interval [0,...
Definition: polynomials.hpp:494
serac::factorial
constexpr SERAC_HOST_DEVICE int factorial(int n)
compute n!
Definition: polynomials.hpp:374
serac::GaussLegendreInterpolation
constexpr SERAC_HOST_DEVICE tensor< T, n > GaussLegendreInterpolation([[maybe_unused]] T x)
Lagrange Interpolating polynomials for nodes at Gauss-Legendre points on the interval [0,...
Definition: polynomials.hpp:558
serac::GaussLegendreNodes
constexpr SERAC_HOST_DEVICE auto GaussLegendreNodes()
The positions of Gauss-Legendre points for different geometries.
Definition: polynomials.hpp:43
serac::tensor
Arbitrary-rank tensor class.
Definition: tensor.hpp:29
tensor.hpp
Implementation of the tensor class used by Functional.
Generated by doxygen 1.9.1