DNDSR 0.2.1
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
BaseFunction.hpp
Go to the documentation of this file.
1#pragma once
2
3#include "DNDS/Defines.hpp"
4#include "Geom/EigenTensor.hpp"
5#include "Geom/PeriodicInfo.hpp"
6#include <array>
7#include <tuple>
8
9namespace DNDS::Geom::Base
10{
11 /// including up to 3 orders or diffs
12 static const int ndiff = 3;
13 static const int ndiffSiz = 20;
14 static const int ndiffSiz2D = 10;
15 static const std::array<int, ndiff + 1> ndiffSizC{1, 4, 10, 20};
16 static const std::array<int, ndiff + 1> ndiffSizC2D{1, 3, 6, 10};
17 template <int dim>
18 constexpr int ndiffSizS(int i)
19 {
20 return dim == 2 ? ndiffSizC2D[i] : ndiffSizC[i];
21 }
22 template <int dim>
23 constexpr int ndiff2order(int rows)
24 {
25 if constexpr (dim == 2)
26 switch (rows)
27 {
28 case 10:
29 return 3;
30 case 6:
31 return 2;
32 case 3:
33 return 1;
34 case 1:
35 return 0;
36 default:
37 return -1;
38 }
39 else // dim ==3
40 switch (rows)
41 {
42 case 20:
43 return 3;
44 case 10:
45 return 2;
46 case 4:
47 return 1;
48 case 1:
49 return 0;
50 default:
51 return -1;
52 }
53 return -1;
54 }
55 static const std::array<std::array<int, 3>, ndiffSiz> diffOperatorOrderList{{
56 //{diffOrderX_0, diffOrderX_1, diffOrder_X2} // indexPlace, diffSeq
57 {{0, 0, 0}}, // 00
58 {{1, 0, 0}}, // 01
59 {{0, 1, 0}}, // 02
60 {{0, 0, 1}}, // 03
61 {{2, 0, 0}}, // 04 00
62 {{0, 2, 0}}, // 05 11
63 {{0, 0, 2}}, // 06 22
64 {{1, 1, 0}}, // 07 01
65 {{0, 1, 1}}, // 08 12
66 {{1, 0, 1}}, // 09 02
67 {{3, 0, 0}}, // 10 000
68 {{0, 3, 0}}, // 11 111
69 {{0, 0, 3}}, // 12 222
70 {{2, 1, 0}}, // 13 001
71 {{1, 2, 0}}, // 14 011
72 {{0, 2, 1}}, // 15 112
73 {{0, 1, 2}}, // 16 122
74 {{1, 0, 2}}, // 17 022
75 {{2, 0, 1}}, // 18 002
76 {{1, 1, 1}}, // 19 012
77 }};
78
79 static const std::array<std::array<int, 3>, ndiffSiz2D> diffOperatorOrderList2D =
80 {{
81 {{0, 0, 0}}, // 00 00
82 {{1, 0, 0}}, // 01 01 0
83 {{0, 1, 0}}, // 02 02 1
84 {{2, 0, 0}}, // 03 04 00
85 {{1, 1, 0}}, // 04 05 01
86 {{0, 2, 0}}, // 05 07 11
87 {{3, 0, 0}}, // 06 10 000
88 {{2, 1, 0}}, // 07 11 001
89 {{1, 2, 0}}, // 08 13 011
90 {{0, 3, 0}}, // 09 16 111
91 }};
92
93 static const std::array<std::array<int, ndiff>, ndiffSiz> diffOperatorDimList{{
94 //{diffOrderX_0, diffOrderX_1, diffOrder_X2} // indexPlace, diffSeq
95 {{}}, // 00
96 {{0}}, // 01
97 {{1}}, // 02
98 {{2}}, // 03
99 {{0, 0}}, // 04 00
100 {{1, 1}}, // 05 11
101 {{2, 2}}, // 06 22
102 {{0, 1}}, // 07 01
103 {{1, 2}}, // 08 12
104 {{0, 2}}, // 09 02
105 {{0, 0, 0}}, // 10 000
106 {{1, 1, 1}}, // 11 111
107 {{2, 2, 2}}, // 12 222
108 {{0, 0, 1}}, // 13 001
109 {{0, 1, 1}}, // 14 011
110 {{1, 1, 2}}, // 15 112
111 {{1, 2, 2}}, // 16 122
112 {{0, 2, 2}}, // 17 022
113 {{0, 0, 2}}, // 18 002
114 {{0, 1, 2}}, // 19 012
115 }};
116
117 static const std::array<std::array<int, ndiff>, ndiffSiz2D> diffOperatorDimList2D =
118 {{
119 {{}}, // 00 00
120 {{0}}, // 01 01 0
121 {{1}}, // 02 02 1
122 {{0, 0}}, // 03 04 00
123 {{0, 1}}, // 04 05 01
124 {{1, 1}}, // 05 07 11
125 {{0, 0, 0}}, // 06 10 000
126 {{0, 0, 1}}, // 07 11 001
127 {{0, 1, 1}}, // 08 13 011
128 {{1, 1, 1}}, // 09 16 111
129 }};
130
131 using t_diffOpIJK2I = std::tuple<
132 int,
133 std::array<int, 3>,
134 std::array<std::array<int, 3>, 3>,
135 std::array<std::array<std::array<int, 3>, 3>, 3>>;
136 template <int dim, int NDiffC>
137 constexpr t_diffOpIJK2I _get_diffOperatorIJK2I(const std::array<std::array<int, 3>, NDiffC> &diffOps)
138 {
139 auto ret = t_diffOpIJK2I();
140 std::get<0>(ret) = 0;
141 auto array3IsSame = [](const std::array<int, 3> &a, const std::array<int, 3> &b)
142 {
143 return (a[0] == b[0]) && (a[1] == b[1]) && (a[2] == b[2]);
144 };
145 auto searchForArray3 = [=](const std::array<int, 3> &a)
146 {
147 int ret = -1;
148 int found = 0;
149 for (int i = 0; i < NDiffC; i++)
150 if (array3IsSame(a, diffOps[i]))
151 ret = i, found++;
152 return ret;
153 };
154
155 for (int d0 = 0; d0 < dim; d0++)
156 {
157 std::array<int, 3> entry{0, 0, 0};
158 entry[d0]++;
159 std::get<1>(ret)[d0] = searchForArray3(entry);
160 }
161 for (int d1 = 0; d1 < dim; d1++)
162 for (int d0 = 0; d0 < dim; d0++)
163 {
164 std::array<int, 3> entry{0, 0, 0};
165 entry[d0]++;
166 entry[d1]++;
167 std::get<2>(ret)[d0][d1] = searchForArray3(entry);
168 }
169 for (int d2 = 0; d2 < dim; d2++)
170 for (int d1 = 0; d1 < dim; d1++)
171 for (int d0 = 0; d0 < dim; d0++)
172 {
173 std::array<int, 3> entry{0, 0, 0};
174 entry[d0]++;
175 entry[d1]++;
176 entry[d2]++;
177 std::get<3>(ret)[d0][d1][d2] = searchForArray3(entry);
178 }
179 return ret;
180 }
181
182 static const t_diffOpIJK2I diffOperatorIJK2I = _get_diffOperatorIJK2I<3, ndiffSiz>(diffOperatorOrderList);
183
184 static const t_diffOpIJK2I diffOperatorIJK2I2D = _get_diffOperatorIJK2I<2, ndiffSiz2D>(diffOperatorOrderList2D);
185
186 template <int dim>
187 constexpr auto &getDiffOperatorIJK2I()
188 {
189 return (dim == 2) ? diffOperatorIJK2I2D : diffOperatorIJK2I;
190 }
191
192 static const std::array<std::array<int, ndiff + 1>, ndiff + 1> dFactorials = {{{1, 0, 0, 0},
193 {1, 1, 0, 0},
194 {1, 2, 2, 0},
195 {1, 3, 6, 6}}};
196
197 static const std::array<int, ndiff * 3 + 1> factorials = {
198 1,
199 1,
200 1 * 2,
201 1 * 2 * 3,
202 1 * 2 * 3 * 4,
203 1 * 2 * 3 * 4 * 5,
204 1 * 2 * 3 * 4 * 5 * 6,
205 1 * 2 * 3 * 4 * 5 * 6 * 7,
206 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8,
207 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9,
208 };
209 static const std::array<int, ndiffSiz2D> diffNCombs2D{
210 1, 1, 1, 1, 2, 1, 1, 3, 3, 1};
211
212 static const std::array<int, ndiffSiz> diffNCombs{
213 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 3, 3, 3, 3, 3, 3, 6};
214
215 inline real iPow(int p, real x)
216 {
217 switch (p)
218 {
219 case 0:
220 return 1.;
221 case 1:
222 return x;
223 case 2:
224 return x * x;
225 case 3:
226 return x * x * x;
227 default:
228 return 1e300;
229 break;
230 }
231 }
232
233 template <int dim, int order>
234 constexpr int PolynomialNDOF() // 2-d specific
235 {
236 switch (dim)
237 {
238 case 2:
239 {
240 switch (order)
241 {
242 case 0:
243 return 1;
244 case 1:
245 return 3;
246 case 2:
247 return 6;
248 case 3:
249 return 10;
250 default:
251 return -1;
252 }
253 }
254 break;
255 case 3:
256 {
257 switch (order)
258 {
259 case 0:
260 return 1;
261 case 1:
262 return 4;
263 case 2:
264 return 10;
265 case 3:
266 return 20;
267 default:
268 return -1;
269 }
270 }
271 break;
272 default:
273 return -1;
274 }
275 }
276
277 inline static real FPolynomial3D(int px, int py, int pz, int dx, int dy, int dz, real x, real y, real z)
278 {
279 int c = dFactorials[px][dx] * dFactorials[py][dy] * dFactorials[pz][dz];
280 return c ? c * iPow(px - dx, x) * iPow(py - dy, y) * iPow(pz - dz, z) : 0.;
281 // return c ? c * std::pow(x, px - dx) * std::pow(y, py - dy) * std::pow(z, pz - dz) : 0.;
282 }
283
284 template <class TDIBJ>
285 void FPolynomialFill2D(TDIBJ &T, real x, real y, real z, real lx, real ly, real lz, int rows, int cols)
286 {
287 using namespace Geom::Base;
288 T.setZero();
289 if (rows == 10 && cols == 10)
290 {
291 T(0, 0) = 1.0;
292 T(0, 1) = x;
293 T(0, 2) = y;
294 T(0, 3) = x * x;
295 T(0, 4) = x * y;
296 T(0, 5) = y * y;
297 T(0, 6) = x * x * x;
298 T(0, 7) = (x * x) * y;
299 T(0, 8) = x * (y * y);
300 T(0, 9) = y * y * y;
301 T(1, 1) = 1.0 / lx;
302 T(1, 3) = (x * 2.0) / lx;
303 T(1, 4) = y / lx;
304 T(1, 6) = ((x * x) * 3.0) / lx;
305 T(1, 7) = (x * y * 2.0) / lx;
306 T(1, 8) = (y * y) / lx;
307 T(2, 2) = 1.0 / ly;
308 T(2, 4) = x / ly;
309 T(2, 5) = (y * 2.0) / ly;
310 T(2, 7) = (x * x) / ly;
311 T(2, 8) = (x * y * 2.0) / ly;
312 T(2, 9) = ((y * y) * 3.0) / ly;
313 T(3, 3) = 1.0 / (lx * lx) * 2.0;
314 T(3, 6) = 1.0 / (lx * lx) * x * 6.0;
315 T(3, 7) = 1.0 / (lx * lx) * y * 2.0;
316 T(4, 4) = 1.0 / (lx * ly);
317 T(4, 7) = (x * 2.0) / (lx * ly);
318 T(4, 8) = (y * 2.0) / (lx * ly);
319 T(5, 5) = 1.0 / (ly * ly) * 2.0;
320 T(5, 8) = 1.0 / (ly * ly) * x * 2.0;
321 T(5, 9) = 1.0 / (ly * ly) * y * 6.0;
322 T(6, 6) = 1.0 / (lx * lx * lx) * 6.0;
323 T(7, 7) = (1.0 / (lx * lx) * 2.0) / ly;
324 T(8, 8) = (1.0 / (ly * ly) * 2.0) / lx;
325 T(9, 9) = 1.0 / (ly * ly * ly) * 6.0;
326 }
327 else if (rows == 3 && cols == 10)
328 {
329 T(0, 0) = 1.0;
330 T(0, 1) = x;
331 T(0, 2) = y;
332 T(0, 3) = x * x;
333 T(0, 4) = x * y;
334 T(0, 5) = y * y;
335 T(0, 6) = x * x * x;
336 T(0, 7) = (x * x) * y;
337 T(0, 8) = x * (y * y);
338 T(0, 9) = y * y * y;
339 T(1, 1) = 1.0 / lx;
340 T(1, 3) = (x * 2.0) / lx;
341 T(1, 4) = y / lx;
342 T(1, 6) = ((x * x) * 3.0) / lx;
343 T(1, 7) = (x * y * 2.0) / lx;
344 T(1, 8) = (y * y) / lx;
345 T(2, 2) = 1.0 / ly;
346 T(2, 4) = x / ly;
347 T(2, 5) = (y * 2.0) / ly;
348 T(2, 7) = (x * x) / ly;
349 T(2, 8) = (x * y * 2.0) / ly;
350 T(2, 9) = ((y * y) * 3.0) / ly;
351 }
352 else if (rows == 1 && cols == 10)
353 {
354 T(0, 0) = 1.0;
355 T(0, 1) = x;
356 T(0, 2) = y;
357 T(0, 3) = x * x;
358 T(0, 4) = x * y;
359 T(0, 5) = y * y;
360 T(0, 6) = x * x * x;
361 T(0, 7) = (x * x) * y;
362 T(0, 8) = x * (y * y);
363 T(0, 9) = y * y * y;
364 }
365 else if (rows == 6 && cols == 6)
366 {
367 T(0, 0) = 1.0;
368 T(0, 1) = x;
369 T(0, 2) = y;
370 T(0, 3) = x * x;
371 T(0, 4) = x * y;
372 T(0, 5) = y * y;
373 T(1, 1) = 1.0 / lx;
374 T(1, 3) = (x * 2.0) / lx;
375 T(1, 4) = y / lx;
376 T(2, 2) = 1.0 / ly;
377 T(2, 4) = x / ly;
378 T(2, 5) = (y * 2.0) / ly;
379 T(3, 3) = 1.0 / (lx * lx) * 2.0;
380 T(4, 4) = 1.0 / (lx * ly);
381 T(5, 5) = 1.0 / (ly * ly) * 2.0;
382 }
383 else if (rows == 3 && cols == 6)
384 {
385 T(0, 0) = 1.0;
386 T(0, 1) = x;
387 T(0, 2) = y;
388 T(0, 3) = x * x;
389 T(0, 4) = x * y;
390 T(0, 5) = y * y;
391 T(1, 1) = 1.0 / lx;
392 T(1, 3) = (x * 2.0) / lx;
393 T(1, 4) = y / lx;
394 T(2, 2) = 1.0 / ly;
395 T(2, 4) = x / ly;
396 T(2, 5) = (y * 2.0) / ly;
397 }
398 else if (rows == 1 && cols == 6)
399 {
400 T(0, 0) = 1.0;
401 T(0, 1) = x;
402 T(0, 2) = y;
403 T(0, 3) = x * x;
404 T(0, 4) = x * y;
405 T(0, 5) = y * y;
406 }
407 else if (rows == 3 && cols == 3)
408 {
409 T(0, 0) = 1.0;
410 T(0, 1) = x;
411 T(0, 2) = y;
412 T(1, 1) = 1.0 / lx;
413 T(2, 2) = 1.0 / ly;
414 }
415 else if (rows == 1 && cols == 3)
416 {
417 T(0, 0) = 1.0;
418 T(0, 1) = x;
419 T(0, 2) = y;
420 }
421 else
422 {
423 for (int idiff = 0; idiff < rows; idiff++)
424 for (int ibase = 0; ibase < cols; ibase++)
425 {
426
427 int px = diffOperatorOrderList2D[ibase][0];
428 int py = diffOperatorOrderList2D[ibase][1];
429 int pz = diffOperatorOrderList2D[ibase][2];
430 int ndx = diffOperatorOrderList2D[idiff][0];
431 int ndy = diffOperatorOrderList2D[idiff][1];
432 int ndz = diffOperatorOrderList2D[idiff][2];
433 T(idiff, ibase) =
434 FPolynomial3D(px, py, pz, ndx, ndy, ndz,
435 x, y, z / 1.) /
436 (iPow(ndx, lx) * iPow(ndy, ly) * iPow(ndz, 1.0));
437 }
438 }
439 }
440
441 template <class TDIBJ>
442 void FPolynomialFill3D(TDIBJ &T, real x, real y, real z, real lx, real ly, real lz, int rows, int cols)
443 {
444 T.setZero();
445 if (rows == 20 && cols == 20)
446 {
447 T(0, 0) = 1.0;
448 T(0, 1) = x;
449 T(0, 2) = y;
450 T(0, 3) = z;
451 T(0, 4) = x * x;
452 T(0, 5) = y * y;
453 T(0, 6) = z * z;
454 T(0, 7) = x * y;
455 T(0, 8) = y * z;
456 T(0, 9) = x * z;
457 T(0, 10) = x * x * x;
458 T(0, 11) = y * y * y;
459 T(0, 12) = z * z * z;
460 T(0, 13) = (x * x) * y;
461 T(0, 14) = x * (y * y);
462 T(0, 15) = (y * y) * z;
463 T(0, 16) = y * (z * z);
464 T(0, 17) = x * (z * z);
465 T(0, 18) = (x * x) * z;
466 T(0, 19) = x * y * z;
467 T(1, 1) = 1.0 / lx;
468 T(1, 4) = (x * 2.0) / lx;
469 T(1, 7) = y / lx;
470 T(1, 9) = z / lx;
471 T(1, 10) = ((x * x) * 3.0) / lx;
472 T(1, 13) = (x * y * 2.0) / lx;
473 T(1, 14) = (y * y) / lx;
474 T(1, 17) = (z * z) / lx;
475 T(1, 18) = (x * z * 2.0) / lx;
476 T(1, 19) = (y * z) / lx;
477 T(2, 2) = 1.0 / ly;
478 T(2, 5) = (y * 2.0) / ly;
479 T(2, 7) = x / ly;
480 T(2, 8) = z / ly;
481 T(2, 11) = ((y * y) * 3.0) / ly;
482 T(2, 13) = (x * x) / ly;
483 T(2, 14) = (x * y * 2.0) / ly;
484 T(2, 15) = (y * z * 2.0) / ly;
485 T(2, 16) = (z * z) / ly;
486 T(2, 19) = (x * z) / ly;
487 T(3, 3) = 1.0 / lz;
488 T(3, 6) = (z * 2.0) / lz;
489 T(3, 8) = y / lz;
490 T(3, 9) = x / lz;
491 T(3, 12) = ((z * z) * 3.0) / lz;
492 T(3, 15) = (y * y) / lz;
493 T(3, 16) = (y * z * 2.0) / lz;
494 T(3, 17) = (x * z * 2.0) / lz;
495 T(3, 18) = (x * x) / lz;
496 T(3, 19) = (x * y) / lz;
497 T(4, 4) = 1.0 / (lx * lx) * 2.0;
498 T(4, 10) = 1.0 / (lx * lx) * x * 6.0;
499 T(4, 13) = 1.0 / (lx * lx) * y * 2.0;
500 T(4, 18) = 1.0 / (lx * lx) * z * 2.0;
501 T(5, 5) = 1.0 / (ly * ly) * 2.0;
502 T(5, 11) = 1.0 / (ly * ly) * y * 6.0;
503 T(5, 14) = 1.0 / (ly * ly) * x * 2.0;
504 T(5, 15) = 1.0 / (ly * ly) * z * 2.0;
505 T(6, 6) = 1.0 / (lz * lz) * 2.0;
506 T(6, 12) = 1.0 / (lz * lz) * z * 6.0;
507 T(6, 16) = 1.0 / (lz * lz) * y * 2.0;
508 T(6, 17) = 1.0 / (lz * lz) * x * 2.0;
509 T(7, 7) = 1.0 / (lx * ly);
510 T(7, 13) = (x * 2.0) / (lx * ly);
511 T(7, 14) = (y * 2.0) / (lx * ly);
512 T(7, 19) = z / (lx * ly);
513 T(8, 8) = 1.0 / (ly * lz);
514 T(8, 15) = (y * 2.0) / (ly * lz);
515 T(8, 16) = (z * 2.0) / (ly * lz);
516 T(8, 19) = x / (ly * lz);
517 T(9, 9) = 1.0 / (lx * lz);
518 T(9, 17) = (z * 2.0) / (lx * lz);
519 T(9, 18) = (x * 2.0) / (lx * lz);
520 T(9, 19) = y / (lx * lz);
521 T(10, 10) = 1.0 / (lx * lx * lx) * 6.0;
522 T(11, 11) = 1.0 / (ly * ly * ly) * 6.0;
523 T(12, 12) = 1.0 / (lz * lz * lz) * 6.0;
524 T(13, 13) = (1.0 / (lx * lx) * 2.0) / ly;
525 T(14, 14) = (1.0 / (ly * ly) * 2.0) / lx;
526 T(15, 15) = (1.0 / (ly * ly) * 2.0) / lz;
527 T(16, 16) = (1.0 / (lz * lz) * 2.0) / ly;
528 T(17, 17) = (1.0 / (lz * lz) * 2.0) / lx;
529 T(18, 18) = (1.0 / (lx * lx) * 2.0) / lz;
530 T(19, 19) = 1.0 / (lx * ly * lz);
531 }
532 else if (rows == 4 && cols == 20)
533 {
534 T(0, 0) = 1.0;
535 T(0, 1) = x;
536 T(0, 2) = y;
537 T(0, 3) = z;
538 T(0, 4) = x * x;
539 T(0, 5) = y * y;
540 T(0, 6) = z * z;
541 T(0, 7) = x * y;
542 T(0, 8) = y * z;
543 T(0, 9) = x * z;
544 T(0, 10) = x * x * x;
545 T(0, 11) = y * y * y;
546 T(0, 12) = z * z * z;
547 T(0, 13) = (x * x) * y;
548 T(0, 14) = x * (y * y);
549 T(0, 15) = (y * y) * z;
550 T(0, 16) = y * (z * z);
551 T(0, 17) = x * (z * z);
552 T(0, 18) = (x * x) * z;
553 T(0, 19) = x * y * z;
554 T(1, 1) = 1.0 / lx;
555 T(1, 4) = (x * 2.0) / lx;
556 T(1, 7) = y / lx;
557 T(1, 9) = z / lx;
558 T(1, 10) = ((x * x) * 3.0) / lx;
559 T(1, 13) = (x * y * 2.0) / lx;
560 T(1, 14) = (y * y) / lx;
561 T(1, 17) = (z * z) / lx;
562 T(1, 18) = (x * z * 2.0) / lx;
563 T(1, 19) = (y * z) / lx;
564 T(2, 2) = 1.0 / ly;
565 T(2, 5) = (y * 2.0) / ly;
566 T(2, 7) = x / ly;
567 T(2, 8) = z / ly;
568 T(2, 11) = ((y * y) * 3.0) / ly;
569 T(2, 13) = (x * x) / ly;
570 T(2, 14) = (x * y * 2.0) / ly;
571 T(2, 15) = (y * z * 2.0) / ly;
572 T(2, 16) = (z * z) / ly;
573 T(2, 19) = (x * z) / ly;
574 T(3, 3) = 1.0 / lz;
575 T(3, 6) = (z * 2.0) / lz;
576 T(3, 8) = y / lz;
577 T(3, 9) = x / lz;
578 T(3, 12) = ((z * z) * 3.0) / lz;
579 T(3, 15) = (y * y) / lz;
580 T(3, 16) = (y * z * 2.0) / lz;
581 T(3, 17) = (x * z * 2.0) / lz;
582 T(3, 18) = (x * x) / lz;
583 T(3, 19) = (x * y) / lz;
584 }
585 else if (rows == 1 && cols == 20)
586 {
587 T(0, 0) = 1.0;
588 T(0, 1) = x;
589 T(0, 2) = y;
590 T(0, 3) = z;
591 T(0, 4) = x * x;
592 T(0, 5) = y * y;
593 T(0, 6) = z * z;
594 T(0, 7) = x * y;
595 T(0, 8) = y * z;
596 T(0, 9) = x * z;
597 T(0, 10) = x * x * x;
598 T(0, 11) = y * y * y;
599 T(0, 12) = z * z * z;
600 T(0, 13) = (x * x) * y;
601 T(0, 14) = x * (y * y);
602 T(0, 15) = (y * y) * z;
603 T(0, 16) = y * (z * z);
604 T(0, 17) = x * (z * z);
605 T(0, 18) = (x * x) * z;
606 T(0, 19) = x * y * z;
607 }
608 else if (rows == 10 && cols == 10)
609 {
610 T(0, 0) = 1.0;
611 T(0, 1) = x;
612 T(0, 2) = y;
613 T(0, 3) = z;
614 T(0, 4) = x * x;
615 T(0, 5) = y * y;
616 T(0, 6) = z * z;
617 T(0, 7) = x * y;
618 T(0, 8) = y * z;
619 T(0, 9) = x * z;
620 T(1, 1) = 1.0 / lx;
621 T(1, 4) = (x * 2.0) / lx;
622 T(1, 7) = y / lx;
623 T(1, 9) = z / lx;
624 T(2, 2) = 1.0 / ly;
625 T(2, 5) = (y * 2.0) / ly;
626 T(2, 7) = x / ly;
627 T(2, 8) = z / ly;
628 T(3, 3) = 1.0 / lz;
629 T(3, 6) = (z * 2.0) / lz;
630 T(3, 8) = y / lz;
631 T(3, 9) = x / lz;
632 T(4, 4) = 1.0 / (lx * lx) * 2.0;
633 T(5, 5) = 1.0 / (ly * ly) * 2.0;
634 T(6, 6) = 1.0 / (lz * lz) * 2.0;
635 T(7, 7) = 1.0 / (lx * ly);
636 T(8, 8) = 1.0 / (ly * lz);
637 T(9, 9) = 1.0 / (lx * lz);
638 }
639 else if (rows == 4 && cols == 10)
640 {
641 T(0, 0) = 1.0;
642 T(0, 1) = x;
643 T(0, 2) = y;
644 T(0, 3) = z;
645 T(0, 4) = x * x;
646 T(0, 5) = y * y;
647 T(0, 6) = z * z;
648 T(0, 7) = x * y;
649 T(0, 8) = y * z;
650 T(0, 9) = x * z;
651 T(1, 1) = 1.0 / lx;
652 T(1, 4) = (x * 2.0) / lx;
653 T(1, 7) = y / lx;
654 T(1, 9) = z / lx;
655 T(2, 2) = 1.0 / ly;
656 T(2, 5) = (y * 2.0) / ly;
657 T(2, 7) = x / ly;
658 T(2, 8) = z / ly;
659 T(3, 3) = 1.0 / lz;
660 T(3, 6) = (z * 2.0) / lz;
661 T(3, 8) = y / lz;
662 T(3, 9) = x / lz;
663 }
664 else if (rows == 1 && cols == 10)
665 {
666 T(0, 0) = 1.0;
667 T(0, 1) = x;
668 T(0, 2) = y;
669 T(0, 3) = z;
670 T(0, 4) = x * x;
671 T(0, 5) = y * y;
672 T(0, 6) = z * z;
673 T(0, 7) = x * y;
674 T(0, 8) = y * z;
675 T(0, 9) = x * z;
676 }
677 else if (rows == 4 && cols == 4)
678 {
679 T(0, 0) = 1.0;
680 T(0, 1) = x;
681 T(0, 2) = y;
682 T(0, 3) = z;
683 T(1, 1) = 1.0 / lx;
684 T(2, 2) = 1.0 / ly;
685 T(3, 3) = 1.0 / lz;
686 }
687 else if (rows == 1 && cols == 4)
688 {
689 T(0, 0) = 1.0;
690 T(0, 1) = x;
691 T(0, 2) = y;
692 T(0, 3) = z;
693 }
694 else
695 {
696 for (int idiff = 0; idiff < rows; idiff++)
697 for (int ibase = 0; ibase < cols; ibase++)
698 {
699 int px = diffOperatorOrderList[ibase][0];
700 int py = diffOperatorOrderList[ibase][1];
701 int pz = diffOperatorOrderList[ibase][2];
702 int ndx = diffOperatorOrderList[idiff][0];
703 int ndy = diffOperatorOrderList[idiff][1];
704 int ndz = diffOperatorOrderList[idiff][2];
705 T(idiff, ibase) =
706 FPolynomial3D(px, py, pz, ndx, ndy, ndz,
707 x, y, z) /
708 (iPow(ndx, lx) * iPow(ndy, ly) * iPow(ndz, lz));
709 }
710 }
711 }
712
713 // #include <unsupported/Eigen/CXX11/TensorSymmetry>
714 template <int dim>
715 inline int GetNDof(int maxOrder)
716 {
717 int maxNDOF = -1;
718 switch (maxOrder)
719 {
720 case 0:
721 maxNDOF = PolynomialNDOF<dim, 0>();
722 break;
723 case 1:
724 maxNDOF = PolynomialNDOF<dim, 1>();
725 break;
726 case 2:
727 maxNDOF = PolynomialNDOF<dim, 2>();
728 break;
729 case 3:
730 maxNDOF = PolynomialNDOF<dim, 3>();
731 break;
732 default:
733 {
734 DNDS_assert_info(false, "maxNDOF invalid");
735 }
736 }
737 return maxNDOF;
738 }
739
740}
Defines.hpp
Core type aliases, constants, and metaprogramming utilities for the DNDS framework.
DNDS_assert_info
#define DNDS_assert_info(expr, info)
Debug-only assertion with an extra std::string info message.
Definition Errors.hpp:117
DNDS::Geom::Base::ndiff2order
constexpr int ndiff2order(int rows)
Definition BaseFunction.hpp:23
DNDS::Geom::Base::t_diffOpIJK2I
std::tuple< int, std::array< int, 3 >, std::array< std::array< int, 3 >, 3 >, std::array< std::array< std::array< int, 3 >, 3 >, 3 > > t_diffOpIJK2I
Definition BaseFunction.hpp:135
DNDS::Geom::Base::ndiffSizS
constexpr int ndiffSizS(int i)
Definition BaseFunction.hpp:18
DNDS::Geom::Base::iPow
real iPow(int p, real x)
Definition BaseFunction.hpp:215
DNDS::Geom::Base::GetNDof
int GetNDof(int maxOrder)
Definition BaseFunction.hpp:715
DNDS::Geom::Base::FPolynomialFill3D
void FPolynomialFill3D(TDIBJ &T, real x, real y, real z, real lx, real ly, real lz, int rows, int cols)
Definition BaseFunction.hpp:442
DNDS::Geom::Base::getDiffOperatorIJK2I
constexpr auto & getDiffOperatorIJK2I()
Definition BaseFunction.hpp:187
DNDS::Geom::Base::FPolynomialFill2D
void FPolynomialFill2D(TDIBJ &T, real x, real y, real z, real lx, real ly, real lz, int rows, int cols)
Definition BaseFunction.hpp:285
DNDS::Geom::Base::_get_diffOperatorIJK2I
constexpr t_diffOpIJK2I _get_diffOperatorIJK2I(const std::array< std::array< int, 3 >, NDiffC > &diffOps)
Definition BaseFunction.hpp:137
DNDS::Geom::Base::PolynomialNDOF
constexpr int PolynomialNDOF()
Definition BaseFunction.hpp:234
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:110
z
tVec z(NCells)
x
tVec x(NCells)
b
tVec b(NCells)
i
res i
Definition test_MeshConnectivity_Interpolate.cpp:205
order
double order
Definition test_ODE.cpp:257
p
const tPoint const tPoint const tPoint & p
Definition test_Reconstruction.cpp:484