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