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 <limits>
7 #include "qualityMeasuresJacobian.h"
8 #include "FuncSpaceData.h"
9 #include "MElement.h"
10 #include "BasisFactory.h"
11 #include "bezierBasis.h"
12 #include "JacobianBasis.h"
13 #include "Numeric.h"
14 #include "fullMatrix.h"
15
16 // For regression tests:
17 #include "GModel.h"
18
19 static const double cTri = 2 / std::sqrt(3);
20 static const double cTet = std::sqrt(2);
21 static const double cPyr = 4 * std::sqrt(2);
22
_computeCoeffLengthVectors(const fullMatrix<double> & mat,fullMatrix<double> & coeff,int type,int numCoeff=-1)23 static void _computeCoeffLengthVectors(const fullMatrix<double> &mat,
24 fullMatrix<double> &coeff, int type,
25 int numCoeff = -1)
26 {
27 int sz1 = numCoeff > -1 ? numCoeff : mat.size1();
28
29 switch(type) {
30 case TYPE_QUA: coeff.resize(sz1, 2); break;
31 case TYPE_TRI: coeff.resize(sz1, 3); break;
32 case TYPE_HEX: coeff.resize(sz1, 3); break;
33 case TYPE_PRI: coeff.resize(sz1, 4); break;
34 case TYPE_TET: coeff.resize(sz1, 6); break;
35 case TYPE_PYR: coeff.resize(sz1, 6); break;
36 default:
37 Msg::Warning("Unknown element type %d for quality computation", type);
38 coeff.resize(0, 0);
39 return;
40 }
41
42 if(type != TYPE_PYR) {
43 for(int i = 0; i < sz1; i++) {
44 coeff(i, 0) = std::sqrt(pow_int(mat(i, 0), 2) + pow_int(mat(i, 1), 2) +
45 pow_int(mat(i, 2), 2));
46 coeff(i, 1) = std::sqrt(pow_int(mat(i, 3), 2) + pow_int(mat(i, 4), 2) +
47 pow_int(mat(i, 5), 2));
48 }
49 if(type == TYPE_TRI) {
50 for(int i = 0; i < sz1; i++) {
51 coeff(i, 2) = std::sqrt(pow_int(mat(i, 3) - mat(i, 0), 2) +
52 pow_int(mat(i, 4) - mat(i, 1), 2) +
53 pow_int(mat(i, 5) - mat(i, 2), 2));
54 }
55 }
56 else if(type != TYPE_QUA) { // if 3D
57 for(int i = 0; i < sz1; i++) {
58 coeff(i, 2) = std::sqrt(pow_int(mat(i, 6), 2) + pow_int(mat(i, 7), 2) +
59 pow_int(mat(i, 8), 2));
60 }
61 }
62 if(type == TYPE_TET || type == TYPE_PRI) {
63 for(int i = 0; i < sz1; i++) {
64 coeff(i, 3) = std::sqrt(pow_int(mat(i, 3) - mat(i, 0), 2) +
65 pow_int(mat(i, 4) - mat(i, 1), 2) +
66 pow_int(mat(i, 5) - mat(i, 2), 2));
67 }
68 }
69 if(type == TYPE_TET) {
70 for(int i = 0; i < sz1; i++) {
71 coeff(i, 4) = std::sqrt(pow_int(mat(i, 6) - mat(i, 0), 2) +
72 pow_int(mat(i, 7) - mat(i, 1), 2) +
73 pow_int(mat(i, 8) - mat(i, 2), 2));
74 coeff(i, 5) = std::sqrt(pow_int(mat(i, 6) - mat(i, 3), 2) +
75 pow_int(mat(i, 7) - mat(i, 4), 2) +
76 pow_int(mat(i, 8) - mat(i, 5), 2));
77 }
78 }
79 }
80 else {
81 for(int i = 0; i < sz1; i++) {
82 coeff(i, 0) =
83 std::sqrt(pow_int(2 * mat(i, 0), 2) + pow_int(2 * mat(i, 1), 2) +
84 pow_int(2 * mat(i, 2), 2));
85 coeff(i, 1) =
86 std::sqrt(pow_int(2 * mat(i, 3), 2) + pow_int(2 * mat(i, 4), 2) +
87 pow_int(2 * mat(i, 5), 2));
88 coeff(i, 2) = std::sqrt(pow_int(mat(i, 6) + mat(i, 0) + mat(i, 3), 2) +
89 pow_int(mat(i, 7) + mat(i, 1) + mat(i, 4), 2) +
90 pow_int(mat(i, 8) + mat(i, 2) + mat(i, 5), 2));
91 coeff(i, 3) = std::sqrt(pow_int(mat(i, 6) - mat(i, 0) + mat(i, 3), 2) +
92 pow_int(mat(i, 7) - mat(i, 1) + mat(i, 4), 2) +
93 pow_int(mat(i, 8) - mat(i, 2) + mat(i, 5), 2));
94 coeff(i, 4) = std::sqrt(pow_int(mat(i, 6) - mat(i, 0) - mat(i, 3), 2) +
95 pow_int(mat(i, 7) - mat(i, 1) - mat(i, 4), 2) +
96 pow_int(mat(i, 8) - mat(i, 2) - mat(i, 5), 2));
97 coeff(i, 5) = std::sqrt(pow_int(mat(i, 6) + mat(i, 0) - mat(i, 3), 2) +
98 pow_int(mat(i, 7) + mat(i, 1) - mat(i, 4), 2) +
99 pow_int(mat(i, 8) + mat(i, 2) - mat(i, 5), 2));
100 }
101 }
102 }
103
_computeIGE(const fullVector<double> & det,const fullMatrix<double> & v,fullVector<double> & ige,int type)104 static void _computeIGE(const fullVector<double> &det,
105 const fullMatrix<double> &v, fullVector<double> &ige,
106 int type)
107 {
108 int sz = std::min(det.size(), v.size1());
109 ige.resize(sz);
110
111 switch(type) {
112 case TYPE_QUA:
113 for(int i = 0; i < sz; i++) { ige(i) = det(i) / v(i, 0) / v(i, 1); }
114 break;
115 case TYPE_HEX:
116 for(int i = 0; i < sz; i++) {
117 ige(i) = det(i) / v(i, 0) / v(i, 1) / v(i, 2);
118 }
119 break;
120 case TYPE_TRI:
121 for(int i = 0; i < sz; i++) {
122 ige(i) = cTri * det(i) *
123 (1 / v(i, 0) / v(i, 1) + 1 / v(i, 0) / v(i, 2) +
124 1 / v(i, 1) / v(i, 2)) /
125 3;
126 }
127 break;
128 case TYPE_PRI:
129 for(int i = 0; i < sz; i++) {
130 ige(i) =
131 cTri * det(i) *
132 (1 / v(i, 0) / v(i, 1) / v(i, 2) + 1 / v(i, 0) / v(i, 3) / v(i, 2) +
133 1 / v(i, 1) / v(i, 3) / v(i, 2)) /
134 3;
135 }
136 break;
137 case TYPE_TET:
138 for(int i = 0; i < sz; i++) {
139 ige(i) =
140 cTet * det(i) *
141 (1 / v(i, 0) / v(i, 5) / v(i, 1) + 1 / v(i, 0) / v(i, 5) / v(i, 2) +
142 1 / v(i, 0) / v(i, 5) / v(i, 3) + 1 / v(i, 0) / v(i, 5) / v(i, 4) +
143 1 / v(i, 1) / v(i, 4) / v(i, 0) + 1 / v(i, 1) / v(i, 4) / v(i, 2) +
144 1 / v(i, 1) / v(i, 4) / v(i, 3) + 1 / v(i, 1) / v(i, 4) / v(i, 5) +
145 1 / v(i, 2) / v(i, 3) / v(i, 0) + 1 / v(i, 2) / v(i, 3) / v(i, 1) +
146 1 / v(i, 2) / v(i, 3) / v(i, 4) + 1 / v(i, 2) / v(i, 3) / v(i, 5)) /
147 12;
148 }
149 break;
150 case TYPE_PYR:
151 for(int i = 0; i < sz; i++) {
152 ige(i) =
153 cPyr * det(i) *
154 (1 / v(i, 0) / v(i, 1) / v(i, 2) + 1 / v(i, 0) / v(i, 1) / v(i, 3) +
155 1 / v(i, 0) / v(i, 1) / v(i, 4) + 1 / v(i, 0) / v(i, 1) / v(i, 5) +
156 1 / v(i, 2) / v(i, 3) / v(i, 4) + 1 / v(i, 2) / v(i, 3) / v(i, 5) +
157 1 / v(i, 4) / v(i, 5) / v(i, 2) + 1 / v(i, 4) / v(i, 5) / v(i, 3)) /
158 8;
159 }
160 break;
161 }
162 }
163
_computeICN(const fullVector<double> & det,const fullMatrix<double> & grad,fullVector<double> & icn,int dim)164 static void _computeICN(const fullVector<double> &det,
165 const fullMatrix<double> &grad, fullVector<double> &icn,
166 int dim)
167 {
168 int sz = std::min(det.size(), grad.size1());
169 icn.resize(sz);
170
171 for(int i = 0; i < sz; i++) {
172 double p = 0;
173 for(int k = 0; k < grad.size2(); ++k) { p += pow_int(grad(i, k), 2); }
174 if(dim == 2)
175 icn(i) = 2 * det(i) / p;
176 else // 3D
177 icn(i) = 3 * std::pow(det(i), 2 / 3.) / p;
178 }
179 }
180
_getQualityFunctionSpace(MElement * el,FuncSpaceData & fsGrad,FuncSpaceData & fsDet,int orderSamplingPoints=0)181 static bool _getQualityFunctionSpace(MElement *el, FuncSpaceData &fsGrad,
182 FuncSpaceData &fsDet,
183 int orderSamplingPoints = 0)
184 {
185 const int type = el->getType();
186
187 if(orderSamplingPoints < 1) { // For computing bounds
188 const int order = el->getPolynomialOrder();
189 const int jacOrder = order * el->getDim();
190
191 switch(type) {
192 case TYPE_TRI:
193 fsGrad = FuncSpaceData(el, order - 1, false);
194 fsDet = FuncSpaceData(el, jacOrder - 2, false);
195 break;
196 case TYPE_TET:
197 fsGrad = FuncSpaceData(el, order - 1, false);
198 fsDet = FuncSpaceData(el, jacOrder - 3, false);
199 break;
200 case TYPE_QUA:
201 case TYPE_HEX:
202 case TYPE_PRI:
203 fsGrad = FuncSpaceData(el, order, false);
204 fsDet = FuncSpaceData(el, jacOrder, false);
205 break;
206 case TYPE_PYR:
207 fsGrad = FuncSpaceData(el, false, order, order - 1, false);
208 fsDet = FuncSpaceData(el, false, jacOrder, jacOrder - 3, false);
209 break;
210 default:
211 Msg::Info("Quality measure not implemented for %s",
212 el->getName().c_str());
213 return false;
214 }
215 }
216 else {
217 const int type = el->getType();
218 switch(type) {
219 case TYPE_TRI:
220 case TYPE_TET:
221 case TYPE_QUA:
222 case TYPE_HEX:
223 case TYPE_PRI:
224 fsGrad = FuncSpaceData(el, orderSamplingPoints, false);
225 fsDet = FuncSpaceData(el, orderSamplingPoints, false);
226 break;
227 case TYPE_PYR:
228 fsGrad = FuncSpaceData(el, true, 1, orderSamplingPoints - 1, false);
229 fsDet = FuncSpaceData(el, true, 1, orderSamplingPoints - 1, false);
230 break;
231 default:
232 Msg::Info("Quality measure not implemented for %s",
233 el->getName().c_str());
234 return false;
235 }
236 }
237 return true;
238 }
239
240 namespace jacobianBasedQuality {
241
minMaxJacobianDeterminant(MElement * el,double & min,double & max,const fullMatrix<double> * normals,bool debug)242 void minMaxJacobianDeterminant(MElement *el, double &min, double &max,
243 const fullMatrix<double> *normals, bool debug)
244 {
245 // Get Jacobian basis
246 const JacobianBasis *jfs = el->getJacobianFuncSpace();
247 if(!jfs) {
248 Msg::Warning("Jacobian function space not implemented for %s",
249 el->getName().c_str());
250 min = 99;
251 max = -99;
252 return;
253 }
254
255 // Sample jacobian determinant
256 fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
257 fullVector<double> coeffLag(jfs->getNumSamplingPnts());
258 el->getNodesCoord(nodesXYZ);
259 jfs->getSignedJacobian(nodesXYZ, coeffLag, normals);
260
261 // Convert into Bezier coeff
262 bezierCoeff::usePools(static_cast<std::size_t>(coeffLag.size()), 0);
263 bezierCoeff *bez = new bezierCoeff(jfs->getFuncSpaceData(), coeffLag, 0);
264
265 // Refine coefficients
266 std::vector<_coeffData *> domains(1, new _coeffDataJac(bez));
267 _subdivideDomains(domains, true, debug);
268
269 // Get extrema
270 min = std::numeric_limits<double>::max();
271 max = -min;
272 for(std::size_t i = 0; i < domains.size(); ++i) {
273 min = std::min(min, domains[i]->minB());
274 max = std::max(max, domains[i]->maxB());
275 domains[i]->deleteBezierCoeff();
276 delete domains[i];
277 }
278 }
279
minIGEMeasure(MElement * el,bool knownValid,bool reversedOk,const fullMatrix<double> * normals,bool debug)280 double minIGEMeasure(MElement *el, bool knownValid, bool reversedOk,
281 const fullMatrix<double> *normals, bool debug)
282 {
283 if(!knownValid) {
284 // Computation of the measure should never be performed to invalid
285 // elements (for which the measure is 0).
286 double jmin, jmax;
287 minMaxJacobianDeterminant(el, jmin, jmax, normals);
288 if((jmin <= 0 && jmax >= 0) || (jmax < 0 && !reversedOk)) return 0;
289 }
290
291 // Get Jacobian and gradient bases
292 const GradientBasis *gradBasis;
293 const JacobianBasis *jacBasis;
294 const int tag = el->getTypeForMSH();
295 FuncSpaceData jacMatSpace, jacDetSpace;
296 if(!_getQualityFunctionSpace(el, jacMatSpace, jacDetSpace)) return 0;
297 gradBasis = BasisFactory::getGradientBasis(tag, jacMatSpace);
298 jacBasis = BasisFactory::getJacobianBasis(tag, jacDetSpace);
299
300 // Sample gradients and jacobian determinant
301 fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
302 fullVector<double> coeffDetLag(jacBasis->getNumSamplingPnts());
303 fullMatrix<double> coeffMatLag(gradBasis->getNumSamplingPoints(), 9);
304 el->getNodesCoord(nodesXYZ);
305 jacBasis->getSignedJacobian(nodesXYZ, coeffDetLag, normals);
306 gradBasis->getAllGradientsFromNodes(nodesXYZ, coeffMatLag);
307 // NB: If coeffDetLag(0) is negative, then all coefficients are negative
308 if(coeffDetLag(0) < 0) coeffDetLag.scale(-1);
309 if(el->getDim() == 2) coeffMatLag.resize(coeffMatLag.size1(), 6, false);
310
311 // Convert into Bezier coeff
312 bezierCoeff::usePools(static_cast<std::size_t>(coeffDetLag.size()),
313 static_cast<std::size_t>(coeffMatLag.size1()) *
314 static_cast<std::size_t>(coeffMatLag.size2()));
315 bezierCoeff *bezDet = new bezierCoeff(jacDetSpace, coeffDetLag, 0);
316 bezierCoeff *bezMat = new bezierCoeff(jacMatSpace, coeffMatLag, 1);
317
318 // Compute measure and refine
319 std::vector<_coeffData *> domains;
320 domains.push_back(new _coeffDataIGE(el->getType(), bezDet, bezMat));
321 _subdivideDomains(domains, false, debug);
322
323 return _getMinAndDeleteDomains(domains);
324 }
325
minICNMeasure(MElement * el,bool knownValid,bool reversedOk,const fullMatrix<double> * normals,bool debug)326 double minICNMeasure(MElement *el, bool knownValid, bool reversedOk,
327 const fullMatrix<double> *normals, bool debug)
328 {
329 if(!knownValid) {
330 // Computation of the measure should never
331 // be performed to invalid elements (for which the measure is 0).
332 double jmin, jmax;
333 minMaxJacobianDeterminant(el, jmin, jmax, normals);
334 if((jmin <= 0 && jmax >= 0) || (jmax < 0 && !reversedOk)) return 0;
335 }
336
337 // Get Jacobian and gradient bases
338 const GradientBasis *gradBasis;
339 const JacobianBasis *jacBasis;
340 const int tag = el->getTypeForMSH();
341 FuncSpaceData jacMatSpace, jacDetSpace;
342 if(!_getQualityFunctionSpace(el, jacMatSpace, jacDetSpace)) return 0;
343 gradBasis = BasisFactory::getGradientBasis(tag, jacMatSpace);
344 jacBasis = BasisFactory::getJacobianBasis(tag, jacDetSpace);
345
346 // Sample gradients and jacobian determinant
347 fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
348 fullVector<double> coeffDetLag(jacBasis->getNumSamplingPnts());
349 fullMatrix<double> coeffMatLag(gradBasis->getNumSamplingPoints(), 9);
350 el->getNodesCoord(nodesXYZ);
351 jacBasis->getSignedIdealJacobian(nodesXYZ, coeffDetLag, normals);
352 gradBasis->getAllIdealGradientsFromNodes(nodesXYZ, coeffMatLag);
353 // NB: If coeffDetLag(0) is negative, then all coefficients are negative
354 if(coeffDetLag(0) < 0) coeffDetLag.scale(-1);
355 if(el->getDim() == 2) coeffMatLag.resize(coeffMatLag.size1(), 6, false);
356
357 // Convert into Bezier coeff
358 bezierCoeff::usePools(static_cast<std::size_t>(coeffDetLag.size()),
359 static_cast<std::size_t>(coeffMatLag.size1()) *
360 static_cast<std::size_t>(coeffMatLag.size2()));
361 bezierCoeff *bezDet = new bezierCoeff(jacDetSpace, coeffDetLag, 0);
362 bezierCoeff *bezMat = new bezierCoeff(jacMatSpace, coeffMatLag, 1);
363
364 // Compute measure and refine
365 std::vector<_coeffData *> domains;
366 domains.push_back(new _coeffDataICN(el->getDim(), bezDet, bezMat));
367 _subdivideDomains(domains, false, debug);
368
369 return _getMinAndDeleteDomains(domains);
370 }
371
sampleJacobianDeterminant(MElement * el,int deg,double & min,double & max,const fullMatrix<double> * normals)372 void sampleJacobianDeterminant(MElement *el, int deg, double &min,
373 double &max, const fullMatrix<double> *normals)
374 {
375 fullVector<double> jac;
376 sampleJacobianDeterminant(el, deg, jac, normals);
377
378 min = std::numeric_limits<double>::max();
379 max = -min;
380 for(int i = 0; i < jac.size(); ++i) {
381 min = std::min(min, jac(i));
382 max = std::max(max, jac(i));
383 }
384 }
385
sampleIGEMeasure(MElement * el,int deg,double & min,double & max)386 void sampleIGEMeasure(MElement *el, int deg, double &min, double &max)
387 {
388 fullVector<double> ige;
389 sampleIGEMeasure(el, deg, ige);
390
391 min = std::numeric_limits<double>::max();
392 max = -min;
393 for(int i = 0; i < ige.size(); ++i) {
394 min = std::min(min, ige(i));
395 max = std::max(max, ige(i));
396 }
397 }
398
sampleICNMeasure(MElement * el,int deg,double & min,double & max)399 void sampleICNMeasure(MElement *el, int deg, double &min, double &max)
400 {
401 fullVector<double> icn;
402 sampleICNMeasure(el, deg, icn);
403
404 min = std::numeric_limits<double>::max();
405 max = -min;
406 for(int i = 0; i < icn.size(); ++i) {
407 min = std::min(min, icn(i));
408 max = std::max(max, icn(i));
409 }
410 }
411
sampleJacobianDeterminant(MElement * el,int deg,fullVector<double> & jac,const fullMatrix<double> * normals)412 void sampleJacobianDeterminant(MElement *el, int deg, fullVector<double> &jac,
413 const fullMatrix<double> *normals)
414 {
415 // Get Jacobian basis
416 const JacobianBasis *jacBasis;
417 FuncSpaceData sampleSpace;
418 if(el->getType() != TYPE_PYR)
419 sampleSpace = FuncSpaceData(el, deg, false);
420 else
421 sampleSpace = FuncSpaceData(TYPE_PYR, true, 1, deg - 1, false);
422 jacBasis = BasisFactory::getJacobianBasis(el->getTypeForMSH(), sampleSpace);
423
424 // Sample jacobian determinant
425 fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
426 el->getNodesCoord(nodesXYZ);
427 jac.resize(jacBasis->getNumSamplingPnts());
428 jacBasis->getSignedJacobian(nodesXYZ, jac, normals);
429 }
430
sampleIGEMeasure(MElement * el,int deg,fullVector<double> & ige)431 void sampleIGEMeasure(MElement *el, int deg, fullVector<double> &ige)
432 {
433 // Get Jacobian and gradient bases
434 const GradientBasis *gradBasis;
435 const JacobianBasis *jacBasis;
436 const int tag = el->getTypeForMSH();
437 FuncSpaceData jacMatSpace, jacDetSpace;
438 if(!_getQualityFunctionSpace(el, jacMatSpace, jacDetSpace, deg)) return;
439 gradBasis = BasisFactory::getGradientBasis(tag, jacMatSpace);
440 jacBasis = BasisFactory::getJacobianBasis(tag, jacDetSpace);
441
442 // Sample gradients and jacobian determinant
443 fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
444 fullVector<double> determinant(jacBasis->getNumSamplingPnts());
445 fullMatrix<double> gradients(gradBasis->getNumSamplingPoints(), 9);
446 el->getNodesCoord(nodesXYZ);
447 jacBasis->getSignedJacobian(nodesXYZ, determinant);
448 gradBasis->getAllGradientsFromNodes(nodesXYZ, gradients);
449 if(el->getDim() == 2) gradients.resize(gradients.size1(), 6, false);
450
451 // Compute measure
452 fullMatrix<double> v;
453 const int type = el->getType();
454 _computeCoeffLengthVectors(gradients, v, type);
455 _computeIGE(determinant, v, ige, type);
456 }
457
sampleICNMeasure(MElement * el,int deg,fullVector<double> & icn)458 void sampleICNMeasure(MElement *el, int deg, fullVector<double> &icn)
459 {
460 // Get Jacobian and gradient bases
461 const GradientBasis *gradBasis;
462 const JacobianBasis *jacBasis;
463 const int tag = el->getTypeForMSH();
464 FuncSpaceData jacMatSpace, jacDetSpace;
465 if(!_getQualityFunctionSpace(el, jacMatSpace, jacDetSpace, deg)) return;
466 gradBasis = BasisFactory::getGradientBasis(tag, jacMatSpace);
467 jacBasis = BasisFactory::getJacobianBasis(tag, jacDetSpace);
468
469 // Sample gradients and jacobian determinant
470 fullMatrix<double> nodesXYZ(el->getNumVertices(), 3);
471 fullVector<double> determinant(jacBasis->getNumSamplingPnts());
472 fullMatrix<double> gradients(gradBasis->getNumSamplingPoints(), 9);
473 el->getNodesCoord(nodesXYZ);
474 jacBasis->getSignedIdealJacobian(nodesXYZ, determinant);
475 gradBasis->getAllIdealGradientsFromNodes(nodesXYZ, gradients);
476
477 // Compute measure
478 _computeICN(determinant, gradients, icn, el->getDim());
479 }
480
481 // Virtual class _coeffData
operator ()(_coeffData * cd1,_coeffData * cd2) const482 bool _lessMinB::operator()(_coeffData *cd1, _coeffData *cd2) const
483 {
484 return cd1->minB() > cd2->minB();
485 }
486
operator ()(_coeffData * cd1,_coeffData * cd2) const487 bool _lessMaxB::operator()(_coeffData *cd1, _coeffData *cd2) const
488 {
489 return cd1->maxB() < cd2->maxB();
490 }
491
492 // Jacobian determinant (for validity of all elements)
_coeffDataJac(const bezierCoeff * coeffs)493 _coeffDataJac::_coeffDataJac(const bezierCoeff *coeffs)
494 : _coeffData(), _coeffs(coeffs)
495 {
496 const bezierCoeff &coeff = (*_coeffs);
497
498 _minL = _maxL = coeff.getCornerCoeff(0);
499 for(int i = 1; i < coeff.getNumCornerCoeff(); i++) {
500 _minL = std::min(_minL, coeff.getCornerCoeff(i));
501 _maxL = std::max(_maxL, coeff.getCornerCoeff(i));
502 }
503 _minB = _maxB = coeff(0);
504 for(int i = 1; i < coeff.getNumCoeff(); i++) {
505 _minB = std::min(_minB, coeff(i));
506 _maxB = std::max(_maxB, coeff(i));
507 }
508 }
509
boundsOk(double minL,double maxL) const510 bool _coeffDataJac::boundsOk(double minL, double maxL) const
511 {
512 double tol = std::max(std::abs(minL), std::abs(maxL)) * 1e-3;
513 return (minL <= 0 || _minB > 0) && (maxL >= 0 || _maxB < 0) &&
514 minL - _minB < tol && _maxB - maxL < tol;
515 // NB: First condition implies minL and minB both positive or both negative
516 }
517
getSubCoeff(std::vector<_coeffData * > & v) const518 void _coeffDataJac::getSubCoeff(std::vector<_coeffData *> &v) const
519 {
520 std::vector<bezierCoeff *> sub;
521 _coeffs->subdivide(sub);
522
523 v.clear();
524 for(std::size_t i = 0; i < sub.size(); i++) {
525 v.push_back(new _coeffDataJac(sub[i]));
526 }
527 }
528
deleteBezierCoeff()529 void _coeffDataJac::deleteBezierCoeff() { delete _coeffs; }
530
531 // IGE measure (Inverse Gradient Error)
_coeffDataIGE(int type,const bezierCoeff * det,const bezierCoeff * mat)532 _coeffDataIGE::_coeffDataIGE(int type, const bezierCoeff *det,
533 const bezierCoeff *mat)
534 : _coeffData(), _coeffDet(det), _coeffMat(mat), _type(type)
535 {
536 _computeAtCorner(_minL, _maxL);
537
538 _minB = 0;
539 if(boundsOk(_minL, _maxL))
540 return;
541 else {
542 _minB = _computeLowerBound();
543 }
544 // computation of _maxB not implemented for now
545 }
546
boundsOk(double minL,double maxL) const547 bool _coeffDataIGE::boundsOk(double minL, double maxL) const
548 {
549 static const double tolmin = 1e-3;
550 static const double tolmax = 1e-2;
551 const double tol = tolmin + (tolmax - tolmin) * std::max(_minB, .0);
552 return minL - _minB < tol;
553 }
554
getSubCoeff(std::vector<_coeffData * > & v) const555 void _coeffDataIGE::getSubCoeff(std::vector<_coeffData *> &v) const
556 {
557 std::vector<bezierCoeff *> subD;
558 std::vector<bezierCoeff *> subM;
559 _coeffDet->subdivide(subD);
560 _coeffMat->subdivide(subM);
561
562 v.clear();
563 for(std::size_t i = 0; i < subD.size(); i++) {
564 v.push_back(new _coeffDataIGE(_type, subD[i], subM[i]));
565 }
566 }
567
deleteBezierCoeff()568 void _coeffDataIGE::deleteBezierCoeff()
569 {
570 delete _coeffDet;
571 delete _coeffMat;
572 }
573
_computeAtCorner(double & min,double & max) const574 void _coeffDataIGE::_computeAtCorner(double &min, double &max) const
575 {
576 fullVector<double> det, ige;
577 fullMatrix<double> mat;
578 _coeffDet->getCornerCoeffs(det);
579 _coeffMat->getCornerCoeffs(mat);
580
581 fullMatrix<double> v;
582 _computeCoeffLengthVectors(mat, v, _type);
583 _computeIGE(det, v, ige, _type);
584
585 min = std::numeric_limits<double>::max();
586 max = -min;
587 for(int i = 0; i < ige.size(); ++i) {
588 min = std::min(min, ige(i));
589 max = std::max(max, ige(i));
590 }
591 }
592
_computeLowerBound() const593 double _coeffDataIGE::_computeLowerBound() const
594 {
595 fullVector<double> det;
596 fullMatrix<double> mat;
597 _coeffDet->setVectorAsProxy(det);
598 _coeffMat->setMatrixAsProxy(mat);
599
600 // Speedup: If one coeff _coeffDet is negative, without bounding
601 // J^2/(a^2+b^2), we would get with certainty a negative lower bound.
602 // Returning 0.
603 for(int i = 0; i < det.size(); ++i) {
604 if(det(i) < 0) { return 0; }
605 }
606
607 fullMatrix<double> v;
608 _computeCoeffLengthVectors(mat, v, _type);
609
610 fullVector<double> prox[6];
611 for(int i = 0; i < v.size2(); ++i) { prox[i].setAsProxy(v, i); }
612
613 const bezierBasisRaiser *raiser = _coeffMat->getBezierBasis()->getRaiser();
614 fullVector<double> coeffDenominator;
615 double result = 0;
616
617 switch(_type) {
618 case TYPE_QUA:
619 raiser->computeCoeff(prox[0], prox[1], coeffDenominator);
620 return _computeBoundRational(det, coeffDenominator, true);
621
622 case TYPE_TRI:
623 raiser->computeCoeff(prox[0], prox[1], coeffDenominator);
624 result += _computeBoundRational(det, coeffDenominator, true);
625 raiser->computeCoeff(prox[0], prox[2], coeffDenominator);
626 result += _computeBoundRational(det, coeffDenominator, true);
627 raiser->computeCoeff(prox[1], prox[2], coeffDenominator);
628 result += _computeBoundRational(det, coeffDenominator, true);
629 // The bound is not sharp, which can lead to a lot of subdivision. This
630 // can be avoided by bounding the total function instead of bounding each
631 // rational function and summing the three bounds. This is done for
632 // TYPE_TET and TYPE_PYR. In order to do that for triangles, it is
633 // needed to implement raising bezier coefficient from different spaces.
634 // If computation of sharp bounds is implemented, change also function
635 // _getMinAndDeleteDomains(..) so that it returns minB instead of some
636 // combination of minB and minL.
637 return cTri * result / 3;
638
639 case TYPE_HEX:
640 raiser->computeCoeff(prox[0], prox[1], prox[2], coeffDenominator);
641 return _computeBoundRational(det, coeffDenominator, true);
642
643 case TYPE_PRI:
644 raiser->computeCoeff(prox[0], prox[1], prox[2], coeffDenominator);
645 result += _computeBoundRational(det, coeffDenominator, true);
646 raiser->computeCoeff(prox[0], prox[3], prox[2], coeffDenominator);
647 result += _computeBoundRational(det, coeffDenominator, true);
648 raiser->computeCoeff(prox[1], prox[3], prox[2], coeffDenominator);
649 result += _computeBoundRational(det, coeffDenominator, true);
650 // Same comment than for TYPE_TRI.
651 return cTri * result / 3;
652
653 case TYPE_TET: {
654 fullVector<double> thirdTerm, coeffNum1, tmp;
655 thirdTerm = prox[1];
656 thirdTerm.axpy(prox[2]);
657 thirdTerm.axpy(prox[3]);
658 thirdTerm.axpy(prox[4]);
659 raiser->computeCoeff(prox[0], prox[5], thirdTerm, coeffNum1);
660 thirdTerm = prox[0];
661 thirdTerm.axpy(prox[2]);
662 thirdTerm.axpy(prox[3]);
663 thirdTerm.axpy(prox[5]);
664 raiser->computeCoeff(prox[1], prox[4], thirdTerm, tmp);
665 coeffNum1.axpy(tmp);
666 thirdTerm = prox[0];
667 thirdTerm.axpy(prox[1]);
668 thirdTerm.axpy(prox[4]);
669 thirdTerm.axpy(prox[5]);
670 raiser->computeCoeff(prox[2], prox[3], thirdTerm, tmp);
671 coeffNum1.axpy(tmp);
672
673 fullVector<double> coeffDen1, coeffDen2;
674 raiser->computeCoeff(prox[0], prox[1], prox[2], coeffDen1);
675 raiser->computeCoeff(prox[3], prox[4], prox[5], coeffDen2);
676
677 fullVector<double> &coeffNumerator = tmp;
678 const bezierBasisRaiser *raiserBis =
679 _coeffDet->getBezierBasis()->getRaiser();
680 raiserBis->computeCoeff(coeffNum1, det, coeffNumerator);
681 raiserBis->computeCoeff(coeffDen1, coeffDen2, coeffDenominator);
682
683 result = _computeBoundRational(coeffNumerator, coeffDenominator, true);
684 return cTet * result / 12;
685 }
686
687 case TYPE_PYR: {
688 fullVector<double> thirdTerm, coeffNum1, tmp;
689 thirdTerm = prox[2];
690 thirdTerm.axpy(prox[3]);
691 thirdTerm.axpy(prox[4]);
692 thirdTerm.axpy(prox[5]);
693 raiser->computeCoeff(prox[0], prox[1], thirdTerm, coeffNum1);
694 thirdTerm = prox[4];
695 thirdTerm.axpy(prox[5]);
696 raiser->computeCoeff(prox[2], prox[3], thirdTerm, tmp);
697 coeffNum1.axpy(tmp);
698 thirdTerm = prox[2];
699 thirdTerm.axpy(prox[3]);
700 raiser->computeCoeff(prox[4], prox[5], thirdTerm, tmp);
701 coeffNum1.axpy(tmp);
702
703 fullVector<double> coeffDen1, coeffDen2;
704 raiser->computeCoeff(prox[0], prox[1], prox[2], coeffDen1);
705 raiser->computeCoeff(prox[3], prox[4], prox[5], coeffDen2);
706
707 fullVector<double> &coeffNumerator = tmp;
708 const bezierBasisRaiser *raiserBis =
709 _coeffDet->getBezierBasis()->getRaiser();
710 raiserBis->computeCoeff(coeffNum1, det, coeffNumerator);
711 raiserBis->computeCoeff(coeffDen1, coeffDen2, coeffDenominator);
712
713 result = _computeBoundRational(coeffNumerator, coeffDenominator, true);
714 return cPyr * result / 8;
715 }
716
717 default: Msg::Info("Unknown element type %d for IGE", _type); return -1;
718 }
719 }
720
721 // ICN measure (Inverse Condition Number)
_coeffDataICN(int dim,const bezierCoeff * det,const bezierCoeff * mat)722 _coeffDataICN::_coeffDataICN(int dim, const bezierCoeff *det,
723 const bezierCoeff *mat)
724 : _coeffData(), _coeffDet(det), _coeffMat(mat), _dim(dim)
725 {
726 _computeAtCorner(_minL, _maxL);
727
728 _minB = 0;
729 if(boundsOk(_minL, _maxL))
730 return;
731 else {
732 _minB = _computeLowerBound();
733 }
734 // _maxB not used for now
735 }
736
boundsOk(double minL,double maxL) const737 bool _coeffDataICN::boundsOk(double minL, double maxL) const
738 {
739 static const double tolmin = 1e-3;
740 static const double tolmax = 1e-2;
741 const double tol = tolmin + (tolmax - tolmin) * std::max(_minB, .0);
742 return minL - _minB < tol;
743 }
744
getSubCoeff(std::vector<_coeffData * > & v) const745 void _coeffDataICN::getSubCoeff(std::vector<_coeffData *> &v) const
746 {
747 std::vector<bezierCoeff *> subD;
748 std::vector<bezierCoeff *> subM;
749 _coeffDet->subdivide(subD);
750 _coeffMat->subdivide(subM);
751
752 v.clear();
753 for(std::size_t i = 0; i < subD.size(); i++) {
754 v.push_back(new _coeffDataICN(_dim, subD[i], subM[i]));
755 }
756 }
757
deleteBezierCoeff()758 void _coeffDataICN::deleteBezierCoeff()
759 {
760 delete _coeffDet;
761 delete _coeffMat;
762 }
763
_computeAtCorner(double & min,double & max) const764 void _coeffDataICN::_computeAtCorner(double &min, double &max) const
765 {
766 fullVector<double> det, icn;
767 fullMatrix<double> mat;
768 _coeffDet->getCornerCoeffs(det);
769 _coeffMat->getCornerCoeffs(mat);
770 _computeICN(det, mat, icn, _dim);
771
772 min = std::numeric_limits<double>::max();
773 max = -min;
774 for(int i = 0; i < icn.size(); i++) {
775 min = std::min(min, icn(i));
776 max = std::max(max, icn(i));
777 }
778 }
779
_computeLowerBound() const780 double _coeffDataICN::_computeLowerBound() const
781 {
782 fullVector<double> det;
783 fullMatrix<double> mat;
784 _coeffDet->setVectorAsProxy(det);
785 _coeffMat->setMatrixAsProxy(mat);
786
787 // Speedup: If one coeff _coeffDet is negative, we would get
788 // a negative lower bound. For now, returning 0.
789 for(int i = 0; i < det.size(); ++i) {
790 if(det(i) < 0) { return 0; }
791 }
792
793 const bezierBasisRaiser *raiser = _coeffMat->getBezierBasis()->getRaiser();
794 if(_dim == 2) {
795 fullVector<double> coeffDenominator;
796 {
797 fullVector<double> prox;
798 for(int k = 0; k < mat.size2(); ++k) {
799 prox.setAsProxy(mat, k);
800 fullVector<double> tmp;
801 raiser->computeCoeff(prox, prox, tmp);
802 if(k == 0) coeffDenominator.resize(tmp.size());
803 coeffDenominator.axpy(tmp, 1);
804 }
805 }
806 return 2 * _computeBoundRational(det, coeffDenominator, true);
807 }
808
809 // 3D element
810 fullVector<double> coeffDenominator;
811 {
812 // P: coefficients of function that bound from above the Frobenius norm of
813 // J (elements of P are automatically set to 0)
814 fullVector<double> P(mat.size1());
815 for(int i = 0; i < mat.size1(); ++i) {
816 for(int k = 0; k < mat.size2(); ++k) { P(i) += mat(i, k) * mat(i, k); }
817 P(i) = std::sqrt(P(i));
818 }
819 raiser->computeCoeff(P, P, P, coeffDenominator);
820 }
821
822 const double boundFraction =
823 _computeBoundRational(det, coeffDenominator, true);
824
825 return 3 * std::pow(boundFraction * boundFraction, 1. / 3);
826 }
827
828 // Miscellaneous
829 template <typename Comp>
_subdivideDomainsMinOrMax(std::vector<_coeffData * > & domains,double & minL,double & maxL,bool debug)830 void _subdivideDomainsMinOrMax(std::vector<_coeffData *> &domains,
831 double &minL, double &maxL, bool debug)
832 {
833 std::vector<_coeffData *> subs;
834 make_heap(domains.begin(), domains.end(), Comp());
835 int k = 0;
836 const int max_subdivision = 1000;
837 while(!domains[0]->boundsOk(minL, maxL) && k + 1 < max_subdivision) {
838 _coeffData *cd = domains[0];
839 pop_heap(domains.begin(), domains.end(), Comp());
840 domains.pop_back();
841 cd->getSubCoeff(subs);
842 cd->deleteBezierCoeff();
843 delete cd;
844
845 for(std::size_t i = 0; i < subs.size(); i++) {
846 minL = std::min(minL, subs[i]->minL());
847 maxL = std::max(maxL, subs[i]->maxL());
848 domains.push_back(subs[i]);
849 push_heap(domains.begin(), domains.end(), Comp());
850 }
851 ++k;
852 }
853 if(debug) { std::cout << "Number of subdivisions: " << k << std::endl; }
854 else if(k == max_subdivision) {
855 Msg::Error("Max subdivision (%d) (size domains %d)", max_subdivision,
856 domains.size());
857 }
858 }
859
_subdivideDomains(std::vector<_coeffData * > & domains,bool alsoMax,bool debug)860 void _subdivideDomains(std::vector<_coeffData *> &domains, bool alsoMax,
861 bool debug)
862 {
863 if(domains.empty()) {
864 Msg::Warning("Empty vector in Bezier subdivision, nothing to do");
865 return;
866 }
867 double minL = domains[0]->minL();
868 double maxL = domains[0]->maxL();
869 for(std::size_t i = 1; i < domains.size(); ++i) {
870 minL = std::min(minL, domains[i]->minL());
871 maxL = std::max(maxL, domains[i]->maxL());
872 }
873
874 _subdivideDomainsMinOrMax<_lessMinB>(domains, minL, maxL, debug);
875 if(alsoMax)
876 _subdivideDomainsMinOrMax<_lessMaxB>(domains, minL, maxL, debug);
877 }
878
_getMinAndDeleteDomains(std::vector<_coeffData * > & domains)879 double _getMinAndDeleteDomains(std::vector<_coeffData *> &domains)
880 {
881 double minB = domains[0]->minB();
882 double minL = domains[0]->minL();
883 domains[0]->deleteBezierCoeff();
884 delete domains[0];
885 for(std::size_t i = 1; i < domains.size(); ++i) {
886 minB = std::min(minB, domains[i]->minB());
887 minL = std::min(minL, domains[i]->minL());
888 domains[i]->deleteBezierCoeff();
889 delete domains[i];
890 }
891 double fact = .5 * (minB + minL);
892 // This is done because, for triangles and prisms, currently, the
893 // computation of bounds is not sharp. It can happen than the IGE measure
894 // is very close to 1 everywhere but that the bound is close to 1 only
895 // after a huge amount of subdivision. In this case, it is better to
896 // return minL instead of minB. The best solution would be to implement
897 // sharp bounds for triangles and prisms, see function
898 // _coeffDataIGE::_computeLowerBound(..). If it is done, change this to
899 // return minB.
900 return fact * minL + (1 - fact) * minB;
901 }
902
_computeBoundRational(const fullVector<double> & numerator,const fullVector<double> & denominator,bool lower,bool positiveDenom)903 double _computeBoundRational(const fullVector<double> &numerator,
904 const fullVector<double> &denominator,
905 bool lower, bool positiveDenom)
906 {
907 if(numerator.size() != denominator.size()) {
908 Msg::Error("In order to compute a bound on a rational function, I need "
909 "vectors of the same size! (%d vs %d)",
910 numerator.size(), denominator.size());
911 return 0;
912 }
913
914 // upper and lower bound of the desired bound:
915 const double inf = std::numeric_limits<double>::max();
916 double upperBound = inf;
917 double lowerBound = -inf;
918
919 if(!positiveDenom) lower = !lower;
920
921 if(lower) {
922 // if lower is true, we seek: bound * den <= num
923 for(int i = 0; i < numerator.size(); ++i) {
924 if(denominator(i) == 0) {
925 if(numerator(i) < 0) return -inf;
926 }
927 else if(denominator(i) > 0) {
928 upperBound = std::min(upperBound, numerator(i) / denominator(i));
929 }
930 else {
931 lowerBound = std::max(lowerBound, numerator(i) / denominator(i));
932 }
933 }
934 if(lowerBound > upperBound)
935 return -inf;
936 else
937 return upperBound;
938 }
939 else {
940 // otherwise, we seek: bound * den >= num
941 for(int i = 0; i < numerator.size(); ++i) {
942 if(denominator(i) == 0) {
943 if(numerator(i) > 0) return inf;
944 }
945 else if(denominator(i) > 0) {
946 lowerBound = std::max(lowerBound, numerator(i) / denominator(i));
947 }
948 else {
949 upperBound = std::min(upperBound, numerator(i) / denominator(i));
950 }
951 }
952 if(lowerBound > upperBound)
953 return inf;
954 else
955 return lowerBound;
956 }
957 }
958
testAllMeasuresAllElements()959 void testAllMeasuresAllElements()
960 {
961 GModel *m = GModel::current();
962 std::set<GEntity *, GEntityPtrFullLessThan> entities;
963 for(auto it = m->firstRegion(); it != m->lastRegion(); it++)
964 entities.insert(*it);
965 for(auto it = m->firstFace(); it != m->lastFace(); it++)
966 entities.insert(*it);
967 for(auto it = m->firstEdge(); it != m->lastEdge(); it++)
968 entities.insert(*it);
969
970 for(auto it = entities.begin(); it != entities.end(); ++it) {
971 unsigned num = (*it)->getNumMeshElements();
972 for(unsigned i = 0; i < num; ++i) {
973 MElement *el = (*it)->getMeshElement(i);
974 testAllMeasures(el);
975 }
976 }
977 }
978
testAllMeasures(MElement * el,const fullMatrix<double> * normals)979 void testAllMeasures(MElement *el, const fullMatrix<double> *normals)
980 {
981 static int orderSampling = 50;
982 static double tol = 1e-5;
983 double minSampled, maxSampled, minAlgo, maxAlgo;
984 std::cout << std::endl;
985 std::cout << "Element #" << el->getNum() << " (type: " << el->getType();
986 std::cout << ", " << el->getTypeForMSH() << ")" << std::endl;
987
988 sampleJacobianDeterminant(el, orderSampling, minSampled, maxSampled,
989 normals);
990 minMaxJacobianDeterminant(el, minAlgo, maxAlgo, normals, true);
991 std::cout << "JAC sampled: " << minSampled << " " << maxSampled;
992 std::cout << " v.s. computed: " << minAlgo << " " << maxAlgo << std::endl;
993 if(minSampled < minAlgo * (1 - tol) || maxSampled > maxAlgo * (1 + tol)) {
994 std::cout << "ERROR sampled measure outside the bounds" << std::endl;
995 return;
996 }
997
998 if(minAlgo <= 0 && maxAlgo >= 0) {
999 std::cout << "GOOD (Invalid)" << std::endl;
1000 return;
1001 }
1002
1003 sampleIGEMeasure(el, orderSampling, minSampled, maxSampled);
1004 minAlgo = minIGEMeasure(el, true, true, normals, true);
1005 std::cout << "IGE sampled: " << minSampled << " " << maxSampled;
1006 std::cout << " v.s. computed: " << minAlgo << " -" << std::endl;
1007 if(minSampled < minAlgo * (1 - tol)) {
1008 std::cout << "ERROR sampled measure smaller than the bound" << std::endl;
1009 return;
1010 }
1011
1012 sampleICNMeasure(el, orderSampling, minSampled, maxSampled);
1013 minAlgo = minICNMeasure(el, true, true, normals, true);
1014 std::cout << "ICN sampled: " << minSampled << " " << maxSampled;
1015 std::cout << " v.s. computed: " << minAlgo << " -" << std::endl;
1016 if(minSampled < minAlgo * (1 - tol)) {
1017 std::cout << "ERROR sampled measure smaller than the bound" << std::endl;
1018 return;
1019 }
1020 std::cout << "GOOD" << std::endl;
1021 }
1022
1023 } // end namespace jacobianBasedQuality
1024