1 /*=========================================================================
2 
3   Program:   Visualization Toolkit
4   Module:    vtkMeshQuality.cxx
5 
6   Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7   All rights reserved.
8   See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9 
10      This software is distributed WITHOUT ANY WARRANTY; without even
11      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12      PURPOSE.  See the above copyright notice for more information.
13 
14   Copyright 2003-2008 Sandia Corporation.
15   Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
16   license for use of this work by or on behalf of the
17   U.S. Government. Redistribution and use in source and binary forms, with
18   or without modification, are permitted provided that this Notice and any
19   statement of authorship are reproduced on all copies.
20 
21   Contact: dcthomp@sandia.gov,pppebay@sandia.gov
22 
23 =========================================================================*/
24 
25 // Hide VTK_DEPRECATED_IN_9_0_0() warnings for this class.
26 #define VTK_DEPRECATION_LEVEL 0
27 
28 #include "vtkMeshQuality.h"
29 #include "vtkCell.h"
30 #include "vtkCellData.h"
31 #include "vtkCellTypes.h"
32 #include "vtkDataArray.h"
33 #include "vtkDataSet.h"
34 #include "vtkDoubleArray.h"
35 #include "vtkInformation.h"
36 #include "vtkInformationVector.h"
37 #include "vtkMath.h"
38 #include "vtkObjectFactory.h"
39 #include "vtkPoints.h"
40 #include "vtkTetra.h"
41 #include "vtkTriangle.h"
42 
43 #include "vtk_verdict.h"
44 
45 vtkStandardNewMacro(vtkMeshQuality);
46 
47 typedef double (*CellQualityType)(vtkCell*);
48 
49 double TetVolume(vtkCell* cell);
50 
51 static const char* QualityMeasureNames[] = { "EdgeRatio", "AspectRatio", "RadiusRatio",
52   "AspectFrobenius", "MedAspectFrobenius", "MaxAspectFrobenius", "MinAngle", "CollapseRatio",
53   "MaxAngle", "Condition", "ScaledJacobian", "Shear", "RelativeSizeSquared", "Shape",
54   "ShapeAndSize", "Distortion", "MaxEdgeRatio", "Skew", "Taper", "Volume", "Stretch", "Diagonal",
55   "Dimension", "Oddy", "ShearAndSize", "Jacobian", "Warpage", "AspectGamma", "Area", "AspectBeta" };
56 
57 double vtkMeshQuality::CurrentTriNormal[3];
58 
PrintSelf(ostream & os,vtkIndent indent)59 void vtkMeshQuality::PrintSelf(ostream& os, vtkIndent indent)
60 {
61   const char onStr[] = "On";
62   const char offStr[] = "Off";
63 
64   this->Superclass::PrintSelf(os, indent);
65 
66   os << indent << "SaveCellQuality:   " << (this->SaveCellQuality ? onStr : offStr) << endl;
67   os << indent << "TriangleQualityMeasure: " << QualityMeasureNames[this->TriangleQualityMeasure]
68      << endl;
69   os << indent << "QuadQualityMeasure: " << QualityMeasureNames[this->QuadQualityMeasure] << endl;
70   os << indent << "TetQualityMeasure: " << QualityMeasureNames[this->TetQualityMeasure] << endl;
71   os << indent << "HexQualityMeasure: " << QualityMeasureNames[this->HexQualityMeasure] << endl;
72   os << indent << "Volume: " << (this->Volume ? onStr : offStr) << endl;
73   os << indent << "CompatibilityMode: " << (this->CompatibilityMode ? onStr : offStr) << endl;
74 }
75 
vtkMeshQuality()76 vtkMeshQuality::vtkMeshQuality()
77 {
78   this->SaveCellQuality = 1; // Default is On
79   this->TriangleQualityMeasure = VTK_QUALITY_ASPECT_RATIO;
80   this->QuadQualityMeasure = VTK_QUALITY_EDGE_RATIO;
81   this->TetQualityMeasure = VTK_QUALITY_ASPECT_RATIO;
82   this->HexQualityMeasure = VTK_QUALITY_MAX_ASPECT_FROBENIUS;
83   this->Volume = 0;
84   this->CompatibilityMode = 0;
85 }
86 
87 vtkMeshQuality::~vtkMeshQuality() = default;
88 
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)89 int vtkMeshQuality::RequestData(vtkInformation* vtkNotUsed(request),
90   vtkInformationVector** inputVector, vtkInformationVector* outputVector)
91 {
92   // get the info objects
93   vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
94   vtkInformation* outInfo = outputVector->GetInformationObject(0);
95 
96   // get the input and output
97   vtkDataSet* in = vtkDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
98   vtkDataSet* out = vtkDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
99 
100   CellQualityType TriangleQuality, QuadQuality, TetQuality, HexQuality;
101   vtkDoubleArray* quality = nullptr;
102   vtkDoubleArray* volume = nullptr;
103   vtkIdType N = in->GetNumberOfCells();
104   double qtrim, qtriM, Eqtri, Eqtri2;
105   double qquam, qquaM, Eqqua, Eqqua2;
106   double qtetm, qtetM, Eqtet, Eqtet2;
107   double qhexm, qhexM, Eqhex, Eqhex2;
108   double q;
109   double V = 0.;
110   vtkIdType ntri = 0;
111   vtkIdType nqua = 0;
112   vtkIdType ntet = 0;
113   vtkIdType nhex = 0;
114   vtkCell* cell;
115   int progressNumer = 0;
116   double progressDenom = 20.;
117 
118   this->CellNormals = in->GetCellData()->GetNormals();
119 
120   if (this->CellNormals)
121     v_set_tri_normal_func(
122       reinterpret_cast<ComputeNormal>(vtkMeshQuality::GetCurrentTriangleNormal));
123   else
124     v_set_tri_normal_func(nullptr);
125 
126   // Initialize the min and max values, std deviations, etc.
127   qtriM = qquaM = qtetM = qhexM = VTK_DOUBLE_MIN;
128   qtrim = qquam = qtetm = qhexm = VTK_DOUBLE_MAX;
129   Eqtri = Eqtri2 = Eqqua = Eqqua2 = Eqtet = Eqtet2 = Eqhex = Eqhex2 = 0.;
130 
131   switch (this->GetTriangleQualityMeasure())
132   {
133     case VTK_QUALITY_AREA:
134       TriangleQuality = TriangleArea;
135       break;
136     case VTK_QUALITY_EDGE_RATIO:
137       TriangleQuality = TriangleEdgeRatio;
138       break;
139     case VTK_QUALITY_ASPECT_RATIO:
140       TriangleQuality = TriangleAspectRatio;
141       break;
142     case VTK_QUALITY_RADIUS_RATIO:
143       TriangleQuality = TriangleRadiusRatio;
144       break;
145     case VTK_QUALITY_ASPECT_FROBENIUS:
146       TriangleQuality = TriangleAspectFrobenius;
147       break;
148     case VTK_QUALITY_MIN_ANGLE:
149       TriangleQuality = TriangleMinAngle;
150       break;
151     case VTK_QUALITY_MAX_ANGLE:
152       TriangleQuality = TriangleMaxAngle;
153       break;
154     case VTK_QUALITY_CONDITION:
155       TriangleQuality = TriangleCondition;
156       break;
157     case VTK_QUALITY_SCALED_JACOBIAN:
158       TriangleQuality = TriangleScaledJacobian;
159       break;
160     case VTK_QUALITY_RELATIVE_SIZE_SQUARED:
161       TriangleQuality = TriangleRelativeSizeSquared;
162       break;
163     case VTK_QUALITY_SHAPE:
164       TriangleQuality = TriangleShape;
165       break;
166     case VTK_QUALITY_SHAPE_AND_SIZE:
167       TriangleQuality = TriangleShapeAndSize;
168       break;
169     case VTK_QUALITY_DISTORTION:
170       TriangleQuality = TriangleDistortion;
171       break;
172     default:
173       vtkWarningMacro("Bad TriangleQualityMeasure (" << this->GetTriangleQualityMeasure()
174                                                      << "), using RadiusRatio instead");
175       TriangleQuality = TriangleRadiusRatio;
176       break;
177   }
178 
179   switch (this->GetQuadQualityMeasure())
180   {
181     case VTK_QUALITY_EDGE_RATIO:
182       QuadQuality = QuadEdgeRatio;
183       break;
184     case VTK_QUALITY_ASPECT_RATIO:
185       QuadQuality = QuadAspectRatio;
186       break;
187     case VTK_QUALITY_RADIUS_RATIO:
188       QuadQuality = QuadRadiusRatio;
189       break;
190     case VTK_QUALITY_MED_ASPECT_FROBENIUS:
191       QuadQuality = QuadMedAspectFrobenius;
192       break;
193     case VTK_QUALITY_MAX_ASPECT_FROBENIUS:
194       QuadQuality = QuadMaxAspectFrobenius;
195       break;
196     case VTK_QUALITY_MIN_ANGLE:
197       QuadQuality = QuadMinAngle;
198       break;
199     case VTK_QUALITY_MAX_EDGE_RATIO:
200       QuadQuality = QuadMaxEdgeRatios;
201       break;
202     case VTK_QUALITY_SKEW:
203       QuadQuality = QuadSkew;
204       break;
205     case VTK_QUALITY_TAPER:
206       QuadQuality = QuadTaper;
207       break;
208     case VTK_QUALITY_WARPAGE:
209       QuadQuality = QuadWarpage;
210       break;
211     case VTK_QUALITY_AREA:
212       QuadQuality = QuadArea;
213       break;
214     case VTK_QUALITY_STRETCH:
215       QuadQuality = QuadStretch;
216       break;
217     // case VTK_QUALITY_MIN_ANGLE:
218     case VTK_QUALITY_MAX_ANGLE:
219       QuadQuality = QuadMaxAngle;
220       break;
221     case VTK_QUALITY_ODDY:
222       QuadQuality = QuadOddy;
223       break;
224     case VTK_QUALITY_CONDITION:
225       QuadQuality = QuadCondition;
226       break;
227     case VTK_QUALITY_JACOBIAN:
228       QuadQuality = QuadJacobian;
229       break;
230     case VTK_QUALITY_SCALED_JACOBIAN:
231       QuadQuality = QuadScaledJacobian;
232       break;
233     case VTK_QUALITY_SHEAR:
234       QuadQuality = QuadShear;
235       break;
236     case VTK_QUALITY_SHAPE:
237       QuadQuality = QuadShape;
238       break;
239     case VTK_QUALITY_RELATIVE_SIZE_SQUARED:
240       QuadQuality = QuadRelativeSizeSquared;
241       break;
242     case VTK_QUALITY_SHAPE_AND_SIZE:
243       QuadQuality = QuadShapeAndSize;
244       break;
245     case VTK_QUALITY_SHEAR_AND_SIZE:
246       QuadQuality = QuadShearAndSize;
247       break;
248     case VTK_QUALITY_DISTORTION:
249       QuadQuality = QuadDistortion;
250       break;
251     default:
252       vtkWarningMacro("Bad QuadQualityMeasure (" << this->GetQuadQualityMeasure()
253                                                  << "), using EdgeRatio instead");
254       QuadQuality = QuadEdgeRatio;
255       break;
256   }
257 
258   switch (this->GetTetQualityMeasure())
259   {
260     case VTK_QUALITY_EDGE_RATIO:
261       TetQuality = TetEdgeRatio;
262       break;
263     case VTK_QUALITY_ASPECT_RATIO:
264       TetQuality = TetAspectRatio;
265       break;
266     case VTK_QUALITY_RADIUS_RATIO:
267       TetQuality = TetRadiusRatio;
268       break;
269     case VTK_QUALITY_ASPECT_FROBENIUS:
270       TetQuality = TetAspectFrobenius;
271       break;
272     case VTK_QUALITY_MIN_ANGLE:
273       TetQuality = TetMinAngle;
274       break;
275     case VTK_QUALITY_COLLAPSE_RATIO:
276       TetQuality = TetCollapseRatio;
277       break;
278     case VTK_QUALITY_ASPECT_BETA:
279       TetQuality = TetAspectBeta;
280       break;
281     case VTK_QUALITY_ASPECT_GAMMA:
282       TetQuality = TetAspectGamma;
283       break;
284     case VTK_QUALITY_VOLUME:
285       TetQuality = TetVolume;
286       break;
287     case VTK_QUALITY_CONDITION:
288       TetQuality = TetCondition;
289       break;
290     case VTK_QUALITY_JACOBIAN:
291       TetQuality = TetJacobian;
292       break;
293     case VTK_QUALITY_SCALED_JACOBIAN:
294       TetQuality = TetScaledJacobian;
295       break;
296     case VTK_QUALITY_SHAPE:
297       TetQuality = TetShape;
298       break;
299     case VTK_QUALITY_RELATIVE_SIZE_SQUARED:
300       TetQuality = TetRelativeSizeSquared;
301       break;
302     case VTK_QUALITY_SHAPE_AND_SIZE:
303       TetQuality = TetShapeandSize;
304       break;
305     case VTK_QUALITY_DISTORTION:
306       TetQuality = TetDistortion;
307       break;
308     default:
309       vtkWarningMacro("Bad TetQualityMeasure (" << this->GetTetQualityMeasure()
310                                                 << "), using RadiusRatio instead");
311       TetQuality = TetRadiusRatio;
312       break;
313   }
314 
315   switch (this->GetHexQualityMeasure())
316   {
317     case VTK_QUALITY_EDGE_RATIO:
318       HexQuality = HexEdgeRatio;
319       break;
320     case VTK_QUALITY_MED_ASPECT_FROBENIUS:
321       HexQuality = HexMedAspectFrobenius;
322       break;
323     case VTK_QUALITY_MAX_ASPECT_FROBENIUS:
324       HexQuality = HexMaxAspectFrobenius;
325       break;
326     case VTK_QUALITY_MAX_EDGE_RATIO:
327       HexQuality = HexMaxEdgeRatio;
328       break;
329     case VTK_QUALITY_SKEW:
330       HexQuality = HexSkew;
331       break;
332     case VTK_QUALITY_TAPER:
333       HexQuality = HexTaper;
334       break;
335     case VTK_QUALITY_VOLUME:
336       HexQuality = HexVolume;
337       break;
338     case VTK_QUALITY_STRETCH:
339       HexQuality = HexStretch;
340       break;
341     case VTK_QUALITY_DIAGONAL:
342       HexQuality = HexDiagonal;
343       break;
344     case VTK_QUALITY_DIMENSION:
345       HexQuality = HexDimension;
346       break;
347     case VTK_QUALITY_ODDY:
348       HexQuality = HexOddy;
349       break;
350     case VTK_QUALITY_CONDITION:
351       HexQuality = HexCondition;
352       break;
353     case VTK_QUALITY_JACOBIAN:
354       HexQuality = HexJacobian;
355       break;
356     case VTK_QUALITY_SCALED_JACOBIAN:
357       HexQuality = HexScaledJacobian;
358       break;
359     case VTK_QUALITY_SHEAR:
360       HexQuality = HexShear;
361       break;
362     case VTK_QUALITY_SHAPE:
363       HexQuality = HexShape;
364       break;
365     case VTK_QUALITY_RELATIVE_SIZE_SQUARED:
366       HexQuality = HexRelativeSizeSquared;
367       break;
368     case VTK_QUALITY_SHAPE_AND_SIZE:
369       HexQuality = HexShapeAndSize;
370       break;
371     case VTK_QUALITY_SHEAR_AND_SIZE:
372       HexQuality = HexShearAndSize;
373       break;
374     case VTK_QUALITY_DISTORTION:
375       HexQuality = HexDistortion;
376       break;
377     default:
378       vtkWarningMacro("Bad HexQualityMeasure (" << this->GetTetQualityMeasure()
379                                                 << "), using MaxAspectFrobenius instead");
380       HexQuality = HexMaxAspectFrobenius;
381       break;
382   }
383 
384   out->ShallowCopy(in);
385 
386   if (this->SaveCellQuality)
387   {
388     quality = vtkDoubleArray::New();
389     if (this->CompatibilityMode)
390     {
391       if (this->Volume)
392       {
393         quality->SetNumberOfComponents(2);
394       }
395       else
396       {
397         quality->SetNumberOfComponents(1);
398       }
399     }
400     else
401     {
402       quality->SetNumberOfComponents(1);
403     }
404     quality->SetNumberOfTuples(N);
405     quality->SetName("Quality");
406     out->GetCellData()->AddArray(quality);
407     out->GetCellData()->SetActiveAttribute("Quality", vtkDataSetAttributes::SCALARS);
408     quality->Delete();
409 
410     if (!this->CompatibilityMode)
411     {
412       if (this->Volume)
413       {
414         volume = vtkDoubleArray::New();
415         volume->SetNumberOfComponents(1);
416         volume->SetNumberOfTuples(N);
417         volume->SetName("Volume");
418         out->GetCellData()->AddArray(volume);
419         volume->Delete();
420       }
421     }
422   }
423 
424   // These measures require the average area/volume for all cells of the same type in the mesh.
425   // Either use the hinted value (computed by a previous vtkMeshQuality filter) or compute it.
426   if (this->GetTriangleQualityMeasure() == VTK_QUALITY_RELATIVE_SIZE_SQUARED ||
427     this->GetTriangleQualityMeasure() == VTK_QUALITY_SHAPE_AND_SIZE ||
428     this->GetQuadQualityMeasure() == VTK_QUALITY_RELATIVE_SIZE_SQUARED ||
429     this->GetQuadQualityMeasure() == VTK_QUALITY_SHAPE_AND_SIZE ||
430     this->GetQuadQualityMeasure() == VTK_QUALITY_SHEAR_AND_SIZE ||
431     this->GetTetQualityMeasure() == VTK_QUALITY_RELATIVE_SIZE_SQUARED ||
432     this->GetTetQualityMeasure() == VTK_QUALITY_SHAPE_AND_SIZE ||
433     this->GetHexQualityMeasure() == VTK_QUALITY_RELATIVE_SIZE_SQUARED ||
434     this->GetHexQualityMeasure() == VTK_QUALITY_SHAPE_AND_SIZE ||
435     this->GetHexQualityMeasure() == VTK_QUALITY_SHEAR_AND_SIZE)
436   {
437     vtkDataArray* triAreaHint = in->GetFieldData()->GetArray("TriArea");
438     vtkDataArray* quadAreaHint = in->GetFieldData()->GetArray("QuadArea");
439     vtkDataArray* tetVolHint = in->GetFieldData()->GetArray("TetVolume");
440     vtkDataArray* hexVolHint = in->GetFieldData()->GetArray("HexVolume");
441 
442     double triAreaTuple[5];
443     double quadAreaTuple[5];
444     double tetVolTuple[5];
445     double hexVolTuple[5];
446 
447     if (triAreaHint && triAreaHint->GetNumberOfTuples() > 0 &&
448       triAreaHint->GetNumberOfComponents() == 5 && quadAreaHint &&
449       quadAreaHint->GetNumberOfTuples() > 0 && quadAreaHint->GetNumberOfComponents() == 5 &&
450       tetVolHint && tetVolHint->GetNumberOfTuples() > 0 &&
451       tetVolHint->GetNumberOfComponents() == 5 && hexVolHint &&
452       hexVolHint->GetNumberOfTuples() > 0 && hexVolHint->GetNumberOfComponents() == 5)
453     {
454       triAreaHint->GetTuple(0, triAreaTuple);
455       quadAreaHint->GetTuple(0, quadAreaTuple);
456       tetVolHint->GetTuple(0, tetVolTuple);
457       hexVolHint->GetTuple(0, hexVolTuple);
458       v_set_tri_size(triAreaTuple[1] / triAreaTuple[4]);
459       v_set_quad_size(quadAreaTuple[1] / quadAreaTuple[4]);
460       v_set_tet_size(tetVolTuple[1] / tetVolTuple[4]);
461       v_set_hex_size(hexVolTuple[1] / hexVolTuple[4]);
462     }
463     else
464     {
465       for (int i = 0; i < 5; ++i)
466       {
467         triAreaTuple[i] = 0;
468         quadAreaTuple[i] = 0;
469         tetVolTuple[i] = 0;
470         hexVolTuple[i] = 0;
471       }
472       for (vtkIdType c = 0; c < N; ++c)
473       {
474         double a, v; // area and volume
475         cell = out->GetCell(c);
476         switch (cell->GetCellType())
477         {
478           case VTK_TRIANGLE:
479             a = TriangleArea(cell);
480             if (a > triAreaTuple[2])
481             {
482               if (triAreaTuple[0] == triAreaTuple[2])
483               { // min == max => min has not been set
484                 triAreaTuple[0] = a;
485               }
486               triAreaTuple[2] = a;
487             }
488             else if (a < triAreaTuple[0])
489             {
490               triAreaTuple[0] = a;
491             }
492             triAreaTuple[1] += a;
493             triAreaTuple[3] += a * a;
494             ntri++;
495             break;
496           case VTK_QUAD:
497             a = QuadArea(cell);
498             if (a > quadAreaTuple[2])
499             {
500               if (quadAreaTuple[0] == quadAreaTuple[2])
501               { // min == max => min has not been set
502                 quadAreaTuple[0] = a;
503               }
504               quadAreaTuple[2] = a;
505             }
506             else if (a < quadAreaTuple[0])
507             {
508               quadAreaTuple[0] = a;
509             }
510             quadAreaTuple[1] += a;
511             quadAreaTuple[3] += a * a;
512             nqua++;
513             break;
514           case VTK_TETRA:
515             v = TetVolume(cell);
516             if (v > tetVolTuple[2])
517             {
518               if (tetVolTuple[0] == tetVolTuple[2])
519               { // min == max => min has not been set
520                 tetVolTuple[0] = v;
521               }
522               tetVolTuple[2] = v;
523             }
524             else if (v < tetVolTuple[0])
525             {
526               tetVolTuple[0] = v;
527             }
528             tetVolTuple[1] += v;
529             tetVolTuple[3] += v * v;
530             ntet++;
531             break;
532           case VTK_HEXAHEDRON:
533             v = HexVolume(cell);
534             if (v > hexVolTuple[2])
535             {
536               if (hexVolTuple[0] == hexVolTuple[2])
537               { // min == max => min has not been set
538                 hexVolTuple[0] = v;
539               }
540               hexVolTuple[2] = v;
541             }
542             else if (v < hexVolTuple[0])
543             {
544               hexVolTuple[0] = v;
545             }
546             hexVolTuple[1] += v;
547             hexVolTuple[3] += v * v;
548             nhex++;
549             break;
550         }
551       }
552       triAreaTuple[4] = ntri;
553       quadAreaTuple[4] = nqua;
554       tetVolTuple[4] = ntet;
555       hexVolTuple[4] = nhex;
556       v_set_tri_size(triAreaTuple[1] / triAreaTuple[4]);
557       v_set_quad_size(quadAreaTuple[1] / quadAreaTuple[4]);
558       v_set_tet_size(tetVolTuple[1] / tetVolTuple[4]);
559       v_set_hex_size(hexVolTuple[1] / hexVolTuple[4]);
560       progressNumer = 20;
561       progressDenom = 40.;
562       ntri = 0;
563       nqua = 0;
564       ntet = 0;
565       nhex = 0;
566 
567       // Save info as field data for downstream filters
568       triAreaHint = vtkDoubleArray::New();
569       triAreaHint->SetName("TriArea");
570       triAreaHint->SetNumberOfComponents(5);
571       triAreaHint->InsertNextTuple(triAreaTuple);
572       out->GetFieldData()->AddArray(triAreaHint);
573       triAreaHint->Delete();
574 
575       quadAreaHint = vtkDoubleArray::New();
576       quadAreaHint->SetName("QuadArea");
577       quadAreaHint->SetNumberOfComponents(5);
578       quadAreaHint->InsertNextTuple(quadAreaTuple);
579       out->GetFieldData()->AddArray(quadAreaHint);
580       quadAreaHint->Delete();
581 
582       tetVolHint = vtkDoubleArray::New();
583       tetVolHint->SetName("TetVolume");
584       tetVolHint->SetNumberOfComponents(5);
585       tetVolHint->InsertNextTuple(tetVolTuple);
586       out->GetFieldData()->AddArray(tetVolHint);
587       tetVolHint->Delete();
588 
589       hexVolHint = vtkDoubleArray::New();
590       hexVolHint->SetName("HexVolume");
591       hexVolHint->SetNumberOfComponents(5);
592       hexVolHint->InsertNextTuple(hexVolTuple);
593       out->GetFieldData()->AddArray(hexVolHint);
594       hexVolHint->Delete();
595     }
596   }
597 
598   int p;
599   vtkIdType c = 0;
600   vtkIdType sz = N / 20 + 1;
601   vtkIdType inner;
602   this->UpdateProgress(progressNumer / progressDenom + 0.01);
603   for (p = 0; p < 20; ++p)
604   {
605     for (inner = 0; (inner < sz && c < N); ++c, ++inner)
606     {
607       cell = out->GetCell(c);
608       V = 0.;
609       switch (cell->GetCellType())
610       {
611         case VTK_TRIANGLE:
612           if (this->CellNormals)
613             this->CellNormals->GetTuple(c, vtkMeshQuality::CurrentTriNormal);
614           q = TriangleQuality(cell);
615           if (q > qtriM)
616           {
617             if (qtrim > qtriM)
618             {
619               qtrim = q;
620             }
621             qtriM = q;
622           }
623           else if (q < qtrim)
624           {
625             qtrim = q;
626           }
627           Eqtri += q;
628           Eqtri2 += q * q;
629           ++ntri;
630           break;
631         case VTK_QUAD:
632           q = QuadQuality(cell);
633           if (q > qquaM)
634           {
635             if (qquam > qquaM)
636             {
637               qquam = q;
638             }
639             qquaM = q;
640           }
641           else if (q < qquam)
642           {
643             qquam = q;
644           }
645           Eqqua += q;
646           Eqqua2 += q * q;
647           ++nqua;
648           break;
649         case VTK_TETRA:
650           q = TetQuality(cell);
651           if (q > qtetM)
652           {
653             if (qtetm > qtetM)
654             {
655               qtetm = q;
656             }
657             qtetM = q;
658           }
659           else if (q < qtetm)
660           {
661             qtetm = q;
662           }
663           Eqtet += q;
664           Eqtet2 += q * q;
665           ++ntet;
666           if (this->Volume)
667           {
668             V = TetVolume(cell);
669             if (!this->CompatibilityMode)
670             {
671               volume->SetTuple1(0, V);
672             }
673           }
674           break;
675         case VTK_HEXAHEDRON:
676           q = HexQuality(cell);
677           if (q > qhexM)
678           {
679             if (qhexm > qhexM)
680             {
681               qhexm = q;
682             }
683             qhexM = q;
684           }
685           else if (q < qhexm)
686           {
687             qhexm = q;
688           }
689           Eqhex += q;
690           Eqhex2 += q * q;
691           ++nhex;
692           break;
693         default:
694           q = 0.;
695       }
696 
697       if (this->SaveCellQuality)
698       {
699         if (this->CompatibilityMode && this->Volume)
700         {
701           quality->SetTuple2(c, V, q);
702         }
703         else
704         {
705           quality->SetTuple1(c, q);
706         }
707       }
708     }
709     this->UpdateProgress(double(p + 1 + progressNumer) / progressDenom);
710   }
711 
712   if (ntri)
713   {
714     Eqtri /= static_cast<double>(ntri);
715     double multFactor = 1. / static_cast<double>(ntri > 1 ? ntri - 1 : ntri);
716     Eqtri2 = multFactor * (Eqtri2 - static_cast<double>(ntri) * Eqtri * Eqtri);
717   }
718   else
719   {
720     qtrim = Eqtri = qtriM = Eqtri2 = 0.;
721   }
722 
723   if (nqua)
724   {
725     Eqqua /= static_cast<double>(nqua);
726     double multFactor = 1. / static_cast<double>(nqua > 1 ? nqua - 1 : nqua);
727     Eqqua2 = multFactor * (Eqqua2 - static_cast<double>(nqua) * Eqqua * Eqqua);
728   }
729   else
730   {
731     qquam = Eqqua = qquaM = Eqqua2 = 0.;
732   }
733 
734   if (ntet)
735   {
736     Eqtet /= static_cast<double>(ntet);
737     double multFactor = 1. / static_cast<double>(ntet > 1 ? ntet - 1 : ntet);
738     Eqtet2 = multFactor * (Eqtet2 - static_cast<double>(ntet) * Eqtet * Eqtet);
739   }
740   else
741   {
742     qtetm = Eqtet = qtetM = Eqtet2 = 0.;
743   }
744 
745   if (nhex)
746   {
747     Eqhex /= static_cast<double>(nhex);
748     double multFactor = 1. / static_cast<double>(nhex > 1 ? nhex - 1 : nhex);
749     Eqhex2 = multFactor * (Eqhex2 - static_cast<double>(nhex) * Eqhex * Eqhex);
750   }
751   else
752   {
753     qhexm = Eqhex = qhexM = Eqhex2 = 0.;
754   }
755 
756   double tuple[5];
757   quality = vtkDoubleArray::New();
758   quality->SetName("Mesh Triangle Quality");
759   quality->SetNumberOfComponents(5);
760   tuple[0] = qtrim;
761   tuple[1] = Eqtri;
762   tuple[2] = qtriM;
763   tuple[3] = Eqtri2;
764   tuple[4] = ntri;
765   quality->InsertNextTuple(tuple);
766   out->GetFieldData()->AddArray(quality);
767   quality->Delete();
768 
769   quality = vtkDoubleArray::New();
770   quality->SetName("Mesh Quadrilateral Quality");
771   quality->SetNumberOfComponents(5);
772   tuple[0] = qquam;
773   tuple[1] = Eqqua;
774   tuple[2] = qquaM;
775   tuple[3] = Eqqua2;
776   tuple[4] = nqua;
777   quality->InsertNextTuple(tuple);
778   out->GetFieldData()->AddArray(quality);
779   quality->Delete();
780 
781   quality = vtkDoubleArray::New();
782   quality->SetName("Mesh Tetrahedron Quality");
783   quality->SetNumberOfComponents(5);
784   tuple[0] = qtetm;
785   tuple[1] = Eqtet;
786   tuple[2] = qtetM;
787   tuple[3] = Eqtet2;
788   tuple[4] = ntet;
789   quality->InsertNextTuple(tuple);
790   out->GetFieldData()->AddArray(quality);
791   quality->Delete();
792 
793   quality = vtkDoubleArray::New();
794   quality->SetName("Mesh Hexahedron Quality");
795   quality->SetNumberOfComponents(5);
796   tuple[0] = qhexm;
797   tuple[1] = Eqhex;
798   tuple[2] = qhexM;
799   tuple[3] = Eqhex2;
800   tuple[4] = nhex;
801   quality->InsertNextTuple(tuple);
802   out->GetFieldData()->AddArray(quality);
803   quality->Delete();
804 
805   return 1;
806 }
807 
GetCurrentTriangleNormal(double point[3],double normal[3])808 int vtkMeshQuality::GetCurrentTriangleNormal(double point[3], double normal[3])
809 {
810   // ignore the location where the normal should be evaluated.
811   (void)point;
812 
813   // copy the cell normal
814   for (int i = 0; i < 3; ++i)
815     normal[i] = vtkMeshQuality::CurrentTriNormal[i];
816   return 1;
817 }
818 
819 // Triangle quality metrics
820 
TriangleArea(vtkCell * cell)821 double vtkMeshQuality::TriangleArea(vtkCell* cell)
822 {
823   double pc[3][3];
824 
825   vtkPoints* p = cell->GetPoints();
826   p->GetPoint(0, pc[0]);
827   p->GetPoint(1, pc[1]);
828   p->GetPoint(2, pc[2]);
829 
830   return v_tri_area(3, pc);
831 }
832 
TriangleEdgeRatio(vtkCell * cell)833 double vtkMeshQuality::TriangleEdgeRatio(vtkCell* cell)
834 {
835   double pc[3][3];
836 
837   vtkPoints* p = cell->GetPoints();
838   p->GetPoint(0, pc[0]);
839   p->GetPoint(1, pc[1]);
840   p->GetPoint(2, pc[2]);
841 
842   return v_tri_edge_ratio(3, pc);
843 }
844 
TriangleAspectRatio(vtkCell * cell)845 double vtkMeshQuality::TriangleAspectRatio(vtkCell* cell)
846 {
847   double pc[3][3];
848 
849   vtkPoints* p = cell->GetPoints();
850   p->GetPoint(0, pc[0]);
851   p->GetPoint(1, pc[1]);
852   p->GetPoint(2, pc[2]);
853 
854   return v_tri_aspect_ratio(3, pc);
855 }
856 
TriangleRadiusRatio(vtkCell * cell)857 double vtkMeshQuality::TriangleRadiusRatio(vtkCell* cell)
858 {
859   double pc[3][3];
860 
861   vtkPoints* p = cell->GetPoints();
862   p->GetPoint(0, pc[0]);
863   p->GetPoint(1, pc[1]);
864   p->GetPoint(2, pc[2]);
865 
866   return v_tri_radius_ratio(3, pc);
867 }
868 
TriangleAspectFrobenius(vtkCell * cell)869 double vtkMeshQuality::TriangleAspectFrobenius(vtkCell* cell)
870 {
871   double pc[3][3];
872 
873   vtkPoints* p = cell->GetPoints();
874   p->GetPoint(0, pc[0]);
875   p->GetPoint(1, pc[1]);
876   p->GetPoint(2, pc[2]);
877 
878   return v_tri_aspect_frobenius(3, pc);
879 }
880 
TriangleMinAngle(vtkCell * cell)881 double vtkMeshQuality::TriangleMinAngle(vtkCell* cell)
882 {
883   double pc[3][3];
884 
885   vtkPoints* p = cell->GetPoints();
886   p->GetPoint(0, pc[0]);
887   p->GetPoint(1, pc[1]);
888   p->GetPoint(2, pc[2]);
889 
890   return v_tri_minimum_angle(3, pc);
891 }
892 
TriangleMaxAngle(vtkCell * cell)893 double vtkMeshQuality::TriangleMaxAngle(vtkCell* cell)
894 {
895   double pc[3][3];
896 
897   vtkPoints* p = cell->GetPoints();
898   p->GetPoint(0, pc[0]);
899   p->GetPoint(1, pc[1]);
900   p->GetPoint(2, pc[2]);
901 
902   return v_tri_maximum_angle(3, pc);
903 }
904 
TriangleCondition(vtkCell * cell)905 double vtkMeshQuality::TriangleCondition(vtkCell* cell)
906 {
907   double pc[3][3];
908 
909   vtkPoints* p = cell->GetPoints();
910   p->GetPoint(0, pc[0]);
911   p->GetPoint(1, pc[1]);
912   p->GetPoint(2, pc[2]);
913 
914   return v_tri_condition(3, pc);
915 }
916 
TriangleScaledJacobian(vtkCell * cell)917 double vtkMeshQuality::TriangleScaledJacobian(vtkCell* cell)
918 {
919   double pc[3][3];
920 
921   vtkPoints* p = cell->GetPoints();
922   p->GetPoint(0, pc[0]);
923   p->GetPoint(1, pc[1]);
924   p->GetPoint(2, pc[2]);
925 
926   return v_tri_scaled_jacobian(3, pc);
927 }
928 
TriangleRelativeSizeSquared(vtkCell * cell)929 double vtkMeshQuality::TriangleRelativeSizeSquared(vtkCell* cell)
930 {
931   double pc[3][3];
932 
933   vtkPoints* p = cell->GetPoints();
934   p->GetPoint(0, pc[0]);
935   p->GetPoint(1, pc[1]);
936   p->GetPoint(2, pc[2]);
937 
938   return v_tri_relative_size_squared(3, pc);
939 }
940 
TriangleShape(vtkCell * cell)941 double vtkMeshQuality::TriangleShape(vtkCell* cell)
942 {
943   double pc[3][3];
944 
945   vtkPoints* p = cell->GetPoints();
946   p->GetPoint(0, pc[0]);
947   p->GetPoint(1, pc[1]);
948   p->GetPoint(2, pc[2]);
949 
950   return v_tri_shape(3, pc);
951 }
952 
TriangleShapeAndSize(vtkCell * cell)953 double vtkMeshQuality::TriangleShapeAndSize(vtkCell* cell)
954 {
955   double pc[3][3];
956 
957   vtkPoints* p = cell->GetPoints();
958   p->GetPoint(0, pc[0]);
959   p->GetPoint(1, pc[1]);
960   p->GetPoint(2, pc[2]);
961 
962   return v_tri_shape_and_size(3, pc);
963 }
964 
TriangleDistortion(vtkCell * cell)965 double vtkMeshQuality::TriangleDistortion(vtkCell* cell)
966 {
967   double pc[3][3];
968 
969   vtkPoints* p = cell->GetPoints();
970   p->GetPoint(0, pc[0]);
971   p->GetPoint(1, pc[1]);
972   p->GetPoint(2, pc[2]);
973 
974   return v_tri_distortion(3, pc);
975 }
976 
977 // Quadrangle quality metrics
978 
QuadEdgeRatio(vtkCell * cell)979 double vtkMeshQuality::QuadEdgeRatio(vtkCell* cell)
980 {
981   double pc[4][3];
982 
983   vtkPoints* p = cell->GetPoints();
984   for (int i = 0; i < 4; ++i)
985     p->GetPoint(i, pc[i]);
986 
987   return v_quad_edge_ratio(4, pc);
988 }
989 
QuadAspectRatio(vtkCell * cell)990 double vtkMeshQuality::QuadAspectRatio(vtkCell* cell)
991 {
992   double pc[4][3];
993 
994   vtkPoints* p = cell->GetPoints();
995   for (int i = 0; i < 4; ++i)
996     p->GetPoint(i, pc[i]);
997 
998   return v_quad_aspect_ratio(4, pc);
999 }
1000 
QuadRadiusRatio(vtkCell * cell)1001 double vtkMeshQuality::QuadRadiusRatio(vtkCell* cell)
1002 {
1003   double pc[4][3];
1004 
1005   vtkPoints* p = cell->GetPoints();
1006   for (int i = 0; i < 4; ++i)
1007     p->GetPoint(i, pc[i]);
1008 
1009   return v_quad_radius_ratio(4, pc);
1010 }
1011 
QuadMedAspectFrobenius(vtkCell * cell)1012 double vtkMeshQuality::QuadMedAspectFrobenius(vtkCell* cell)
1013 {
1014   double pc[4][3];
1015 
1016   vtkPoints* p = cell->GetPoints();
1017   for (int i = 0; i < 4; ++i)
1018     p->GetPoint(i, pc[i]);
1019 
1020   return v_quad_med_aspect_frobenius(4, pc);
1021 }
1022 
QuadMaxAspectFrobenius(vtkCell * cell)1023 double vtkMeshQuality::QuadMaxAspectFrobenius(vtkCell* cell)
1024 {
1025   double pc[4][3];
1026 
1027   vtkPoints* p = cell->GetPoints();
1028   for (int i = 0; i < 4; ++i)
1029     p->GetPoint(i, pc[i]);
1030 
1031   return v_quad_max_aspect_frobenius(4, pc);
1032 }
1033 
QuadMinAngle(vtkCell * cell)1034 double vtkMeshQuality::QuadMinAngle(vtkCell* cell)
1035 {
1036   double pc[4][3];
1037 
1038   vtkPoints* p = cell->GetPoints();
1039   for (int i = 0; i < 4; ++i)
1040     p->GetPoint(i, pc[i]);
1041 
1042   return v_quad_minimum_angle(4, pc);
1043 }
1044 
QuadMaxEdgeRatios(vtkCell * cell)1045 double vtkMeshQuality::QuadMaxEdgeRatios(vtkCell* cell)
1046 {
1047   double pc[4][3];
1048 
1049   vtkPoints* p = cell->GetPoints();
1050   for (int i = 0; i < 4; ++i)
1051     p->GetPoint(i, pc[i]);
1052 
1053   return v_quad_max_edge_ratio(4, pc);
1054 }
1055 
QuadSkew(vtkCell * cell)1056 double vtkMeshQuality::QuadSkew(vtkCell* cell)
1057 {
1058   double pc[4][3];
1059 
1060   vtkPoints* p = cell->GetPoints();
1061   for (int i = 0; i < 4; ++i)
1062     p->GetPoint(i, pc[i]);
1063 
1064   return v_quad_skew(4, pc);
1065 }
1066 
QuadTaper(vtkCell * cell)1067 double vtkMeshQuality::QuadTaper(vtkCell* cell)
1068 {
1069   double pc[4][3];
1070 
1071   vtkPoints* p = cell->GetPoints();
1072   for (int i = 0; i < 4; ++i)
1073     p->GetPoint(i, pc[i]);
1074 
1075   return v_quad_taper(4, pc);
1076 }
1077 
QuadWarpage(vtkCell * cell)1078 double vtkMeshQuality::QuadWarpage(vtkCell* cell)
1079 {
1080   double pc[4][3];
1081 
1082   vtkPoints* p = cell->GetPoints();
1083   for (int i = 0; i < 4; ++i)
1084     p->GetPoint(i, pc[i]);
1085 
1086   return v_quad_warpage(4, pc);
1087 }
1088 
QuadArea(vtkCell * cell)1089 double vtkMeshQuality::QuadArea(vtkCell* cell)
1090 {
1091   double pc[4][3];
1092 
1093   vtkPoints* p = cell->GetPoints();
1094   for (int i = 0; i < 4; ++i)
1095     p->GetPoint(i, pc[i]);
1096 
1097   return v_quad_area(4, pc);
1098 }
1099 
QuadStretch(vtkCell * cell)1100 double vtkMeshQuality::QuadStretch(vtkCell* cell)
1101 {
1102   double pc[4][3];
1103 
1104   vtkPoints* p = cell->GetPoints();
1105   for (int i = 0; i < 4; ++i)
1106     p->GetPoint(i, pc[i]);
1107 
1108   return v_quad_stretch(4, pc);
1109 }
1110 
1111 #if 0
1112 // FIXME
1113 double vtkMeshQuality::QuadMinAngle( vtkCell* cell )
1114 {
1115   double pc[4][3];
1116 
1117   vtkPoints *p = cell->GetPoints();
1118   for ( int i = 0; i < 4; ++i )
1119     p->GetPoint( i, pc[i] );
1120 
1121   return v_quad_minimum_angle( 4, pc );
1122 }
1123 #endif // 0
1124 
QuadMaxAngle(vtkCell * cell)1125 double vtkMeshQuality::QuadMaxAngle(vtkCell* cell)
1126 {
1127   double pc[4][3];
1128 
1129   vtkPoints* p = cell->GetPoints();
1130   for (int i = 0; i < 4; ++i)
1131     p->GetPoint(i, pc[i]);
1132 
1133   return v_quad_maximum_angle(4, pc);
1134 }
1135 
QuadOddy(vtkCell * cell)1136 double vtkMeshQuality::QuadOddy(vtkCell* cell)
1137 {
1138   double pc[4][3];
1139 
1140   vtkPoints* p = cell->GetPoints();
1141   for (int i = 0; i < 4; ++i)
1142     p->GetPoint(i, pc[i]);
1143 
1144   return v_quad_oddy(4, pc);
1145 }
1146 
QuadCondition(vtkCell * cell)1147 double vtkMeshQuality::QuadCondition(vtkCell* cell)
1148 {
1149   double pc[4][3];
1150 
1151   vtkPoints* p = cell->GetPoints();
1152   for (int i = 0; i < 4; ++i)
1153     p->GetPoint(i, pc[i]);
1154 
1155   return v_quad_condition(4, pc);
1156 }
1157 
QuadJacobian(vtkCell * cell)1158 double vtkMeshQuality::QuadJacobian(vtkCell* cell)
1159 {
1160   double pc[4][3];
1161 
1162   vtkPoints* p = cell->GetPoints();
1163   for (int i = 0; i < 4; ++i)
1164     p->GetPoint(i, pc[i]);
1165 
1166   return v_quad_jacobian(4, pc);
1167 }
1168 
QuadScaledJacobian(vtkCell * cell)1169 double vtkMeshQuality::QuadScaledJacobian(vtkCell* cell)
1170 {
1171   double pc[4][3];
1172 
1173   vtkPoints* p = cell->GetPoints();
1174   for (int i = 0; i < 4; ++i)
1175     p->GetPoint(i, pc[i]);
1176 
1177   return v_quad_scaled_jacobian(4, pc);
1178 }
1179 
QuadShear(vtkCell * cell)1180 double vtkMeshQuality::QuadShear(vtkCell* cell)
1181 {
1182   double pc[4][3];
1183 
1184   vtkPoints* p = cell->GetPoints();
1185   for (int i = 0; i < 4; ++i)
1186     p->GetPoint(i, pc[i]);
1187 
1188   return v_quad_shear(4, pc);
1189 }
1190 
QuadShape(vtkCell * cell)1191 double vtkMeshQuality::QuadShape(vtkCell* cell)
1192 {
1193   double pc[4][3];
1194 
1195   vtkPoints* p = cell->GetPoints();
1196   for (int i = 0; i < 4; ++i)
1197     p->GetPoint(i, pc[i]);
1198 
1199   return v_quad_shape(4, pc);
1200 }
1201 
QuadRelativeSizeSquared(vtkCell * cell)1202 double vtkMeshQuality::QuadRelativeSizeSquared(vtkCell* cell)
1203 {
1204   double pc[4][3];
1205 
1206   vtkPoints* p = cell->GetPoints();
1207   for (int i = 0; i < 4; ++i)
1208     p->GetPoint(i, pc[i]);
1209 
1210   return v_quad_relative_size_squared(4, pc);
1211 }
1212 
QuadShapeAndSize(vtkCell * cell)1213 double vtkMeshQuality::QuadShapeAndSize(vtkCell* cell)
1214 {
1215   double pc[4][3];
1216 
1217   vtkPoints* p = cell->GetPoints();
1218   for (int i = 0; i < 4; ++i)
1219     p->GetPoint(i, pc[i]);
1220 
1221   return v_quad_shape_and_size(4, pc);
1222 }
1223 
QuadShearAndSize(vtkCell * cell)1224 double vtkMeshQuality::QuadShearAndSize(vtkCell* cell)
1225 {
1226   double pc[4][3];
1227 
1228   vtkPoints* p = cell->GetPoints();
1229   for (int i = 0; i < 4; ++i)
1230     p->GetPoint(i, pc[i]);
1231 
1232   return v_quad_shear_and_size(4, pc);
1233 }
1234 
QuadDistortion(vtkCell * cell)1235 double vtkMeshQuality::QuadDistortion(vtkCell* cell)
1236 {
1237   double pc[4][3];
1238 
1239   vtkPoints* p = cell->GetPoints();
1240   for (int i = 0; i < 4; ++i)
1241     p->GetPoint(i, pc[i]);
1242 
1243   return v_quad_distortion(4, pc);
1244 }
1245 
1246 // Volume of a tetrahedron, for compatibility with the original vtkMeshQuality
1247 
TetVolume(vtkCell * cell)1248 double TetVolume(vtkCell* cell)
1249 {
1250   double x0[3];
1251   double x1[3];
1252   double x2[3];
1253   double x3[3];
1254   cell->Points->GetPoint(0, x0);
1255   cell->Points->GetPoint(1, x1);
1256   cell->Points->GetPoint(2, x2);
1257   cell->Points->GetPoint(3, x3);
1258   return vtkTetra::ComputeVolume(x0, x1, x2, x3);
1259 }
1260 
1261 // Tetrahedral quality metrics
1262 
TetEdgeRatio(vtkCell * cell)1263 double vtkMeshQuality::TetEdgeRatio(vtkCell* cell)
1264 {
1265   double pc[4][3];
1266 
1267   vtkPoints* p = cell->GetPoints();
1268   for (int i = 0; i < 4; ++i)
1269     p->GetPoint(i, pc[i]);
1270 
1271   return v_tet_edge_ratio(4, pc);
1272 }
1273 
TetAspectRatio(vtkCell * cell)1274 double vtkMeshQuality::TetAspectRatio(vtkCell* cell)
1275 {
1276   double pc[4][3];
1277 
1278   vtkPoints* p = cell->GetPoints();
1279   for (int i = 0; i < 4; ++i)
1280     p->GetPoint(i, pc[i]);
1281 
1282   return v_tet_aspect_ratio(4, pc);
1283 }
1284 
TetRadiusRatio(vtkCell * cell)1285 double vtkMeshQuality::TetRadiusRatio(vtkCell* cell)
1286 {
1287   double pc[4][3];
1288 
1289   vtkPoints* p = cell->GetPoints();
1290   for (int i = 0; i < 4; ++i)
1291     p->GetPoint(i, pc[i]);
1292 
1293   return v_tet_radius_ratio(4, pc);
1294 }
1295 
TetAspectBeta(vtkCell * cell)1296 double vtkMeshQuality::TetAspectBeta(vtkCell* cell)
1297 {
1298   double pc[4][3];
1299 
1300   vtkPoints* p = cell->GetPoints();
1301   for (int i = 0; i < 4; ++i)
1302     p->GetPoint(i, pc[i]);
1303 
1304   return v_tet_aspect_beta(4, pc);
1305 }
1306 
TetAspectFrobenius(vtkCell * cell)1307 double vtkMeshQuality::TetAspectFrobenius(vtkCell* cell)
1308 {
1309   double pc[4][3];
1310 
1311   vtkPoints* p = cell->GetPoints();
1312   for (int i = 0; i < 4; ++i)
1313     p->GetPoint(i, pc[i]);
1314 
1315   return v_tet_aspect_frobenius(4, pc);
1316 }
1317 
TetMinAngle(vtkCell * cell)1318 double vtkMeshQuality::TetMinAngle(vtkCell* cell)
1319 {
1320   double pc[4][3];
1321 
1322   vtkPoints* p = cell->GetPoints();
1323   for (int i = 0; i < 4; ++i)
1324     p->GetPoint(i, pc[i]);
1325 
1326   return v_tet_minimum_angle(4, pc);
1327 }
1328 
TetCollapseRatio(vtkCell * cell)1329 double vtkMeshQuality::TetCollapseRatio(vtkCell* cell)
1330 {
1331   double pc[4][3];
1332 
1333   vtkPoints* p = cell->GetPoints();
1334   for (int i = 0; i < 4; ++i)
1335     p->GetPoint(i, pc[i]);
1336 
1337   return v_tet_collapse_ratio(4, pc);
1338 }
1339 
TetAspectGamma(vtkCell * cell)1340 double vtkMeshQuality::TetAspectGamma(vtkCell* cell)
1341 {
1342   double pc[4][3];
1343 
1344   vtkPoints* p = cell->GetPoints();
1345   for (int i = 0; i < 4; ++i)
1346     p->GetPoint(i, pc[i]);
1347 
1348   return v_tet_aspect_gamma(4, pc);
1349 }
1350 
TetVolume(vtkCell * cell)1351 double vtkMeshQuality::TetVolume(vtkCell* cell)
1352 {
1353   double pc[4][3];
1354 
1355   vtkPoints* p = cell->GetPoints();
1356   for (int i = 0; i < 4; ++i)
1357     p->GetPoint(i, pc[i]);
1358 
1359   return v_tet_volume(4, pc);
1360 }
1361 
TetCondition(vtkCell * cell)1362 double vtkMeshQuality::TetCondition(vtkCell* cell)
1363 {
1364   double pc[4][3];
1365 
1366   vtkPoints* p = cell->GetPoints();
1367   for (int i = 0; i < 4; ++i)
1368     p->GetPoint(i, pc[i]);
1369 
1370   return v_tet_condition(4, pc);
1371 }
1372 
TetJacobian(vtkCell * cell)1373 double vtkMeshQuality::TetJacobian(vtkCell* cell)
1374 {
1375   double pc[4][3];
1376 
1377   vtkPoints* p = cell->GetPoints();
1378   for (int i = 0; i < 4; ++i)
1379     p->GetPoint(i, pc[i]);
1380 
1381   return v_tet_jacobian(4, pc);
1382 }
1383 
TetScaledJacobian(vtkCell * cell)1384 double vtkMeshQuality::TetScaledJacobian(vtkCell* cell)
1385 {
1386   double pc[4][3];
1387 
1388   vtkPoints* p = cell->GetPoints();
1389   for (int i = 0; i < 4; ++i)
1390     p->GetPoint(i, pc[i]);
1391 
1392   return v_tet_scaled_jacobian(4, pc);
1393 }
1394 
TetShape(vtkCell * cell)1395 double vtkMeshQuality::TetShape(vtkCell* cell)
1396 {
1397   double pc[4][3];
1398 
1399   vtkPoints* p = cell->GetPoints();
1400   for (int i = 0; i < 4; ++i)
1401     p->GetPoint(i, pc[i]);
1402 
1403   return v_tet_shape(4, pc);
1404 }
1405 
TetRelativeSizeSquared(vtkCell * cell)1406 double vtkMeshQuality::TetRelativeSizeSquared(vtkCell* cell)
1407 {
1408   double pc[4][3];
1409 
1410   vtkPoints* p = cell->GetPoints();
1411   for (int i = 0; i < 4; ++i)
1412     p->GetPoint(i, pc[i]);
1413 
1414   return v_tet_relative_size_squared(4, pc);
1415 }
1416 
TetShapeandSize(vtkCell * cell)1417 double vtkMeshQuality::TetShapeandSize(vtkCell* cell)
1418 {
1419   double pc[4][3];
1420 
1421   vtkPoints* p = cell->GetPoints();
1422   for (int i = 0; i < 4; ++i)
1423     p->GetPoint(i, pc[i]);
1424 
1425   return v_tet_shape_and_size(4, pc);
1426 }
1427 
TetDistortion(vtkCell * cell)1428 double vtkMeshQuality::TetDistortion(vtkCell* cell)
1429 {
1430   double pc[4][3];
1431 
1432   vtkPoints* p = cell->GetPoints();
1433   for (int i = 0; i < 4; ++i)
1434     p->GetPoint(i, pc[i]);
1435 
1436   return v_tet_distortion(4, pc);
1437 }
1438 
1439 // Hexahedral quality metrics
1440 
HexEdgeRatio(vtkCell * cell)1441 double vtkMeshQuality::HexEdgeRatio(vtkCell* cell)
1442 {
1443   double pc[8][3];
1444   vtkPoints* p = cell->GetPoints();
1445 
1446   for (int i = 0; i < 8; ++i)
1447     p->GetPoint(i, pc[i]);
1448 
1449   return v_hex_edge_ratio(8, pc);
1450 }
1451 
HexMedAspectFrobenius(vtkCell * cell)1452 double vtkMeshQuality::HexMedAspectFrobenius(vtkCell* cell)
1453 {
1454   double pc[8][3];
1455 
1456   vtkPoints* p = cell->GetPoints();
1457   for (int i = 0; i < 8; ++i)
1458     p->GetPoint(i, pc[i]);
1459 
1460   return v_hex_med_aspect_frobenius(8, pc);
1461 }
1462 
HexMaxAspectFrobenius(vtkCell * cell)1463 double vtkMeshQuality::HexMaxAspectFrobenius(vtkCell* cell)
1464 {
1465   double pc[8][3];
1466 
1467   vtkPoints* p = cell->GetPoints();
1468   for (int i = 0; i < 8; ++i)
1469     p->GetPoint(i, pc[i]);
1470 
1471   return v_hex_max_aspect_frobenius(8, pc);
1472 }
1473 
HexMaxEdgeRatio(vtkCell * cell)1474 double vtkMeshQuality::HexMaxEdgeRatio(vtkCell* cell)
1475 {
1476   double pc[8][3];
1477 
1478   vtkPoints* p = cell->GetPoints();
1479   for (int i = 0; i < 8; ++i)
1480     p->GetPoint(i, pc[i]);
1481 
1482   return v_hex_max_edge_ratio(8, pc);
1483 }
1484 
HexSkew(vtkCell * cell)1485 double vtkMeshQuality::HexSkew(vtkCell* cell)
1486 {
1487   double pc[8][3];
1488   vtkPoints* p = cell->GetPoints();
1489 
1490   for (int i = 0; i < 8; ++i)
1491     p->GetPoint(i, pc[i]);
1492 
1493   return v_hex_skew(8, pc);
1494 }
1495 
HexTaper(vtkCell * cell)1496 double vtkMeshQuality::HexTaper(vtkCell* cell)
1497 {
1498   double pc[8][3];
1499   vtkPoints* p = cell->GetPoints();
1500 
1501   for (int i = 0; i < 8; ++i)
1502     p->GetPoint(i, pc[i]);
1503 
1504   return v_hex_taper(8, pc);
1505 }
1506 
HexVolume(vtkCell * cell)1507 double vtkMeshQuality::HexVolume(vtkCell* cell)
1508 {
1509   double pc[8][3];
1510   vtkPoints* p = cell->GetPoints();
1511 
1512   for (int i = 0; i < 8; ++i)
1513     p->GetPoint(i, pc[i]);
1514 
1515   return v_hex_volume(8, pc);
1516 }
1517 
HexStretch(vtkCell * cell)1518 double vtkMeshQuality::HexStretch(vtkCell* cell)
1519 {
1520   double pc[8][3];
1521   vtkPoints* p = cell->GetPoints();
1522 
1523   for (int i = 0; i < 8; ++i)
1524     p->GetPoint(i, pc[i]);
1525 
1526   return v_hex_stretch(8, pc);
1527 }
1528 
HexDiagonal(vtkCell * cell)1529 double vtkMeshQuality::HexDiagonal(vtkCell* cell)
1530 {
1531   double pc[8][3];
1532   vtkPoints* p = cell->GetPoints();
1533 
1534   for (int i = 0; i < 8; ++i)
1535     p->GetPoint(i, pc[i]);
1536 
1537   return v_hex_diagonal(8, pc);
1538 }
1539 
HexDimension(vtkCell * cell)1540 double vtkMeshQuality::HexDimension(vtkCell* cell)
1541 {
1542   double pc[8][3];
1543   vtkPoints* p = cell->GetPoints();
1544 
1545   for (int i = 0; i < 8; ++i)
1546     p->GetPoint(i, pc[i]);
1547 
1548   return v_hex_dimension(8, pc);
1549 }
1550 
HexOddy(vtkCell * cell)1551 double vtkMeshQuality::HexOddy(vtkCell* cell)
1552 {
1553   double pc[8][3];
1554   vtkPoints* p = cell->GetPoints();
1555 
1556   for (int i = 0; i < 8; ++i)
1557     p->GetPoint(i, pc[i]);
1558 
1559   return v_hex_oddy(8, pc);
1560 }
1561 
HexCondition(vtkCell * cell)1562 double vtkMeshQuality::HexCondition(vtkCell* cell)
1563 {
1564   double pc[8][3];
1565   vtkPoints* p = cell->GetPoints();
1566 
1567   for (int i = 0; i < 8; ++i)
1568     p->GetPoint(i, pc[i]);
1569 
1570   return v_hex_condition(8, pc);
1571 }
1572 
HexJacobian(vtkCell * cell)1573 double vtkMeshQuality::HexJacobian(vtkCell* cell)
1574 {
1575   double pc[8][3];
1576   vtkPoints* p = cell->GetPoints();
1577 
1578   for (int i = 0; i < 8; ++i)
1579     p->GetPoint(i, pc[i]);
1580 
1581   return v_hex_jacobian(8, pc);
1582 }
1583 
HexScaledJacobian(vtkCell * cell)1584 double vtkMeshQuality::HexScaledJacobian(vtkCell* cell)
1585 {
1586   double pc[8][3];
1587   vtkPoints* p = cell->GetPoints();
1588 
1589   for (int i = 0; i < 8; ++i)
1590     p->GetPoint(i, pc[i]);
1591 
1592   return v_hex_scaled_jacobian(8, pc);
1593 }
1594 
HexShear(vtkCell * cell)1595 double vtkMeshQuality::HexShear(vtkCell* cell)
1596 {
1597   double pc[8][3];
1598   vtkPoints* p = cell->GetPoints();
1599 
1600   for (int i = 0; i < 8; ++i)
1601     p->GetPoint(i, pc[i]);
1602 
1603   return v_hex_shear(8, pc);
1604 }
1605 
HexShape(vtkCell * cell)1606 double vtkMeshQuality::HexShape(vtkCell* cell)
1607 {
1608   double pc[8][3];
1609   vtkPoints* p = cell->GetPoints();
1610 
1611   for (int i = 0; i < 8; ++i)
1612     p->GetPoint(i, pc[i]);
1613 
1614   return v_hex_shape(8, pc);
1615 }
1616 
HexRelativeSizeSquared(vtkCell * cell)1617 double vtkMeshQuality::HexRelativeSizeSquared(vtkCell* cell)
1618 {
1619   double pc[8][3];
1620   vtkPoints* p = cell->GetPoints();
1621 
1622   for (int i = 0; i < 8; ++i)
1623     p->GetPoint(i, pc[i]);
1624 
1625   return v_hex_relative_size_squared(8, pc);
1626 }
1627 
HexShapeAndSize(vtkCell * cell)1628 double vtkMeshQuality::HexShapeAndSize(vtkCell* cell)
1629 {
1630   double pc[8][3];
1631   vtkPoints* p = cell->GetPoints();
1632 
1633   for (int i = 0; i < 8; ++i)
1634     p->GetPoint(i, pc[i]);
1635 
1636   return v_hex_shape_and_size(8, pc);
1637 }
1638 
HexShearAndSize(vtkCell * cell)1639 double vtkMeshQuality::HexShearAndSize(vtkCell* cell)
1640 {
1641   double pc[8][3];
1642   vtkPoints* p = cell->GetPoints();
1643 
1644   for (int i = 0; i < 8; ++i)
1645     p->GetPoint(i, pc[i]);
1646 
1647   return v_hex_shear_and_size(8, pc);
1648 }
1649 
HexDistortion(vtkCell * cell)1650 double vtkMeshQuality::HexDistortion(vtkCell* cell)
1651 {
1652   double pc[8][3];
1653   vtkPoints* p = cell->GetPoints();
1654 
1655   for (int i = 0; i < 8; ++i)
1656     p->GetPoint(i, pc[i]);
1657 
1658   return v_hex_distortion(8, pc);
1659 }
1660