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