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