gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
OrthogonalPoly.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include <stdexcept>
7 #include "OrthogonalPoly.h"
8 
9 double OrthogonalPoly::EvalLobatto(int order, double x)
10 {
11  double L = 0;
12  double xsquare = pow(x, 2);
13  switch(order) {
14  case(0): L = 0.5 * (1 - 1 * x); return L;
15  case(1): L = 0.5 * (1 + x); return L;
16  case(2):
17  L = (-1 + pow(x, 2));
18  L = L * 0.5 * pow(3. / 2., 0.5);
19  return L;
20  case(3):
21  L = x * (-1 + xsquare);
22  L = L * 0.5 * pow(5. / 2., 0.5);
23  return L;
24  case(4):
25  L = 1 + xsquare * (-6 + 5. * xsquare);
26  L = L * 1. / 8. * pow(7. / 2., 0.5);
27  return L;
28  case(5):
29  L = x * (3 + xsquare * (-10 + 7. * xsquare));
30  L = L * 3. / 8. * pow(2., -0.5);
31  return L;
32  case(6):
33  L = -1 + xsquare * (15 + xsquare * (-35 + 21 * xsquare));
34  L = L * 1. / 16. * pow(11. / 2., 0.5);
35  return L;
36  case(7):
37  L = x * (-5 + xsquare * (35 + xsquare * (-63 + 33 * xsquare)));
38  L = L * 1. / 16. * pow(13. / 2., 0.5);
39  return L;
40  case(8):
41  L =
42  5 + xsquare * (-140 + xsquare * (630 + xsquare * (-924 + 429 * xsquare)));
43  L = L * 1. / 128. * pow(15. / 2., 0.5);
44  return L;
45  case(9):
46  L = x *
47  (35 + xsquare *
48  (-420 + xsquare * (1386 + xsquare * (-1716 + 715 * xsquare))));
49  L = L * 1. / 128. * pow(17. / 2., 0.5);
50  return L;
51  case(10):
52  L = -7 +
53  xsquare *
54  (315 +
55  xsquare *
56  (-2310 + xsquare * (6006 + xsquare * (-6435 + 2431 * xsquare))));
57  L = L * 1. / 256. * pow(19. / 2., 0.5);
58  return L;
59  case(11):
60  L =
61  x *
62  (-63 +
63  xsquare *
64  (1155 +
65  xsquare *
66  (-6006 + xsquare * (12870 + xsquare * (-12155 + 4199 * xsquare)))));
67  L = L * 1. / 256. * pow(21. / 2., 0.5);
68  return L;
69  case(12):
70  L = 21 +
71  xsquare *
72  (-1386 +
73  xsquare *
74  (15015 +
75  xsquare *
76  (-60060 +
77  xsquare * (109395 + xsquare * (-92378 + 29393 * xsquare)))));
78  L = L * 1. / 1024. * pow(23. / 2., 0.5);
79  return L;
80  case(13):
81  L =
82  x *
83  (231 +
84  xsquare *
85  (-6006 +
86  xsquare *
87  (45045 +
88  xsquare *
89  (-145860 +
90  xsquare * (230945 + xsquare * (-176358 + 52003 * xsquare))))));
91  L = L * 5. / 1024. * pow(2., -0.5);
92  return L;
93  case(14):
94  L =
95  -33 +
96  xsquare *
97  (3003 +
98  xsquare *
99  (-45045 +
100  xsquare *
101  (255255 +
102  xsquare * (-692835 +
103  xsquare * (969969 + xsquare * (-676039 +
104  185725 * xsquare))))));
105  L = L * 3. / 2048. * pow(3. / 2., 0.5);
106  return L;
107  case(15):
108  L = x *
109  (-429 +
110  xsquare *
111  (15015 +
112  xsquare *
113  (-153153 +
114  xsquare *
115  (692835 +
116  xsquare *
117  (-1616615 +
118  xsquare *
119  (2028117 + xsquare * (-1300075 + 334305 * xsquare)))))));
120  L = L * 1. / 2048. * pow(29. / 2., 0.5);
121  return L;
122 
123  default: throw std::runtime_error("Lobatto functions are written for orders =< 15");
124  }
125 }
126 
127 double OrthogonalPoly::EvalDLobatto(int order, double x)
128 {
129  double dL = 0;
130  double xsquare = pow(x, 2);
131  switch(order) {
132  case(0): dL = -0.5; return dL;
133  case(1): dL = 0.5; return dL;
134  case(2):
135  dL = 2 * x;
136  dL = dL * 0.5 * pow(3. / 2., 0.5);
137  return dL;
138  case(3):
139  dL = -1 + 3 * xsquare;
140  dL = dL * 0.5 * pow(5. / 2., 0.5);
141  return dL;
142  case(4):
143  dL = x * (-12 + 20 * xsquare);
144  dL = dL * 1. / 8. * pow(7. / 2., 0.5);
145  return dL;
146  case(5):
147  dL = 3 + xsquare * (-30 + 35. * xsquare);
148  dL = dL * 3. / 8. * pow(2., -0.5);
149  return dL;
150  case(6):
151  dL = x * (30 + xsquare * (-140 + 126 * xsquare));
152  dL = dL * 1. / 16. * pow(11. / 2., 0.5);
153  return dL;
154  case(7):
155  dL = -5 + xsquare * (105 + xsquare * (-315 + 231 * xsquare));
156  dL = dL * 1. / 16. * pow(13. / 2., 0.5);
157  return dL;
158  case(8):
159  dL = x * (-280 + xsquare * (2520 + xsquare * (-5544 + 3432 * xsquare)));
160  dL = dL * 1. / 128. * pow(15. / 2., 0.5);
161  return dL;
162  case(9):
163  dL =
164  35 + xsquare *
165  (-1260 + xsquare * (6930 + xsquare * (-12012 + 6435 * xsquare)));
166  dL = dL * 1. / 128. * pow(17. / 2., 0.5);
167  return dL;
168  case(10):
169  dL = x *
170  (630 +
171  xsquare *
172  (-9240 + xsquare * (36036 + xsquare * (-51480 + 24310 * xsquare))));
173  dL = dL * 1. / 256. * pow(19. / 2., 0.5);
174  return dL;
175  case(11):
176  dL =
177  -63 +
178  xsquare *
179  (3465 +
180  xsquare * (-30030 +
181  xsquare * (90090 + xsquare * (-109395 + 46189 * xsquare))));
182  dL = dL * 1. / 256. * pow(21. / 2., 0.5);
183  return dL;
184  case(12):
185  dL =
186  x * (-2772 +
187  xsquare *
188  (60060 +
189  xsquare *
190  (-360360 +
191  xsquare * (875160 + xsquare * (-923780 + 352716 * xsquare)))));
192  dL = dL * 1. / 1024. * pow(23. / 2., 0.5);
193  return dL;
194  case(13):
195  dL =
196  231 +
197  xsquare *
198  (-18018 +
199  xsquare *
200  (225225 +
201  xsquare *
202  (-1021020 +
203  xsquare * (2078505 + xsquare * (-1939938 + 676039 * xsquare)))));
204  dL = dL * 5. / 1024. * pow(2., -0.5);
205  return dL;
206  case(14):
207  dL =
208  x *
209  (6006 +
210  xsquare *
211  (-180180 +
212  xsquare *
213  (1531530 +
214  xsquare * (-5542680 +
215  xsquare * (9699690 + xsquare * (-8112468 +
216  2600150 * xsquare))))));
217  dL = dL * 3. / 2048. * pow(3. / 2., 0.5);
218  return dL;
219  case(15):
220  dL = -429 +
221  xsquare *
222  (45045 +
223  xsquare *
224  (-765765 +
225  xsquare *
226  (4849845 +
227  xsquare *
228  (-14549535 +
229  xsquare * (22309287 +
230  xsquare * (-16900975 + 5014575 * xsquare))))));
231  dL = dL * 1. / 2048. * pow(29. / 2., 0.5);
232  return dL;
233 
234  default: throw std::runtime_error("Lobatto functions are written for orders =< 15");
235  }
236 }
237 
238 double OrthogonalPoly::EvalKernelFunction(int order, double x)
239 {
240  double phi = 0;
241  double xsquare = pow(x, 2);
242  switch(order) {
243  case(0): phi = -pow(6, 0.5); return phi;
244  case(1): phi = -x * pow(10, 0.5); return phi;
245  case(2):
246  phi = 1 - 5 * xsquare;
247  phi = phi * 0.5 * pow(7. / 2., 0.5);
248  return phi;
249  case(3):
250  phi = x * (3 - 7 * xsquare);
251  phi = phi * 3. / 2. * pow(2, -0.5);
252  return phi;
253  case(4):
254  phi = -1 + xsquare * (14 - 21 * xsquare);
255  phi = phi * 1. / 4. * pow(11. / 2., 0.5);
256  return phi;
257  case(5):
258  phi = x * (-5 + xsquare * (30 - 33 * xsquare));
259  phi = phi * 1. / 4. * pow(13. / 2., 0.5);
260  return phi;
261  case(6):
262  phi = 5 + xsquare * (-135 + xsquare * (495 - 429 * xsquare));
263  phi = phi * 1. / 32. * pow(15. / 2., 0.5);
264  return phi;
265  case(7):
266  phi = x * (35 + xsquare * (-385 + xsquare * (1001 - 715 * xsquare)));
267  phi = phi * 1. / 32. * pow(17. / 2., 0.5);
268  return phi;
269  case(8):
270  phi = -7 + xsquare *
271  (308 + xsquare * (-2002 + xsquare * (4004 - 2431 * xsquare)));
272  phi = phi * 1. / 64. * pow(19. / 2., 0.5);
273  return phi;
274  case(9):
275  phi =
276  x *
277  (-63 + xsquare *
278  (1092 + xsquare * (-4914 + xsquare * (7956 - 4199 * xsquare))));
279  phi = phi * 1. / 64. * pow(21. / 2., 0.5);
280  return phi;
281  case(10):
282  phi =
283  21 +
284  xsquare *
285  (-1365 +
286  xsquare *
287  (13650 + xsquare * (-46410 + xsquare * (62985 - 29393 * xsquare))));
288  phi = phi * 1. / 256. * pow(23. / 2., 0.5);
289  return phi;
290  case(11):
291  phi =
292  x * (231 +
293  xsquare *
294  (-5775 +
295  xsquare *
296  (39270 +
297  xsquare * (-106590 + xsquare * (124355 - 52003 * xsquare)))));
298  phi = phi * 5. / 256. * pow(2, -0.5);
299  return phi;
300  case(12):
301  phi =
302  -33 +
303  xsquare *
304  (2970 +
305  xsquare *
306  (-42075 +
307  xsquare *
308  (213180 +
309  xsquare * (-479655 + xsquare * (490314 - 185725 * xsquare)))));
310  phi = phi * 3. / 512. * pow(3. / 2., 0.5);
311  return phi;
312  case(13):
313  phi =
314  x *
315  (-429 +
316  xsquare *
317  (14586 +
318  xsquare *
319  (-138567 +
320  xsquare *
321  (554268 + xsquare * (-1062347 +
322  xsquare * (965770 - 334305 * xsquare))))));
323  phi = phi * 1. / 512. * pow(29. / 2., 0.5);
324  return phi;
325  default: throw std::runtime_error("Lobatto functions are written for orders =< 15");
326  }
327 }
328 
329 double OrthogonalPoly::EvalDKernelFunction(int order, double x)
330 {
331  double dphi = 0;
332  double xsquare = pow(x, 2);
333  switch(order) {
334  case(0): dphi = 0; return dphi;
335  case(1): dphi = -pow(10, 0.5); return dphi;
336  case(2):
337  dphi = -10 * x;
338  ;
339  dphi = dphi * 0.5 * pow(7. / 2., 0.5);
340  return dphi;
341  case(3):
342  dphi = 3 - 21 * xsquare;
343  dphi = dphi * 3. / 2. * pow(2, -0.5);
344  return dphi;
345  case(4):
346  dphi = x * (28 - 84 * xsquare);
347  dphi = dphi * 1. / 4. * pow(11. / 2., 0.5);
348  return dphi;
349  case(5):
350  dphi = -5 + xsquare * (90 - 165 * xsquare);
351  dphi = dphi * 1. / 4. * pow(13. / 2., 0.5);
352  return dphi;
353  case(6):
354  dphi = x * (-270 + xsquare * (1980 - 2574 * xsquare));
355  dphi = dphi * 1. / 32. * pow(15. / 2., 0.5);
356  return dphi;
357  case(7):
358  dphi = 35 + xsquare * (-1155 + xsquare * (5005 - 5005 * xsquare));
359  dphi = dphi * 1. / 32. * pow(17. / 2., 0.5);
360  return dphi;
361  case(8):
362  dphi = x * (616 + xsquare * (-8008 + xsquare * (24024 - 19448 * xsquare)));
363  dphi = dphi * 1. / 64. * pow(19. / 2., 0.5);
364  return dphi;
365  case(9):
366  dphi =
367  -63 + xsquare *
368  (3276 + xsquare * (-24570 + xsquare * (55692 - 37791 * xsquare)));
369  dphi = dphi * 1. / 64. * pow(21. / 2., 0.5);
370  return dphi;
371  case(10):
372  dphi =
373  x *
374  (-2730 +
375  xsquare *
376  (54600 + xsquare * (-278460 + xsquare * (503880 - 293930 * xsquare))));
377 
378  dphi = dphi * 1. / 256. * pow(23. / 2., 0.5);
379  return dphi;
380  case(11):
381  dphi = 231 +
382  xsquare *
383  (-17325 +
384  xsquare *
385  (196350 +
386  xsquare * (-746130 + xsquare * (1119195 - 572033 * xsquare))));
387  dphi = dphi * 5. / 256. * pow(2, -0.5);
388  return dphi;
389  case(12):
390  dphi =
391  x * (5940 +
392  xsquare *
393  (-168300 +
394  xsquare * (1279080 +
395  xsquare * (-3837240 +
396  xsquare * (4903140 - 2228700 * xsquare)))));
397  dphi = dphi * 3. / 512. * pow(3. / 2., 0.5);
398  return dphi;
399  case(13):
400  dphi =
401  -429 +
402  xsquare *
403  (43758 +
404  xsquare *
405  (-692835 +
406  xsquare * (3879876 +
407  xsquare * (-9561123 +
408  xsquare * (10623470 - 4345965 * xsquare)))));
409  dphi = dphi * 1. / 512. * pow(29. / 2., 0.5);
410  return dphi;
411  default: throw std::runtime_error("Lobatto functions are written for orders =< 15");
412  }
413 }
414 
415 double OrthogonalPoly::EvalLegendre(int order, double x)
416 {
417  double L = 0;
418  double xsquare = pow(x, 2);
419  switch(order) {
420  case(0): L = 1; return L;
421  case(1): L = x; return L;
422  case(2): L = 3. / 2. * xsquare - 1. / 2.; return L;
423  case(3): L = 0.5 * x * (5 * xsquare - 3); return L;
424  case(4):
425  L = (3 + xsquare * (35 * xsquare - 30));
426  L = 1. / 8. * L;
427  return L;
428  case(5):
429  L = x * (xsquare * (63 * xsquare - 70) + 15);
430  L = 1. / 8. * L;
431  return L;
432  case(6):
433  L = ((231 * xsquare - 315) * xsquare + 105) * xsquare - 5;
434  L = 1. / 16. * L;
435  return L;
436  case(7):
437  L = x * (((429 * xsquare - 693) * xsquare + 315) * xsquare - 35);
438  L = 1. / 16. * L;
439  return L;
440  case(8):
441  L =
442  (((6435 * xsquare - 12012) * xsquare + 6930) * xsquare - 1260) * xsquare +
443  35;
444  L = 1. / 128. * L;
445  return L;
446  case(9):
447  L = ((((12155 * xsquare - 25740) * xsquare + 18018) * xsquare - 4620) *
448  xsquare +
449  315) *
450  x;
451  L = 1. / 128. * L;
452  return L;
453  case(10):
454  L = ((((46189 * xsquare - 109395) * xsquare + 90090) * xsquare - 30030) *
455  xsquare +
456  3465) *
457  xsquare -
458  63;
459  L = 1. / 256. * L;
460  return L;
461  default: throw std::runtime_error("Legendre functions are written for orders =< 10");
462  }
463 }
464 
465 double OrthogonalPoly::EvalDLegendre(int order, double x)
466 {
467  double dL = 0;
468  double xsquare = pow(x, 2);
469  switch(order) {
470  case(0): dL = 0; return dL;
471  case(1): dL = 1; return dL;
472  case(2): dL = 3 * x; return dL;
473  case(3): dL = 0.5 * (15 * xsquare - 3); return dL;
474  case(4):
475  dL = x * (140 * xsquare - 60);
476  dL = 1. / 8. * dL;
477  return dL;
478  case(5):
479  dL = 15 + xsquare * (315 * xsquare - 210);
480  dL = 1. / 8. * dL;
481  return dL;
482  case(6):
483  dL = x * (210 + xsquare * (1386 * xsquare - 1260));
484  dL = 1. / 16. * dL;
485  return dL;
486  case(7):
487  dL = ((xsquare * 3003 - 3465) * xsquare + 945) * xsquare - 35;
488  dL = 1. / 16. * dL;
489  return dL;
490  case(8):
491  dL = x * (((51480 * xsquare - 72072) * xsquare + 27720) * xsquare - 2520);
492  dL = 1. / 128. * dL;
493  return dL;
494  case(9):
495  dL = 315 +
496  xsquare * (-13860 +
497  xsquare * (90090 + xsquare * (-180180 + 109395 * xsquare)));
498  dL = 1. / 128. * dL;
499  return dL;
500  case(10):
501  dL = x * (6930 +
502  xsquare *
503  (-120120 +
504  xsquare * (540540 + xsquare * (-875160 + 461890 * xsquare))));
505  dL = 1. / 256. * dL;
506  return dL;
507  default: throw std::runtime_error("Legendre functions are written for orders =< 10");
508  }
509 }
OrthogonalPoly::EvalDLegendre
double EvalDLegendre(int order, double x)
Definition: OrthogonalPoly.cpp:465
OrthogonalPoly::EvalLegendre
double EvalLegendre(int order, double x)
Definition: OrthogonalPoly.cpp:415
OrthogonalPoly::EvalKernelFunction
double EvalKernelFunction(int order, double x)
Definition: OrthogonalPoly.cpp:238
OrthogonalPoly::EvalLobatto
double EvalLobatto(int order, double x)
Definition: OrthogonalPoly.cpp:9
OrthogonalPoly::EvalDLobatto
double EvalDLobatto(int order, double x)
Definition: OrthogonalPoly.cpp:127
OrthogonalPoly.h
OrthogonalPoly::EvalDKernelFunction
double EvalDKernelFunction(int order, double x)
Definition: OrthogonalPoly.cpp:329