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