-
Notifications
You must be signed in to change notification settings - Fork 2
/
CondNumBasis.cpp
613 lines (567 loc) · 25.6 KB
/
CondNumBasis.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
// Gmsh - Copyright (C) 1997-2018 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues
#include "CondNumBasis.h"
#include "GmshDefines.h"
#include "GmshMessage.h"
#include <vector>
#include "polynomialBasis.h"
#include "pyramidalBasis.h"
#include "pointsGenerators.h"
#include "BasisFactory.h"
#include "Numeric.h"
namespace {
// Compute the determinant of a 3x3 matrix
inline double calcDet3x3(double M11, double M12, double M13, double M21,
double M22, double M23, double M31, double M32,
double M33)
{
return M11 * (M22 * M33 - M23 * M32) - M12 * (M21 * M33 - M23 * M31) +
M13 * (M21 * M32 - M22 * M31);
}
// Compute the squared Frobenius norm of the inverse of a matrix
template <bool sign>
inline double calcInvCondNum2D(double dxdX, double dxdY, double dydX,
double dydY, double dzdX, double dzdY,
double nx, double ny, double nz)
{
const double dxdXSq = dxdX * dxdX, dydXSq = dydX * dydX,
dzdXSq = dzdX * dzdX;
const double dxdYSq = dxdY * dxdY, dydYSq = dydY * dydY,
dzdYSq = dzdY * dzdY;
const double Dx = dxdXSq - dxdYSq, Dy = dydXSq - dydYSq;
const double Cx = dxdX * dxdY, Cy = dydX * dydY;
const double S1 = dzdYSq * dzdYSq;
const double S2 = (dzdXSq - Dy - Dx) * dzdYSq;
const double S3 = (Cy + Cx) * dzdX * dzdY;
const double S4 = dzdXSq * dzdXSq;
const double S5 = (Dy + Dx) * dzdXSq;
const double S6 = Cx * Cy;
const double S7 = dydXSq * dydXSq;
const double S8 = Dy * dydXSq;
const double S9 = dxdXSq * dxdXSq;
const double S10 = Dx * dxdXSq;
const double S11 = Dy * Dy;
const double S12 = Dx * Dy;
const double S13 = Dx * Dx;
const double S = 2. * (S2 + S5 + S12) + 4. * (S7 - S8 + S9 - S10) +
8. * (S3 + S6) + S1 + S4 + S11 + S13;
const double N = dxdXSq + dxdYSq + dydXSq + dydYSq + dzdXSq + dzdYSq;
const double sqrtS = (S > 0.0) ? sqrt(S) : 0.0;
const double sigma1Sq = 0.5 * (N + sqrtS), sigma2Sq = 0.5 * (N - sqrtS);
const double iCN = 2. * sqrt(sigma1Sq * sigma2Sq) / (sigma1Sq + sigma2Sq);
if(sign) {
const double lnx = dydX * dzdY - dzdX * dydY,
lny = dzdX * dxdY -
dxdX * dzdY, // Local normal from mapping gradients
lnz = dxdX * dydY - dydX * dxdY;
const double dp = lnx * nx + lny * ny +
lnz * nz; // Dot product to determine element validity
return (dp >= 0.) ? iCN : -iCN;
}
else
return iCN;
// return std::min(sqrt(sigma1Sq), sqrt(sigma2Sq)) /
// std::max(sqrt(sigma1Sq), sqrt(sigma2Sq));
}
// Compute the squared Frobenius norm of the inverse of a matrix
template <bool sign>
inline double calcInvCondNum3D(double J11, double J12, double J13, double J21,
double J22, double J23, double J31, double J32,
double J33)
{
const double D = calcDet3x3(J11, J12, J13, J21, J22, J23, J31, J32, J33);
if(D == 0.) return 0.;
const double I11 = J22 * J33 - J23 * J32, I12 = J13 * J32 - J12 * J33,
I13 = J12 * J23 - J13 * J22, I21 = J23 * J31 - J21 * J33,
I22 = J11 * J33 - J13 * J31, I23 = J13 * J21 - J11 * J23,
I31 = J21 * J32 - J22 * J31, I32 = J12 * J31 - J11 * J32,
I33 = J11 * J22 - J12 * J21;
const double nSqJ = J11 * J11 + J12 * J12 + J13 * J13 + J21 * J21 +
J22 * J22 + J23 * J23 + J31 * J31 + J32 * J32 +
J33 * J33;
const double nSqDInvJ = I11 * I11 + I12 * I12 + I13 * I13 + I21 * I21 +
I22 * I22 + I23 * I23 + I31 * I31 + I32 * I32 +
I33 * I33;
if(sign)
return 3. * D / sqrt(nSqJ * nSqDInvJ);
else
return 3. * std::fabs(D) / sqrt(nSqJ * nSqDInvJ);
}
// Compute condition number and its gradients
// w.r.t. node positions, at one location in a 2D element
template <bool sign>
inline void calcGradInvCondNum2D(double dxdX, double dxdY, double dydX,
double dydY, double dzdX, double dzdY,
double nx, double ny, double nz, int i,
int numMapNodes,
const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY,
fullMatrix<double> &IDI)
{
const double EpsDegen = 1.e-6;
bool posJac = true;
if(sign) {
const double lnx = dydX * dzdY - dzdX * dydY,
lny = dzdX * dxdY -
dxdX * dzdY, // Local normal from mapping gradients
lnz = dxdX * dydY - dydX * dxdY;
const double dp = lnx * nx + lny * ny +
lnz * nz; // Dot product to determine element validity
posJac = (dp >= 0.);
}
const double dxdXSq = dxdX * dxdX, dydXSq = dydX * dydX,
dzdXSq = dzdX * dzdX;
const double dxdYSq = dxdY * dxdY, dydYSq = dydY * dydY,
dzdYSq = dzdY * dzdY;
const double Dx = dxdXSq - dxdYSq, Dy = dydXSq - dydYSq;
const double Cx = dxdX * dxdY, Cy = dydX * dydY;
const double S1 = dzdYSq * dzdYSq;
const double S2 = (dzdXSq - Dy - Dx) * dzdYSq;
const double S3 = (Cy + Cx) * dzdX * dzdY;
const double S4 = dzdXSq * dzdXSq;
const double S5 = (Dy + Dx) * dzdXSq;
const double S6 = Cx * Cy;
const double S7 = dydXSq * dydXSq;
const double S8 = Dy * dydXSq;
const double S9 = dxdXSq * dxdXSq;
const double S10 = Dx * dxdXSq;
const double S11 = Dy * Dy;
const double S12 = Dx * Dy;
const double S13 = Dx * Dx;
const double S = 2. * (S2 + S5 + S12) + 4. * (S7 - S8 + S9 - S10) +
8. * (S3 + S6) + S1 + S4 + S11 + S13;
if(S == 0.) { // S == 0. -> Ideal element
for(int j = 0; j < 3 * numMapNodes; j++) IDI(i, j) = 0.;
IDI(i, 3 * numMapNodes) = posJac ? 1. : -1.;
return;
}
const double N = dxdXSq + dxdYSq + dydXSq + dydYSq + dzdXSq + dzdYSq;
const double sqrtS = sqrt(S), invSqrtS = 1. / sqrtS;
const double sigma1Sq = 0.5 * (N + sqrtS), sigma2Sq = 0.5 * (N - sqrtS);
const bool degen =
(sigma2Sq < EpsDegen * sigma1Sq); // Check for degenerate element
const double sum = sigma1Sq + sigma2Sq, invSum = 1. / sum;
const double prod = sigma1Sq * sigma2Sq;
const double sqrtProd = sqrt(prod);
const double halfICN = sqrtProd * invSum;
IDI(i, 3 * numMapNodes) = posJac ? 2. * halfICN : -2. * halfICN;
if(degen) { // Degenerate element: special formula for gradients
const double nnXx = dzdX * ny - dydX * nz, nnXy = dxdX * nz - dzdX * nx,
nnXz = dydX * nx - dxdX * ny;
const double nnYx = dzdY * ny - dydY * nz, nnYy = dxdY * nz - dzdY * nx,
nnYz = dydY * nx - dxdY * ny;
const double fact = 2. / N;
for(int j = 0; j < numMapNodes; j++) {
const double &dPhidX = gSMatX(i, j);
const double &dPhidY = gSMatY(i, j);
IDI(i, j) = fact * (dPhidY * nnXx - dPhidX * nnYx);
IDI(i, j + numMapNodes) = fact * (dPhidY * nnXy - dPhidX * nnYy);
IDI(i, j + 2 * numMapNodes) = fact * (dPhidY * nnXz - dPhidX * nnYz);
}
return;
}
const double invSqrtProd = 1. / sqrtProd;
for(int j = 0; j < numMapNodes; j++) {
const double &dPhidX = gSMatX(i, j);
const double &dPhidY = gSMatY(i, j);
const double ddxdXSqdxj = 2. * dPhidX * dxdX,
ddxdYSqdxj = 2. * dPhidY * dxdY;
const double dDxdxj = ddxdXSqdxj - ddxdYSqdxj;
const double dCxdxj = dPhidX * dxdY + dxdX * dPhidY;
const double dS2dxj = -dDxdxj * dzdYSq;
const double dS3dxj = dCxdxj * dzdX * dzdY;
const double dS5dxj = dDxdxj * dzdXSq;
const double dS6dxj = dCxdxj * Cy;
const double dS9dxj = 2. * ddxdXSqdxj * dxdXSq;
const double dS10dxj = dDxdxj * dxdXSq + Dx * ddxdXSqdxj;
const double dS12dxj = dDxdxj * Dy;
const double dS13dxj = 2. * dDxdxj * Dx;
const double dSdxj = 2. * (dS2dxj + dS5dxj + dS12dxj) +
4. * (dS9dxj - dS10dxj) + 8. * (dS3dxj + dS6dxj) +
dS13dxj;
const double dNdxj = ddxdXSqdxj + ddxdYSqdxj;
const double dsqrtSdxj = 0.5 * dSdxj * invSqrtS;
const double dsigma1Sqdxj = 0.5 * (dNdxj + dsqrtSdxj),
dsigma2Sqdxj = 0.5 * (dNdxj - dsqrtSdxj);
const double dSumdxj = dsigma1Sqdxj + dsigma2Sqdxj;
const double dProddxj = dsigma1Sqdxj * sigma2Sq + sigma1Sq * dsigma2Sqdxj;
const double diCNdxj =
(dProddxj * sum - 2. * prod * dSumdxj) * invSum * invSum * invSqrtProd;
IDI(i, j) = posJac ? diCNdxj : -diCNdxj;
const double ddydXSqdyj = 2. * dPhidX * dydX,
ddydYSqdyj = 2. * dPhidY * dydY;
const double dDydyj = ddydXSqdyj - ddydYSqdyj;
const double dCydyj = dPhidX * dydY + dydX * dPhidY;
const double dS2dyj = -dDydyj * dzdYSq;
const double dS3dyj = dCydyj * dzdX * dzdY;
const double dS5dyj = dDydyj * dzdXSq;
const double dS6dyj = Cx * dCydyj;
const double dS7dyj = 2. * ddydXSqdyj * dydXSq;
const double dS8dyj = dDydyj * dydXSq + Dy * ddydXSqdyj;
const double dS11dyj = 2. * dDydyj * Dy;
const double dS12dyj = Dx * dDydyj;
const double dSdyj = 2. * (dS2dyj + dS5dyj + dS12dyj) +
4. * (dS7dyj - dS8dyj) + 8. * (dS3dyj + dS6dyj) +
dS11dyj;
const double dNdyj = ddydXSqdyj + ddydYSqdyj;
const double dsqrtSdyj = 0.5 * dSdyj * invSqrtS;
const double dsigma1Sqdyj = 0.5 * (dNdyj + dsqrtSdyj),
dsigma2Sqdyj = 0.5 * (dNdyj - dsqrtSdyj);
const double dSumdyj = dsigma1Sqdyj + dsigma2Sqdyj;
const double dProddyj = dsigma1Sqdyj * sigma2Sq + sigma1Sq * dsigma2Sqdyj;
const double diCNdyj =
(dProddyj * sum - 2. * prod * dSumdyj) * invSum * invSum * invSqrtProd;
IDI(i, j + numMapNodes) = posJac ? diCNdyj : -diCNdyj;
const double ddzdXSqdzj = 2. * dPhidX * dzdX,
ddzdYSqdzj = 2. * dPhidY * dzdY;
const double dS1dzj = 2. * ddzdYSqdzj * dzdYSq;
const double dS2dzj = (dzdXSq - Dy - Dx) * ddzdYSqdzj;
const double dS3dzj = (Cy + Cx) * (ddzdXSqdzj * dzdY + dzdX * ddzdYSqdzj);
const double dS4dzj = 2. * ddzdXSqdzj * dzdXSq;
const double dS5dzj = (Dy + Dx) * ddzdXSqdzj;
const double dSdzj =
2. * (dS2dzj + dS5dzj) + 8. * dS3dzj + dS1dzj + dS4dzj;
const double dNdzj = ddzdXSqdzj + ddzdYSqdzj;
const double dsqrtSdzj = 0.5 * dSdzj * invSqrtS;
const double dsigma1Sqdzj = 0.5 * (dNdzj + dsqrtSdzj),
dsigma2Sqdzj = 0.5 * (dNdzj - dsqrtSdzj);
const double dSumdzj = dsigma1Sqdzj + dsigma2Sqdzj;
const double dProddzj = dsigma1Sqdzj * sigma2Sq + sigma1Sq * dsigma2Sqdzj;
const double diCNdzj =
(dProddzj * sum - 2. * prod * dSumdzj) * invSum * invSum * invSqrtProd;
IDI(i, j + 2 * numMapNodes) = posJac ? diCNdzj : -diCNdzj;
}
}
// Compute condition number and its gradients
// w.r.t. node positions, at one location in a 3D element
template <bool sign>
inline void calcGradInvCondNum3D(
double dxdX, double dxdY, double dxdZ, double dydX, double dydY,
double dydZ, double dzdX, double dzdY, double dzdZ, int i, int numMapNodes,
const fullMatrix<double> &gSMatX, const fullMatrix<double> &gSMatY,
const fullMatrix<double> &gSMatZ, fullMatrix<double> &IDI)
{
const double normJSq = dxdX * dxdX + dxdY * dxdY + dxdZ * dxdZ +
dydX * dydX + dydY * dydY + dydZ * dydZ +
dzdX * dzdX + dzdY * dzdY + dzdZ * dzdZ;
const double I11 = dydY * dzdZ - dydZ * dzdY,
I12 = dxdZ * dzdY - dxdY * dzdZ,
I13 = dxdY * dydZ - dxdZ * dydY,
I21 = dydZ * dzdX - dydX * dzdZ,
I22 = dxdX * dzdZ - dxdZ * dzdX,
I23 = dxdZ * dydX - dxdX * dydZ,
I31 = dydX * dzdY - dydY * dzdX,
I32 = dxdY * dzdX - dxdX * dzdY,
I33 = dxdX * dydY - dxdY * dydX;
const double normISq = I11 * I11 + I12 * I12 + I13 * I13 + I21 * I21 +
I22 * I22 + I23 * I23 + I31 * I31 + I32 * I32 +
I33 * I33;
const double invProd = 1. / (normJSq * normISq),
invSqrtProd = sqrt(invProd);
const double D =
calcDet3x3(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY, dzdZ);
const bool reverse = (!sign && (D < 0.));
const double sICN = 3. * D * invSqrtProd;
IDI(i, 3 * numMapNodes) = reverse ? -sICN : sICN;
for(int j = 0; j < numMapNodes; j++) {
const double &dPhidX = gSMatX(i, j);
const double &dPhidY = gSMatY(i, j);
const double &dPhidZ = gSMatZ(i, j);
const double dNormJSqdxj =
2. * (dPhidX * dxdX + dPhidY * dxdY + dPhidZ * dxdZ);
const double dNormISqdxj = 2. * ((dPhidZ * dzdY - dPhidY * dzdZ) * I12 +
(dPhidY * dydZ - dPhidZ * dydY) * I13 +
(dPhidX * dzdZ - dPhidZ * dzdX) * I22 +
(dPhidZ * dydX - dPhidX * dydZ) * I23 +
(dPhidY * dzdX - dPhidX * dzdY) * I32 +
(dPhidX * dydY - dPhidY * dydX) * I33);
const double dProddxj = dNormJSqdxj * normISq + dNormISqdxj * normJSq;
const double dDdxj = dPhidX * dydY * dzdZ + dzdX * dPhidY * dydZ +
dydX * dzdY * dPhidZ - dzdX * dydY * dPhidZ -
dPhidX * dzdY * dydZ - dydX * dPhidY * dzdZ;
const double dsICNdxj =
3. * (dDdxj * invSqrtProd - 0.5 * D * dProddxj * invProd * invSqrtProd);
IDI(i, j) = reverse ? -dsICNdxj : dsICNdxj;
const double dNormJSqdyj =
2. * (dPhidX * dydX + dPhidY * dydY + dPhidZ * dydZ);
const double dNormISqdyj = 2. * ((dPhidY * dzdZ - dPhidZ * dzdY) * I11 +
(dxdY * dPhidZ - dxdZ * dPhidY) * I13 +
(dPhidZ * dzdX - dPhidX * dzdZ) * I21 +
(dxdZ * dPhidX - dxdX * dPhidZ) * I23 +
(dPhidX * dzdY - dPhidY * dzdX) * I31 +
(dxdX * dPhidY - dxdY * dPhidX) * I33);
const double dProddyj = dNormJSqdyj * normISq + dNormISqdyj * normJSq;
const double dDdyj = dxdX * dPhidY * dzdZ + dzdX * dxdY * dPhidZ +
dPhidX * dzdY * dxdZ - dzdX * dPhidY * dxdZ -
dxdX * dzdY * dPhidZ - dPhidX * dxdY * dzdZ;
const double dsICNdyj =
3. * (dDdyj * invSqrtProd - 0.5 * D * dProddyj * invProd * invSqrtProd);
IDI(i, j + numMapNodes) = reverse ? -dsICNdyj : dsICNdyj;
const double dNormJSqdzj =
2. * (dPhidX * dzdX + dPhidY * dzdY + dPhidZ * dzdZ);
const double dNormISqdzj = 2. * ((dydY * dPhidZ - dydZ * dPhidY) * I11 +
(dxdZ * dPhidY - dxdY * dPhidZ) * I12 +
(dydZ * dPhidX - dydX * dPhidZ) * I21 +
(dxdX * dPhidZ - dxdZ * dPhidX) * I22 +
(dydX * dPhidY - dydY * dPhidX) * I31 +
(dxdY * dPhidX - dxdX * dPhidY) * I32);
const double dProddzj = dNormJSqdzj * normISq + dNormISqdzj * normJSq;
const double dDdzj = dxdX * dydY * dPhidZ + dPhidX * dxdY * dydZ +
dydX * dPhidY * dxdZ - dPhidX * dydY * dxdZ -
dxdX * dPhidY * dydZ - dydX * dxdY * dPhidZ;
const double dsICNdzj =
3. * (dDdzj * invSqrtProd - 0.5 * D * dProddzj * invProd * invSqrtProd);
IDI(i, j + 2 * numMapNodes) = reverse ? -dsICNdzj : dsICNdzj;
}
}
} // namespace
CondNumBasis::CondNumBasis(int tag, int cnOrder)
: _tag(tag), _dim(ElementType::getDimension(tag)),
_condNumOrder(cnOrder >= 0 ? cnOrder : condNumOrder(tag))
{
if(ElementType::getParentType(tag) == TYPE_TRIH) {
_nCondNumNodes = 1;
_nMapNodes = 4;
_nPrimMapNodes = 4;
return;
}
const bool serendip = false;
FuncSpaceData data =
(ElementType::getParentType(tag) == TYPE_PYR) ?
FuncSpaceData(true, tag, true, 1, _condNumOrder - 1, &serendip) :
FuncSpaceData(true, tag, _condNumOrder, &serendip);
fullMatrix<double> lagPoints; // Sampling points
gmshGeneratePoints(data, lagPoints);
_nCondNumNodes = lagPoints.size1();
_nMapNodes = BasisFactory::getNodalBasis(tag)->getNumShapeFunctions();
// Store shape function gradients of mapping at condition number nodes
_gradBasis = BasisFactory::getGradientBasis(data);
// Compute shape function gradients of primary mapping at barycenter,
// in order to compute normal to straight element
const int parentType = ElementType::getParentType(tag);
const int primMapType = ElementType::getType(parentType, 1, false);
const nodalBasis *primMapBasis = BasisFactory::getNodalBasis(primMapType);
_nPrimMapNodes = primMapBasis->getNumShapeFunctions();
double xBar = 0., yBar = 0., zBar = 0.;
double barycenter[3] = {0., 0., 0.};
for(int i = 0; i < _nPrimMapNodes; i++) {
for(int j = 0; j < primMapBasis->points.size2(); ++j) {
barycenter[j] += primMapBasis->points(i, j);
}
}
barycenter[0] /= _nPrimMapNodes;
barycenter[1] /= _nPrimMapNodes;
barycenter[2] /= _nPrimMapNodes;
double(*barDPsi)[3] = new double[_nPrimMapNodes][3];
primMapBasis->df(xBar, yBar, zBar, barDPsi);
// TODO: Make primGradShape from ideal element
primGradShapeBarycenterX.resize(_nPrimMapNodes);
primGradShapeBarycenterY.resize(_nPrimMapNodes);
primGradShapeBarycenterZ.resize(_nPrimMapNodes);
for(int j = 0; j < _nPrimMapNodes; j++) {
primGradShapeBarycenterX(j) = barDPsi[j][0];
primGradShapeBarycenterY(j) = barDPsi[j][1];
primGradShapeBarycenterZ(j) = barDPsi[j][2];
}
delete[] barDPsi;
}
int CondNumBasis::condNumOrder(int tag)
{
const int parentType = ElementType::getParentType(tag);
const int order = ElementType::getOrder(tag);
return condNumOrder(parentType, order);
}
int CondNumBasis::condNumOrder(int parentType, int order)
{
switch(parentType) {
case TYPE_PNT: return 0;
case TYPE_LIN: return order - 1;
case TYPE_TRI: return (order == 1) ? 0 : order;
case TYPE_QUA: return order;
case TYPE_TET: return (order == 1) ? 0 : order;
case TYPE_PRI: return order;
case TYPE_HEX: return order;
case TYPE_PYR: return order;
case TYPE_TRIH: return 0;
default:
Msg::Error("Unknown element type %d, return order 0", parentType);
return 0;
}
}
// Calculate the inverse condition number in Frobenius norm for one element,
// with normal vectors to straight element for regularization. Evaluation points
// depend on the given matrices for shape function gradients.
template <bool sign>
inline void CondNumBasis::getInvCondNumGeneral(
int nCondNumNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
fullVector<double> &condNum) const
{
switch(_dim) {
case 0: {
for(int i = 0; i < nCondNumNodes; i++) condNum(i) = 1.;
break;
}
case 1: {
Msg::Warning("Inverse condition number not implemented in 1D");
condNum.setAll(0.);
break;
}
case 2: {
fullMatrix<double> dxyzdX(nCondNumNodes, 3), dxyzdY(nCondNumNodes, 3);
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
for(int i = 0; i < nCondNumNodes; i++) {
const double &dxdX = dxyzdX(i, 0), &dydX = dxyzdX(i, 1),
&dzdX = dxyzdX(i, 2);
const double &dxdY = dxyzdY(i, 0), &dydY = dxyzdY(i, 1),
&dzdY = dxyzdY(i, 2);
const double &nx = normals(0, 0), &ny = normals(0, 1),
&nz = normals(0, 2);
condNum(i) =
calcInvCondNum2D<sign>(dxdX, dxdY, dydX, dydY, dzdX, dzdY, nx, ny, nz);
}
break;
}
case 3: {
if(ElementType::getParentType(_tag) == TYPE_TRIH) {
for(int i = 0; i < nCondNumNodes; i++) condNum(i) = 1.;
break;
}
fullMatrix<double> dxyzdX(nCondNumNodes, 3), dxyzdY(nCondNumNodes, 3),
dxyzdZ(nCondNumNodes, 3);
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
gSMatZ.mult(nodesXYZ, dxyzdZ);
for(int i = 0; i < nCondNumNodes; i++) {
const double &dxdX = dxyzdX(i, 0), &dydX = dxyzdX(i, 1),
&dzdX = dxyzdX(i, 2);
const double &dxdY = dxyzdY(i, 0), &dydY = dxyzdY(i, 1),
&dzdY = dxyzdY(i, 2);
const double &dxdZ = dxyzdZ(i, 0), &dydZ = dxyzdZ(i, 1),
&dzdZ = dxyzdZ(i, 2);
condNum(i) = calcInvCondNum3D<sign>(dxdX, dxdY, dxdZ, dydX, dydY, dydZ,
dzdX, dzdY, dzdZ);
}
break;
}
}
}
void CondNumBasis::getInvCondNumGeneral(int nCondNumNodes,
const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY,
const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ,
fullVector<double> &invCond) const
{
fullMatrix<double> dumNormals;
getInvCondNumGeneral<false>(nCondNumNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ,
dumNormals, invCond);
}
void CondNumBasis::getSignedInvCondNumGeneral(
int nCondNumNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
fullVector<double> &invCond) const
{
getInvCondNumGeneral<true>(nCondNumNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ,
normals, invCond);
}
// Calculate the inverse condition number in Frobenius norm and its gradients
// w.r.t. node position, with normal vectors to straight element for
// regularization. Evaluation points depend on the given matrices for shape
// function gradients.
template <bool sign>
inline void CondNumBasis::getInvCondNumAndGradientsGeneral(
int nCondNumNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
fullMatrix<double> &IDI) const
{
fullMatrix<double> JDJ(nCondNumNodes, 3 * _nMapNodes + 1);
switch(_dim) {
case 0: {
for(int i = 0; i < nCondNumNodes; i++) {
for(int j = 0; j < _nMapNodes; j++) {
IDI(i, j) = 0.;
IDI(i, j + 1 * _nMapNodes) = 0.;
IDI(i, j + 2 * _nMapNodes) = 0.;
}
IDI(i, 3 * _nMapNodes) = 1.;
}
break;
}
case 1: {
Msg::Warning("Inverse condition number not implemented in 1D");
IDI.setAll(0.);
break;
}
case 2: {
fullMatrix<double> dxyzdX(nCondNumNodes, 3), dxyzdY(nCondNumNodes, 3);
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
for(int i = 0; i < nCondNumNodes; i++) {
const double &dxdX = dxyzdX(i, 0), &dydX = dxyzdX(i, 1),
&dzdX = dxyzdX(i, 2);
const double &dxdY = dxyzdY(i, 0), &dydY = dxyzdY(i, 1),
&dzdY = dxyzdY(i, 2);
const double &nx = normals(0, 0), &ny = normals(0, 1),
&nz = normals(0, 2);
calcGradInvCondNum2D<sign>(dxdX, dxdY, dydX, dydY, dzdX, dzdY, nx, ny, nz,
i, _nMapNodes, gSMatX, gSMatY, IDI);
}
break;
}
case 3: {
if(ElementType::getParentType(_tag) == TYPE_TRIH) {
for(int i = 0; i < nCondNumNodes; i++) {
for(int j = 0; j < _nMapNodes; j++) {
IDI(i, j) = 0.;
IDI(i, j + 1 * _nMapNodes) = 0.;
IDI(i, j + 2 * _nMapNodes) = 0.;
}
IDI(i, 3 * _nMapNodes) = 1.;
}
break;
}
fullMatrix<double> dxyzdX(nCondNumNodes, 3), dxyzdY(nCondNumNodes, 3),
dxyzdZ(nCondNumNodes, 3);
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
gSMatZ.mult(nodesXYZ, dxyzdZ);
for(int i = 0; i < nCondNumNodes; i++) {
const double &dxdX = dxyzdX(i, 0), &dydX = dxyzdX(i, 1),
&dzdX = dxyzdX(i, 2);
const double &dxdY = dxyzdY(i, 0), &dydY = dxyzdY(i, 1),
&dzdY = dxyzdY(i, 2);
const double &dxdZ = dxyzdZ(i, 0), &dydZ = dxyzdZ(i, 1),
&dzdZ = dxyzdZ(i, 2);
calcGradInvCondNum3D<sign>(dxdX, dxdY, dxdZ, dydX, dydY, dydZ, dzdX, dzdY,
dzdZ, i, _nMapNodes, gSMatX, gSMatY, gSMatZ,
IDI);
}
break;
}
}
}
void CondNumBasis::getInvCondNumAndGradientsGeneral(
int nCondNumNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullMatrix<double> &IDI) const
{
fullMatrix<double> dumNormals;
getInvCondNumAndGradientsGeneral<false>(nCondNumNodes, gSMatX, gSMatY, gSMatZ,
nodesXYZ, dumNormals, IDI);
}
void CondNumBasis::getSignedInvCondNumAndGradientsGeneral(
int nCondNumNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
fullMatrix<double> &IDI) const
{
getInvCondNumAndGradientsGeneral<true>(nCondNumNodes, gSMatX, gSMatY, gSMatZ,
nodesXYZ, normals, IDI);
}