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