DNDSR 0.1.0.dev1+gcd065ad
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 int dFactorials[ndiff + 1][ndiff + 1] = {
193 {1, 0, 0, 0},
194 {1, 1, 0, 0},
195 {1, 2, 2, 0},
196 {1, 3, 6, 6}};
197
198 static const std::array<int, ndiff * 3 + 1> factorials = {
199 1,
200 1,
201 1 * 2,
202 1 * 2 * 3,
203 1 * 2 * 3 * 4,
204 1 * 2 * 3 * 4 * 5,
205 1 * 2 * 3 * 4 * 5 * 6,
206 1 * 2 * 3 * 4 * 5 * 6 * 7,
207 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8,
208 1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9,
209 };
210 static const std::array<int, ndiffSiz2D> diffNCombs2D{
211 1, 1, 1, 1, 2, 1, 1, 3, 3, 1};
212
213 static const std::array<int, ndiffSiz> diffNCombs{
214 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 3, 3, 3, 3, 3, 3, 6};
215
216 inline real iPow(int p, real x)
217 {
218 switch (p)
219 {
220 case 0:
221 return 1.;
222 case 1:
223 return x;
224 case 2:
225 return x * x;
226 case 3:
227 return x * x * x;
228 default:
229 return 1e300;
230 break;
231 }
232 }
233
234 template <int dim, int order>
235 inline constexpr int PolynomialNDOF() // 2-d specific
236 {
237 switch (dim)
238 {
239 case 2:
240 {
241 switch (order)
242 {
243 case 0:
244 return 1;
245 case 1:
246 return 3;
247 case 2:
248 return 6;
249 case 3:
250 return 10;
251 default:
252 return -1;
253 }
254 }
255 break;
256 case 3:
257 {
258 switch (order)
259 {
260 case 0:
261 return 1;
262 case 1:
263 return 4;
264 case 2:
265 return 10;
266 case 3:
267 return 20;
268 default:
269 return -1;
270 }
271 }
272 break;
273 default:
274 return -1;
275 }
276 }
277
278 inline static real FPolynomial3D(int px, int py, int pz, int dx, int dy, int dz, real x, real y, real z)
279 {
280 int c = dFactorials[px][dx] * dFactorials[py][dy] * dFactorials[pz][dz];
281 return c ? c * iPow(px - dx, x) * iPow(py - dy, y) * iPow(pz - dz, z) : 0.;
282 // return c ? c * std::pow(x, px - dx) * std::pow(y, py - dy) * std::pow(z, pz - dz) : 0.;
283 }
284
285 template <class TDIBJ>
286 void FPolynomialFill2D(TDIBJ &T, real x, real y, real z, real lx, real ly, real lz, int rows, int cols)
287 {
288 using namespace Geom::Base;
289 T.setZero();
290 if (rows == 10 && cols == 10)
291 {
292 T(0, 0) = 1.0;
293 T(0, 1) = x;
294 T(0, 2) = y;
295 T(0, 3) = x * x;
296 T(0, 4) = x * y;
297 T(0, 5) = y * y;
298 T(0, 6) = x * x * x;
299 T(0, 7) = (x * x) * y;
300 T(0, 8) = x * (y * y);
301 T(0, 9) = y * y * y;
302 T(1, 1) = 1.0 / lx;
303 T(1, 3) = (x * 2.0) / lx;
304 T(1, 4) = y / lx;
305 T(1, 6) = ((x * x) * 3.0) / lx;
306 T(1, 7) = (x * y * 2.0) / lx;
307 T(1, 8) = (y * y) / lx;
308 T(2, 2) = 1.0 / ly;
309 T(2, 4) = x / ly;
310 T(2, 5) = (y * 2.0) / ly;
311 T(2, 7) = (x * x) / ly;
312 T(2, 8) = (x * y * 2.0) / ly;
313 T(2, 9) = ((y * y) * 3.0) / ly;
314 T(3, 3) = 1.0 / (lx * lx) * 2.0;
315 T(3, 6) = 1.0 / (lx * lx) * x * 6.0;
316 T(3, 7) = 1.0 / (lx * lx) * y * 2.0;
317 T(4, 4) = 1.0 / (lx * ly);
318 T(4, 7) = (x * 2.0) / (lx * ly);
319 T(4, 8) = (y * 2.0) / (lx * ly);
320 T(5, 5) = 1.0 / (ly * ly) * 2.0;
321 T(5, 8) = 1.0 / (ly * ly) * x * 2.0;
322 T(5, 9) = 1.0 / (ly * ly) * y * 6.0;
323 T(6, 6) = 1.0 / (lx * lx * lx) * 6.0;
324 T(7, 7) = (1.0 / (lx * lx) * 2.0) / ly;
325 T(8, 8) = (1.0 / (ly * ly) * 2.0) / lx;
326 T(9, 9) = 1.0 / (ly * ly * ly) * 6.0;
327 }
328 else if (rows == 3 && cols == 10)
329 {
330 T(0, 0) = 1.0;
331 T(0, 1) = x;
332 T(0, 2) = y;
333 T(0, 3) = x * x;
334 T(0, 4) = x * y;
335 T(0, 5) = y * y;
336 T(0, 6) = x * x * x;
337 T(0, 7) = (x * x) * y;
338 T(0, 8) = x * (y * y);
339 T(0, 9) = y * y * y;
340 T(1, 1) = 1.0 / lx;
341 T(1, 3) = (x * 2.0) / lx;
342 T(1, 4) = y / lx;
343 T(1, 6) = ((x * x) * 3.0) / lx;
344 T(1, 7) = (x * y * 2.0) / lx;
345 T(1, 8) = (y * y) / lx;
346 T(2, 2) = 1.0 / ly;
347 T(2, 4) = x / ly;
348 T(2, 5) = (y * 2.0) / ly;
349 T(2, 7) = (x * x) / ly;
350 T(2, 8) = (x * y * 2.0) / ly;
351 T(2, 9) = ((y * y) * 3.0) / ly;
352 }
353 else if (rows == 1 && cols == 10)
354 {
355 T(0, 0) = 1.0;
356 T(0, 1) = x;
357 T(0, 2) = y;
358 T(0, 3) = x * x;
359 T(0, 4) = x * y;
360 T(0, 5) = y * y;
361 T(0, 6) = x * x * x;
362 T(0, 7) = (x * x) * y;
363 T(0, 8) = x * (y * y);
364 T(0, 9) = y * y * y;
365 }
366 else if (rows == 6 && cols == 6)
367 {
368 T(0, 0) = 1.0;
369 T(0, 1) = x;
370 T(0, 2) = y;
371 T(0, 3) = x * x;
372 T(0, 4) = x * y;
373 T(0, 5) = y * y;
374 T(1, 1) = 1.0 / lx;
375 T(1, 3) = (x * 2.0) / lx;
376 T(1, 4) = y / lx;
377 T(2, 2) = 1.0 / ly;
378 T(2, 4) = x / ly;
379 T(2, 5) = (y * 2.0) / ly;
380 T(3, 3) = 1.0 / (lx * lx) * 2.0;
381 T(4, 4) = 1.0 / (lx * ly);
382 T(5, 5) = 1.0 / (ly * ly) * 2.0;
383 }
384 else if (rows == 3 && cols == 6)
385 {
386 T(0, 0) = 1.0;
387 T(0, 1) = x;
388 T(0, 2) = y;
389 T(0, 3) = x * x;
390 T(0, 4) = x * y;
391 T(0, 5) = y * y;
392 T(1, 1) = 1.0 / lx;
393 T(1, 3) = (x * 2.0) / lx;
394 T(1, 4) = y / lx;
395 T(2, 2) = 1.0 / ly;
396 T(2, 4) = x / ly;
397 T(2, 5) = (y * 2.0) / ly;
398 }
399 else if (rows == 1 && cols == 6)
400 {
401 T(0, 0) = 1.0;
402 T(0, 1) = x;
403 T(0, 2) = y;
404 T(0, 3) = x * x;
405 T(0, 4) = x * y;
406 T(0, 5) = y * y;
407 }
408 else if (rows == 3 && cols == 3)
409 {
410 T(0, 0) = 1.0;
411 T(0, 1) = x;
412 T(0, 2) = y;
413 T(1, 1) = 1.0 / lx;
414 T(2, 2) = 1.0 / ly;
415 }
416 else if (rows == 1 && cols == 3)
417 {
418 T(0, 0) = 1.0;
419 T(0, 1) = x;
420 T(0, 2) = y;
421 }
422 else
423 {
424 for (int idiff = 0; idiff < rows; idiff++)
425 for (int ibase = 0; ibase < cols; ibase++)
426 {
427
428 int px = diffOperatorOrderList2D[ibase][0];
429 int py = diffOperatorOrderList2D[ibase][1];
430 int pz = diffOperatorOrderList2D[ibase][2];
431 int ndx = diffOperatorOrderList2D[idiff][0];
432 int ndy = diffOperatorOrderList2D[idiff][1];
433 int ndz = diffOperatorOrderList2D[idiff][2];
434 T(idiff, ibase) =
435 FPolynomial3D(px, py, pz, ndx, ndy, ndz,
436 x, y, z / 1.) /
437 (iPow(ndx, lx) * iPow(ndy, ly) * iPow(ndz, 1.0));
438 }
439 }
440 }
441
442 template <class TDIBJ>
443 void FPolynomialFill3D(TDIBJ &T, real x, real y, real z, real lx, real ly, real lz, int rows, int cols)
444 {
445 T.setZero();
446 if (rows == 20 && cols == 20)
447 {
448 T(0, 0) = 1.0;
449 T(0, 1) = x;
450 T(0, 2) = y;
451 T(0, 3) = z;
452 T(0, 4) = x * x;
453 T(0, 5) = y * y;
454 T(0, 6) = z * z;
455 T(0, 7) = x * y;
456 T(0, 8) = y * z;
457 T(0, 9) = x * z;
458 T(0, 10) = x * x * x;
459 T(0, 11) = y * y * y;
460 T(0, 12) = z * z * z;
461 T(0, 13) = (x * x) * y;
462 T(0, 14) = x * (y * y);
463 T(0, 15) = (y * y) * z;
464 T(0, 16) = y * (z * z);
465 T(0, 17) = x * (z * z);
466 T(0, 18) = (x * x) * z;
467 T(0, 19) = x * y * z;
468 T(1, 1) = 1.0 / lx;
469 T(1, 4) = (x * 2.0) / lx;
470 T(1, 7) = y / lx;
471 T(1, 9) = z / lx;
472 T(1, 10) = ((x * x) * 3.0) / lx;
473 T(1, 13) = (x * y * 2.0) / lx;
474 T(1, 14) = (y * y) / lx;
475 T(1, 17) = (z * z) / lx;
476 T(1, 18) = (x * z * 2.0) / lx;
477 T(1, 19) = (y * z) / lx;
478 T(2, 2) = 1.0 / ly;
479 T(2, 5) = (y * 2.0) / ly;
480 T(2, 7) = x / ly;
481 T(2, 8) = z / ly;
482 T(2, 11) = ((y * y) * 3.0) / ly;
483 T(2, 13) = (x * x) / ly;
484 T(2, 14) = (x * y * 2.0) / ly;
485 T(2, 15) = (y * z * 2.0) / ly;
486 T(2, 16) = (z * z) / ly;
487 T(2, 19) = (x * z) / ly;
488 T(3, 3) = 1.0 / lz;
489 T(3, 6) = (z * 2.0) / lz;
490 T(3, 8) = y / lz;
491 T(3, 9) = x / lz;
492 T(3, 12) = ((z * z) * 3.0) / lz;
493 T(3, 15) = (y * y) / lz;
494 T(3, 16) = (y * z * 2.0) / lz;
495 T(3, 17) = (x * z * 2.0) / lz;
496 T(3, 18) = (x * x) / lz;
497 T(3, 19) = (x * y) / lz;
498 T(4, 4) = 1.0 / (lx * lx) * 2.0;
499 T(4, 10) = 1.0 / (lx * lx) * x * 6.0;
500 T(4, 13) = 1.0 / (lx * lx) * y * 2.0;
501 T(4, 18) = 1.0 / (lx * lx) * z * 2.0;
502 T(5, 5) = 1.0 / (ly * ly) * 2.0;
503 T(5, 11) = 1.0 / (ly * ly) * y * 6.0;
504 T(5, 14) = 1.0 / (ly * ly) * x * 2.0;
505 T(5, 15) = 1.0 / (ly * ly) * z * 2.0;
506 T(6, 6) = 1.0 / (lz * lz) * 2.0;
507 T(6, 12) = 1.0 / (lz * lz) * z * 6.0;
508 T(6, 16) = 1.0 / (lz * lz) * y * 2.0;
509 T(6, 17) = 1.0 / (lz * lz) * x * 2.0;
510 T(7, 7) = 1.0 / (lx * ly);
511 T(7, 13) = (x * 2.0) / (lx * ly);
512 T(7, 14) = (y * 2.0) / (lx * ly);
513 T(7, 19) = z / (lx * ly);
514 T(8, 8) = 1.0 / (ly * lz);
515 T(8, 15) = (y * 2.0) / (ly * lz);
516 T(8, 16) = (z * 2.0) / (ly * lz);
517 T(8, 19) = x / (ly * lz);
518 T(9, 9) = 1.0 / (lx * lz);
519 T(9, 17) = (z * 2.0) / (lx * lz);
520 T(9, 18) = (x * 2.0) / (lx * lz);
521 T(9, 19) = y / (lx * lz);
522 T(10, 10) = 1.0 / (lx * lx * lx) * 6.0;
523 T(11, 11) = 1.0 / (ly * ly * ly) * 6.0;
524 T(12, 12) = 1.0 / (lz * lz * lz) * 6.0;
525 T(13, 13) = (1.0 / (lx * lx) * 2.0) / ly;
526 T(14, 14) = (1.0 / (ly * ly) * 2.0) / lx;
527 T(15, 15) = (1.0 / (ly * ly) * 2.0) / lz;
528 T(16, 16) = (1.0 / (lz * lz) * 2.0) / ly;
529 T(17, 17) = (1.0 / (lz * lz) * 2.0) / lx;
530 T(18, 18) = (1.0 / (lx * lx) * 2.0) / lz;
531 T(19, 19) = 1.0 / (lx * ly * lz);
532 }
533 else if (rows == 4 && cols == 20)
534 {
535 T(0, 0) = 1.0;
536 T(0, 1) = x;
537 T(0, 2) = y;
538 T(0, 3) = z;
539 T(0, 4) = x * x;
540 T(0, 5) = y * y;
541 T(0, 6) = z * z;
542 T(0, 7) = x * y;
543 T(0, 8) = y * z;
544 T(0, 9) = x * z;
545 T(0, 10) = x * x * x;
546 T(0, 11) = y * y * y;
547 T(0, 12) = z * z * z;
548 T(0, 13) = (x * x) * y;
549 T(0, 14) = x * (y * y);
550 T(0, 15) = (y * y) * z;
551 T(0, 16) = y * (z * z);
552 T(0, 17) = x * (z * z);
553 T(0, 18) = (x * x) * z;
554 T(0, 19) = x * y * z;
555 T(1, 1) = 1.0 / lx;
556 T(1, 4) = (x * 2.0) / lx;
557 T(1, 7) = y / lx;
558 T(1, 9) = z / lx;
559 T(1, 10) = ((x * x) * 3.0) / lx;
560 T(1, 13) = (x * y * 2.0) / lx;
561 T(1, 14) = (y * y) / lx;
562 T(1, 17) = (z * z) / lx;
563 T(1, 18) = (x * z * 2.0) / lx;
564 T(1, 19) = (y * z) / lx;
565 T(2, 2) = 1.0 / ly;
566 T(2, 5) = (y * 2.0) / ly;
567 T(2, 7) = x / ly;
568 T(2, 8) = z / ly;
569 T(2, 11) = ((y * y) * 3.0) / ly;
570 T(2, 13) = (x * x) / ly;
571 T(2, 14) = (x * y * 2.0) / ly;
572 T(2, 15) = (y * z * 2.0) / ly;
573 T(2, 16) = (z * z) / ly;
574 T(2, 19) = (x * z) / ly;
575 T(3, 3) = 1.0 / lz;
576 T(3, 6) = (z * 2.0) / lz;
577 T(3, 8) = y / lz;
578 T(3, 9) = x / lz;
579 T(3, 12) = ((z * z) * 3.0) / lz;
580 T(3, 15) = (y * y) / lz;
581 T(3, 16) = (y * z * 2.0) / lz;
582 T(3, 17) = (x * z * 2.0) / lz;
583 T(3, 18) = (x * x) / lz;
584 T(3, 19) = (x * y) / lz;
585 }
586 else if (rows == 1 && cols == 20)
587 {
588 T(0, 0) = 1.0;
589 T(0, 1) = x;
590 T(0, 2) = y;
591 T(0, 3) = z;
592 T(0, 4) = x * x;
593 T(0, 5) = y * y;
594 T(0, 6) = z * z;
595 T(0, 7) = x * y;
596 T(0, 8) = y * z;
597 T(0, 9) = x * z;
598 T(0, 10) = x * x * x;
599 T(0, 11) = y * y * y;
600 T(0, 12) = z * z * z;
601 T(0, 13) = (x * x) * y;
602 T(0, 14) = x * (y * y);
603 T(0, 15) = (y * y) * z;
604 T(0, 16) = y * (z * z);
605 T(0, 17) = x * (z * z);
606 T(0, 18) = (x * x) * z;
607 T(0, 19) = x * y * z;
608 }
609 else if (rows == 10 && cols == 10)
610 {
611 T(0, 0) = 1.0;
612 T(0, 1) = x;
613 T(0, 2) = y;
614 T(0, 3) = z;
615 T(0, 4) = x * x;
616 T(0, 5) = y * y;
617 T(0, 6) = z * z;
618 T(0, 7) = x * y;
619 T(0, 8) = y * z;
620 T(0, 9) = x * z;
621 T(1, 1) = 1.0 / lx;
622 T(1, 4) = (x * 2.0) / lx;
623 T(1, 7) = y / lx;
624 T(1, 9) = z / lx;
625 T(2, 2) = 1.0 / ly;
626 T(2, 5) = (y * 2.0) / ly;
627 T(2, 7) = x / ly;
628 T(2, 8) = z / ly;
629 T(3, 3) = 1.0 / lz;
630 T(3, 6) = (z * 2.0) / lz;
631 T(3, 8) = y / lz;
632 T(3, 9) = x / lz;
633 T(4, 4) = 1.0 / (lx * lx) * 2.0;
634 T(5, 5) = 1.0 / (ly * ly) * 2.0;
635 T(6, 6) = 1.0 / (lz * lz) * 2.0;
636 T(7, 7) = 1.0 / (lx * ly);
637 T(8, 8) = 1.0 / (ly * lz);
638 T(9, 9) = 1.0 / (lx * lz);
639 }
640 else if (rows == 4 && cols == 10)
641 {
642 T(0, 0) = 1.0;
643 T(0, 1) = x;
644 T(0, 2) = y;
645 T(0, 3) = z;
646 T(0, 4) = x * x;
647 T(0, 5) = y * y;
648 T(0, 6) = z * z;
649 T(0, 7) = x * y;
650 T(0, 8) = y * z;
651 T(0, 9) = x * z;
652 T(1, 1) = 1.0 / lx;
653 T(1, 4) = (x * 2.0) / lx;
654 T(1, 7) = y / lx;
655 T(1, 9) = z / lx;
656 T(2, 2) = 1.0 / ly;
657 T(2, 5) = (y * 2.0) / ly;
658 T(2, 7) = x / ly;
659 T(2, 8) = z / ly;
660 T(3, 3) = 1.0 / lz;
661 T(3, 6) = (z * 2.0) / lz;
662 T(3, 8) = y / lz;
663 T(3, 9) = x / lz;
664 }
665 else if (rows == 1 && cols == 10)
666 {
667 T(0, 0) = 1.0;
668 T(0, 1) = x;
669 T(0, 2) = y;
670 T(0, 3) = z;
671 T(0, 4) = x * x;
672 T(0, 5) = y * y;
673 T(0, 6) = z * z;
674 T(0, 7) = x * y;
675 T(0, 8) = y * z;
676 T(0, 9) = x * z;
677 }
678 else if (rows == 4 && cols == 4)
679 {
680 T(0, 0) = 1.0;
681 T(0, 1) = x;
682 T(0, 2) = y;
683 T(0, 3) = z;
684 T(1, 1) = 1.0 / lx;
685 T(2, 2) = 1.0 / ly;
686 T(3, 3) = 1.0 / lz;
687 }
688 else if (rows == 1 && cols == 4)
689 {
690 T(0, 0) = 1.0;
691 T(0, 1) = x;
692 T(0, 2) = y;
693 T(0, 3) = z;
694 }
695 else
696 {
697 for (int idiff = 0; idiff < rows; idiff++)
698 for (int ibase = 0; ibase < cols; ibase++)
699 {
700 int px = diffOperatorOrderList[ibase][0];
701 int py = diffOperatorOrderList[ibase][1];
702 int pz = diffOperatorOrderList[ibase][2];
703 int ndx = diffOperatorOrderList[idiff][0];
704 int ndy = diffOperatorOrderList[idiff][1];
705 int ndz = diffOperatorOrderList[idiff][2];
706 T(idiff, ibase) =
707 FPolynomial3D(px, py, pz, ndx, ndy, ndz,
708 x, y, z) /
709 (iPow(ndx, lx) * iPow(ndy, ly) * iPow(ndz, lz));
710 }
711 }
712 }
713
714 // #include <unsupported/Eigen/CXX11/TensorSymmetry>
715 template <int dim>
716 inline int GetNDof(int maxOrder)
717 {
718 int maxNDOF = -1;
719 switch (maxOrder)
720 {
721 case 0:
722 maxNDOF = PolynomialNDOF<dim, 0>();
723 break;
724 case 1:
725 maxNDOF = PolynomialNDOF<dim, 1>();
726 break;
727 case 2:
728 maxNDOF = PolynomialNDOF<dim, 2>();
729 break;
730 case 3:
731 maxNDOF = PolynomialNDOF<dim, 3>();
732 break;
733 default:
734 {
735 DNDS_assert_info(false, "maxNDOF invalid");
736 }
737 }
738 return maxNDOF;
739 }
740
741}
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:113
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:216
DNDS::Geom::Base::GetNDof
int GetNDof(int maxOrder)
Definition BaseFunction.hpp:716
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:443
DNDS::Geom::Base::getDiffOperatorIJK2I
constexpr auto & getDiffOperatorIJK2I()
Definition BaseFunction.hpp:187
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::FPolynomialFill2D
void FPolynomialFill2D(TDIBJ &T, real x, real y, real z, real lx, real ly, real lz, int rows, int cols)
Definition BaseFunction.hpp:286
DNDS::Geom::Base::PolynomialNDOF
constexpr int PolynomialNDOF()
Definition BaseFunction.hpp:235
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
z
tVec z(NCells)
x
tVec x(NCells)
b
tVec b(NCells)
order
double order
Definition test_ODE.cpp:257