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