DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
DiffTensors.hpp
Go to the documentation of this file.
1#pragma once
2
3#include "BaseFunction.hpp"
4
5namespace DNDS::Geom::Base
6{
7 template <int dim = 3, int rank = 0, int powV = 1, class VLe, class VRi>
8 real NormSymDiffOrderTensorV(VLe &&Le, VRi &&Ri)
9 {
10 real ret = 0;
11 if constexpr (dim == 3)
12 {
13 if constexpr (powV == 1)
14 {
15 if constexpr (rank == 0)
16 {
17 ret += Le(0, 0) * Ri(0, 0) * diffNCombs[0];
18 }
19 else if constexpr (rank == 1)
20 {
21 ret += Le(0, 0) * Ri(0, 0) * diffNCombs[0 + 1];
22 ret += Le(1, 0) * Ri(1, 0) * diffNCombs[1 + 1];
23 ret += Le(2, 0) * Ri(2, 0) * diffNCombs[2 + 1];
24 }
25 else if constexpr (rank == 2)
26 {
27 ret += Le(0, 0) * Ri(0, 0) * diffNCombs[0 + 4];
28 ret += Le(1, 0) * Ri(1, 0) * diffNCombs[1 + 4];
29 ret += Le(2, 0) * Ri(2, 0) * diffNCombs[2 + 4];
30 ret += Le(3, 0) * Ri(3, 0) * diffNCombs[3 + 4];
31 ret += Le(4, 0) * Ri(4, 0) * diffNCombs[4 + 4];
32 ret += Le(5, 0) * Ri(5, 0) * diffNCombs[5 + 4];
33 }
34 else if constexpr (rank == 3)
35 {
36 ret += Le(0, 0) * Ri(0, 0) * diffNCombs[0 + 10];
37 ret += Le(1, 0) * Ri(1, 0) * diffNCombs[1 + 10];
38 ret += Le(2, 0) * Ri(2, 0) * diffNCombs[2 + 10];
39 ret += Le(3, 0) * Ri(3, 0) * diffNCombs[3 + 10];
40 ret += Le(4, 0) * Ri(4, 0) * diffNCombs[4 + 10];
41 ret += Le(5, 0) * Ri(5, 0) * diffNCombs[5 + 10];
42 ret += Le(6, 0) * Ri(6, 0) * diffNCombs[6 + 10];
43 ret += Le(7, 0) * Ri(7, 0) * diffNCombs[7 + 10];
44 ret += Le(8, 0) * Ri(8, 0) * diffNCombs[8 + 10];
45 ret += Le(9, 0) * Ri(9, 0) * diffNCombs[9 + 10];
46 }
47 else
48 {
49 DNDS_assert(false);
50 }
51 }
52 else
53 {
54 if constexpr (rank == 0)
55 {
56 ret += Le(0, 0) * Ri(0, 0) * std::pow(diffNCombs[0], powV);
57 }
58 else if constexpr (rank == 1)
59 {
60 ret += Le(0, 0) * Ri(0, 0) * std::pow(diffNCombs[0 + 1], powV);
61 ret += Le(1, 0) * Ri(1, 0) * std::pow(diffNCombs[1 + 1], powV);
62 ret += Le(2, 0) * Ri(2, 0) * std::pow(diffNCombs[2 + 1], powV);
63 }
64 else if constexpr (rank == 2)
65 {
66 ret += Le(0, 0) * Ri(0, 0) * std::pow(diffNCombs[0 + 4], powV);
67 ret += Le(1, 0) * Ri(1, 0) * std::pow(diffNCombs[1 + 4], powV);
68 ret += Le(2, 0) * Ri(2, 0) * std::pow(diffNCombs[2 + 4], powV);
69 ret += Le(3, 0) * Ri(3, 0) * std::pow(diffNCombs[3 + 4], powV);
70 ret += Le(4, 0) * Ri(4, 0) * std::pow(diffNCombs[4 + 4], powV);
71 ret += Le(5, 0) * Ri(5, 0) * std::pow(diffNCombs[5 + 4], powV);
72 }
73 else if constexpr (rank == 3)
74 {
75 ret += Le(0, 0) * Ri(0, 0) * std::pow(diffNCombs[0 + 10], powV);
76 ret += Le(1, 0) * Ri(1, 0) * std::pow(diffNCombs[1 + 10], powV);
77 ret += Le(2, 0) * Ri(2, 0) * std::pow(diffNCombs[2 + 10], powV);
78 ret += Le(3, 0) * Ri(3, 0) * std::pow(diffNCombs[3 + 10], powV);
79 ret += Le(4, 0) * Ri(4, 0) * std::pow(diffNCombs[4 + 10], powV);
80 ret += Le(5, 0) * Ri(5, 0) * std::pow(diffNCombs[5 + 10], powV);
81 ret += Le(6, 0) * Ri(6, 0) * std::pow(diffNCombs[6 + 10], powV);
82 ret += Le(7, 0) * Ri(7, 0) * std::pow(diffNCombs[7 + 10], powV);
83 ret += Le(8, 0) * Ri(8, 0) * std::pow(diffNCombs[8 + 10], powV);
84 ret += Le(9, 0) * Ri(9, 0) * std::pow(diffNCombs[9 + 10], powV);
85 }
86 else
87 {
88 DNDS_assert(false);
89 }
90 }
91 }
92 else
93 {
94 if constexpr (powV == 1)
95 {
96 if constexpr (rank == 0)
97 {
98 ret += Le(0) * Ri(0) * diffNCombs2D[0];
99 }
100 else if constexpr (rank == 1)
101 {
102 ret += Le(0) * Ri(0) * diffNCombs2D[0 + 1];
103 ret += Le(1) * Ri(1) * diffNCombs2D[1 + 1];
104 }
105 else if constexpr (rank == 2)
106 {
107 ret += Le(0) * Ri(0) * diffNCombs2D[0 + 3];
108 ret += Le(1) * Ri(1) * diffNCombs2D[1 + 3];
109 ret += Le(2) * Ri(2) * diffNCombs2D[2 + 3];
110 }
111 else if constexpr (rank == 3)
112 {
113 ret += Le(0) * Ri(0) * diffNCombs2D[0 + 6];
114 ret += Le(1) * Ri(1) * diffNCombs2D[1 + 6];
115 ret += Le(2) * Ri(2) * diffNCombs2D[2 + 6];
116 ret += Le(3) * Ri(3) * diffNCombs2D[3 + 6];
117 }
118 else
119 {
120 DNDS_assert(false);
121 }
122 }
123 else
124 {
125 if constexpr (rank == 0)
126 {
127 ret += Le(0) * Ri(0) * std::pow(diffNCombs2D[0], powV);
128 }
129 else if constexpr (rank == 1)
130 {
131 ret += Le(0) * Ri(0) * std::pow(diffNCombs2D[0 + 1], powV);
132 ret += Le(1) * Ri(1) * std::pow(diffNCombs2D[1 + 1], powV);
133 }
134 else if constexpr (rank == 2)
135 {
136 ret += Le(0) * Ri(0) * std::pow(diffNCombs2D[0 + 3], powV);
137 ret += Le(1) * Ri(1) * std::pow(diffNCombs2D[1 + 3], powV);
138 ret += Le(2) * Ri(2) * std::pow(diffNCombs2D[2 + 3], powV);
139 }
140 else if constexpr (rank == 3)
141 {
142 ret += Le(0) * Ri(0) * std::pow(diffNCombs2D[0 + 6], powV);
143 ret += Le(1) * Ri(1) * std::pow(diffNCombs2D[1 + 6], powV);
144 ret += Le(2) * Ri(2) * std::pow(diffNCombs2D[2 + 6], powV);
145 ret += Le(3) * Ri(3) * std::pow(diffNCombs2D[3 + 6], powV);
146 }
147 else
148 {
149 DNDS_assert(false);
150 }
151 }
152 }
153 return ret;
154 }
155
156 template <int dim, int rank, class VLe, class Trans>
157 void TransSymDiffOrderTensorV(VLe &&Le, Trans &&trans)
158 {
159 if constexpr (dim == 3)
160 {
161 if constexpr (rank == 0)
162 {
163 }
164 else if constexpr (rank == 1)
165 {
166 Le = trans * Le;
167 }
168 else if constexpr (rank == 2)
169 {
170 /*
171 00
172 11
173 22
174 01
175 12
176 02*/
177 Eigen::Matrix3d symTensorR2{
178 {Le(0), Le(3), Le(5)},
179 {Le(3), Le(1), Le(4)},
180 {Le(5), Le(4), Le(2)}};
181 symTensorR2 = trans * symTensorR2 * trans.transpose();
182 Le(0) = symTensorR2(0, 0), Le(1) = symTensorR2(1, 1), Le(2) = symTensorR2(2, 2);
183 Le(3) = symTensorR2(0, 1), Le(4) = symTensorR2(1, 2), Le(5) = symTensorR2(0, 2);
184 }
185 else if constexpr (rank == 3)
186 {
187
188 /*
189 000
190 111
191 222
192 001
193 011
194 112
195 122
196 022
197 002
198 012*/
199 DNDS::ETensor::ETensorR3<real, 3, 3, 3> symTensorR3;
200 symTensorR3(0, 0, 0) = Le(0);
201 symTensorR3(1, 1, 1) = Le(1);
202 symTensorR3(2, 2, 2) = Le(2);
203
204 symTensorR3(0, 0, 1) = symTensorR3(0, 1, 0) = symTensorR3(1, 0, 0) = Le(3);
205 symTensorR3(0, 1, 1) = symTensorR3(1, 1, 0) = symTensorR3(1, 0, 1) = Le(4);
206 symTensorR3(1, 1, 2) = symTensorR3(1, 2, 1) = symTensorR3(2, 1, 1) = Le(5);
207 symTensorR3(1, 2, 2) = symTensorR3(2, 2, 1) = symTensorR3(2, 1, 2) = Le(6);
208 symTensorR3(0, 2, 2) = symTensorR3(2, 2, 0) = symTensorR3(2, 0, 2) = Le(7);
209 symTensorR3(0, 0, 2) = symTensorR3(0, 2, 0) = symTensorR3(2, 0, 0) = Le(8);
210
211 symTensorR3(0, 1, 2) = symTensorR3(1, 2, 0) = symTensorR3(2, 0, 1) =
212 symTensorR3(0, 2, 1) = symTensorR3(2, 1, 0) = symTensorR3(1, 0, 2) = Le(9);
213
214 symTensorR3.MatTransform0(trans.transpose());
215 symTensorR3.MatTransform1(trans.transpose());
216 symTensorR3.MatTransform2(trans.transpose());
217 Le(0) = symTensorR3(0, 0, 0);
218 Le(1) = symTensorR3(1, 1, 1);
219 Le(2) = symTensorR3(2, 2, 2);
220 Le(3) = symTensorR3(0, 0, 1);
221 Le(4) = symTensorR3(0, 1, 1);
222 Le(5) = symTensorR3(1, 1, 2);
223 Le(6) = symTensorR3(1, 2, 2);
224 Le(7) = symTensorR3(0, 2, 2);
225 Le(8) = symTensorR3(0, 0, 2);
226 Le(9) = symTensorR3(0, 1, 2);
227 }
228 else
229 {
230 DNDS_assert(false);
231 }
232 }
233 else // 2-d tensor
234 {
235 if constexpr (rank == 0)
236 {
237 }
238 else if constexpr (rank == 1)
239 {
240 Le = trans * Le;
241 }
242 else if constexpr (rank == 2)
243 {
244 Eigen::Matrix2d symTensorR2{{Le(0), Le(1)},
245 {Le(1), Le(2)}};
246 symTensorR2 = trans * symTensorR2 * trans.transpose();
247 Le(0) = symTensorR2(0, 0), Le(1) = symTensorR2(0, 1), Le(2) = symTensorR2(1, 1);
248 }
249 else if constexpr (rank == 3)
250 {
251
252 DNDS::ETensor::ETensorR3<real, 2, 2, 2> symTensorR3;
253 symTensorR3(0, 0, 0) = Le(0);
254 symTensorR3(0, 0, 1) = symTensorR3(0, 1, 0) = symTensorR3(1, 0, 0) = Le(1);
255 symTensorR3(0, 1, 1) = symTensorR3(1, 0, 1) = symTensorR3(1, 1, 0) = Le(2);
256 symTensorR3(1, 1, 1) = Le(3);
257 symTensorR3.MatTransform0(trans.transpose());
258 symTensorR3.MatTransform1(trans.transpose());
259 symTensorR3.MatTransform2(trans.transpose());
260 Le(0) = symTensorR3(0, 0, 0);
261 Le(1) = symTensorR3(0, 0, 1);
262 Le(2) = symTensorR3(0, 1, 1);
263 Le(3) = symTensorR3(1, 1, 1);
264 }
265 else
266 {
267 DNDS_assert(false);
268 }
269 }
270 }
271
272 /**
273 * @brief convert partial derivatives (diffs) from D_dxi to D_dx
274 *
275 * @tparam dim
276 * @tparam TMat
277 * @param mat standard concatenation of diffs, mat_{i,j}, i=0: D^0_xi, i=1: D^1_{xi_0}, i=2, D^1_{xi_1}... see BaseFunction.hpp
278 * @param dXijdXi dXi_j / d_X_i, Xi: old coord, X: new coord
279 */
280 template <int dim, class TMat>
281 inline void ConvertDiffsLinMap(TMat &&mat, const Geom::tGPoint &dXijdXi)
282 {
283 int rows = mat.rows();
284 if constexpr (dim == 2)
285 switch (rows)
286 {
287 case 10:
288 for (int iB = 0; iB < mat.cols(); iB++)
289 {
290 TransSymDiffOrderTensorV<2, 3>(
291 mat(Eigen::seq(Eigen::fix<6>, Eigen::fix<9>), iB),
292 dXijdXi({0, 1}, {0, 1}));
293 }
294 case 6:
295 for (int iB = 0; iB < mat.cols(); iB++)
296 TransSymDiffOrderTensorV<2, 2>(
297 mat(Eigen::seq(Eigen::fix<3>, Eigen::fix<5>), iB),
298 dXijdXi({0, 1}, {0, 1}));
299
300 case 3:
301 mat({1, 2}, EigenAll) = dXijdXi({0, 1}, {0, 1}) * mat({1, 2}, EigenAll);
302 case 1:
303 break;
304
305 default:
306 std::cerr << mat.rows() << std::endl;
307 DNDS_assert(false);
308 break;
309 }
310 else // dim ==3
311 switch (rows)
312 {
313 case 20:
314 for (int iB = 0; iB < mat.cols(); iB++)
315 {
316 TransSymDiffOrderTensorV<3, 3>(
317 mat(Eigen::seq(Eigen::fix<10>, Eigen::fix<19>), iB),
318 dXijdXi);
319 }
320 case 10:
321 for (int iB = 0; iB < mat.cols(); iB++)
322 TransSymDiffOrderTensorV<3, 2>(
323 mat(Eigen::seq(Eigen::fix<4>, Eigen::fix<9>), iB),
324 dXijdXi);
325
326 case 4:
327 mat({1, 2, 3}, EigenAll) = dXijdXi * mat({1, 2, 3}, EigenAll);
328 case 1:
329 break;
330
331 default:
332 std::cerr << mat.rows() << std::endl;
333 DNDS_assert(false);
334 break;
335 }
336 }
337
338 template <int dim, class TMat, class TDiBjB>
339 inline void ConvertDiffsFullMap(TMat &&mat, TDiBjB &&DxiDx)
340 {
341 int rows = mat.rows();
342 int nBase = mat.cols();
343
344 int cmaxDiffOrder = ndiff2order<dim>(rows);
345 if (cmaxDiffOrder < 0)
346 {
347 std::cerr << mat.rows() << std::endl;
348 DNDS_assert(false);
349 }
350
351 if (cmaxDiffOrder == 0)
352 return;
353 static const auto &diffOperatorIJK2IcDim = getDiffOperatorIJK2I<dim>();
354 MatrixXR out;
355 out.resize(rows, nBase);
356
357 static const int nDiffSizC1 = ndiffSizS<dim>(1);
358 static const auto seq0t1 = Eigen::seq(Eigen::fix<1>, Eigen::fix<nDiffSizC1 - 1>);
359 for (int iB = 0; iB < nBase; iB++)
360 for (int ii = 1; ii < nDiffSizC1; ii++)
361 {
362 int i = dim == 2 ? diffOperatorDimList2D[ii][0] : diffOperatorDimList[ii][0];
363 real &ccv = out(ii, iB);
364 ccv = 0;
365 for (int m = 0; m < dim; m++)
366 ccv += DxiDx(std::get<1>(diffOperatorIJK2IcDim)[i], m) *
367 mat(std::get<1>(diffOperatorIJK2IcDim)[m], iB);
368 }
369 mat(seq0t1, EigenAll) = out(seq0t1, EigenAll);
370 if (cmaxDiffOrder == 1)
371 return;
372
373 static const int nDiffSizC2 = ndiffSizS<dim>(2);
374 static const auto seq1t2 = Eigen::seq(Eigen::fix<nDiffSizC1>, Eigen::fix<nDiffSizC2 - 1>);
375 for (int iB = 0; iB < nBase; iB++)
376 for (int ii = nDiffSizC1; ii < nDiffSizC2; ii++)
377 {
378 int i = dim == 2 ? diffOperatorDimList2D[ii][0] : diffOperatorDimList[ii][0];
379 int j = dim == 2 ? diffOperatorDimList2D[ii][1] : diffOperatorDimList[ii][1];
380 real &ccv = out(ii, iB);
381 ccv = 0;
382 for (int m = 0; m < dim; m++)
383 ccv += DxiDx(std::get<2>(diffOperatorIJK2IcDim)[i][j], m) *
384 mat(std::get<1>(diffOperatorIJK2IcDim)[m], iB);
385 for (int m = 0; m < dim; m++)
386 for (int n = 0; n < dim; n++)
387 ccv += DxiDx(std::get<1>(diffOperatorIJK2IcDim)[i], m) *
388 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[j], n) *
389 mat(std::get<2>(diffOperatorIJK2IcDim)[m][n], iB);
390 }
391 mat(seq1t2, EigenAll) = out(seq1t2, EigenAll);
392 if (cmaxDiffOrder == 2)
393 return;
394
395 static const int nDiffSizC3 = ndiffSizS<dim>(3);
396 static const auto seq2t3 = Eigen::seq(Eigen::fix<nDiffSizC2>, Eigen::fix<nDiffSizC3 - 1>);
397 for (int iB = 0; iB < nBase; iB++)
398 for (int ii = nDiffSizC2; ii < nDiffSizC3; ii++)
399 {
400 int i = dim == 2 ? diffOperatorDimList2D[ii][0] : diffOperatorDimList[ii][0];
401 int j = dim == 2 ? diffOperatorDimList2D[ii][1] : diffOperatorDimList[ii][1];
402 int k = dim == 2 ? diffOperatorDimList2D[ii][2] : diffOperatorDimList[ii][2];
403 real &ccv = out(ii, iB);
404 ccv = 0;
405 for (int m = 0; m < dim; m++)
406 ccv += DxiDx(std::get<3>(diffOperatorIJK2IcDim)[i][j][k], m) *
407 mat(std::get<1>(diffOperatorIJK2IcDim)[m], iB);
408 for (int m = 0; m < dim; m++)
409 for (int n = 0; n < dim; n++)
410 ccv +=
411 DxiDx(std::get<2>(diffOperatorIJK2IcDim)[i][j], m) *
412 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[k], n) *
413 mat(std::get<2>(diffOperatorIJK2IcDim)[m][n], iB) +
414 DxiDx(std::get<2>(diffOperatorIJK2IcDim)[i][k], m) *
415 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[j], n) *
416 mat(std::get<2>(diffOperatorIJK2IcDim)[m][n], iB) +
417 DxiDx(std::get<2>(diffOperatorIJK2IcDim)[j][k], m) *
418 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[i], n) *
419 mat(std::get<2>(diffOperatorIJK2IcDim)[m][n], iB);
420 for (int m = 0; m < dim; m++)
421 for (int n = 0; n < dim; n++)
422 for (int p = 0; p < dim; p++)
423 ccv +=
424 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[i], m) *
425 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[j], n) *
426 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[k], p) *
427 mat(std::get<3>(diffOperatorIJK2IcDim)[m][n][p], iB);
428 }
429 mat(seq2t3, EigenAll) = out(seq2t3, EigenAll);
430 if (cmaxDiffOrder == 3)
431 return;
432 }
433
434 class CFVPeriodicity : public Geom::Periodicity
435 {
436 public:
437 using tBase = Geom::Periodicity;
438 using tBase::tBase;
439
440 CFVPeriodicity(const tBase &vBase) : tBase(vBase) {} // copy from tBase
441
442 template <int dim, class TU>
443 DNDS_DEVICE_CALLABLE void TransDiValueInplace(TU &u, Geom::t_index id)
444 {
445 using namespace Geom;
446 DNDS_assert(FaceIDIsPeriodic(id));
447 t_index i{0};
448 if (FaceIDIsPeriodicDonor(id))
449 i = -id - 3;
450 else
451 i = -id;
452 ConvertDiffsLinMap<dim>(u, rotation[i].map());
453 }
454
455 template <int dim, class TU>
456 DNDS_DEVICE_CALLABLE void TransDiValueBackInplace(TU &u, Geom::t_index id)
457 {
458 using namespace Geom;
459 DNDS_assert(FaceIDIsPeriodic(id));
460 t_index i{0};
461 if (FaceIDIsPeriodicDonor(id))
462 i = -id - 3;
463 else
464 i = -id;
465 ConvertDiffsLinMap<dim>(u, rotation[i].map().transpose());
466 }
467 };
468
469 /**
470 * \brief
471 * invert compact form of all diffs of shape functions
472 * \arg DxDxi NDiff x 3 shape where NDiff is number of derivatives determined by `maxDiff` which is max derivative order
473 * \arg DxiDx NDiff x 3 shape
474 */
475 template <int dim, int maxDiff, class TDiBjA, class TDiBjB>
476 void DxDxi2DxiDx(TDiBjA &&DxDxi, TDiBjB &&DxiDx)
477 {
478 static const int nDiffSizCC = ndiffSizS<dim>(maxDiff);
479 if (maxDiff <= 0) // diff 0 not altered, DxDxi(0, :) are x values and DxiDx(0, :) are xi values (supposedly)
480 return;
481 static const int nDiffSizC1 = ndiffSizS<dim>(1);
482 Eigen::Matrix<real, 3, 3> Dx_j_Dxi_i;
483 if constexpr (dim == 2)
484 {
485 Dx_j_Dxi_i.setZero();
486 Dx_j_Dxi_i(2, 2) = 1;
487 Dx_j_Dxi_i({0, 1}, {0, 1}) = DxDxi({1, 2}, {0, 1});
488 }
489 else
490 Dx_j_Dxi_i = DxDxi({1, 2, 3}, {0, 1, 2});
491 auto lu_cur = Dx_j_Dxi_i.fullPivLu();
492 Eigen::Matrix<real, 3, 3> Dxi_j_Dx_i = lu_cur.inverse(); // Dxi_j_Dx_i Dx_k_Dxi_j = delta_ik
493
494 if constexpr (dim == 2)
495 DxiDx({1, 2}, {0, 1, 2}) = Dxi_j_Dx_i({0, 1}, {0, 1, 2});
496 else
497 DxiDx({1, 2, 3}, {0, 1, 2}) = Dxi_j_Dx_i;
498
499 if (maxDiff == 1)
500 return;
501 static const auto &diffOperatorIJK2IcDim = getDiffOperatorIJK2I<dim>();
502
503 static const int nDiffSizC2 = ndiffSizS<dim>(2);
504 static const auto seq1t2 = Eigen::seq(Eigen::fix<nDiffSizC1>, Eigen::fix<nDiffSizC2 - 1>); // 2D: 3 4 5
505 Eigen::Matrix<real, nDiffSizCC, 3> M_RHS;
506 M_RHS(seq1t2, EigenAll).setZero();
507 for (int j = 0; j < dim; j++)
508 // for (int i = 0; i < dim; i++)
509 // for (int k = 0; k < dim; k++)
510 for (int ii = nDiffSizC1; ii < nDiffSizC2; ii++)
511 {
512 int i = dim == 2 ? diffOperatorDimList2D[ii][0] : diffOperatorDimList[ii][0];
513 int k = dim == 2 ? diffOperatorDimList2D[ii][1] : diffOperatorDimList[ii][1];
514 real &ccv = M_RHS(ii, j);
515 ccv = 0;
516 for (int m = 0; m < dim; m++)
517 for (int n = 0; n < dim; n++)
518 ccv -= DxDxi(std::get<2>(diffOperatorIJK2IcDim)[m][n], j) *
519 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[i], m) *
520 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[k], n);
521 }
522 DxiDx(seq1t2, EigenAll) = M_RHS(seq1t2, EigenAll) /* ikj */ * Dxi_j_Dx_i;
523
524 if (maxDiff == 2)
525 return;
526 static const int nDiffSizC3 = ndiffSizS<dim>(3);
527 static const auto seq2t3 = Eigen::seq(Eigen::fix<nDiffSizC2>, Eigen::fix<nDiffSizC3 - 1>); // 2D: 6 7 8 9
528 M_RHS(seq2t3, EigenAll).setZero();
529 for (int j = 0; j < dim; j++)
530 // for (int i = 0; i < dim; i++)
531 // for (int k = 0; k < dim; k++)
532 // for (int l = 0; l < dim; l++)
533 for (int ii = nDiffSizC2; ii < nDiffSizC3; ii++)
534 {
535 int i = dim == 2 ? diffOperatorDimList2D[ii][0] : diffOperatorDimList[ii][0];
536 int k = dim == 2 ? diffOperatorDimList2D[ii][1] : diffOperatorDimList[ii][1];
537 int l = dim == 2 ? diffOperatorDimList2D[ii][2] : diffOperatorDimList[ii][2];
538 real &ccv = M_RHS(ii, j);
539 ccv = 0;
540 for (int m = 0; m < dim; m++)
541 for (int n = 0; n < dim; n++)
542 ccv -= DxDxi(std::get<2>(diffOperatorIJK2IcDim)[m][n], j) *
543 DxiDx(std::get<2>(diffOperatorIJK2IcDim)[i][k], m) *
544 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[l], n) +
545 DxDxi(std::get<2>(diffOperatorIJK2IcDim)[m][n], j) *
546 DxiDx(std::get<2>(diffOperatorIJK2IcDim)[k][l], m) *
547 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[i], n) +
548 DxDxi(std::get<2>(diffOperatorIJK2IcDim)[m][n], j) *
549 DxiDx(std::get<2>(diffOperatorIJK2IcDim)[i][l], m) *
550 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[k], n);
551 for (int m = 0; m < dim; m++)
552 for (int n = 0; n < dim; n++)
553 for (int p = 0; p < dim; p++)
554 ccv -= DxDxi(std::get<3>(diffOperatorIJK2IcDim)[m][n][p], j) *
555 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[i], m) *
556 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[k], n) *
557 DxiDx(std::get<1>(diffOperatorIJK2IcDim)[l], p);
558 }
559 DxiDx(seq2t3, EigenAll) = M_RHS(seq2t3, EigenAll) /* ikj */ * Dxi_j_Dx_i;
560 if (maxDiff == 3)
561 return;
562 DNDS_assert_info(false, "maxDiff is too large!");
563 }
564}
DNDS_DEVICE_CALLABLE
#define DNDS_DEVICE_CALLABLE
Definition Defines.hpp:76
DNDS_assert_info
#define DNDS_assert_info(expr, info)
Debug-only assertion with an extra std::string info message.
Definition Errors.hpp:113
DNDS_assert
#define DNDS_assert(expr)
Debug-only assertion (compiled out when DNDS_NDEBUG is defined). Prints the expression + file/line + ...
Definition Errors.hpp:108
DNDS::ETensor::ETensorR3::MatTransform0
void MatTransform0(const Tmat &Rmat)
Definition EigenTensor.hpp:60
DNDS::ETensor::ETensorR3::MatTransform1
void MatTransform1(const Tmat &Rmat)
Definition EigenTensor.hpp:71
DNDS::ETensor::ETensorR3::MatTransform2
void MatTransform2(const Tmat &Rmat)
Definition EigenTensor.hpp:82
DNDS::Geom::Base::CFVPeriodicity::TransDiValueBackInplace
DNDS_DEVICE_CALLABLE void TransDiValueBackInplace(TU &u, Geom::t_index id)
Definition DiffTensors.hpp:456
DNDS::Geom::Base::CFVPeriodicity::TransDiValueInplace
DNDS_DEVICE_CALLABLE void TransDiValueInplace(TU &u, Geom::t_index id)
Definition DiffTensors.hpp:443
DNDS::Geom::Base::CFVPeriodicity::CFVPeriodicity
CFVPeriodicity(const tBase &vBase)
Definition DiffTensors.hpp:440
DNDS::Geom::Base::ConvertDiffsLinMap
void ConvertDiffsLinMap(TMat &&mat, const Geom::tGPoint &dXijdXi)
convert partial derivatives (diffs) from D_dxi to D_dx
Definition DiffTensors.hpp:281
DNDS::Geom::Base::NormSymDiffOrderTensorV
real NormSymDiffOrderTensorV(VLe &&Le, VRi &&Ri)
Definition DiffTensors.hpp:8
DNDS::Geom::Base::DxDxi2DxiDx
void DxDxi2DxiDx(TDiBjA &&DxDxi, TDiBjB &&DxiDx)
invert compact form of all diffs of shape functions
Definition DiffTensors.hpp:476
DNDS::Geom::Base::ConvertDiffsFullMap
void ConvertDiffsFullMap(TMat &&mat, TDiBjB &&DxiDx)
Definition DiffTensors.hpp:339
DNDS::Geom::Base::TransSymDiffOrderTensorV
void TransSymDiffOrderTensorV(VLe &&Le, Trans &&trans)
Definition DiffTensors.hpp:157
DNDS::Geom::t_index
int32_t t_index
Definition Geometric.hpp:6
DNDS::Geom::FaceIDIsPeriodic
DNDS_DEVICE_CALLABLE bool FaceIDIsPeriodic(t_index id)
Definition BoundaryCondition.hpp:95
DNDS::Geom::FaceIDIsPeriodicDonor
DNDS_DEVICE_CALLABLE bool FaceIDIsPeriodicDonor(t_index id)
Definition BoundaryCondition.hpp:130
DNDS::Geom::tGPoint
Eigen::Matrix3d tGPoint
Definition Geometric.hpp:11
DNDS::MatrixXR
Eigen::Matrix< real, Eigen::Dynamic, Eigen::Dynamic > MatrixXR
Column-major dynamic Eigen matrix of reals (default layout).
Definition Defines.hpp:205
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
DNDS::Geom::Periodicity::rotation
std::array< tGPointPortable, 4 > rotation
Definition PeriodicInfo.hpp:178
trans
ArrayTransformer< DNDS::real, DynamicSize > trans
Definition test_ArrayTransformer.cpp:169
n
Eigen::Vector3d n(1.0, 0.0, 0.0)