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