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