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