Salome HOME
[bos #23982] EDF 22984 - aspect ratio of hexa
[modules/smesh.git] / src / Controls / SMESH_Controls.cxx
1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 // Lesser General Public License for more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22
23 #include "SMESH_ControlsDef.hxx"
24
25 #include "SMDS_BallElement.hxx"
26 #include "SMDS_FacePosition.hxx"
27 #include "SMDS_Iterator.hxx"
28 #include "SMDS_Mesh.hxx"
29 #include "SMDS_MeshElement.hxx"
30 #include "SMDS_MeshNode.hxx"
31 #include "SMDS_VolumeTool.hxx"
32 #include "SMESHDS_GroupBase.hxx"
33 #include "SMESHDS_GroupOnFilter.hxx"
34 #include "SMESHDS_Mesh.hxx"
35 #include "SMESH_MeshAlgos.hxx"
36 #include "SMESH_OctreeNode.hxx"
37
38 #include <GEOMUtils.hxx>
39 #include <Basics_Utils.hxx>
40
41 #include <BRepAdaptor_Surface.hxx>
42 #include <BRepBndLib.hxx>
43 #include <BRepBuilderAPI_Copy.hxx>
44 #include <BRepClass3d_SolidClassifier.hxx>
45 #include <BRepClass_FaceClassifier.hxx>
46 #include <BRep_Tool.hxx>
47 #include <GeomLib_IsPlanarSurface.hxx>
48 #include <Geom_CylindricalSurface.hxx>
49 #include <Geom_Plane.hxx>
50 #include <Geom_Surface.hxx>
51 #include <NCollection_Map.hxx>
52 #include <Precision.hxx>
53 #include <ShapeAnalysis_Surface.hxx>
54 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
55 #include <TColStd_MapOfInteger.hxx>
56 #include <TColStd_SequenceOfAsciiString.hxx>
57 #include <TColgp_Array1OfXYZ.hxx>
58 #include <TopAbs.hxx>
59 #include <TopExp.hxx>
60 #include <TopoDS.hxx>
61 #include <TopoDS_Edge.hxx>
62 #include <TopoDS_Face.hxx>
63 #include <TopoDS_Iterator.hxx>
64 #include <TopoDS_Shape.hxx>
65 #include <TopoDS_Vertex.hxx>
66 #include <gp_Ax3.hxx>
67 #include <gp_Cylinder.hxx>
68 #include <gp_Dir.hxx>
69 #include <gp_Pln.hxx>
70 #include <gp_Pnt.hxx>
71 #include <gp_Vec.hxx>
72 #include <gp_XYZ.hxx>
73
74 #include <vtkMeshQuality.h>
75
76 #include <set>
77 #include <limits>
78
79 /*
80                             AUXILIARY METHODS
81 */
82
83 namespace {
84
85   const double theEps = 1e-100;
86   const double theInf = 1e+100;
87
88   inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode )
89   {
90     return gp_XYZ(aNode->X(), aNode->Y(), aNode->Z() );
91   }
92
93   inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
94   {
95     gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
96
97     return v1.Magnitude() < gp::Resolution() ||
98       v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
99   }
100
101   inline double getCos2( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
102   {
103     gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
104     double dot = v1 * v2, len1 = v1.SquareMagnitude(), len2 = v2.SquareMagnitude();
105
106     return ( dot < 0 || len1 < gp::Resolution() || len2 < gp::Resolution() ? -1 :
107              dot * dot / len1 / len2 );
108   }
109
110   inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
111   {
112     gp_Vec aVec1( P2 - P1 );
113     gp_Vec aVec2( P3 - P1 );
114     return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
115   }
116
117   inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
118   {
119     return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
120   }
121
122
123
124   inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
125   {
126     double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
127     return aDist;
128   }
129
130   int getNbMultiConnection( const SMDS_Mesh* theMesh, const int theId )
131   {
132     if ( theMesh == 0 )
133       return 0;
134
135     const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
136     if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge/* || anEdge->NbNodes() != 2 */)
137       return 0;
138
139     // for each pair of nodes in anEdge (there are 2 pairs in a quadratic edge)
140     // count elements containing both nodes of the pair.
141     // Note that there may be such cases for a quadratic edge (a horizontal line):
142     //
143     //  Case 1          Case 2
144     //  |     |      |        |      |
145     //  |     |      |        |      |
146     //  +-----+------+  +-----+------+
147     //  |            |  |            |
148     //  |            |  |            |
149     // result should be 2 in both cases
150     //
151     int aResult0 = 0, aResult1 = 0;
152     // last node, it is a medium one in a quadratic edge
153     const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 );
154     const SMDS_MeshNode*    aNode0 = anEdge->GetNode( 0 );
155     const SMDS_MeshNode*    aNode1 = anEdge->GetNode( 1 );
156     if ( aNode1 == aLastNode ) aNode1 = 0;
157
158     SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator();
159     while( anElemIter->more() ) {
160       const SMDS_MeshElement* anElem = anElemIter->next();
161       if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
162         SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
163         while ( anIter->more() ) {
164           if ( const SMDS_MeshElement* anElemNode = anIter->next() ) {
165             if ( anElemNode == aNode0 ) {
166               aResult0++;
167               if ( !aNode1 ) break; // not a quadratic edge
168             }
169             else if ( anElemNode == aNode1 )
170               aResult1++;
171           }
172         }
173       }
174     }
175     int aResult = std::max ( aResult0, aResult1 );
176
177     return aResult;
178   }
179
180   gp_XYZ getNormale( const SMDS_MeshFace* theFace, bool* ok=0 )
181   {
182     int aNbNode = theFace->NbNodes();
183
184     gp_XYZ q1 = gpXYZ( theFace->GetNode(1)) - gpXYZ( theFace->GetNode(0));
185     gp_XYZ q2 = gpXYZ( theFace->GetNode(2)) - gpXYZ( theFace->GetNode(0));
186     gp_XYZ n  = q1 ^ q2;
187     if ( aNbNode > 3 ) {
188       gp_XYZ q3 = gpXYZ( theFace->GetNode(3)) - gpXYZ( theFace->GetNode(0));
189       n += q2 ^ q3;
190     }
191     double len = n.Modulus();
192     bool zeroLen = ( len <= std::numeric_limits<double>::min());
193     if ( !zeroLen )
194       n /= len;
195
196     if (ok) *ok = !zeroLen;
197
198     return n;
199   }
200 }
201
202
203
204 using namespace SMESH::Controls;
205
206 /*
207  *                               FUNCTORS
208  */
209
210 //================================================================================
211 /*
212   Class       : NumericalFunctor
213   Description : Base class for numerical functors
214 */
215 //================================================================================
216
217 NumericalFunctor::NumericalFunctor():
218   myMesh(NULL)
219 {
220   myPrecision = -1;
221 }
222
223 void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh )
224 {
225   myMesh = theMesh;
226 }
227
228 bool NumericalFunctor::GetPoints(const int       theId,
229                                  TSequenceOfXYZ& theRes ) const
230 {
231   theRes.clear();
232
233   if ( myMesh == 0 )
234     return false;
235
236   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
237   if ( !IsApplicable( anElem ))
238     return false;
239
240   return GetPoints( anElem, theRes );
241 }
242
243 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
244                                  TSequenceOfXYZ&         theRes )
245 {
246   theRes.clear();
247
248   if ( anElem == 0 )
249     return false;
250
251   theRes.reserve( anElem->NbNodes() );
252   theRes.setElement( anElem );
253
254   // Get nodes of the element
255   SMDS_NodeIteratorPtr anIter= anElem->interlacedNodesIterator();
256   if ( anIter ) {
257     SMESH_NodeXYZ p;
258     while( anIter->more() ) {
259       if ( p.Set( anIter->next() ))
260         theRes.push_back( p );
261     }
262   }
263
264   return true;
265 }
266
267 long  NumericalFunctor::GetPrecision() const
268 {
269   return myPrecision;
270 }
271
272 void  NumericalFunctor::SetPrecision( const long thePrecision )
273 {
274   myPrecision = thePrecision;
275   myPrecisionValue = pow( 10., (double)( myPrecision ) );
276 }
277
278 double NumericalFunctor::GetValue( long theId )
279 {
280   double aVal = 0;
281
282   myCurrElement = myMesh->FindElement( theId );
283
284   TSequenceOfXYZ P;
285   if ( GetPoints( theId, P )) // elem type is checked here
286     aVal = Round( GetValue( P ));
287
288   return aVal;
289 }
290
291 double NumericalFunctor::Round( const double & aVal )
292 {
293   return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal;
294 }
295
296 //================================================================================
297 /*!
298  * \brief Return true if a value can be computed for a given element.
299  *        Some NumericalFunctor's are meaningful for elements of a certain
300  *        geometry only.
301  */
302 //================================================================================
303
304 bool NumericalFunctor::IsApplicable( const SMDS_MeshElement* element ) const
305 {
306   return element && element->GetType() == this->GetType();
307 }
308
309 bool NumericalFunctor::IsApplicable( long theElementId ) const
310 {
311   return IsApplicable( myMesh->FindElement( theElementId ));
312 }
313
314 //================================================================================
315 /*!
316  * \brief Return histogram of functor values
317  *  \param nbIntervals - number of intervals
318  *  \param nbEvents - number of mesh elements having values within i-th interval
319  *  \param funValues - boundaries of intervals
320  *  \param elements - elements to check vulue of; empty list means "of all"
321  *  \param minmax - boundaries of diapason of values to divide into intervals
322  */
323 //================================================================================
324
325 void NumericalFunctor::GetHistogram(int                     nbIntervals,
326                                     std::vector<int>&       nbEvents,
327                                     std::vector<double>&    funValues,
328                                     const std::vector<int>& elements,
329                                     const double*           minmax,
330                                     const bool              isLogarithmic)
331 {
332   if ( nbIntervals < 1 ||
333        !myMesh ||
334        !myMesh->GetMeshInfo().NbElements( GetType() ))
335     return;
336   nbEvents.resize( nbIntervals, 0 );
337   funValues.resize( nbIntervals+1 );
338
339   // get all values sorted
340   std::multiset< double > values;
341   if ( elements.empty() )
342   {
343     SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator( GetType() );
344     while ( elemIt->more() )
345       values.insert( GetValue( elemIt->next()->GetID() ));
346   }
347   else
348   {
349     std::vector<int>::const_iterator id = elements.begin();
350     for ( ; id != elements.end(); ++id )
351       values.insert( GetValue( *id ));
352   }
353
354   if ( minmax )
355   {
356     funValues[0] = minmax[0];
357     funValues[nbIntervals] = minmax[1];
358   }
359   else
360   {
361     funValues[0] = *values.begin();
362     funValues[nbIntervals] = *values.rbegin();
363   }
364   // case nbIntervals == 1
365   if ( nbIntervals == 1 )
366   {
367     nbEvents[0] = values.size();
368     return;
369   }
370   // case of 1 value
371   if (funValues.front() == funValues.back())
372   {
373     nbEvents.resize( 1 );
374     nbEvents[0] = values.size();
375     funValues[1] = funValues.back();
376     funValues.resize( 2 );
377   }
378   // generic case
379   std::multiset< double >::iterator min = values.begin(), max;
380   for ( int i = 0; i < nbIntervals; ++i )
381   {
382     // find end value of i-th interval
383     double r = (i+1) / double(nbIntervals);
384     if (isLogarithmic && funValues.front() > 1e-07 && funValues.back() > 1e-07) {
385       double logmin = log10(funValues.front());
386       double lval = logmin + r * (log10(funValues.back()) - logmin);
387       funValues[i+1] = pow(10.0, lval);
388     }
389     else {
390       funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
391     }
392
393     // count values in the i-th interval if there are any
394     if ( min != values.end() && *min <= funValues[i+1] )
395     {
396       // find the first value out of the interval
397       max = values.upper_bound( funValues[i+1] ); // max is greater than funValues[i+1], or end()
398       nbEvents[i] = std::distance( min, max );
399       min = max;
400     }
401   }
402   // add values larger than minmax[1]
403   nbEvents.back() += std::distance( min, values.end() );
404 }
405
406 //=======================================================================
407 /*
408   Class       : Volume
409   Description : Functor calculating volume of a 3D element
410 */
411 //================================================================================
412
413 double Volume::GetValue( long theElementId )
414 {
415   if ( theElementId && myMesh ) {
416     SMDS_VolumeTool aVolumeTool;
417     if ( aVolumeTool.Set( myMesh->FindElement( theElementId )))
418       return aVolumeTool.GetSize();
419   }
420   return 0;
421 }
422
423 double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const
424 {
425   return Value;
426 }
427
428 SMDSAbs_ElementType Volume::GetType() const
429 {
430   return SMDSAbs_Volume;
431 }
432
433 //=======================================================================
434 /*
435   Class       : MaxElementLength2D
436   Description : Functor calculating maximum length of 2D element
437 */
438 //================================================================================
439
440 double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P )
441 {
442   if(P.size() == 0)
443     return 0.;
444   double aVal = 0;
445   int len = P.size();
446   if( len == 3 ) { // triangles
447     double L1 = getDistance(P( 1 ),P( 2 ));
448     double L2 = getDistance(P( 2 ),P( 3 ));
449     double L3 = getDistance(P( 3 ),P( 1 ));
450     aVal = Max(L1,Max(L2,L3));
451   }
452   else if( len == 4 ) { // quadrangles
453     double L1 = getDistance(P( 1 ),P( 2 ));
454     double L2 = getDistance(P( 2 ),P( 3 ));
455     double L3 = getDistance(P( 3 ),P( 4 ));
456     double L4 = getDistance(P( 4 ),P( 1 ));
457     double D1 = getDistance(P( 1 ),P( 3 ));
458     double D2 = getDistance(P( 2 ),P( 4 ));
459     aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
460   }
461   else if( len == 6 ) { // quadratic triangles
462     double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
463     double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
464     double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
465     aVal = Max(L1,Max(L2,L3));
466   }
467   else if( len == 8 || len == 9 ) { // quadratic quadrangles
468     double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
469     double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
470     double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
471     double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
472     double D1 = getDistance(P( 1 ),P( 5 ));
473     double D2 = getDistance(P( 3 ),P( 7 ));
474     aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
475   }
476   // Diagonals are undefined for concave polygons
477   // else if ( P.getElementEntity() == SMDSEntity_Quad_Polygon && P.size() > 2 ) // quad polygon
478   // {
479   //   // sides
480   //   aVal = getDistance( P( 1 ), P( P.size() )) + getDistance( P( P.size() ), P( P.size()-1 ));
481   //   for ( size_t i = 1; i < P.size()-1; i += 2 )
482   //   {
483   //     double L = getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 ));
484   //     aVal = Max( aVal, L );
485   //   }
486   //   // diagonals
487   //   for ( int i = P.size()-5; i > 0; i -= 2 )
488   //     for ( int j = i + 4; j < P.size() + i - 2; i += 2 )
489   //     {
490   //       double D = getDistance( P( i ), P( j ));
491   //       aVal = Max( aVal, D );
492   //     }
493   // }
494   // { // polygons
495
496   // }
497
498   if( myPrecision >= 0 )
499   {
500     double prec = pow( 10., (double)myPrecision );
501     aVal = floor( aVal * prec + 0.5 ) / prec;
502   }
503   return aVal;
504 }
505
506 double MaxElementLength2D::GetValue( long theElementId )
507 {
508   TSequenceOfXYZ P;
509   return GetPoints( theElementId, P ) ? GetValue(P) : 0.0;
510 }
511
512 double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const
513 {
514   return Value;
515 }
516
517 SMDSAbs_ElementType MaxElementLength2D::GetType() const
518 {
519   return SMDSAbs_Face;
520 }
521
522 //=======================================================================
523 /*
524   Class       : MaxElementLength3D
525   Description : Functor calculating maximum length of 3D element
526 */
527 //================================================================================
528
529 double MaxElementLength3D::GetValue( long theElementId )
530 {
531   TSequenceOfXYZ P;
532   if( GetPoints( theElementId, P ) ) {
533     double aVal = 0;
534     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
535     SMDSAbs_EntityType      aType = aElem->GetEntityType();
536     int len = P.size();
537     switch ( aType ) {
538     case SMDSEntity_Tetra: { // tetras
539       double L1 = getDistance(P( 1 ),P( 2 ));
540       double L2 = getDistance(P( 2 ),P( 3 ));
541       double L3 = getDistance(P( 3 ),P( 1 ));
542       double L4 = getDistance(P( 1 ),P( 4 ));
543       double L5 = getDistance(P( 2 ),P( 4 ));
544       double L6 = getDistance(P( 3 ),P( 4 ));
545       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
546       break;
547     }
548     case SMDSEntity_Pyramid: { // pyramids
549       double L1 = getDistance(P( 1 ),P( 2 ));
550       double L2 = getDistance(P( 2 ),P( 3 ));
551       double L3 = getDistance(P( 3 ),P( 4 ));
552       double L4 = getDistance(P( 4 ),P( 1 ));
553       double L5 = getDistance(P( 1 ),P( 5 ));
554       double L6 = getDistance(P( 2 ),P( 5 ));
555       double L7 = getDistance(P( 3 ),P( 5 ));
556       double L8 = getDistance(P( 4 ),P( 5 ));
557       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
558       aVal = Max(aVal,Max(L7,L8));
559       break;
560     }
561     case SMDSEntity_Penta: { // pentas
562       double L1 = getDistance(P( 1 ),P( 2 ));
563       double L2 = getDistance(P( 2 ),P( 3 ));
564       double L3 = getDistance(P( 3 ),P( 1 ));
565       double L4 = getDistance(P( 4 ),P( 5 ));
566       double L5 = getDistance(P( 5 ),P( 6 ));
567       double L6 = getDistance(P( 6 ),P( 4 ));
568       double L7 = getDistance(P( 1 ),P( 4 ));
569       double L8 = getDistance(P( 2 ),P( 5 ));
570       double L9 = getDistance(P( 3 ),P( 6 ));
571       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
572       aVal = Max(aVal,Max(Max(L7,L8),L9));
573       break;
574     }
575     case SMDSEntity_Hexa: { // hexas
576       double L1 = getDistance(P( 1 ),P( 2 ));
577       double L2 = getDistance(P( 2 ),P( 3 ));
578       double L3 = getDistance(P( 3 ),P( 4 ));
579       double L4 = getDistance(P( 4 ),P( 1 ));
580       double L5 = getDistance(P( 5 ),P( 6 ));
581       double L6 = getDistance(P( 6 ),P( 7 ));
582       double L7 = getDistance(P( 7 ),P( 8 ));
583       double L8 = getDistance(P( 8 ),P( 5 ));
584       double L9 = getDistance(P( 1 ),P( 5 ));
585       double L10= getDistance(P( 2 ),P( 6 ));
586       double L11= getDistance(P( 3 ),P( 7 ));
587       double L12= getDistance(P( 4 ),P( 8 ));
588       double D1 = getDistance(P( 1 ),P( 7 ));
589       double D2 = getDistance(P( 2 ),P( 8 ));
590       double D3 = getDistance(P( 3 ),P( 5 ));
591       double D4 = getDistance(P( 4 ),P( 6 ));
592       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
593       aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
594       aVal = Max(aVal,Max(L11,L12));
595       aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
596       break;
597     }
598     case SMDSEntity_Hexagonal_Prism: { // hexagonal prism
599       for ( int i1 = 1; i1 < 12; ++i1 )
600         for ( int i2 = i1+1; i1 <= 12; ++i1 )
601           aVal = Max( aVal, getDistance(P( i1 ),P( i2 )));
602       break;
603     }
604     case SMDSEntity_Quad_Tetra: { // quadratic tetras
605       double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
606       double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
607       double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
608       double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
609       double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
610       double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
611       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
612       break;
613     }
614     case SMDSEntity_Quad_Pyramid: { // quadratic pyramids
615       double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
616       double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
617       double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
618       double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
619       double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
620       double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
621       double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
622       double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
623       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
624       aVal = Max(aVal,Max(L7,L8));
625       break;
626     }
627     case SMDSEntity_Quad_Penta:
628     case SMDSEntity_BiQuad_Penta: { // quadratic pentas
629       double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
630       double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
631       double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
632       double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
633       double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
634       double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
635       double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
636       double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
637       double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
638       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
639       aVal = Max(aVal,Max(Max(L7,L8),L9));
640       break;
641     }
642     case SMDSEntity_Quad_Hexa:
643     case SMDSEntity_TriQuad_Hexa: { // quadratic hexas
644       double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
645       double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
646       double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
647       double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
648       double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
649       double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
650       double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
651       double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
652       double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
653       double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
654       double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
655       double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
656       double D1 = getDistance(P( 1 ),P( 7 ));
657       double D2 = getDistance(P( 2 ),P( 8 ));
658       double D3 = getDistance(P( 3 ),P( 5 ));
659       double D4 = getDistance(P( 4 ),P( 6 ));
660       aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
661       aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
662       aVal = Max(aVal,Max(L11,L12));
663       aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
664       break;
665     }
666     case SMDSEntity_Quad_Polyhedra:
667     case SMDSEntity_Polyhedra: { // polys
668       // get the maximum distance between all pairs of nodes
669       for( int i = 1; i <= len; i++ ) {
670         for( int j = 1; j <= len; j++ ) {
671           if( j > i ) { // optimization of the loop
672             double D = getDistance( P(i), P(j) );
673             aVal = Max( aVal, D );
674           }
675         }
676       }
677       break;
678     }
679     case SMDSEntity_Node:
680     case SMDSEntity_0D:
681     case SMDSEntity_Edge:
682     case SMDSEntity_Quad_Edge:
683     case SMDSEntity_Triangle:
684     case SMDSEntity_Quad_Triangle:
685     case SMDSEntity_BiQuad_Triangle:
686     case SMDSEntity_Quadrangle:
687     case SMDSEntity_Quad_Quadrangle:
688     case SMDSEntity_BiQuad_Quadrangle:
689     case SMDSEntity_Polygon:
690     case SMDSEntity_Quad_Polygon:
691     case SMDSEntity_Ball:
692     case SMDSEntity_Last: return 0;
693     } // switch ( aType )
694
695     if( myPrecision >= 0 )
696     {
697       double prec = pow( 10., (double)myPrecision );
698       aVal = floor( aVal * prec + 0.5 ) / prec;
699     }
700     return aVal;
701   }
702   return 0.;
703 }
704
705 double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const
706 {
707   return Value;
708 }
709
710 SMDSAbs_ElementType MaxElementLength3D::GetType() const
711 {
712   return SMDSAbs_Volume;
713 }
714
715 //=======================================================================
716 /*
717   Class       : MinimumAngle
718   Description : Functor for calculation of minimum angle
719 */
720 //================================================================================
721
722 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
723 {
724   if ( P.size() < 3 )
725     return 0.;
726
727   double aMaxCos2;
728
729   aMaxCos2 = getCos2( P( P.size() ), P( 1 ), P( 2 ));
730   aMaxCos2 = Max( aMaxCos2, getCos2( P( P.size()-1 ), P( P.size() ), P( 1 )));
731
732   for ( size_t i = 2; i < P.size(); i++ )
733   {
734     double A0 = getCos2( P( i-1 ), P( i ), P( i+1 ) );
735     aMaxCos2 = Max( aMaxCos2, A0 );
736   }
737   if ( aMaxCos2 < 0 )
738     return 0; // all nodes coincide
739
740   double cos = sqrt( aMaxCos2 );
741   if ( cos >=  1 ) return 0;
742   return acos( cos ) * 180.0 / M_PI;
743 }
744
745 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
746 {
747   //const double aBestAngle = PI / nbNodes;
748   const double aBestAngle = 180.0 - ( 360.0 / double(nbNodes) );
749   return ( fabs( aBestAngle - Value ));
750 }
751
752 SMDSAbs_ElementType MinimumAngle::GetType() const
753 {
754   return SMDSAbs_Face;
755 }
756
757
758 //================================================================================
759 /*
760   Class       : AspectRatio
761   Description : Functor for calculating aspect ratio
762 */
763 //================================================================================
764
765 double AspectRatio::GetValue( long theId )
766 {
767   double aVal = 0;
768   myCurrElement = myMesh->FindElement( theId );
769   TSequenceOfXYZ P;
770   if ( GetPoints( myCurrElement, P ))
771     aVal = Round( GetValue( P ));
772   return aVal;
773 }
774
775 double AspectRatio::GetValue( const TSequenceOfXYZ& P )
776 {
777   // According to "Mesh quality control" by Nadir Bouhamau referring to
778   // Pascal Jean Frey and Paul-Louis George. Maillages, applications aux elements finis.
779   // Hermes Science publications, Paris 1999 ISBN 2-7462-0024-4
780   // PAL10872
781
782   int nbNodes = P.size();
783
784   if ( nbNodes < 3 )
785     return 0;
786
787   // Compute aspect ratio
788
789   if ( nbNodes == 3 ) {
790     // Compute lengths of the sides
791     double aLen1 = getDistance( P( 1 ), P( 2 ));
792     double aLen2 = getDistance( P( 2 ), P( 3 ));
793     double aLen3 = getDistance( P( 3 ), P( 1 ));
794     // Q = alfa * h * p / S, where
795     //
796     // alfa = sqrt( 3 ) / 6
797     // h - length of the longest edge
798     // p - half perimeter
799     // S - triangle surface
800     const double     alfa = sqrt( 3. ) / 6.;
801     double         maxLen = Max( aLen1, Max( aLen2, aLen3 ));
802     double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.;
803     double         anArea = getArea( P( 1 ), P( 2 ), P( 3 ));
804     if ( anArea <= theEps  )
805       return theInf;
806     return alfa * maxLen * half_perimeter / anArea;
807   }
808   else if ( nbNodes == 6 ) { // quadratic triangles
809     // Compute lengths of the sides
810     double aLen1 = getDistance( P( 1 ), P( 3 ));
811     double aLen2 = getDistance( P( 3 ), P( 5 ));
812     double aLen3 = getDistance( P( 5 ), P( 1 ));
813     // algo same as for the linear triangle
814     const double     alfa = sqrt( 3. ) / 6.;
815     double         maxLen = Max( aLen1, Max( aLen2, aLen3 ));
816     double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.;
817     double         anArea = getArea( P( 1 ), P( 3 ), P( 5 ));
818     if ( anArea <= theEps )
819       return theInf;
820     return alfa * maxLen * half_perimeter / anArea;
821   }
822   else if( nbNodes == 4 ) { // quadrangle
823     // Compute lengths of the sides
824     double aLen[4];
825     aLen[0] = getDistance( P(1), P(2) );
826     aLen[1] = getDistance( P(2), P(3) );
827     aLen[2] = getDistance( P(3), P(4) );
828     aLen[3] = getDistance( P(4), P(1) );
829     // Compute lengths of the diagonals
830     double aDia[2];
831     aDia[0] = getDistance( P(1), P(3) );
832     aDia[1] = getDistance( P(2), P(4) );
833     // Compute areas of all triangles which can be built
834     // taking three nodes of the quadrangle
835     double anArea[4];
836     anArea[0] = getArea( P(1), P(2), P(3) );
837     anArea[1] = getArea( P(1), P(2), P(4) );
838     anArea[2] = getArea( P(1), P(3), P(4) );
839     anArea[3] = getArea( P(2), P(3), P(4) );
840     // Q = alpha * L * C1 / C2, where
841     //
842     // alpha = sqrt( 1/32 )
843     // L = max( L1, L2, L3, L4, D1, D2 )
844     // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 )
845     // C2 = min( S1, S2, S3, S4 )
846     // Li - lengths of the edges
847     // Di - lengths of the diagonals
848     // Si - areas of the triangles
849     const double alpha = sqrt( 1 / 32. );
850     double L = Max( aLen[ 0 ],
851                     Max( aLen[ 1 ],
852                          Max( aLen[ 2 ],
853                               Max( aLen[ 3 ],
854                                    Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
855     double C1 = sqrt( aLen[0] * aLen[0] +
856                       aLen[1] * aLen[1] +
857                       aLen[2] * aLen[2] +
858                       aLen[3] * aLen[3] );
859     double C2 = Min( anArea[ 0 ],
860                      Min( anArea[ 1 ],
861                           Min( anArea[ 2 ], anArea[ 3 ] ) ) );
862     if ( C2 <= theEps )
863       return theInf;
864     return alpha * L * C1 / C2;
865   }
866   else if( nbNodes == 8 || nbNodes == 9 ) { // nbNodes==8 - quadratic quadrangle
867     // Compute lengths of the sides
868     double aLen[4];
869     aLen[0] = getDistance( P(1), P(3) );
870     aLen[1] = getDistance( P(3), P(5) );
871     aLen[2] = getDistance( P(5), P(7) );
872     aLen[3] = getDistance( P(7), P(1) );
873     // Compute lengths of the diagonals
874     double aDia[2];
875     aDia[0] = getDistance( P(1), P(5) );
876     aDia[1] = getDistance( P(3), P(7) );
877     // Compute areas of all triangles which can be built
878     // taking three nodes of the quadrangle
879     double anArea[4];
880     anArea[0] = getArea( P(1), P(3), P(5) );
881     anArea[1] = getArea( P(1), P(3), P(7) );
882     anArea[2] = getArea( P(1), P(5), P(7) );
883     anArea[3] = getArea( P(3), P(5), P(7) );
884     // Q = alpha * L * C1 / C2, where
885     //
886     // alpha = sqrt( 1/32 )
887     // L = max( L1, L2, L3, L4, D1, D2 )
888     // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 )
889     // C2 = min( S1, S2, S3, S4 )
890     // Li - lengths of the edges
891     // Di - lengths of the diagonals
892     // Si - areas of the triangles
893     const double alpha = sqrt( 1 / 32. );
894     double L = Max( aLen[ 0 ],
895                     Max( aLen[ 1 ],
896                          Max( aLen[ 2 ],
897                               Max( aLen[ 3 ],
898                                    Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
899     double C1 = sqrt( aLen[0] * aLen[0] +
900                       aLen[1] * aLen[1] +
901                       aLen[2] * aLen[2] +
902                       aLen[3] * aLen[3] );
903     double C2 = Min( anArea[ 0 ],
904                      Min( anArea[ 1 ],
905                           Min( anArea[ 2 ], anArea[ 3 ] ) ) );
906     if ( C2 <= theEps )
907       return theInf;
908     return alpha * L * C1 / C2;
909   }
910   return 0;
911 }
912
913 bool AspectRatio::IsApplicable( const SMDS_MeshElement* element ) const
914 {
915   return ( NumericalFunctor::IsApplicable( element ) && !element->IsPoly() );
916 }
917
918 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
919 {
920   // the aspect ratio is in the range [1.0,infinity]
921   // < 1.0 = very bad, zero area
922   // 1.0 = good
923   // infinity = bad
924   return ( Value < 0.9 ) ? 1000 : Value / 1000.;
925 }
926
927 SMDSAbs_ElementType AspectRatio::GetType() const
928 {
929   return SMDSAbs_Face;
930 }
931
932
933 //================================================================================
934 /*
935   Class       : AspectRatio3D
936   Description : Functor for calculating aspect ratio
937 */
938 //================================================================================
939
940 namespace{
941
942   inline double getHalfPerimeter(double theTria[3]){
943     return (theTria[0] + theTria[1] + theTria[2])/2.0;
944   }
945
946   inline double getArea(double theHalfPerim, double theTria[3]){
947     return sqrt(theHalfPerim*
948                 (theHalfPerim-theTria[0])*
949                 (theHalfPerim-theTria[1])*
950                 (theHalfPerim-theTria[2]));
951   }
952
953   inline double getVolume(double theLen[6]){
954     double a2 = theLen[0]*theLen[0];
955     double b2 = theLen[1]*theLen[1];
956     double c2 = theLen[2]*theLen[2];
957     double d2 = theLen[3]*theLen[3];
958     double e2 = theLen[4]*theLen[4];
959     double f2 = theLen[5]*theLen[5];
960     double P = 4.0*a2*b2*d2;
961     double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
962     double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
963     return sqrt(P-Q+R)/12.0;
964   }
965
966   inline double getVolume2(double theLen[6]){
967     double a2 = theLen[0]*theLen[0];
968     double b2 = theLen[1]*theLen[1];
969     double c2 = theLen[2]*theLen[2];
970     double d2 = theLen[3]*theLen[3];
971     double e2 = theLen[4]*theLen[4];
972     double f2 = theLen[5]*theLen[5];
973
974     double P = a2*e2*(b2+c2+d2+f2-a2-e2);
975     double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
976     double R = c2*d2*(a2+b2+e2+f2-c2-d2);
977     double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
978
979     return sqrt(P+Q+R-S)/12.0;
980   }
981
982   inline double getVolume(const TSequenceOfXYZ& P){
983     gp_Vec aVec1( P( 2 ) - P( 1 ) );
984     gp_Vec aVec2( P( 3 ) - P( 1 ) );
985     gp_Vec aVec3( P( 4 ) - P( 1 ) );
986     gp_Vec anAreaVec( aVec1 ^ aVec2 );
987     return fabs(aVec3 * anAreaVec) / 6.0;
988   }
989
990   inline double getMaxHeight(double theLen[6])
991   {
992     double aHeight = std::max(theLen[0],theLen[1]);
993     aHeight = std::max(aHeight,theLen[2]);
994     aHeight = std::max(aHeight,theLen[3]);
995     aHeight = std::max(aHeight,theLen[4]);
996     aHeight = std::max(aHeight,theLen[5]);
997     return aHeight;
998   }
999
1000   //================================================================================
1001   /*!
1002    * \brief Standard quality of a tetrahedron but not normalized
1003    */
1004   //================================================================================
1005
1006   double tetQualityByHomardMethod( const gp_XYZ & p1,
1007                                    const gp_XYZ & p2,
1008                                    const gp_XYZ & p3,
1009                                    const gp_XYZ & p4 )
1010   {
1011     gp_XYZ edgeVec[6];
1012     edgeVec[0] = ( p1 - p2 );
1013     edgeVec[1] = ( p2 - p3 );
1014     edgeVec[2] = ( p3 - p1 );
1015     edgeVec[3] = ( p4 - p1 );
1016     edgeVec[4] = ( p4 - p2 );
1017     edgeVec[5] = ( p4 - p3 );
1018
1019     double maxEdgeLen2            = edgeVec[0].SquareModulus();
1020     maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[1].SquareModulus() );
1021     maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[2].SquareModulus() );
1022     maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[3].SquareModulus() );
1023     maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[4].SquareModulus() );
1024     maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[5].SquareModulus() );
1025     double maxEdgeLen = Sqrt( maxEdgeLen2 );
1026
1027     gp_XYZ cross01 = edgeVec[0] ^ edgeVec[1];
1028     double sumArea = ( cross01                 ).Modulus(); // actually double area
1029     sumArea       += ( edgeVec[0] ^ edgeVec[3] ).Modulus();
1030     sumArea       += ( edgeVec[1] ^ edgeVec[4] ).Modulus();
1031     sumArea       += ( edgeVec[2] ^ edgeVec[5] ).Modulus();
1032
1033     double sixVolume = Abs( cross01 * edgeVec[4] ); // 6 * volume
1034     double quality   = maxEdgeLen * sumArea / sixVolume; // not normalized!!!
1035     return quality;
1036   }
1037
1038   //================================================================================
1039   /*!
1040    * \brief HOMARD method of hexahedron quality
1041    * 1. Decompose the hexa into 24 tetra: each face is splitted into 4 triangles by
1042    *    adding the diagonals and every triangle is connected to the center of the hexa.
1043    * 2. Compute the quality of every tetra with the same formula as for the standard quality,
1044    *    except that the factor for the normalization is not the same because the final goal
1045    *    is to have a quality equal to 1 for a perfect cube. So the formula is:
1046    *    qual = max(lengthes of 6 edges) * (sum of surfaces of 4 faces) / (7.6569*6*volume)
1047    * 3. The quality of the hexa is the highest value of the qualities of the 24 tetra
1048    */
1049   //================================================================================
1050
1051   double hexQualityByHomardMethod( const TSequenceOfXYZ& P )
1052   {
1053     gp_XYZ quadCenter[6];
1054     quadCenter[0] = ( P(1) + P(2) + P(3) + P(4) ) / 4.;
1055     quadCenter[1] = ( P(5) + P(6) + P(7) + P(8) ) / 4.;
1056     quadCenter[2] = ( P(1) + P(2) + P(6) + P(5) ) / 4.;
1057     quadCenter[3] = ( P(2) + P(3) + P(7) + P(6) ) / 4.;
1058     quadCenter[4] = ( P(3) + P(4) + P(8) + P(7) ) / 4.;
1059     quadCenter[5] = ( P(1) + P(4) + P(8) + P(5) ) / 4.;
1060
1061     gp_XYZ hexCenter = ( P(1) + P(2) + P(3) + P(4) + P(5) + P(6) + P(7) + P(8) ) / 8.;
1062
1063     // quad 1 ( 1 2 3 4 )
1064     double quality =        tetQualityByHomardMethod( P(1), P(2), quadCenter[0], hexCenter );
1065     quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[0], hexCenter ));
1066     quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[0], hexCenter ));
1067     quality = Max( quality, tetQualityByHomardMethod( P(4), P(1), quadCenter[0], hexCenter ));
1068     // quad 2 ( 5 6 7 8 )
1069     quality = Max( quality, tetQualityByHomardMethod( P(5), P(6), quadCenter[1], hexCenter ));
1070     quality = Max( quality, tetQualityByHomardMethod( P(6), P(7), quadCenter[1], hexCenter ));
1071     quality = Max( quality, tetQualityByHomardMethod( P(7), P(8), quadCenter[1], hexCenter ));
1072     quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[1], hexCenter ));
1073     // quad 3 ( 1 2 6 5 )
1074     quality = Max( quality, tetQualityByHomardMethod( P(1), P(2), quadCenter[2], hexCenter ));
1075     quality = Max( quality, tetQualityByHomardMethod( P(2), P(6), quadCenter[2], hexCenter ));
1076     quality = Max( quality, tetQualityByHomardMethod( P(6), P(5), quadCenter[2], hexCenter ));
1077     quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[2], hexCenter ));
1078     // quad 4 ( 2 3 7 6 )
1079     quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[3], hexCenter ));
1080     quality = Max( quality, tetQualityByHomardMethod( P(3), P(7), quadCenter[3], hexCenter ));
1081     quality = Max( quality, tetQualityByHomardMethod( P(7), P(6), quadCenter[3], hexCenter ));
1082     quality = Max( quality, tetQualityByHomardMethod( P(6), P(2), quadCenter[3], hexCenter ));
1083     // quad 5 ( 3 4 8 7 )
1084     quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[4], hexCenter ));
1085     quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[4], hexCenter ));
1086     quality = Max( quality, tetQualityByHomardMethod( P(8), P(7), quadCenter[4], hexCenter ));
1087     quality = Max( quality, tetQualityByHomardMethod( P(7), P(3), quadCenter[4], hexCenter ));
1088     // quad 6 ( 1 4 8 5 )
1089     quality = Max( quality, tetQualityByHomardMethod( P(1), P(4), quadCenter[5], hexCenter ));
1090     quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[5], hexCenter ));
1091     quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[5], hexCenter ));
1092     quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[5], hexCenter ));
1093
1094     return quality / 7.65685424949;
1095   }
1096 }
1097
1098 double AspectRatio3D::GetValue( long theId )
1099 {
1100   double aVal = 0;
1101   myCurrElement = myMesh->FindElement( theId );
1102   if ( myCurrElement && myCurrElement->GetVtkType() == VTK_TETRA )
1103   {
1104     // Action from CoTech | ACTION 31.3:
1105     // EURIWARE BO: Homogenize the formulas used to calculate the Controls in SMESH to fit with
1106     // those of ParaView. The library used by ParaView for those calculations can be reused in SMESH.
1107     vtkUnstructuredGrid* grid = const_cast<SMDS_Mesh*>( myMesh )->GetGrid();
1108     if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->GetVtkID() ))
1109       aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell ));
1110   }
1111   else
1112   {
1113     TSequenceOfXYZ P;
1114     if ( GetPoints( myCurrElement, P ))
1115       aVal = Round( GetValue( P ));
1116   }
1117   return aVal;
1118 }
1119
1120 bool AspectRatio3D::IsApplicable( const SMDS_MeshElement* element ) const
1121 {
1122   return ( NumericalFunctor::IsApplicable( element ) && !element->IsPoly() );
1123 }
1124
1125 double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
1126 {
1127   double aQuality = 0.0;
1128   if(myCurrElement->IsPoly()) return aQuality;
1129
1130   int nbNodes = P.size();
1131
1132   if( myCurrElement->IsQuadratic() ) {
1133     if     (nbNodes==10) nbNodes=4; // quadratic tetrahedron
1134     else if(nbNodes==13) nbNodes=5; // quadratic pyramid
1135     else if(nbNodes==15) nbNodes=6; // quadratic pentahedron
1136     else if(nbNodes==20) nbNodes=8; // quadratic hexahedron
1137     else if(nbNodes==27) nbNodes=8; // tri-quadratic hexahedron
1138     else return aQuality;
1139   }
1140
1141   switch(nbNodes) {
1142   case 4:{
1143     double aLen[6] = {
1144       getDistance(P( 1 ),P( 2 )), // a
1145       getDistance(P( 2 ),P( 3 )), // b
1146       getDistance(P( 3 ),P( 1 )), // c
1147       getDistance(P( 2 ),P( 4 )), // d
1148       getDistance(P( 3 ),P( 4 )), // e
1149       getDistance(P( 1 ),P( 4 ))  // f
1150     };
1151     double aTria[4][3] = {
1152       {aLen[0],aLen[1],aLen[2]}, // abc
1153       {aLen[0],aLen[3],aLen[5]}, // adf
1154       {aLen[1],aLen[3],aLen[4]}, // bde
1155       {aLen[2],aLen[4],aLen[5]}  // cef
1156     };
1157     double aSumArea = 0.0;
1158     double aHalfPerimeter = getHalfPerimeter(aTria[0]);
1159     double anArea = getArea(aHalfPerimeter,aTria[0]);
1160     aSumArea += anArea;
1161     aHalfPerimeter = getHalfPerimeter(aTria[1]);
1162     anArea = getArea(aHalfPerimeter,aTria[1]);
1163     aSumArea += anArea;
1164     aHalfPerimeter = getHalfPerimeter(aTria[2]);
1165     anArea = getArea(aHalfPerimeter,aTria[2]);
1166     aSumArea += anArea;
1167     aHalfPerimeter = getHalfPerimeter(aTria[3]);
1168     anArea = getArea(aHalfPerimeter,aTria[3]);
1169     aSumArea += anArea;
1170     double aVolume = getVolume(P);
1171     //double aVolume = getVolume(aLen);
1172     double aHeight = getMaxHeight(aLen);
1173     static double aCoeff = sqrt(2.0)/12.0;
1174     if ( aVolume > DBL_MIN )
1175       aQuality = aCoeff*aHeight*aSumArea/aVolume;
1176     break;
1177   }
1178   case 5:{
1179     {
1180       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
1181       aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
1182     }
1183     {
1184       gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
1185       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1186     }
1187     {
1188       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
1189       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1190     }
1191     {
1192       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
1193       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1194     }
1195     break;
1196   }
1197   case 6:{
1198     {
1199       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
1200       aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
1201     }
1202     {
1203       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
1204       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1205     }
1206     {
1207       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
1208       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1209     }
1210     {
1211       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
1212       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1213     }
1214     {
1215       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
1216       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1217     }
1218     {
1219       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
1220       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1221     }
1222     break;
1223   }
1224   case 8:{
1225
1226     return hexQualityByHomardMethod( P ); // bos #23982
1227
1228
1229     {
1230       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
1231       aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
1232     }
1233     {
1234       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
1235       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1236     }
1237     {
1238       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
1239       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1240     }
1241     {
1242       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
1243       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1244     }
1245     {
1246       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
1247       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1248     }
1249     {
1250       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
1251       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1252     }
1253     {
1254       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
1255       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1256     }
1257     {
1258       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
1259       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1260     }
1261     {
1262       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
1263       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1264     }
1265     {
1266       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
1267       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1268     }
1269     {
1270       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
1271       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1272     }
1273     {
1274       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
1275       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1276     }
1277     {
1278       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
1279       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1280     }
1281     {
1282       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
1283       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1284     }
1285     {
1286       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
1287       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1288     }
1289     {
1290       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
1291       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1292     }
1293     {
1294       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
1295       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1296     }
1297     {
1298       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
1299       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1300     }
1301     {
1302       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
1303       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1304     }
1305     {
1306       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
1307       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1308     }
1309     {
1310       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
1311       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1312     }
1313     {
1314       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1315       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1316     }
1317     {
1318       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
1319       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1320     }
1321     {
1322       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
1323       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1324     }
1325     {
1326       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1327       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1328     }
1329     {
1330       gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
1331       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1332     }
1333     {
1334       gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
1335       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1336     }
1337     {
1338       gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
1339       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1340     }
1341     {
1342       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
1343       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1344     }
1345     {
1346       gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
1347       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1348     }
1349     {
1350       gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
1351       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1352     }
1353     {
1354       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
1355       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1356     }
1357     {
1358       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
1359       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1360     }
1361     break;
1362   }
1363   case 12:
1364     {
1365       gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )};
1366       aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8]));
1367     }
1368     {
1369       gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )};
1370       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1371     }
1372     {
1373       gp_XYZ aXYZ[8] = {P( 3 ),P( 4 ),P( 6 ),P( 1 ),P( 9 ),P( 10 ),P( 12 ),P( 7 )};
1374       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1375     }
1376     break;
1377   } // switch(nbNodes)
1378
1379   if ( nbNodes > 4 ) {
1380     // evaluate aspect ratio of quadrangle faces
1381     AspectRatio aspect2D;
1382     SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes );
1383     int nbFaces = SMDS_VolumeTool::NbFaces( type );
1384     TSequenceOfXYZ points(4);
1385     for ( int i = 0; i < nbFaces; ++i ) { // loop on faces of a volume
1386       if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 )
1387         continue;
1388       const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true );
1389       for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadrangle face
1390         points( p + 1 ) = P( pInd[ p ] + 1 );
1391       aQuality = std::max( aQuality, aspect2D.GetValue( points ));
1392     }
1393   }
1394   return aQuality;
1395 }
1396
1397 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
1398 {
1399   // the aspect ratio is in the range [1.0,infinity]
1400   // 1.0 = good
1401   // infinity = bad
1402   return Value / 1000.;
1403 }
1404
1405 SMDSAbs_ElementType AspectRatio3D::GetType() const
1406 {
1407   return SMDSAbs_Volume;
1408 }
1409
1410
1411 //================================================================================
1412 /*
1413   Class       : Warping
1414   Description : Functor for calculating warping
1415 */
1416 //================================================================================
1417
1418 bool Warping::IsApplicable( const SMDS_MeshElement* element ) const
1419 {
1420   return NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4;
1421 }
1422
1423 double Warping::GetValue( const TSequenceOfXYZ& P )
1424 {
1425   if ( P.size() != 4 )
1426     return 0;
1427
1428   gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4.;
1429
1430   double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
1431   double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
1432   double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
1433   double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
1434
1435   double val = Max( Max( A1, A2 ), Max( A3, A4 ) );
1436
1437   const double eps = 0.1; // val is in degrees
1438
1439   return val < eps ? 0. : val;
1440 }
1441
1442 double Warping::ComputeA( const gp_XYZ& thePnt1,
1443                           const gp_XYZ& thePnt2,
1444                           const gp_XYZ& thePnt3,
1445                           const gp_XYZ& theG ) const
1446 {
1447   double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
1448   double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
1449   double L = Min( aLen1, aLen2 ) * 0.5;
1450   if ( L < theEps )
1451     return theInf;
1452
1453   gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG;
1454   gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG;
1455   gp_XYZ N  = GI.Crossed( GJ );
1456
1457   if ( N.Modulus() < gp::Resolution() )
1458     return M_PI / 2;
1459
1460   N.Normalize();
1461
1462   double H = ( thePnt2 - theG ).Dot( N );
1463   return asin( fabs( H / L ) ) * 180. / M_PI;
1464 }
1465
1466 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
1467 {
1468   // the warp is in the range [0.0,PI/2]
1469   // 0.0 = good (no warp)
1470   // PI/2 = bad  (face pliee)
1471   return Value;
1472 }
1473
1474 SMDSAbs_ElementType Warping::GetType() const
1475 {
1476   return SMDSAbs_Face;
1477 }
1478
1479
1480 //================================================================================
1481 /*
1482   Class       : Taper
1483   Description : Functor for calculating taper
1484 */
1485 //================================================================================
1486
1487 bool Taper::IsApplicable( const SMDS_MeshElement* element ) const
1488 {
1489   return ( NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4 );
1490 }
1491
1492 double Taper::GetValue( const TSequenceOfXYZ& P )
1493 {
1494   if ( P.size() != 4 )
1495     return 0.;
1496
1497   // Compute taper
1498   double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) );
1499   double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) );
1500   double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) );
1501   double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) );
1502
1503   double JA = 0.25 * ( J1 + J2 + J3 + J4 );
1504   if ( JA <= theEps )
1505     return theInf;
1506
1507   double T1 = fabs( ( J1 - JA ) / JA );
1508   double T2 = fabs( ( J2 - JA ) / JA );
1509   double T3 = fabs( ( J3 - JA ) / JA );
1510   double T4 = fabs( ( J4 - JA ) / JA );
1511
1512   double val = Max( Max( T1, T2 ), Max( T3, T4 ) );
1513
1514   const double eps = 0.01;
1515
1516   return val < eps ? 0. : val;
1517 }
1518
1519 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
1520 {
1521   // the taper is in the range [0.0,1.0]
1522   // 0.0 = good (no taper)
1523   // 1.0 = bad  (les cotes opposes sont allignes)
1524   return Value;
1525 }
1526
1527 SMDSAbs_ElementType Taper::GetType() const
1528 {
1529   return SMDSAbs_Face;
1530 }
1531
1532 //================================================================================
1533 /*
1534   Class       : Skew
1535   Description : Functor for calculating skew in degrees
1536 */
1537 //================================================================================
1538
1539 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
1540 {
1541   gp_XYZ p12 = ( p2 + p1 ) / 2.;
1542   gp_XYZ p23 = ( p3 + p2 ) / 2.;
1543   gp_XYZ p31 = ( p3 + p1 ) / 2.;
1544
1545   gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
1546
1547   return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 );
1548 }
1549
1550 bool Skew::IsApplicable( const SMDS_MeshElement* element ) const
1551 {
1552   return ( NumericalFunctor::IsApplicable( element ) && element->NbNodes() <= 4 );
1553 }
1554
1555 double Skew::GetValue( const TSequenceOfXYZ& P )
1556 {
1557   if ( P.size() != 3 && P.size() != 4 )
1558     return 0.;
1559
1560   // Compute skew
1561   const double PI2 = M_PI / 2.;
1562   if ( P.size() == 3 )
1563   {
1564     double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
1565     double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
1566     double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
1567
1568     return Max( A0, Max( A1, A2 ) ) * 180. / M_PI;
1569   }
1570   else
1571   {
1572     gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2.;
1573     gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2.;
1574     gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2.;
1575     gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2.;
1576
1577     gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
1578     double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
1579       ? 0. : fabs( PI2 - v1.Angle( v2 ) );
1580
1581     double val = A * 180. / M_PI;
1582
1583     const double eps = 0.1; // val is in degrees
1584
1585     return val < eps ? 0. : val;
1586   }
1587 }
1588
1589 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
1590 {
1591   // the skew is in the range [0.0,PI/2].
1592   // 0.0 = good
1593   // PI/2 = bad
1594   return Value;
1595 }
1596
1597 SMDSAbs_ElementType Skew::GetType() const
1598 {
1599   return SMDSAbs_Face;
1600 }
1601
1602
1603 //================================================================================
1604 /*
1605   Class       : Area
1606   Description : Functor for calculating area
1607 */
1608 //================================================================================
1609
1610 double Area::GetValue( const TSequenceOfXYZ& P )
1611 {
1612   double val = 0.0;
1613   if ( P.size() > 2 )
1614   {
1615     gp_Vec aVec1( P(2) - P(1) );
1616     gp_Vec aVec2( P(3) - P(1) );
1617     gp_Vec SumVec = aVec1 ^ aVec2;
1618
1619     for (size_t i=4; i<=P.size(); i++)
1620     {
1621       gp_Vec aVec1( P(i-1) - P(1) );
1622       gp_Vec aVec2( P(i  ) - P(1) );
1623       gp_Vec tmp = aVec1 ^ aVec2;
1624       SumVec.Add(tmp);
1625     }
1626     val = SumVec.Magnitude() * 0.5;
1627   }
1628   return val;
1629 }
1630
1631 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
1632 {
1633   // meaningless as it is not a quality control functor
1634   return Value;
1635 }
1636
1637 SMDSAbs_ElementType Area::GetType() const
1638 {
1639   return SMDSAbs_Face;
1640 }
1641
1642 //================================================================================
1643 /*
1644   Class       : Length
1645   Description : Functor for calculating length of edge
1646 */
1647 //================================================================================
1648
1649 double Length::GetValue( const TSequenceOfXYZ& P )
1650 {
1651   switch ( P.size() ) {
1652   case 2:  return getDistance( P( 1 ), P( 2 ) );
1653   case 3:  return getDistance( P( 1 ), P( 2 ) ) + getDistance( P( 2 ), P( 3 ) );
1654   default: return 0.;
1655   }
1656 }
1657
1658 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
1659 {
1660   // meaningless as it is not quality control functor
1661   return Value;
1662 }
1663
1664 SMDSAbs_ElementType Length::GetType() const
1665 {
1666   return SMDSAbs_Edge;
1667 }
1668
1669 //================================================================================
1670 /*
1671   Class       : Length3D
1672   Description : Functor for calculating minimal length of element edge
1673 */
1674 //================================================================================
1675
1676 Length3D::Length3D():
1677   Length2D ( SMDSAbs_Volume )
1678 {
1679 }
1680
1681 //================================================================================
1682 /*
1683   Class       : Length2D
1684   Description : Functor for calculating minimal length of element edge
1685 */
1686 //================================================================================
1687
1688 Length2D::Length2D( SMDSAbs_ElementType type ):
1689   myType ( type )
1690 {
1691 }
1692
1693 bool Length2D::IsApplicable( const SMDS_MeshElement* element ) const
1694 {
1695   return ( NumericalFunctor::IsApplicable( element ) &&
1696            element->GetEntityType() != SMDSEntity_Polyhedra );
1697 }
1698
1699 double Length2D::GetValue( const TSequenceOfXYZ& P )
1700 {
1701   double aVal = 0;
1702   int len = P.size();
1703   SMDSAbs_EntityType aType = P.getElementEntity();
1704
1705   switch (aType) {
1706   case SMDSEntity_Edge:
1707     if (len == 2)
1708       aVal = getDistance( P( 1 ), P( 2 ) );
1709     break;
1710   case SMDSEntity_Quad_Edge:
1711     if (len == 3) // quadratic edge
1712       aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
1713     break;
1714   case SMDSEntity_Triangle:
1715     if (len == 3){ // triangles
1716       double L1 = getDistance(P( 1 ),P( 2 ));
1717       double L2 = getDistance(P( 2 ),P( 3 ));
1718       double L3 = getDistance(P( 3 ),P( 1 ));
1719       aVal = Min(L1,Min(L2,L3));
1720     }
1721     break;
1722   case SMDSEntity_Quadrangle:
1723     if (len == 4){ // quadrangles
1724       double L1 = getDistance(P( 1 ),P( 2 ));
1725       double L2 = getDistance(P( 2 ),P( 3 ));
1726       double L3 = getDistance(P( 3 ),P( 4 ));
1727       double L4 = getDistance(P( 4 ),P( 1 ));
1728       aVal = Min(Min(L1,L2),Min(L3,L4));
1729     }
1730     break;
1731   case SMDSEntity_Quad_Triangle:
1732   case SMDSEntity_BiQuad_Triangle:
1733     if (len >= 6){ // quadratic triangles
1734       double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1735       double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1736       double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
1737       aVal = Min(L1,Min(L2,L3));
1738     }
1739     break;
1740   case SMDSEntity_Quad_Quadrangle:
1741   case SMDSEntity_BiQuad_Quadrangle:
1742     if (len >= 8){ // quadratic quadrangles
1743       double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1744       double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1745       double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
1746       double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
1747       aVal = Min(Min(L1,L2),Min(L3,L4));
1748     }
1749     break;
1750   case SMDSEntity_Tetra:
1751     if (len == 4){ // tetrahedra
1752       double L1 = getDistance(P( 1 ),P( 2 ));
1753       double L2 = getDistance(P( 2 ),P( 3 ));
1754       double L3 = getDistance(P( 3 ),P( 1 ));
1755       double L4 = getDistance(P( 1 ),P( 4 ));
1756       double L5 = getDistance(P( 2 ),P( 4 ));
1757       double L6 = getDistance(P( 3 ),P( 4 ));
1758       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1759     }
1760     break;
1761   case SMDSEntity_Pyramid:
1762     if (len == 5){ // pyramid
1763       double L1 = getDistance(P( 1 ),P( 2 ));
1764       double L2 = getDistance(P( 2 ),P( 3 ));
1765       double L3 = getDistance(P( 3 ),P( 4 ));
1766       double L4 = getDistance(P( 4 ),P( 1 ));
1767       double L5 = getDistance(P( 1 ),P( 5 ));
1768       double L6 = getDistance(P( 2 ),P( 5 ));
1769       double L7 = getDistance(P( 3 ),P( 5 ));
1770       double L8 = getDistance(P( 4 ),P( 5 ));
1771
1772       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1773       aVal = Min(aVal,Min(L7,L8));
1774     }
1775     break;
1776   case SMDSEntity_Penta:
1777     if (len == 6) { // pentahedron
1778       double L1 = getDistance(P( 1 ),P( 2 ));
1779       double L2 = getDistance(P( 2 ),P( 3 ));
1780       double L3 = getDistance(P( 3 ),P( 1 ));
1781       double L4 = getDistance(P( 4 ),P( 5 ));
1782       double L5 = getDistance(P( 5 ),P( 6 ));
1783       double L6 = getDistance(P( 6 ),P( 4 ));
1784       double L7 = getDistance(P( 1 ),P( 4 ));
1785       double L8 = getDistance(P( 2 ),P( 5 ));
1786       double L9 = getDistance(P( 3 ),P( 6 ));
1787
1788       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1789       aVal = Min(aVal,Min(Min(L7,L8),L9));
1790     }
1791     break;
1792   case SMDSEntity_Hexa:
1793     if (len == 8){ // hexahedron
1794       double L1 = getDistance(P( 1 ),P( 2 ));
1795       double L2 = getDistance(P( 2 ),P( 3 ));
1796       double L3 = getDistance(P( 3 ),P( 4 ));
1797       double L4 = getDistance(P( 4 ),P( 1 ));
1798       double L5 = getDistance(P( 5 ),P( 6 ));
1799       double L6 = getDistance(P( 6 ),P( 7 ));
1800       double L7 = getDistance(P( 7 ),P( 8 ));
1801       double L8 = getDistance(P( 8 ),P( 5 ));
1802       double L9 = getDistance(P( 1 ),P( 5 ));
1803       double L10= getDistance(P( 2 ),P( 6 ));
1804       double L11= getDistance(P( 3 ),P( 7 ));
1805       double L12= getDistance(P( 4 ),P( 8 ));
1806
1807       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1808       aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
1809       aVal = Min(aVal,Min(L11,L12));
1810     }
1811     break;
1812   case SMDSEntity_Quad_Tetra:
1813     if (len == 10){ // quadratic tetrahedron
1814       double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
1815       double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
1816       double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
1817       double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1818       double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
1819       double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
1820       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1821     }
1822     break;
1823   case SMDSEntity_Quad_Pyramid:
1824     if (len == 13){ // quadratic pyramid
1825       double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
1826       double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
1827       double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1828       double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1829       double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1830       double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
1831       double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
1832       double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
1833       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1834       aVal = Min(aVal,Min(L7,L8));
1835     }
1836     break;
1837   case SMDSEntity_Quad_Penta:
1838   case SMDSEntity_BiQuad_Penta:
1839     if (len >= 15){ // quadratic pentahedron
1840       double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
1841       double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
1842       double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1843       double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1844       double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
1845       double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
1846       double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
1847       double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
1848       double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
1849       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1850       aVal = Min(aVal,Min(Min(L7,L8),L9));
1851     }
1852     break;
1853   case SMDSEntity_Quad_Hexa:
1854   case SMDSEntity_TriQuad_Hexa:
1855     if (len >= 20) { // quadratic hexahedron
1856       double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
1857       double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
1858       double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
1859       double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
1860       double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
1861       double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
1862       double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
1863       double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
1864       double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
1865       double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
1866       double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
1867       double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
1868       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1869       aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
1870       aVal = Min(aVal,Min(L11,L12));
1871     }
1872     break;
1873   case SMDSEntity_Polygon:
1874     if ( len > 1 ) {
1875       aVal = getDistance( P(1), P( P.size() ));
1876       for ( size_t i = 1; i < P.size(); ++i )
1877         aVal = Min( aVal, getDistance( P( i ), P( i+1 )));
1878     }
1879     break;
1880   case SMDSEntity_Quad_Polygon:
1881     if ( len > 2 ) {
1882       aVal = getDistance( P(1), P( P.size() )) + getDistance( P(P.size()), P( P.size()-1 ));
1883       for ( size_t i = 1; i < P.size()-1; i += 2 )
1884         aVal = Min( aVal, getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 )));
1885     }
1886     break;
1887   case SMDSEntity_Hexagonal_Prism:
1888     if (len == 12) { // hexagonal prism
1889       double L1 = getDistance(P( 1 ),P( 2 ));
1890       double L2 = getDistance(P( 2 ),P( 3 ));
1891       double L3 = getDistance(P( 3 ),P( 4 ));
1892       double L4 = getDistance(P( 4 ),P( 5 ));
1893       double L5 = getDistance(P( 5 ),P( 6 ));
1894       double L6 = getDistance(P( 6 ),P( 1 ));
1895
1896       double L7 = getDistance(P( 7 ), P( 8 ));
1897       double L8 = getDistance(P( 8 ), P( 9 ));
1898       double L9 = getDistance(P( 9 ), P( 10 ));
1899       double L10= getDistance(P( 10 ),P( 11 ));
1900       double L11= getDistance(P( 11 ),P( 12 ));
1901       double L12= getDistance(P( 12 ),P( 7 ));
1902
1903       double L13 = getDistance(P( 1 ),P( 7 ));
1904       double L14 = getDistance(P( 2 ),P( 8 ));
1905       double L15 = getDistance(P( 3 ),P( 9 ));
1906       double L16 = getDistance(P( 4 ),P( 10 ));
1907       double L17 = getDistance(P( 5 ),P( 11 ));
1908       double L18 = getDistance(P( 6 ),P( 12 ));
1909       aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1910       aVal = Min(aVal, Min(Min(Min(L7,L8),Min(L9,L10)),Min(L11,L12)));
1911       aVal = Min(aVal, Min(Min(Min(L13,L14),Min(L15,L16)),Min(L17,L18)));
1912     }
1913     break;
1914   case SMDSEntity_Polyhedra:
1915   {
1916   }
1917   break;
1918   default:
1919     return 0;
1920   }
1921
1922   if (aVal < 0 ) {
1923     return 0.;
1924   }
1925
1926   if ( myPrecision >= 0 )
1927   {
1928     double prec = pow( 10., (double)( myPrecision ) );
1929     aVal = floor( aVal * prec + 0.5 ) / prec;
1930   }
1931
1932   return aVal;
1933 }
1934
1935 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1936 {
1937   // meaningless as it is not a quality control functor
1938   return Value;
1939 }
1940
1941 SMDSAbs_ElementType Length2D::GetType() const
1942 {
1943   return myType;
1944 }
1945
1946 Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
1947   myLength(theLength)
1948 {
1949   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
1950   if(thePntId1 > thePntId2){
1951     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
1952   }
1953 }
1954
1955 bool Length2D::Value::operator<(const Length2D::Value& x) const
1956 {
1957   if(myPntId[0] < x.myPntId[0]) return true;
1958   if(myPntId[0] == x.myPntId[0])
1959     if(myPntId[1] < x.myPntId[1]) return true;
1960   return false;
1961 }
1962
1963 void Length2D::GetValues(TValues& theValues)
1964 {
1965   if ( myType == SMDSAbs_Face )
1966   {
1967     for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
1968     {
1969       const SMDS_MeshFace* anElem = anIter->next();
1970       if ( anElem->IsQuadratic() )
1971       {
1972         // use special nodes iterator
1973         SMDS_NodeIteratorPtr anIter = anElem->interlacedNodesIterator();
1974         long aNodeId[4] = { 0,0,0,0 };
1975         gp_Pnt P[4];
1976
1977         double aLength = 0;
1978         if ( anIter->more() )
1979         {
1980           const SMDS_MeshNode* aNode = anIter->next();
1981           P[0] = P[1] = SMESH_NodeXYZ( aNode );
1982           aNodeId[0] = aNodeId[1] = aNode->GetID();
1983           aLength = 0;
1984         }
1985         for ( ; anIter->more(); )
1986         {
1987           const SMDS_MeshNode* N1 = anIter->next();
1988           P[2] = SMESH_NodeXYZ( N1 );
1989           aNodeId[2] = N1->GetID();
1990           aLength = P[1].Distance(P[2]);
1991           if(!anIter->more()) break;
1992           const SMDS_MeshNode* N2 = anIter->next();
1993           P[3] = SMESH_NodeXYZ( N2 );
1994           aNodeId[3] = N2->GetID();
1995           aLength += P[2].Distance(P[3]);
1996           Value aValue1(aLength,aNodeId[1],aNodeId[2]);
1997           Value aValue2(aLength,aNodeId[2],aNodeId[3]);
1998           P[1] = P[3];
1999           aNodeId[1] = aNodeId[3];
2000           theValues.insert(aValue1);
2001           theValues.insert(aValue2);
2002         }
2003         aLength += P[2].Distance(P[0]);
2004         Value aValue1(aLength,aNodeId[1],aNodeId[2]);
2005         Value aValue2(aLength,aNodeId[2],aNodeId[0]);
2006         theValues.insert(aValue1);
2007         theValues.insert(aValue2);
2008       }
2009       else {
2010         SMDS_NodeIteratorPtr aNodesIter = anElem->nodeIterator();
2011         long aNodeId[2] = {0,0};
2012         gp_Pnt P[3];
2013
2014         double aLength;
2015         const SMDS_MeshElement* aNode;
2016         if ( aNodesIter->more())
2017         {
2018           aNode = aNodesIter->next();
2019           P[0] = P[1] = SMESH_NodeXYZ( aNode );
2020           aNodeId[0] = aNodeId[1] = aNode->GetID();
2021           aLength = 0;
2022         }
2023         for( ; aNodesIter->more(); )
2024         {
2025           aNode = aNodesIter->next();
2026           long anId = aNode->GetID();
2027
2028           P[2] = SMESH_NodeXYZ( aNode );
2029
2030           aLength = P[1].Distance(P[2]);
2031
2032           Value aValue(aLength,aNodeId[1],anId);
2033           aNodeId[1] = anId;
2034           P[1] = P[2];
2035           theValues.insert(aValue);
2036         }
2037
2038         aLength = P[0].Distance(P[1]);
2039
2040         Value aValue(aLength,aNodeId[0],aNodeId[1]);
2041         theValues.insert(aValue);
2042       }
2043     }
2044   }
2045   else
2046   {
2047     // not implemented
2048   }
2049 }
2050
2051 //================================================================================
2052 /*
2053   Class       : Deflection2D
2054   Description : computes distance between a face center and an underlying surface
2055 */
2056 //================================================================================
2057
2058 double Deflection2D::GetValue( const TSequenceOfXYZ& P )
2059 {
2060   if ( myMesh && P.getElement() )
2061   {
2062     // get underlying surface
2063     if ( myShapeIndex != P.getElement()->getshapeId() )
2064     {
2065       mySurface.Nullify();
2066       myShapeIndex = P.getElement()->getshapeId();
2067       const TopoDS_Shape& S =
2068         static_cast< const SMESHDS_Mesh* >( myMesh )->IndexToShape( myShapeIndex );
2069       if ( !S.IsNull() && S.ShapeType() == TopAbs_FACE )
2070       {
2071         mySurface = new ShapeAnalysis_Surface( BRep_Tool::Surface( TopoDS::Face( S )));
2072
2073         GeomLib_IsPlanarSurface isPlaneCheck( mySurface->Surface() );
2074         if ( isPlaneCheck.IsPlanar() )
2075           myPlane.reset( new gp_Pln( isPlaneCheck.Plan() ));
2076         else
2077           myPlane.reset();
2078       }
2079     }
2080     // project gravity center to the surface
2081     if ( !mySurface.IsNull() )
2082     {
2083       gp_XYZ gc(0,0,0);
2084       gp_XY  uv(0,0);
2085       int nbUV = 0;
2086       for ( size_t i = 0; i < P.size(); ++i )
2087       {
2088         gc += P(i+1);
2089
2090         if ( SMDS_FacePositionPtr fPos = P.getElement()->GetNode( i )->GetPosition() )
2091         {
2092           uv.ChangeCoord(1) += fPos->GetUParameter();
2093           uv.ChangeCoord(2) += fPos->GetVParameter();
2094           ++nbUV;
2095         }
2096       }
2097       gc /= P.size();
2098       if ( nbUV ) uv /= nbUV;
2099
2100       double maxLen = MaxElementLength2D().GetValue( P );
2101       double    tol = 1e-3 * maxLen;
2102       double dist;
2103       if ( myPlane )
2104       {
2105         dist = myPlane->Distance( gc );
2106         if ( dist < tol )
2107           dist = 0;
2108       }
2109       else
2110       {
2111         if ( uv.X() != 0 && uv.Y() != 0 ) // faster way
2112           mySurface->NextValueOfUV( uv, gc, tol, 0.5 * maxLen );
2113         else
2114           mySurface->ValueOfUV( gc, tol );
2115         dist = mySurface->Gap();
2116       }
2117       return Round( dist );
2118     }
2119   }
2120   return 0;
2121 }
2122
2123 void Deflection2D::SetMesh( const SMDS_Mesh* theMesh )
2124 {
2125   NumericalFunctor::SetMesh( dynamic_cast<const SMESHDS_Mesh* >( theMesh ));
2126   myShapeIndex = -100;
2127   myPlane.reset();
2128 }
2129
2130 SMDSAbs_ElementType Deflection2D::GetType() const
2131 {
2132   return SMDSAbs_Face;
2133 }
2134
2135 double Deflection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
2136 {
2137   // meaningless as it is not quality control functor
2138   return Value;
2139 }
2140
2141 //================================================================================
2142 /*
2143   Class       : MultiConnection
2144   Description : Functor for calculating number of faces conneted to the edge
2145 */
2146 //================================================================================
2147
2148 double MultiConnection::GetValue( const TSequenceOfXYZ& /*P*/ )
2149 {
2150   return 0;
2151 }
2152 double MultiConnection::GetValue( long theId )
2153 {
2154   return getNbMultiConnection( myMesh, theId );
2155 }
2156
2157 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
2158 {
2159   // meaningless as it is not quality control functor
2160   return Value;
2161 }
2162
2163 SMDSAbs_ElementType MultiConnection::GetType() const
2164 {
2165   return SMDSAbs_Edge;
2166 }
2167
2168 //================================================================================
2169 /*
2170   Class       : MultiConnection2D
2171   Description : Functor for calculating number of faces conneted to the edge
2172 */
2173 //================================================================================
2174
2175 double MultiConnection2D::GetValue( const TSequenceOfXYZ& /*P*/ )
2176 {
2177   return 0;
2178 }
2179
2180 double MultiConnection2D::GetValue( long theElementId )
2181 {
2182   int aResult = 0;
2183
2184   const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId);
2185   SMDSAbs_ElementType aType = aFaceElem->GetType();
2186
2187   switch (aType) {
2188   case SMDSAbs_Face:
2189     {
2190       int i = 0, len = aFaceElem->NbNodes();
2191       SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
2192       if (!anIter) break;
2193
2194       const SMDS_MeshNode *aNode, *aNode0 = 0;
2195       TColStd_MapOfInteger aMap, aMapPrev;
2196
2197       for (i = 0; i <= len; i++) {
2198         aMapPrev = aMap;
2199         aMap.Clear();
2200
2201         int aNb = 0;
2202         if (anIter->more()) {
2203           aNode = (SMDS_MeshNode*)anIter->next();
2204         } else {
2205           if (i == len)
2206             aNode = aNode0;
2207           else
2208             break;
2209         }
2210         if (!aNode) break;
2211         if (i == 0) aNode0 = aNode;
2212
2213         SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
2214         while (anElemIter->more()) {
2215           const SMDS_MeshElement* anElem = anElemIter->next();
2216           if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
2217             int anId = anElem->GetID();
2218
2219             aMap.Add(anId);
2220             if (aMapPrev.Contains(anId)) {
2221               aNb++;
2222             }
2223           }
2224         }
2225         aResult = Max(aResult, aNb);
2226       }
2227     }
2228     break;
2229   default:
2230     aResult = 0;
2231   }
2232
2233   return aResult;
2234 }
2235
2236 double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
2237 {
2238   // meaningless as it is not quality control functor
2239   return Value;
2240 }
2241
2242 SMDSAbs_ElementType MultiConnection2D::GetType() const
2243 {
2244   return SMDSAbs_Face;
2245 }
2246
2247 MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
2248 {
2249   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
2250   if(thePntId1 > thePntId2){
2251     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
2252   }
2253 }
2254
2255 bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const
2256 {
2257   if(myPntId[0] < x.myPntId[0]) return true;
2258   if(myPntId[0] == x.myPntId[0])
2259     if(myPntId[1] < x.myPntId[1]) return true;
2260   return false;
2261 }
2262
2263 void MultiConnection2D::GetValues(MValues& theValues)
2264 {
2265   if ( !myMesh ) return;
2266   for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
2267   {
2268     const SMDS_MeshFace*     anElem = anIter->next();
2269     SMDS_NodeIteratorPtr aNodesIter = anElem->interlacedNodesIterator();
2270
2271     const SMDS_MeshNode* aNode1 = anElem->GetNode( anElem->NbNodes() - 1 );
2272     const SMDS_MeshNode* aNode2;
2273     for ( ; aNodesIter->more(); )
2274     {
2275       aNode2 = aNodesIter->next();
2276
2277       Value aValue ( aNode1->GetID(), aNode2->GetID() );
2278       MValues::iterator aItr = theValues.insert( std::make_pair( aValue, 0 )).first;
2279       aItr->second++;
2280       aNode1 = aNode2;
2281     }
2282   }
2283   return;
2284 }
2285
2286 //================================================================================
2287 /*
2288   Class       : BallDiameter
2289   Description : Functor returning diameter of a ball element
2290 */
2291 //================================================================================
2292
2293 double BallDiameter::GetValue( long theId )
2294 {
2295   double diameter = 0;
2296
2297   if ( const SMDS_BallElement* ball =
2298        myMesh->DownCast< SMDS_BallElement >( myMesh->FindElement( theId )))
2299   {
2300     diameter = ball->GetDiameter();
2301   }
2302   return diameter;
2303 }
2304
2305 double BallDiameter::GetBadRate( double Value, int /*nbNodes*/ ) const
2306 {
2307   // meaningless as it is not a quality control functor
2308   return Value;
2309 }
2310
2311 SMDSAbs_ElementType BallDiameter::GetType() const
2312 {
2313   return SMDSAbs_Ball;
2314 }
2315
2316 //================================================================================
2317 /*
2318   Class       : NodeConnectivityNumber
2319   Description : Functor returning number of elements connected to a node
2320 */
2321 //================================================================================
2322
2323 double NodeConnectivityNumber::GetValue( long theId )
2324 {
2325   double nb = 0;
2326
2327   if ( const SMDS_MeshNode* node = myMesh->FindNode( theId ))
2328   {
2329     SMDSAbs_ElementType type;
2330     if ( myMesh->NbVolumes() > 0 )
2331       type = SMDSAbs_Volume;
2332     else if ( myMesh->NbFaces() > 0 )
2333       type = SMDSAbs_Face;
2334     else if ( myMesh->NbEdges() > 0 )
2335       type = SMDSAbs_Edge;
2336     else
2337       return 0;
2338     nb = node->NbInverseElements( type );
2339   }
2340   return nb;
2341 }
2342
2343 double NodeConnectivityNumber::GetBadRate( double Value, int /*nbNodes*/ ) const
2344 {
2345   return Value;
2346 }
2347
2348 SMDSAbs_ElementType NodeConnectivityNumber::GetType() const
2349 {
2350   return SMDSAbs_Node;
2351 }
2352
2353 /*
2354                             PREDICATES
2355 */
2356
2357 //================================================================================
2358 /*
2359   Class       : BadOrientedVolume
2360   Description : Predicate bad oriented volumes
2361 */
2362 //================================================================================
2363
2364 BadOrientedVolume::BadOrientedVolume()
2365 {
2366   myMesh = 0;
2367 }
2368
2369 void BadOrientedVolume::SetMesh( const SMDS_Mesh* theMesh )
2370 {
2371   myMesh = theMesh;
2372 }
2373
2374 bool BadOrientedVolume::IsSatisfy( long theId )
2375 {
2376   if ( myMesh == 0 )
2377     return false;
2378
2379   SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
2380
2381   bool isOk = true;
2382   if ( vTool.IsPoly() )
2383   {
2384     isOk = true;
2385     for ( int i = 0; i < vTool.NbFaces() && isOk; ++i )
2386       isOk = vTool.IsFaceExternal( i );
2387   }
2388   else
2389   {
2390     isOk = vTool.IsForward();
2391   }
2392   return !isOk;
2393 }
2394
2395 SMDSAbs_ElementType BadOrientedVolume::GetType() const
2396 {
2397   return SMDSAbs_Volume;
2398 }
2399
2400 /*
2401   Class       : BareBorderVolume
2402 */
2403
2404 bool BareBorderVolume::IsSatisfy(long theElementId )
2405 {
2406   SMDS_VolumeTool  myTool;
2407   if ( myTool.Set( myMesh->FindElement(theElementId)))
2408   {
2409     for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
2410       if ( myTool.IsFreeFace( iF ))
2411       {
2412         const SMDS_MeshNode** n = myTool.GetFaceNodes(iF);
2413         std::vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF));
2414         if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*Nomedium=*/false))
2415           return true;
2416       }
2417   }
2418   return false;
2419 }
2420
2421 //================================================================================
2422 /*
2423   Class       : BareBorderFace
2424 */
2425 //================================================================================
2426
2427 bool BareBorderFace::IsSatisfy(long theElementId )
2428 {
2429   bool ok = false;
2430   if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
2431   {
2432     if ( face->GetType() == SMDSAbs_Face )
2433     {
2434       int nbN = face->NbCornerNodes();
2435       for ( int i = 0; i < nbN && !ok; ++i )
2436       {
2437         // check if a link is shared by another face
2438         const SMDS_MeshNode* n1 = face->GetNode( i );
2439         const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
2440         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
2441         bool isShared = false;
2442         while ( !isShared && fIt->more() )
2443         {
2444           const SMDS_MeshElement* f = fIt->next();
2445           isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
2446         }
2447         if ( !isShared )
2448         {
2449           const int iQuad = face->IsQuadratic();
2450           myLinkNodes.resize( 2 + iQuad);
2451           myLinkNodes[0] = n1;
2452           myLinkNodes[1] = n2;
2453           if ( iQuad )
2454             myLinkNodes[2] = face->GetNode( i+nbN );
2455           ok = !myMesh->FindElement( myLinkNodes, SMDSAbs_Edge, /*noMedium=*/false);
2456         }
2457       }
2458     }
2459   }
2460   return ok;
2461 }
2462
2463 //================================================================================
2464 /*
2465   Class       : OverConstrainedVolume
2466 */
2467 //================================================================================
2468
2469 bool OverConstrainedVolume::IsSatisfy(long theElementId )
2470 {
2471   // An element is over-constrained if it has N-1 free borders where
2472   // N is the number of edges/faces for a 2D/3D element.
2473   SMDS_VolumeTool  myTool;
2474   if ( myTool.Set( myMesh->FindElement(theElementId)))
2475   {
2476     int nbSharedFaces = 0;
2477     for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
2478       if ( !myTool.IsFreeFace( iF ) && ++nbSharedFaces > 1 )
2479         break;
2480     return ( nbSharedFaces == 1 );
2481   }
2482   return false;
2483 }
2484
2485 //================================================================================
2486 /*
2487   Class       : OverConstrainedFace
2488 */
2489 //================================================================================
2490
2491 bool OverConstrainedFace::IsSatisfy(long theElementId )
2492 {
2493   // An element is over-constrained if it has N-1 free borders where
2494   // N is the number of edges/faces for a 2D/3D element.
2495   if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
2496     if ( face->GetType() == SMDSAbs_Face )
2497     {
2498       int nbSharedBorders = 0;
2499       int nbN = face->NbCornerNodes();
2500       for ( int i = 0; i < nbN; ++i )
2501       {
2502         // check if a link is shared by another face
2503         const SMDS_MeshNode* n1 = face->GetNode( i );
2504         const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
2505         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
2506         bool isShared = false;
2507         while ( !isShared && fIt->more() )
2508         {
2509           const SMDS_MeshElement* f = fIt->next();
2510           isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
2511         }
2512         if ( isShared && ++nbSharedBorders > 1 )
2513           break;
2514       }
2515       return ( nbSharedBorders == 1 );
2516     }
2517   return false;
2518 }
2519
2520 //================================================================================
2521 /*
2522   Class       : CoincidentNodes
2523   Description : Predicate of Coincident nodes
2524 */
2525 //================================================================================
2526
2527 CoincidentNodes::CoincidentNodes()
2528 {
2529   myToler = 1e-5;
2530 }
2531
2532 bool CoincidentNodes::IsSatisfy( long theElementId )
2533 {
2534   return myCoincidentIDs.Contains( theElementId );
2535 }
2536
2537 SMDSAbs_ElementType CoincidentNodes::GetType() const
2538 {
2539   return SMDSAbs_Node;
2540 }
2541
2542 void CoincidentNodes::SetTolerance( const double theToler )
2543 {
2544   if ( myToler != theToler )
2545   {
2546     SetMesh(0);
2547     myToler = theToler;
2548   }
2549 }
2550
2551 void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
2552 {
2553   myMeshModifTracer.SetMesh( theMesh );
2554   if ( myMeshModifTracer.IsMeshModified() )
2555   {
2556     TIDSortedNodeSet nodesToCheck;
2557     SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
2558     while ( nIt->more() )
2559       nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
2560
2561     std::list< std::list< const SMDS_MeshNode*> > nodeGroups;
2562     SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler );
2563
2564     myCoincidentIDs.Clear();
2565     std::list< std::list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin();
2566     for ( ; groupIt != nodeGroups.end(); ++groupIt )
2567     {
2568       std::list< const SMDS_MeshNode*>& coincNodes = *groupIt;
2569       std::list< const SMDS_MeshNode*>::iterator n = coincNodes.begin();
2570       for ( ; n != coincNodes.end(); ++n )
2571         myCoincidentIDs.Add( (*n)->GetID() );
2572     }
2573   }
2574 }
2575
2576 //================================================================================
2577 /*
2578   Class       : CoincidentElements
2579   Description : Predicate of Coincident Elements
2580   Note        : This class is suitable only for visualization of Coincident Elements
2581 */
2582 //================================================================================
2583
2584 CoincidentElements::CoincidentElements()
2585 {
2586   myMesh = 0;
2587 }
2588
2589 void CoincidentElements::SetMesh( const SMDS_Mesh* theMesh )
2590 {
2591   myMesh = theMesh;
2592 }
2593
2594 bool CoincidentElements::IsSatisfy( long theElementId )
2595 {
2596   if ( !myMesh ) return false;
2597
2598   if ( const SMDS_MeshElement* e = myMesh->FindElement( theElementId ))
2599   {
2600     if ( e->GetType() != GetType() ) return false;
2601     std::set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() );
2602     const int nbNodes = e->NbNodes();
2603     SMDS_ElemIteratorPtr invIt = (*elemNodes.begin())->GetInverseElementIterator( GetType() );
2604     while ( invIt->more() )
2605     {
2606       const SMDS_MeshElement* e2 = invIt->next();
2607       if ( e2 == e || e2->NbNodes() != nbNodes ) continue;
2608
2609       bool sameNodes = true;
2610       for ( size_t i = 0; i < elemNodes.size() && sameNodes; ++i )
2611         sameNodes = ( elemNodes.count( e2->GetNode( i )));
2612       if ( sameNodes )
2613         return true;
2614     }
2615   }
2616   return false;
2617 }
2618
2619 SMDSAbs_ElementType CoincidentElements1D::GetType() const
2620 {
2621   return SMDSAbs_Edge;
2622 }
2623 SMDSAbs_ElementType CoincidentElements2D::GetType() const
2624 {
2625   return SMDSAbs_Face;
2626 }
2627 SMDSAbs_ElementType CoincidentElements3D::GetType() const
2628 {
2629   return SMDSAbs_Volume;
2630 }
2631
2632
2633 //================================================================================
2634 /*
2635   Class       : FreeBorders
2636   Description : Predicate for free borders
2637 */
2638 //================================================================================
2639
2640 FreeBorders::FreeBorders()
2641 {
2642   myMesh = 0;
2643 }
2644
2645 void FreeBorders::SetMesh( const SMDS_Mesh* theMesh )
2646 {
2647   myMesh = theMesh;
2648 }
2649
2650 bool FreeBorders::IsSatisfy( long theId )
2651 {
2652   return getNbMultiConnection( myMesh, theId ) == 1;
2653 }
2654
2655 SMDSAbs_ElementType FreeBorders::GetType() const
2656 {
2657   return SMDSAbs_Edge;
2658 }
2659
2660
2661 //================================================================================
2662 /*
2663   Class       : FreeEdges
2664   Description : Predicate for free Edges
2665 */
2666 //================================================================================
2667
2668 FreeEdges::FreeEdges()
2669 {
2670   myMesh = 0;
2671 }
2672
2673 void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
2674 {
2675   myMesh = theMesh;
2676 }
2677
2678 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId  )
2679 {
2680   SMDS_ElemIteratorPtr anElemIter = theNodes[ 0 ]->GetInverseElementIterator(SMDSAbs_Face);
2681   while( anElemIter->more() )
2682   {
2683     if ( const SMDS_MeshElement* anElem = anElemIter->next())
2684     {
2685       const int anId = anElem->GetID();
2686       if ( anId != theFaceId && anElem->GetNodeIndex( theNodes[1] ) >= 0 )
2687         return false;
2688     }
2689   }
2690   return true;
2691 }
2692
2693 bool FreeEdges::IsSatisfy( long theId )
2694 {
2695   if ( myMesh == 0 )
2696     return false;
2697
2698   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2699   if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
2700     return false;
2701
2702   SMDS_NodeIteratorPtr anIter = aFace->interlacedNodesIterator();
2703   if ( !anIter )
2704     return false;
2705
2706   int i = 0, nbNodes = aFace->NbNodes();
2707   std::vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
2708   while( anIter->more() )
2709     if ( ! ( aNodes[ i++ ] = anIter->next() ))
2710       return false;
2711   aNodes[ nbNodes ] = aNodes[ 0 ];
2712
2713   for ( i = 0; i < nbNodes; i++ )
2714     if ( IsFreeEdge( &aNodes[ i ], theId ) )
2715       return true;
2716
2717   return false;
2718 }
2719
2720 SMDSAbs_ElementType FreeEdges::GetType() const
2721 {
2722   return SMDSAbs_Face;
2723 }
2724
2725 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
2726   myElemId(theElemId)
2727 {
2728   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
2729   if(thePntId1 > thePntId2){
2730     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
2731   }
2732 }
2733
2734 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
2735   if(myPntId[0] < x.myPntId[0]) return true;
2736   if(myPntId[0] == x.myPntId[0])
2737     if(myPntId[1] < x.myPntId[1]) return true;
2738   return false;
2739 }
2740
2741 inline void UpdateBorders(const FreeEdges::Border& theBorder,
2742                           FreeEdges::TBorders& theRegistry,
2743                           FreeEdges::TBorders& theContainer)
2744 {
2745   if(theRegistry.find(theBorder) == theRegistry.end()){
2746     theRegistry.insert(theBorder);
2747     theContainer.insert(theBorder);
2748   }else{
2749     theContainer.erase(theBorder);
2750   }
2751 }
2752
2753 void FreeEdges::GetBoreders(TBorders& theBorders)
2754 {
2755   TBorders aRegistry;
2756   for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
2757   {
2758     const SMDS_MeshFace* anElem = anIter->next();
2759     long anElemId = anElem->GetID();
2760     SMDS_NodeIteratorPtr aNodesIter = anElem->interlacedNodesIterator();
2761     if ( !aNodesIter->more() ) continue;
2762     long aNodeId[2] = {0,0};
2763     aNodeId[0] = anElem->GetNode( anElem->NbNodes()-1 )->GetID();
2764     for ( ; aNodesIter->more(); )
2765     {
2766       aNodeId[1] = aNodesIter->next()->GetID();
2767       Border aBorder( anElemId, aNodeId[0], aNodeId[1] );
2768       UpdateBorders( aBorder, aRegistry, theBorders );
2769       aNodeId[0] = aNodeId[1];
2770     }
2771   }
2772 }
2773
2774 //================================================================================
2775 /*
2776   Class       : FreeNodes
2777   Description : Predicate for free nodes
2778 */
2779 //================================================================================
2780
2781 FreeNodes::FreeNodes()
2782 {
2783   myMesh = 0;
2784 }
2785
2786 void FreeNodes::SetMesh( const SMDS_Mesh* theMesh )
2787 {
2788   myMesh = theMesh;
2789 }
2790
2791 bool FreeNodes::IsSatisfy( long theNodeId )
2792 {
2793   const SMDS_MeshNode* aNode = myMesh->FindNode( theNodeId );
2794   if (!aNode)
2795     return false;
2796
2797   return (aNode->NbInverseElements() < 1);
2798 }
2799
2800 SMDSAbs_ElementType FreeNodes::GetType() const
2801 {
2802   return SMDSAbs_Node;
2803 }
2804
2805
2806 //================================================================================
2807 /*
2808   Class       : FreeFaces
2809   Description : Predicate for free faces
2810 */
2811 //================================================================================
2812
2813 FreeFaces::FreeFaces()
2814 {
2815   myMesh = 0;
2816 }
2817
2818 void FreeFaces::SetMesh( const SMDS_Mesh* theMesh )
2819 {
2820   myMesh = theMesh;
2821 }
2822
2823 bool FreeFaces::IsSatisfy( long theId )
2824 {
2825   if (!myMesh) return false;
2826   // check that faces nodes refers to less than two common volumes
2827   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2828   if ( !aFace || aFace->GetType() != SMDSAbs_Face )
2829     return false;
2830
2831   int nbNode = aFace->NbNodes();
2832
2833   // collect volumes to check that number of volumes with count equal nbNode not less than 2
2834   typedef std::map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
2835   typedef std::map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
2836   TMapOfVolume mapOfVol;
2837
2838   SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator();
2839   while ( nodeItr->more() )
2840   {
2841     const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
2842     if ( !aNode ) continue;
2843     SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume);
2844     while ( volItr->more() )
2845     {
2846       SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
2847       TItrMapOfVolume    itr = mapOfVol.insert( std::make_pair( aVol, 0 )).first;
2848       (*itr).second++;
2849     }
2850   }
2851   int nbVol = 0;
2852   TItrMapOfVolume volItr = mapOfVol.begin();
2853   TItrMapOfVolume volEnd = mapOfVol.end();
2854   for ( ; volItr != volEnd; ++volItr )
2855     if ( (*volItr).second >= nbNode )
2856       nbVol++;
2857   // face is not free if number of volumes constructed on their nodes more than one
2858   return (nbVol < 2);
2859 }
2860
2861 SMDSAbs_ElementType FreeFaces::GetType() const
2862 {
2863   return SMDSAbs_Face;
2864 }
2865
2866 //================================================================================
2867 /*
2868   Class       : LinearOrQuadratic
2869   Description : Predicate to verify whether a mesh element is linear
2870 */
2871 //================================================================================
2872
2873 LinearOrQuadratic::LinearOrQuadratic()
2874 {
2875   myMesh = 0;
2876 }
2877
2878 void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh )
2879 {
2880   myMesh = theMesh;
2881 }
2882
2883 bool LinearOrQuadratic::IsSatisfy( long theId )
2884 {
2885   if (!myMesh) return false;
2886   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2887   if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) )
2888     return false;
2889   return (!anElem->IsQuadratic());
2890 }
2891
2892 void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType )
2893 {
2894   myType = theType;
2895 }
2896
2897 SMDSAbs_ElementType LinearOrQuadratic::GetType() const
2898 {
2899   return myType;
2900 }
2901
2902 //================================================================================
2903 /*
2904   Class       : GroupColor
2905   Description : Functor for check color of group to which mesh element belongs to
2906 */
2907 //================================================================================
2908
2909 GroupColor::GroupColor()
2910 {
2911 }
2912
2913 bool GroupColor::IsSatisfy( long theId )
2914 {
2915   return myIDs.count( theId );
2916 }
2917
2918 void GroupColor::SetType( SMDSAbs_ElementType theType )
2919 {
2920   myType = theType;
2921 }
2922
2923 SMDSAbs_ElementType GroupColor::GetType() const
2924 {
2925   return myType;
2926 }
2927
2928 static bool isEqual( const Quantity_Color& theColor1,
2929                      const Quantity_Color& theColor2 )
2930 {
2931   // tolerance to compare colors
2932   const double tol = 5*1e-3;
2933   return ( fabs( theColor1.Red()   - theColor2.Red() )   < tol &&
2934            fabs( theColor1.Green() - theColor2.Green() ) < tol &&
2935            fabs( theColor1.Blue()  - theColor2.Blue() )  < tol );
2936 }
2937
2938 void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
2939 {
2940   myIDs.clear();
2941
2942   const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
2943   if ( !aMesh )
2944     return;
2945
2946   int nbGrp = aMesh->GetNbGroups();
2947   if ( !nbGrp )
2948     return;
2949
2950   // iterates on groups and find necessary elements ids
2951   const std::set<SMESHDS_GroupBase*>&       aGroups = aMesh->GetGroups();
2952   std::set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
2953   for (; GrIt != aGroups.end(); GrIt++)
2954   {
2955     SMESHDS_GroupBase* aGrp = (*GrIt);
2956     if ( !aGrp )
2957       continue;
2958     // check type and color of group
2959     if ( !isEqual( myColor, aGrp->GetColor() ))
2960       continue;
2961
2962     // IPAL52867 (prevent infinite recursion via GroupOnFilter)
2963     if ( SMESHDS_GroupOnFilter * gof = dynamic_cast< SMESHDS_GroupOnFilter* >( aGrp ))
2964       if ( gof->GetPredicate().get() == this )
2965         continue;
2966
2967     SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType();
2968     if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) {
2969       // add elements IDS into control
2970       int aSize = aGrp->Extent();
2971       for (int i = 0; i < aSize; i++)
2972         myIDs.insert( aGrp->GetID(i+1) );
2973     }
2974   }
2975 }
2976
2977 void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
2978 {
2979   Kernel_Utils::Localizer loc;
2980   TCollection_AsciiString aStr = theStr;
2981   aStr.RemoveAll( ' ' );
2982   aStr.RemoveAll( '\t' );
2983   for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) )
2984     aStr.Remove( aPos, 2 );
2985   Standard_Real clr[3];
2986   clr[0] = clr[1] = clr[2] = 0.;
2987   for ( int i = 0; i < 3; i++ ) {
2988     TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 );
2989     if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() )
2990       clr[i] = tmpStr.RealValue();
2991   }
2992   myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB );
2993 }
2994
2995 //=======================================================================
2996 // name    : GetRangeStr
2997 // Purpose : Get range as a string.
2998 //           Example: "1,2,3,50-60,63,67,70-"
2999 //=======================================================================
3000
3001 void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
3002 {
3003   theResStr.Clear();
3004   theResStr += TCollection_AsciiString( myColor.Red() );
3005   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() );
3006   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
3007 }
3008
3009 //================================================================================
3010 /*
3011   Class       : ElemGeomType
3012   Description : Predicate to check element geometry type
3013 */
3014 //================================================================================
3015
3016 ElemGeomType::ElemGeomType()
3017 {
3018   myMesh = 0;
3019   myType = SMDSAbs_All;
3020   myGeomType = SMDSGeom_TRIANGLE;
3021 }
3022
3023 void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh )
3024 {
3025   myMesh = theMesh;
3026 }
3027
3028 bool ElemGeomType::IsSatisfy( long theId )
3029 {
3030   if (!myMesh) return false;
3031   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3032   if ( !anElem )
3033     return false;
3034   const SMDSAbs_ElementType anElemType = anElem->GetType();
3035   if ( myType != SMDSAbs_All && anElemType != myType )
3036     return false;
3037   bool isOk = ( anElem->GetGeomType() == myGeomType );
3038   return isOk;
3039 }
3040
3041 void ElemGeomType::SetType( SMDSAbs_ElementType theType )
3042 {
3043   myType = theType;
3044 }
3045
3046 SMDSAbs_ElementType ElemGeomType::GetType() const
3047 {
3048   return myType;
3049 }
3050
3051 void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType )
3052 {
3053   myGeomType = theType;
3054 }
3055
3056 SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
3057 {
3058   return myGeomType;
3059 }
3060
3061 //================================================================================
3062 /*
3063   Class       : ElemEntityType
3064   Description : Predicate to check element entity type
3065 */
3066 //================================================================================
3067
3068 ElemEntityType::ElemEntityType():
3069   myMesh( 0 ),
3070   myType( SMDSAbs_All ),
3071   myEntityType( SMDSEntity_0D )
3072 {
3073 }
3074
3075 void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh )
3076 {
3077   myMesh = theMesh;
3078 }
3079
3080 bool ElemEntityType::IsSatisfy( long theId )
3081 {
3082   if ( !myMesh ) return false;
3083   if ( myType == SMDSAbs_Node )
3084     return myMesh->FindNode( theId );
3085   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3086   return ( anElem &&
3087            myEntityType == anElem->GetEntityType() );
3088 }
3089
3090 void ElemEntityType::SetType( SMDSAbs_ElementType theType )
3091 {
3092   myType = theType;
3093 }
3094
3095 SMDSAbs_ElementType ElemEntityType::GetType() const
3096 {
3097   return myType;
3098 }
3099
3100 void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType )
3101 {
3102   myEntityType = theEntityType;
3103 }
3104
3105 SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const
3106 {
3107   return myEntityType;
3108 }
3109
3110 //================================================================================
3111 /*!
3112  * \brief Class ConnectedElements
3113  */
3114 //================================================================================
3115
3116 ConnectedElements::ConnectedElements():
3117   myNodeID(0), myType( SMDSAbs_All ), myOkIDsReady( false ) {}
3118
3119 SMDSAbs_ElementType ConnectedElements::GetType() const
3120 { return myType; }
3121
3122 int ConnectedElements::GetNode() const
3123 { return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ
3124
3125 std::vector<double> ConnectedElements::GetPoint() const
3126 { return myXYZ; }
3127
3128 void ConnectedElements::clearOkIDs()
3129 { myOkIDsReady = false; myOkIDs.clear(); }
3130
3131 void ConnectedElements::SetType( SMDSAbs_ElementType theType )
3132 {
3133   if ( myType != theType || myMeshModifTracer.IsMeshModified() )
3134     clearOkIDs();
3135   myType = theType;
3136 }
3137
3138 void ConnectedElements::SetMesh( const SMDS_Mesh* theMesh )
3139 {
3140   myMeshModifTracer.SetMesh( theMesh );
3141   if ( myMeshModifTracer.IsMeshModified() )
3142   {
3143     clearOkIDs();
3144     if ( !myXYZ.empty() )
3145       SetPoint( myXYZ[0], myXYZ[1], myXYZ[2] ); // find a node near myXYZ it in a new mesh
3146   }
3147 }
3148
3149 void ConnectedElements::SetNode( int nodeID )
3150 {
3151   myNodeID = nodeID;
3152   myXYZ.clear();
3153
3154   bool isSameDomain = false;
3155   if ( myOkIDsReady && myMeshModifTracer.GetMesh() && !myMeshModifTracer.IsMeshModified() )
3156     if ( const SMDS_MeshNode* n = myMeshModifTracer.GetMesh()->FindNode( myNodeID ))
3157     {
3158       SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( myType );
3159       while ( !isSameDomain && eIt->more() )
3160         isSameDomain = IsSatisfy( eIt->next()->GetID() );
3161     }
3162   if ( !isSameDomain )
3163     clearOkIDs();
3164 }
3165
3166 void ConnectedElements::SetPoint( double x, double y, double z )
3167 {
3168   myXYZ.resize(3);
3169   myXYZ[0] = x;
3170   myXYZ[1] = y;
3171   myXYZ[2] = z;
3172   myNodeID = 0;
3173
3174   bool isSameDomain = false;
3175
3176   // find myNodeID by myXYZ if possible
3177   if ( myMeshModifTracer.GetMesh() )
3178   {
3179     SMESHUtils::Deleter<SMESH_ElementSearcher> searcher
3180       ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() ));
3181
3182     std::vector< const SMDS_MeshElement* > foundElems;
3183     searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems );
3184
3185     if ( !foundElems.empty() )
3186     {
3187       myNodeID = foundElems[0]->GetNode(0)->GetID();
3188       if ( myOkIDsReady && !myMeshModifTracer.IsMeshModified() )
3189         isSameDomain = IsSatisfy( foundElems[0]->GetID() );
3190     }
3191   }
3192   if ( !isSameDomain )
3193     clearOkIDs();
3194 }
3195
3196 bool ConnectedElements::IsSatisfy( long theElementId )
3197 {
3198   // Here we do NOT check if the mesh has changed, we do it in Set...() only!!!
3199
3200   if ( !myOkIDsReady )
3201   {