Back to index

salome-smesh  6.5.0
SMESH_Controls.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 #include "SMESH_ControlsDef.hxx"
00024 
00025 #include <set>
00026 #include <limits>
00027 
00028 #include <BRepAdaptor_Surface.hxx>
00029 #include <BRepClass_FaceClassifier.hxx>
00030 #include <BRep_Tool.hxx>
00031 
00032 #include <TopAbs.hxx>
00033 #include <TopoDS.hxx>
00034 #include <TopoDS_Edge.hxx>
00035 #include <TopoDS_Face.hxx>
00036 #include <TopoDS_Shape.hxx>
00037 #include <TopoDS_Vertex.hxx>
00038 #include <TopoDS_Iterator.hxx>
00039 
00040 #include <Geom_CylindricalSurface.hxx>
00041 #include <Geom_Plane.hxx>
00042 #include <Geom_Surface.hxx>
00043 
00044 #include <Precision.hxx>
00045 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
00046 #include <TColStd_MapOfInteger.hxx>
00047 #include <TColStd_SequenceOfAsciiString.hxx>
00048 #include <TColgp_Array1OfXYZ.hxx>
00049 
00050 #include <gp_Ax3.hxx>
00051 #include <gp_Cylinder.hxx>
00052 #include <gp_Dir.hxx>
00053 #include <gp_Pln.hxx>
00054 #include <gp_Pnt.hxx>
00055 #include <gp_Vec.hxx>
00056 #include <gp_XYZ.hxx>
00057 
00058 #include "SMDS_Mesh.hxx"
00059 #include "SMDS_Iterator.hxx"
00060 #include "SMDS_MeshElement.hxx"
00061 #include "SMDS_MeshNode.hxx"
00062 #include "SMDS_VolumeTool.hxx"
00063 #include "SMDS_QuadraticFaceOfNodes.hxx"
00064 #include "SMDS_QuadraticEdge.hxx"
00065 
00066 #include "SMESHDS_Mesh.hxx"
00067 #include "SMESHDS_GroupBase.hxx"
00068 
00069 #include "SMESH_OctreeNode.hxx"
00070 
00071 #include <vtkMeshQuality.h>
00072 
00073 /*
00074                             AUXILIARY METHODS
00075 */
00076 
00077 namespace{
00078 
00079   inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode )
00080   {
00081     return gp_XYZ(aNode->X(), aNode->Y(), aNode->Z() );
00082   }
00083 
00084   inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
00085   {
00086     gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
00087 
00088     return v1.Magnitude() < gp::Resolution() ||
00089       v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
00090   }
00091 
00092   inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
00093   {
00094     gp_Vec aVec1( P2 - P1 );
00095     gp_Vec aVec2( P3 - P1 );
00096     return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
00097   }
00098 
00099   inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
00100   {
00101     return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
00102   }
00103 
00104 
00105 
00106   inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
00107   {
00108     double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
00109     return aDist;
00110   }
00111 
00112   int getNbMultiConnection( const SMDS_Mesh* theMesh, const int theId )
00113   {
00114     if ( theMesh == 0 )
00115       return 0;
00116 
00117     const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
00118     if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge/* || anEdge->NbNodes() != 2 */)
00119       return 0;
00120 
00121     // for each pair of nodes in anEdge (there are 2 pairs in a quadratic edge)
00122     // count elements containing both nodes of the pair.
00123     // Note that there may be such cases for a quadratic edge (a horizontal line):
00124     //
00125     //  Case 1          Case 2
00126     //  |     |      |        |      |
00127     //  |     |      |        |      |
00128     //  +-----+------+  +-----+------+ 
00129     //  |            |  |            |
00130     //  |            |  |            |
00131     // result sould be 2 in both cases
00132     //
00133     int aResult0 = 0, aResult1 = 0;
00134      // last node, it is a medium one in a quadratic edge
00135     const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 );
00136     const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 );
00137     const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 );
00138     if ( aNode1 == aLastNode ) aNode1 = 0;
00139 
00140     SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator();
00141     while( anElemIter->more() ) {
00142       const SMDS_MeshElement* anElem = anElemIter->next();
00143       if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
00144         SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
00145         while ( anIter->more() ) {
00146           if ( const SMDS_MeshElement* anElemNode = anIter->next() ) {
00147             if ( anElemNode == aNode0 ) {
00148               aResult0++;
00149               if ( !aNode1 ) break; // not a quadratic edge
00150             }
00151             else if ( anElemNode == aNode1 )
00152               aResult1++;
00153           }
00154         }
00155       }
00156     }
00157     int aResult = std::max ( aResult0, aResult1 );
00158 
00159 //     TColStd_MapOfInteger aMap;
00160 
00161 //     SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
00162 //     if ( anIter != 0 ) {
00163 //       while( anIter->more() ) {
00164 //      const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
00165 //      if ( aNode == 0 )
00166 //        return 0;
00167 //      SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
00168 //      while( anElemIter->more() ) {
00169 //        const SMDS_MeshElement* anElem = anElemIter->next();
00170 //        if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
00171 //          int anId = anElem->GetID();
00172 
00173 //          if ( anIter->more() )              // i.e. first node
00174 //            aMap.Add( anId );
00175 //          else if ( aMap.Contains( anId ) )
00176 //            aResult++;
00177 //        }
00178 //      }
00179 //       }
00180 //     }
00181 
00182     return aResult;
00183   }
00184 
00185   gp_XYZ getNormale( const SMDS_MeshFace* theFace, bool* ok=0 )
00186   {
00187     int aNbNode = theFace->NbNodes();
00188 
00189     gp_XYZ q1 = gpXYZ( theFace->GetNode(1)) - gpXYZ( theFace->GetNode(0));
00190     gp_XYZ q2 = gpXYZ( theFace->GetNode(2)) - gpXYZ( theFace->GetNode(0));
00191     gp_XYZ n  = q1 ^ q2;
00192     if ( aNbNode > 3 ) {
00193       gp_XYZ q3 = gpXYZ( theFace->GetNode(3)) - gpXYZ( theFace->GetNode(0));
00194       n += q2 ^ q3;
00195     }
00196     double len = n.Modulus();
00197     bool zeroLen = ( len <= numeric_limits<double>::min());
00198     if ( !zeroLen )
00199       n /= len;
00200 
00201     if (ok) *ok = !zeroLen;
00202 
00203     return n;
00204   }
00205 }
00206 
00207 
00208 
00209 using namespace SMESH::Controls;
00210 
00211 /*
00212  *                               FUNCTORS
00213  */
00214 
00215 /*
00216   Class       : NumericalFunctor
00217   Description : Base class for numerical functors
00218 */
00219 NumericalFunctor::NumericalFunctor():
00220   myMesh(NULL)
00221 {
00222   myPrecision = -1;
00223 }
00224 
00225 void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh )
00226 {
00227   myMesh = theMesh;
00228 }
00229 
00230 bool NumericalFunctor::GetPoints(const int theId,
00231                                  TSequenceOfXYZ& theRes ) const
00232 {
00233   theRes.clear();
00234 
00235   if ( myMesh == 0 )
00236     return false;
00237 
00238   return GetPoints( myMesh->FindElement( theId ), theRes );
00239 }
00240 
00241 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
00242                                  TSequenceOfXYZ&         theRes )
00243 {
00244   theRes.clear();
00245 
00246   if ( anElem == 0)
00247     return false;
00248 
00249   theRes.reserve( anElem->NbNodes() );
00250 
00251   // Get nodes of the element
00252   SMDS_ElemIteratorPtr anIter;
00253 
00254   if ( anElem->IsQuadratic() ) {
00255     switch ( anElem->GetType() ) {
00256     case SMDSAbs_Edge:
00257       anIter = dynamic_cast<const SMDS_VtkEdge*>
00258         (anElem)->interlacedNodesElemIterator();
00259       break;
00260     case SMDSAbs_Face:
00261       anIter = dynamic_cast<const SMDS_VtkFace*>
00262         (anElem)->interlacedNodesElemIterator();
00263       break;
00264     default:
00265       anIter = anElem->nodesIterator();
00266       //return false;
00267     }
00268   }
00269   else {
00270     anIter = anElem->nodesIterator();
00271   }
00272 
00273   if ( anIter ) {
00274     while( anIter->more() ) {
00275       if ( const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( anIter->next() ))
00276         theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
00277     }
00278   }
00279 
00280   return true;
00281 }
00282 
00283 long  NumericalFunctor::GetPrecision() const
00284 {
00285   return myPrecision;
00286 }
00287 
00288 void  NumericalFunctor::SetPrecision( const long thePrecision )
00289 {
00290   myPrecision = thePrecision;
00291   myPrecisionValue = pow( 10., (double)( myPrecision ) );
00292 }
00293 
00294 double NumericalFunctor::GetValue( long theId )
00295 {
00296   double aVal = 0;
00297 
00298   myCurrElement = myMesh->FindElement( theId );
00299 
00300   TSequenceOfXYZ P;
00301   if ( GetPoints( theId, P ))
00302     aVal = Round( GetValue( P ));
00303 
00304   return aVal;
00305 }
00306 
00307 double NumericalFunctor::Round( const double & aVal )
00308 {
00309   return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal;
00310 }
00311 
00312 //================================================================================
00321 //================================================================================
00322 
00323 void NumericalFunctor::GetHistogram(int                  nbIntervals,
00324                                     std::vector<int>&    nbEvents,
00325                                     std::vector<double>& funValues,
00326                                     const vector<int>&   elements,
00327                                     const double*        minmax)
00328 {
00329   if ( nbIntervals < 1 ||
00330        !myMesh ||
00331        !myMesh->GetMeshInfo().NbElements( GetType() ))
00332     return;
00333   nbEvents.resize( nbIntervals, 0 );
00334   funValues.resize( nbIntervals+1 );
00335 
00336   // get all values sorted
00337   std::multiset< double > values;
00338   if ( elements.empty() )
00339   {
00340     SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator(GetType());
00341     while ( elemIt->more() )
00342       values.insert( GetValue( elemIt->next()->GetID() ));
00343   }
00344   else
00345   {
00346     vector<int>::const_iterator id = elements.begin();
00347     for ( ; id != elements.end(); ++id )
00348       values.insert( GetValue( *id ));
00349   }
00350 
00351   if ( minmax )
00352   {
00353     funValues[0] = minmax[0];
00354     funValues[nbIntervals] = minmax[1];
00355   }
00356   else
00357   {
00358     funValues[0] = *values.begin();
00359     funValues[nbIntervals] = *values.rbegin();
00360   }
00361   // case nbIntervals == 1
00362   if ( nbIntervals == 1 )
00363   {
00364     nbEvents[0] = values.size();
00365     return;
00366   }
00367   // case of 1 value
00368   if (funValues.front() == funValues.back())
00369   {
00370     nbEvents.resize( 1 );
00371     nbEvents[0] = values.size();
00372     funValues[1] = funValues.back();
00373     funValues.resize( 2 );
00374   }
00375   // generic case
00376   std::multiset< double >::iterator min = values.begin(), max;
00377   for ( int i = 0; i < nbIntervals; ++i )
00378   {
00379     // find end value of i-th interval
00380     double r = (i+1) / double( nbIntervals );
00381     funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
00382 
00383     // count values in the i-th interval if there are any
00384     if ( min != values.end() && *min <= funValues[i+1] )
00385     {
00386       // find the first value out of the interval
00387       max = values.upper_bound( funValues[i+1] ); // max is greater than funValues[i+1], or end()
00388       nbEvents[i] = std::distance( min, max );
00389       min = max;
00390     }
00391   }
00392   // add values larger than minmax[1]
00393   nbEvents.back() += std::distance( min, values.end() );
00394 }
00395 
00396 //=======================================================================
00397 //function : GetValue
00398 //purpose  : 
00399 //=======================================================================
00400 
00401 double Volume::GetValue( long theElementId )
00402 {
00403   if ( theElementId && myMesh ) {
00404     SMDS_VolumeTool aVolumeTool;
00405     if ( aVolumeTool.Set( myMesh->FindElement( theElementId )))
00406       return aVolumeTool.GetSize();
00407   }
00408   return 0;
00409 }
00410 
00411 //=======================================================================
00412 //function : GetBadRate
00413 //purpose  : meaningless as it is not quality control functor
00414 //=======================================================================
00415 
00416 double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const
00417 {
00418   return Value;
00419 }
00420 
00421 //=======================================================================
00422 //function : GetType
00423 //purpose  : 
00424 //=======================================================================
00425 
00426 SMDSAbs_ElementType Volume::GetType() const
00427 {
00428   return SMDSAbs_Volume;
00429 }
00430 
00431 
00432 /*
00433   Class       : MaxElementLength2D
00434   Description : Functor calculating maximum length of 2D element
00435 */
00436 
00437 double MaxElementLength2D::GetValue( long theElementId )
00438 {
00439   TSequenceOfXYZ P;
00440   if( GetPoints( theElementId, P ) ) {
00441     double aVal = 0;
00442     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
00443     SMDSAbs_ElementType aType = aElem->GetType();
00444     int len = P.size();
00445     switch( aType ) {
00446     case SMDSAbs_Face:
00447       if( len == 3 ) { // triangles
00448         double L1 = getDistance(P( 1 ),P( 2 ));
00449         double L2 = getDistance(P( 2 ),P( 3 ));
00450         double L3 = getDistance(P( 3 ),P( 1 ));
00451         aVal = Max(L1,Max(L2,L3));
00452         break;
00453       }
00454       else if( len == 4 ) { // quadrangles
00455         double L1 = getDistance(P( 1 ),P( 2 ));
00456         double L2 = getDistance(P( 2 ),P( 3 ));
00457         double L3 = getDistance(P( 3 ),P( 4 ));
00458         double L4 = getDistance(P( 4 ),P( 1 ));
00459         double D1 = getDistance(P( 1 ),P( 3 ));
00460         double D2 = getDistance(P( 2 ),P( 4 ));
00461         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
00462         break;
00463       }
00464       else if( len == 6 ) { // quadratic triangles
00465         double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
00466         double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
00467         double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
00468         aVal = Max(L1,Max(L2,L3));
00469         break;
00470       }
00471       else if( len == 8 || len == 9 ) { // quadratic quadrangles
00472         double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
00473         double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
00474         double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
00475         double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
00476         double D1 = getDistance(P( 1 ),P( 5 ));
00477         double D2 = getDistance(P( 3 ),P( 7 ));
00478         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
00479         break;
00480       }
00481     }
00482 
00483     if( myPrecision >= 0 )
00484     {
00485       double prec = pow( 10., (double)myPrecision );
00486       aVal = floor( aVal * prec + 0.5 ) / prec;
00487     }
00488     return aVal;
00489   }
00490   return 0.;
00491 }
00492 
00493 double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const
00494 {
00495   return Value;
00496 }
00497 
00498 SMDSAbs_ElementType MaxElementLength2D::GetType() const
00499 {
00500   return SMDSAbs_Face;
00501 }
00502 
00503 /*
00504   Class       : MaxElementLength3D
00505   Description : Functor calculating maximum length of 3D element
00506 */
00507 
00508 double MaxElementLength3D::GetValue( long theElementId )
00509 {
00510   TSequenceOfXYZ P;
00511   if( GetPoints( theElementId, P ) ) {
00512     double aVal = 0;
00513     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
00514     SMDSAbs_ElementType aType = aElem->GetType();
00515     int len = P.size();
00516     switch( aType ) {
00517     case SMDSAbs_Volume:
00518       if( len == 4 ) { // tetras
00519         double L1 = getDistance(P( 1 ),P( 2 ));
00520         double L2 = getDistance(P( 2 ),P( 3 ));
00521         double L3 = getDistance(P( 3 ),P( 1 ));
00522         double L4 = getDistance(P( 1 ),P( 4 ));
00523         double L5 = getDistance(P( 2 ),P( 4 ));
00524         double L6 = getDistance(P( 3 ),P( 4 ));
00525         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
00526         break;
00527       }
00528       else if( len == 5 ) { // pyramids
00529         double L1 = getDistance(P( 1 ),P( 2 ));
00530         double L2 = getDistance(P( 2 ),P( 3 ));
00531         double L3 = getDistance(P( 3 ),P( 4 ));
00532         double L4 = getDistance(P( 4 ),P( 1 ));
00533         double L5 = getDistance(P( 1 ),P( 5 ));
00534         double L6 = getDistance(P( 2 ),P( 5 ));
00535         double L7 = getDistance(P( 3 ),P( 5 ));
00536         double L8 = getDistance(P( 4 ),P( 5 ));
00537         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
00538         aVal = Max(aVal,Max(L7,L8));
00539         break;
00540       }
00541       else if( len == 6 ) { // pentas
00542         double L1 = getDistance(P( 1 ),P( 2 ));
00543         double L2 = getDistance(P( 2 ),P( 3 ));
00544         double L3 = getDistance(P( 3 ),P( 1 ));
00545         double L4 = getDistance(P( 4 ),P( 5 ));
00546         double L5 = getDistance(P( 5 ),P( 6 ));
00547         double L6 = getDistance(P( 6 ),P( 4 ));
00548         double L7 = getDistance(P( 1 ),P( 4 ));
00549         double L8 = getDistance(P( 2 ),P( 5 ));
00550         double L9 = getDistance(P( 3 ),P( 6 ));
00551         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
00552         aVal = Max(aVal,Max(Max(L7,L8),L9));
00553         break;
00554       }
00555       else if( len == 8 ) { // hexas
00556         double L1 = getDistance(P( 1 ),P( 2 ));
00557         double L2 = getDistance(P( 2 ),P( 3 ));
00558         double L3 = getDistance(P( 3 ),P( 4 ));
00559         double L4 = getDistance(P( 4 ),P( 1 ));
00560         double L5 = getDistance(P( 5 ),P( 6 ));
00561         double L6 = getDistance(P( 6 ),P( 7 ));
00562         double L7 = getDistance(P( 7 ),P( 8 ));
00563         double L8 = getDistance(P( 8 ),P( 5 ));
00564         double L9 = getDistance(P( 1 ),P( 5 ));
00565         double L10= getDistance(P( 2 ),P( 6 ));
00566         double L11= getDistance(P( 3 ),P( 7 ));
00567         double L12= getDistance(P( 4 ),P( 8 ));
00568         double D1 = getDistance(P( 1 ),P( 7 ));
00569         double D2 = getDistance(P( 2 ),P( 8 ));
00570         double D3 = getDistance(P( 3 ),P( 5 ));
00571         double D4 = getDistance(P( 4 ),P( 6 ));
00572         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
00573         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
00574         aVal = Max(aVal,Max(L11,L12));
00575         aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
00576         break;
00577       }
00578       else if( len == 12 ) { // hexagonal prism
00579         for ( int i1 = 1; i1 < 12; ++i1 )
00580           for ( int i2 = i1+1; i1 <= 12; ++i1 )
00581             aVal = Max( aVal, getDistance(P( i1 ),P( i2 )));
00582         break;
00583       }
00584       else if( len == 10 ) { // quadratic tetras
00585         double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
00586         double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
00587         double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
00588         double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
00589         double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
00590         double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
00591         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
00592         break;
00593       }
00594       else if( len == 13 ) { // quadratic pyramids
00595         double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
00596         double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
00597         double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
00598         double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
00599         double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
00600         double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
00601         double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
00602         double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
00603         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
00604         aVal = Max(aVal,Max(L7,L8));
00605         break;
00606       }
00607       else if( len == 15 ) { // quadratic pentas
00608         double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
00609         double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
00610         double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
00611         double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
00612         double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
00613         double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
00614         double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
00615         double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
00616         double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
00617         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
00618         aVal = Max(aVal,Max(Max(L7,L8),L9));
00619         break;
00620       }
00621       else if( len == 20 || len == 27 ) { // quadratic hexas
00622         double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
00623         double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
00624         double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
00625         double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
00626         double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
00627         double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
00628         double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
00629         double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
00630         double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
00631         double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
00632         double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
00633         double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
00634         double D1 = getDistance(P( 1 ),P( 7 ));
00635         double D2 = getDistance(P( 2 ),P( 8 ));
00636         double D3 = getDistance(P( 3 ),P( 5 ));
00637         double D4 = getDistance(P( 4 ),P( 6 ));
00638         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
00639         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
00640         aVal = Max(aVal,Max(L11,L12));
00641         aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
00642         break;
00643       }
00644       else if( len > 1 && aElem->IsPoly() ) { // polys
00645         // get the maximum distance between all pairs of nodes
00646         for( int i = 1; i <= len; i++ ) {
00647           for( int j = 1; j <= len; j++ ) {
00648             if( j > i ) { // optimization of the loop
00649               double D = getDistance( P(i), P(j) );
00650               aVal = Max( aVal, D );
00651             }
00652           }
00653         }
00654       }
00655     }
00656 
00657     if( myPrecision >= 0 )
00658     {
00659       double prec = pow( 10., (double)myPrecision );
00660       aVal = floor( aVal * prec + 0.5 ) / prec;
00661     }
00662     return aVal;
00663   }
00664   return 0.;
00665 }
00666 
00667 double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const
00668 {
00669   return Value;
00670 }
00671 
00672 SMDSAbs_ElementType MaxElementLength3D::GetType() const
00673 {
00674   return SMDSAbs_Volume;
00675 }
00676 
00677 
00678 /*
00679   Class       : MinimumAngle
00680   Description : Functor for calculation of minimum angle
00681 */
00682 
00683 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
00684 {
00685   double aMin;
00686 
00687   if (P.size() <3)
00688     return 0.;
00689 
00690   aMin = getAngle(P( P.size() ), P( 1 ), P( 2 ));
00691   aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 )));
00692 
00693   for (int i=2; i<P.size();i++){
00694       double A0 = getAngle( P( i-1 ), P( i ), P( i+1 ) );
00695     aMin = Min(aMin,A0);
00696   }
00697 
00698   return aMin * 180.0 / M_PI;
00699 }
00700 
00701 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
00702 {
00703   //const double aBestAngle = PI / nbNodes;
00704   const double aBestAngle = 180.0 - ( 360.0 / double(nbNodes) );
00705   return ( fabs( aBestAngle - Value ));
00706 }
00707 
00708 SMDSAbs_ElementType MinimumAngle::GetType() const
00709 {
00710   return SMDSAbs_Face;
00711 }
00712 
00713 
00714 /*
00715   Class       : AspectRatio
00716   Description : Functor for calculating aspect ratio
00717 */
00718 double AspectRatio::GetValue( const TSequenceOfXYZ& P )
00719 {
00720   // According to "Mesh quality control" by Nadir Bouhamau referring to
00721   // Pascal Jean Frey and Paul-Louis George. Maillages, applications aux elements finis.
00722   // Hermes Science publications, Paris 1999 ISBN 2-7462-0024-4
00723   // PAL10872
00724 
00725   int nbNodes = P.size();
00726 
00727   if ( nbNodes < 3 )
00728     return 0;
00729 
00730   // Compute aspect ratio
00731 
00732   if ( nbNodes == 3 ) {
00733     // Compute lengths of the sides
00734     std::vector< double > aLen (nbNodes);
00735     for ( int i = 0; i < nbNodes - 1; i++ )
00736       aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
00737     aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
00738     // Q = alfa * h * p / S, where
00739     //
00740     // alfa = sqrt( 3 ) / 6
00741     // h - length of the longest edge
00742     // p - half perimeter
00743     // S - triangle surface
00744     const double alfa = sqrt( 3. ) / 6.;
00745     double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
00746     double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
00747     double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
00748     if ( anArea <= Precision::Confusion() )
00749       return 0.;
00750     return alfa * maxLen * half_perimeter / anArea;
00751   }
00752   else if ( nbNodes == 6 ) { // quadratic triangles
00753     // Compute lengths of the sides
00754     std::vector< double > aLen (3);
00755     aLen[0] = getDistance( P(1), P(3) );
00756     aLen[1] = getDistance( P(3), P(5) );
00757     aLen[2] = getDistance( P(5), P(1) );
00758     // Q = alfa * h * p / S, where
00759     //
00760     // alfa = sqrt( 3 ) / 6
00761     // h - length of the longest edge
00762     // p - half perimeter
00763     // S - triangle surface
00764     const double alfa = sqrt( 3. ) / 6.;
00765     double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
00766     double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
00767     double anArea = getArea( P(1), P(3), P(5) );
00768     if ( anArea <= Precision::Confusion() )
00769       return 0.;
00770     return alfa * maxLen * half_perimeter / anArea;
00771   }
00772   else if( nbNodes == 4 ) { // quadrangle
00773     // Compute lengths of the sides
00774     std::vector< double > aLen (4);
00775     aLen[0] = getDistance( P(1), P(2) );
00776     aLen[1] = getDistance( P(2), P(3) );
00777     aLen[2] = getDistance( P(3), P(4) );
00778     aLen[3] = getDistance( P(4), P(1) );
00779     // Compute lengths of the diagonals
00780     std::vector< double > aDia (2);
00781     aDia[0] = getDistance( P(1), P(3) );
00782     aDia[1] = getDistance( P(2), P(4) );
00783     // Compute areas of all triangles which can be built
00784     // taking three nodes of the quadrangle
00785     std::vector< double > anArea (4);
00786     anArea[0] = getArea( P(1), P(2), P(3) );
00787     anArea[1] = getArea( P(1), P(2), P(4) );
00788     anArea[2] = getArea( P(1), P(3), P(4) );
00789     anArea[3] = getArea( P(2), P(3), P(4) );
00790     // Q = alpha * L * C1 / C2, where
00791     //
00792     // alpha = sqrt( 1/32 )
00793     // L = max( L1, L2, L3, L4, D1, D2 )
00794     // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
00795     // C2 = min( S1, S2, S3, S4 )
00796     // Li - lengths of the edges
00797     // Di - lengths of the diagonals
00798     // Si - areas of the triangles
00799     const double alpha = sqrt( 1 / 32. );
00800     double L = Max( aLen[ 0 ],
00801                  Max( aLen[ 1 ],
00802                    Max( aLen[ 2 ],
00803                      Max( aLen[ 3 ],
00804                        Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
00805     double C1 = sqrt( ( aLen[0] * aLen[0] +
00806                         aLen[1] * aLen[1] +
00807                         aLen[2] * aLen[2] +
00808                         aLen[3] * aLen[3] ) / 4. );
00809     double C2 = Min( anArea[ 0 ],
00810                   Min( anArea[ 1 ],
00811                     Min( anArea[ 2 ], anArea[ 3 ] ) ) );
00812     if ( C2 <= Precision::Confusion() )
00813       return 0.;
00814     return alpha * L * C1 / C2;
00815   }
00816   else if( nbNodes == 8 || nbNodes == 9 ) { // nbNodes==8 - quadratic quadrangle
00817     // Compute lengths of the sides
00818     std::vector< double > aLen (4);
00819     aLen[0] = getDistance( P(1), P(3) );
00820     aLen[1] = getDistance( P(3), P(5) );
00821     aLen[2] = getDistance( P(5), P(7) );
00822     aLen[3] = getDistance( P(7), P(1) );
00823     // Compute lengths of the diagonals
00824     std::vector< double > aDia (2);
00825     aDia[0] = getDistance( P(1), P(5) );
00826     aDia[1] = getDistance( P(3), P(7) );
00827     // Compute areas of all triangles which can be built
00828     // taking three nodes of the quadrangle
00829     std::vector< double > anArea (4);
00830     anArea[0] = getArea( P(1), P(3), P(5) );
00831     anArea[1] = getArea( P(1), P(3), P(7) );
00832     anArea[2] = getArea( P(1), P(5), P(7) );
00833     anArea[3] = getArea( P(3), P(5), P(7) );
00834     // Q = alpha * L * C1 / C2, where
00835     //
00836     // alpha = sqrt( 1/32 )
00837     // L = max( L1, L2, L3, L4, D1, D2 )
00838     // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
00839     // C2 = min( S1, S2, S3, S4 )
00840     // Li - lengths of the edges
00841     // Di - lengths of the diagonals
00842     // Si - areas of the triangles
00843     const double alpha = sqrt( 1 / 32. );
00844     double L = Max( aLen[ 0 ],
00845                  Max( aLen[ 1 ],
00846                    Max( aLen[ 2 ],
00847                      Max( aLen[ 3 ],
00848                        Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
00849     double C1 = sqrt( ( aLen[0] * aLen[0] +
00850                         aLen[1] * aLen[1] +
00851                         aLen[2] * aLen[2] +
00852                         aLen[3] * aLen[3] ) / 4. );
00853     double C2 = Min( anArea[ 0 ],
00854                   Min( anArea[ 1 ],
00855                     Min( anArea[ 2 ], anArea[ 3 ] ) ) );
00856     if ( C2 <= Precision::Confusion() )
00857       return 0.;
00858     return alpha * L * C1 / C2;
00859   }
00860   return 0;
00861 }
00862 
00863 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
00864 {
00865   // the aspect ratio is in the range [1.0,infinity]
00866   // < 1.0 = very bad, zero area
00867   // 1.0 = good
00868   // infinity = bad
00869   return ( Value < 0.9 ) ? 1000 : Value / 1000.;
00870 }
00871 
00872 SMDSAbs_ElementType AspectRatio::GetType() const
00873 {
00874   return SMDSAbs_Face;
00875 }
00876 
00877 
00878 /*
00879   Class       : AspectRatio3D
00880   Description : Functor for calculating aspect ratio
00881 */
00882 namespace{
00883 
00884   inline double getHalfPerimeter(double theTria[3]){
00885     return (theTria[0] + theTria[1] + theTria[2])/2.0;
00886   }
00887 
00888   inline double getArea(double theHalfPerim, double theTria[3]){
00889     return sqrt(theHalfPerim*
00890                 (theHalfPerim-theTria[0])*
00891                 (theHalfPerim-theTria[1])*
00892                 (theHalfPerim-theTria[2]));
00893   }
00894 
00895   inline double getVolume(double theLen[6]){
00896     double a2 = theLen[0]*theLen[0];
00897     double b2 = theLen[1]*theLen[1];
00898     double c2 = theLen[2]*theLen[2];
00899     double d2 = theLen[3]*theLen[3];
00900     double e2 = theLen[4]*theLen[4];
00901     double f2 = theLen[5]*theLen[5];
00902     double P = 4.0*a2*b2*d2;
00903     double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
00904     double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
00905     return sqrt(P-Q+R)/12.0;
00906   }
00907 
00908   inline double getVolume2(double theLen[6]){
00909     double a2 = theLen[0]*theLen[0];
00910     double b2 = theLen[1]*theLen[1];
00911     double c2 = theLen[2]*theLen[2];
00912     double d2 = theLen[3]*theLen[3];
00913     double e2 = theLen[4]*theLen[4];
00914     double f2 = theLen[5]*theLen[5];
00915 
00916     double P = a2*e2*(b2+c2+d2+f2-a2-e2);
00917     double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
00918     double R = c2*d2*(a2+b2+e2+f2-c2-d2);
00919     double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
00920 
00921     return sqrt(P+Q+R-S)/12.0;
00922   }
00923 
00924   inline double getVolume(const TSequenceOfXYZ& P){
00925     gp_Vec aVec1( P( 2 ) - P( 1 ) );
00926     gp_Vec aVec2( P( 3 ) - P( 1 ) );
00927     gp_Vec aVec3( P( 4 ) - P( 1 ) );
00928     gp_Vec anAreaVec( aVec1 ^ aVec2 );
00929     return fabs(aVec3 * anAreaVec) / 6.0;
00930   }
00931 
00932   inline double getMaxHeight(double theLen[6])
00933   {
00934     double aHeight = std::max(theLen[0],theLen[1]);
00935     aHeight = std::max(aHeight,theLen[2]);
00936     aHeight = std::max(aHeight,theLen[3]);
00937     aHeight = std::max(aHeight,theLen[4]);
00938     aHeight = std::max(aHeight,theLen[5]);
00939     return aHeight;
00940   }
00941 
00942 }
00943 
00944 double AspectRatio3D::GetValue( long theId )
00945 {
00946   double aVal = 0;
00947   myCurrElement = myMesh->FindElement( theId );
00948   if ( myCurrElement && myCurrElement->GetVtkType() == VTK_TETRA )
00949   {
00950     // Action from CoTech | ACTION 31.3:
00951     // EURIWARE BO: Homogenize the formulas used to calculate the Controls in SMESH to fit with
00952     // those of ParaView. The library used by ParaView for those calculations can be reused in SMESH.
00953     vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid();
00954     if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() ))
00955       aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell ));
00956   }
00957   else
00958   {
00959     TSequenceOfXYZ P;
00960     if ( GetPoints( myCurrElement, P ))
00961       aVal = Round( GetValue( P ));
00962   }
00963   return aVal;
00964 }
00965 
00966 double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
00967 {
00968   double aQuality = 0.0;
00969   if(myCurrElement->IsPoly()) return aQuality;
00970 
00971   int nbNodes = P.size();
00972 
00973   if(myCurrElement->IsQuadratic()) {
00974     if(nbNodes==10) nbNodes=4; // quadratic tetrahedron
00975     else if(nbNodes==13) nbNodes=5; // quadratic pyramid
00976     else if(nbNodes==15) nbNodes=6; // quadratic pentahedron
00977     else if(nbNodes==20) nbNodes=8; // quadratic hexahedron
00978     else if(nbNodes==27) nbNodes=8; // quadratic hexahedron
00979     else return aQuality;
00980   }
00981 
00982   switch(nbNodes) {
00983   case 4:{
00984     double aLen[6] = {
00985       getDistance(P( 1 ),P( 2 )), // a
00986       getDistance(P( 2 ),P( 3 )), // b
00987       getDistance(P( 3 ),P( 1 )), // c
00988       getDistance(P( 2 ),P( 4 )), // d
00989       getDistance(P( 3 ),P( 4 )), // e
00990       getDistance(P( 1 ),P( 4 ))  // f
00991     };
00992     double aTria[4][3] = {
00993       {aLen[0],aLen[1],aLen[2]}, // abc
00994       {aLen[0],aLen[3],aLen[5]}, // adf
00995       {aLen[1],aLen[3],aLen[4]}, // bde
00996       {aLen[2],aLen[4],aLen[5]}  // cef
00997     };
00998     double aSumArea = 0.0;
00999     double aHalfPerimeter = getHalfPerimeter(aTria[0]);
01000     double anArea = getArea(aHalfPerimeter,aTria[0]);
01001     aSumArea += anArea;
01002     aHalfPerimeter = getHalfPerimeter(aTria[1]);
01003     anArea = getArea(aHalfPerimeter,aTria[1]);
01004     aSumArea += anArea;
01005     aHalfPerimeter = getHalfPerimeter(aTria[2]);
01006     anArea = getArea(aHalfPerimeter,aTria[2]);
01007     aSumArea += anArea;
01008     aHalfPerimeter = getHalfPerimeter(aTria[3]);
01009     anArea = getArea(aHalfPerimeter,aTria[3]);
01010     aSumArea += anArea;
01011     double aVolume = getVolume(P);
01012     //double aVolume = getVolume(aLen);
01013     double aHeight = getMaxHeight(aLen);
01014     static double aCoeff = sqrt(2.0)/12.0;
01015     if ( aVolume > DBL_MIN )
01016       aQuality = aCoeff*aHeight*aSumArea/aVolume;
01017     break;
01018   }
01019   case 5:{
01020     {
01021       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
01022       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01023     }
01024     {
01025       gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
01026       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01027     }
01028     {
01029       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
01030       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01031     }
01032     {
01033       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
01034       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01035     }
01036     break;
01037   }
01038   case 6:{
01039     {
01040       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
01041       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01042     }
01043     {
01044       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
01045       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01046     }
01047     {
01048       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
01049       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01050     }
01051     {
01052       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
01053       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01054     }
01055     {
01056       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
01057       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01058     }
01059     {
01060       gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
01061       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01062     }
01063     break;
01064   }
01065   case 8:{
01066     {
01067       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
01068       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01069     }
01070     {
01071       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
01072       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01073     }
01074     {
01075       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
01076       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01077     }
01078     {
01079       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
01080       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01081     }
01082     {
01083       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
01084       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01085     }
01086     {
01087       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
01088       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01089     }
01090     {
01091       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
01092       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01093     }
01094     {
01095       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
01096       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01097     }
01098     {
01099       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
01100       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01101     }
01102     {
01103       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
01104       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01105     }
01106     {
01107       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
01108       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01109     }
01110     {
01111       gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
01112       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01113     }
01114     {
01115       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
01116       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01117     }
01118     {
01119       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
01120       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01121     }
01122     {
01123       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
01124       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01125     }
01126     {
01127       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
01128       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01129     }
01130     {
01131       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
01132       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01133     }
01134     {
01135       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
01136       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01137     }
01138     {
01139       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
01140       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01141     }
01142     {
01143       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
01144       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01145     }
01146     {
01147       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
01148       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01149     }
01150     {
01151       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
01152       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01153     }
01154     {
01155       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
01156       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01157     }
01158     {
01159       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
01160       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01161     }
01162     {
01163       gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
01164       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01165     }
01166     {
01167       gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
01168       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01169     }
01170     {
01171       gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
01172       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01173     }
01174     {
01175       gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
01176       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01177     }
01178     {
01179       gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
01180       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01181     }
01182     {
01183       gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
01184       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01185     }
01186     {
01187       gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
01188       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01189     }
01190     {
01191       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
01192       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01193     }
01194     {
01195       gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
01196       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
01197     }
01198     break;
01199   }
01200   case 12:
01201     {
01202       gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )};
01203       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
01204     }
01205     {
01206       gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )};
01207       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
01208     }
01209     {
01210       gp_XYZ aXYZ[8] = {P( 3 ),P( 4 ),P( 6 ),P( 1 ),P( 9 ),P( 10 ),P( 12 ),P( 7 )};
01211       aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
01212     }
01213     break;
01214   } // switch(nbNodes)
01215 
01216   if ( nbNodes > 4 ) {
01217     // avaluate aspect ratio of quadranle faces
01218     AspectRatio aspect2D;
01219     SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes );
01220     int nbFaces = SMDS_VolumeTool::NbFaces( type );
01221     TSequenceOfXYZ points(4);
01222     for ( int i = 0; i < nbFaces; ++i ) { // loop on faces of a volume
01223       if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 )
01224         continue;
01225       const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true );
01226       for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face
01227         points( p + 1 ) = P( pInd[ p ] + 1 );
01228       aQuality = std::max( aQuality, aspect2D.GetValue( points ));
01229     }
01230   }
01231   return aQuality;
01232 }
01233 
01234 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
01235 {
01236   // the aspect ratio is in the range [1.0,infinity]
01237   // 1.0 = good
01238   // infinity = bad
01239   return Value / 1000.;
01240 }
01241 
01242 SMDSAbs_ElementType AspectRatio3D::GetType() const
01243 {
01244   return SMDSAbs_Volume;
01245 }
01246 
01247 
01248 /*
01249   Class       : Warping
01250   Description : Functor for calculating warping
01251 */
01252 double Warping::GetValue( const TSequenceOfXYZ& P )
01253 {
01254   if ( P.size() != 4 )
01255     return 0;
01256 
01257   gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4.;
01258 
01259   double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
01260   double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
01261   double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
01262   double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
01263 
01264   return Max( Max( A1, A2 ), Max( A3, A4 ) );
01265 }
01266 
01267 double Warping::ComputeA( const gp_XYZ& thePnt1,
01268                           const gp_XYZ& thePnt2,
01269                           const gp_XYZ& thePnt3,
01270                           const gp_XYZ& theG ) const
01271 {
01272   double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
01273   double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
01274   double L = Min( aLen1, aLen2 ) * 0.5;
01275   if ( L < Precision::Confusion())
01276     return 0.;
01277 
01278   gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG;
01279   gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG;
01280   gp_XYZ N  = GI.Crossed( GJ );
01281 
01282   if ( N.Modulus() < gp::Resolution() )
01283     return M_PI / 2;
01284 
01285   N.Normalize();
01286 
01287   double H = ( thePnt2 - theG ).Dot( N );
01288   return asin( fabs( H / L ) ) * 180. / M_PI;
01289 }
01290 
01291 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
01292 {
01293   // the warp is in the range [0.0,PI/2]
01294   // 0.0 = good (no warp)
01295   // PI/2 = bad  (face pliee)
01296   return Value;
01297 }
01298 
01299 SMDSAbs_ElementType Warping::GetType() const
01300 {
01301   return SMDSAbs_Face;
01302 }
01303 
01304 
01305 /*
01306   Class       : Taper
01307   Description : Functor for calculating taper
01308 */
01309 double Taper::GetValue( const TSequenceOfXYZ& P )
01310 {
01311   if ( P.size() != 4 )
01312     return 0.;
01313 
01314   // Compute taper
01315   double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2.;
01316   double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2.;
01317   double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2.;
01318   double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2.;
01319 
01320   double JA = 0.25 * ( J1 + J2 + J3 + J4 );
01321   if ( JA <= Precision::Confusion() )
01322     return 0.;
01323 
01324   double T1 = fabs( ( J1 - JA ) / JA );
01325   double T2 = fabs( ( J2 - JA ) / JA );
01326   double T3 = fabs( ( J3 - JA ) / JA );
01327   double T4 = fabs( ( J4 - JA ) / JA );
01328 
01329   return Max( Max( T1, T2 ), Max( T3, T4 ) );
01330 }
01331 
01332 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
01333 {
01334   // the taper is in the range [0.0,1.0]
01335   // 0.0  = good (no taper)
01336   // 1.0 = bad  (les cotes opposes sont allignes)
01337   return Value;
01338 }
01339 
01340 SMDSAbs_ElementType Taper::GetType() const
01341 {
01342   return SMDSAbs_Face;
01343 }
01344 
01345 
01346 /*
01347   Class       : Skew
01348   Description : Functor for calculating skew in degrees
01349 */
01350 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
01351 {
01352   gp_XYZ p12 = ( p2 + p1 ) / 2.;
01353   gp_XYZ p23 = ( p3 + p2 ) / 2.;
01354   gp_XYZ p31 = ( p3 + p1 ) / 2.;
01355 
01356   gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
01357 
01358   return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 );
01359 }
01360 
01361 double Skew::GetValue( const TSequenceOfXYZ& P )
01362 {
01363   if ( P.size() != 3 && P.size() != 4 )
01364     return 0.;
01365 
01366   // Compute skew
01367   static double PI2 = M_PI / 2.;
01368   if ( P.size() == 3 )
01369   {
01370     double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
01371     double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
01372     double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
01373 
01374     return Max( A0, Max( A1, A2 ) ) * 180. / M_PI;
01375   }
01376   else
01377   {
01378     gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2.;
01379     gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2.;
01380     gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2.;
01381     gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2.;
01382 
01383     gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
01384     double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
01385       ? 0. : fabs( PI2 - v1.Angle( v2 ) );
01386 
01387     //BUG SWP12743
01388     if ( A < Precision::Angular() )
01389       return 0.;
01390 
01391     return A * 180. / M_PI;
01392   }
01393 }
01394 
01395 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
01396 {
01397   // the skew is in the range [0.0,PI/2].
01398   // 0.0 = good
01399   // PI/2 = bad
01400   return Value;
01401 }
01402 
01403 SMDSAbs_ElementType Skew::GetType() const
01404 {
01405   return SMDSAbs_Face;
01406 }
01407 
01408 
01409 /*
01410   Class       : Area
01411   Description : Functor for calculating area
01412 */
01413 double Area::GetValue( const TSequenceOfXYZ& P )
01414 {
01415   double val = 0.0;
01416   if ( P.size() > 2 ) {
01417     gp_Vec aVec1( P(2) - P(1) );
01418     gp_Vec aVec2( P(3) - P(1) );
01419     gp_Vec SumVec = aVec1 ^ aVec2;
01420     for (int i=4; i<=P.size(); i++) {
01421       gp_Vec aVec1( P(i-1) - P(1) );
01422       gp_Vec aVec2( P(i) - P(1) );
01423       gp_Vec tmp = aVec1 ^ aVec2;
01424       SumVec.Add(tmp);
01425     }
01426     val = SumVec.Magnitude() * 0.5;
01427   }
01428   return val;
01429 }
01430 
01431 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
01432 {
01433   // meaningless as it is not a quality control functor
01434   return Value;
01435 }
01436 
01437 SMDSAbs_ElementType Area::GetType() const
01438 {
01439   return SMDSAbs_Face;
01440 }
01441 
01442 
01443 /*
01444   Class       : Length
01445   Description : Functor for calculating length of edge
01446 */
01447 double Length::GetValue( const TSequenceOfXYZ& P )
01448 {
01449   switch ( P.size() ) {
01450   case 2:  return getDistance( P( 1 ), P( 2 ) );
01451   case 3:  return getDistance( P( 1 ), P( 2 ) ) + getDistance( P( 2 ), P( 3 ) );
01452   default: return 0.;
01453   }
01454 }
01455 
01456 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
01457 {
01458   // meaningless as it is not quality control functor
01459   return Value;
01460 }
01461 
01462 SMDSAbs_ElementType Length::GetType() const
01463 {
01464   return SMDSAbs_Edge;
01465 }
01466 
01467 /*
01468   Class       : Length2D
01469   Description : Functor for calculating length of edge
01470 */
01471 
01472 double Length2D::GetValue( long theElementId)
01473 {
01474   TSequenceOfXYZ P;
01475 
01476   //cout<<"Length2D::GetValue"<<endl;
01477   if (GetPoints(theElementId,P)){
01478     //for(int jj=1; jj<=P.size(); jj++)
01479     //  cout<<"jj="<<jj<<" P("<<P(jj).X()<<","<<P(jj).Y()<<","<<P(jj).Z()<<")"<<endl;
01480 
01481     double  aVal;// = GetValue( P );
01482     const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
01483     SMDSAbs_ElementType aType = aElem->GetType();
01484 
01485     int len = P.size();
01486 
01487     switch (aType){
01488     case SMDSAbs_All:
01489     case SMDSAbs_Node:
01490     case SMDSAbs_Edge:
01491       if (len == 2){
01492         aVal = getDistance( P( 1 ), P( 2 ) );
01493         break;
01494       }
01495       else if (len == 3){ // quadratic edge
01496         aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
01497         break;
01498       }
01499     case SMDSAbs_Face:
01500       if (len == 3){ // triangles
01501         double L1 = getDistance(P( 1 ),P( 2 ));
01502         double L2 = getDistance(P( 2 ),P( 3 ));
01503         double L3 = getDistance(P( 3 ),P( 1 ));
01504         aVal = Max(L1,Max(L2,L3));
01505         break;
01506       }
01507       else if (len == 4){ // quadrangles
01508         double L1 = getDistance(P( 1 ),P( 2 ));
01509         double L2 = getDistance(P( 2 ),P( 3 ));
01510         double L3 = getDistance(P( 3 ),P( 4 ));
01511         double L4 = getDistance(P( 4 ),P( 1 ));
01512         aVal = Max(Max(L1,L2),Max(L3,L4));
01513         break;
01514       }
01515       if (len == 6){ // quadratic triangles
01516         double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
01517         double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
01518         double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
01519         aVal = Max(L1,Max(L2,L3));
01520         //cout<<"L1="<<L1<<" L2="<<L2<<"L3="<<L3<<" aVal="<<aVal<<endl;
01521         break;
01522       }
01523       else if (len == 8){ // quadratic quadrangles
01524         double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
01525         double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
01526         double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
01527         double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
01528         aVal = Max(Max(L1,L2),Max(L3,L4));
01529         break;
01530       }
01531     case SMDSAbs_Volume:
01532       if (len == 4){ // tetraidrs
01533         double L1 = getDistance(P( 1 ),P( 2 ));
01534         double L2 = getDistance(P( 2 ),P( 3 ));
01535         double L3 = getDistance(P( 3 ),P( 1 ));
01536         double L4 = getDistance(P( 1 ),P( 4 ));
01537         double L5 = getDistance(P( 2 ),P( 4 ));
01538         double L6 = getDistance(P( 3 ),P( 4 ));
01539         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
01540         break;
01541       }
01542       else if (len == 5){ // piramids
01543         double L1 = getDistance(P( 1 ),P( 2 ));
01544         double L2 = getDistance(P( 2 ),P( 3 ));
01545         double L3 = getDistance(P( 3 ),P( 4 ));
01546         double L4 = getDistance(P( 4 ),P( 1 ));
01547         double L5 = getDistance(P( 1 ),P( 5 ));
01548         double L6 = getDistance(P( 2 ),P( 5 ));
01549         double L7 = getDistance(P( 3 ),P( 5 ));
01550         double L8 = getDistance(P( 4 ),P( 5 ));
01551 
01552         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
01553         aVal = Max(aVal,Max(L7,L8));
01554         break;
01555       }
01556       else if (len == 6){ // pentaidres
01557         double L1 = getDistance(P( 1 ),P( 2 ));
01558         double L2 = getDistance(P( 2 ),P( 3 ));
01559         double L3 = getDistance(P( 3 ),P( 1 ));
01560         double L4 = getDistance(P( 4 ),P( 5 ));
01561         double L5 = getDistance(P( 5 ),P( 6 ));
01562         double L6 = getDistance(P( 6 ),P( 4 ));
01563         double L7 = getDistance(P( 1 ),P( 4 ));
01564         double L8 = getDistance(P( 2 ),P( 5 ));
01565         double L9 = getDistance(P( 3 ),P( 6 ));
01566 
01567         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
01568         aVal = Max(aVal,Max(Max(L7,L8),L9));
01569         break;
01570       }
01571       else if (len == 8){ // hexaider
01572         double L1 = getDistance(P( 1 ),P( 2 ));
01573         double L2 = getDistance(P( 2 ),P( 3 ));
01574         double L3 = getDistance(P( 3 ),P( 4 ));
01575         double L4 = getDistance(P( 4 ),P( 1 ));
01576         double L5 = getDistance(P( 5 ),P( 6 ));
01577         double L6 = getDistance(P( 6 ),P( 7 ));
01578         double L7 = getDistance(P( 7 ),P( 8 ));
01579         double L8 = getDistance(P( 8 ),P( 5 ));
01580         double L9 = getDistance(P( 1 ),P( 5 ));
01581         double L10= getDistance(P( 2 ),P( 6 ));
01582         double L11= getDistance(P( 3 ),P( 7 ));
01583         double L12= getDistance(P( 4 ),P( 8 ));
01584 
01585         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
01586         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
01587         aVal = Max(aVal,Max(L11,L12));
01588         break;
01589 
01590       }
01591 
01592       if (len == 10){ // quadratic tetraidrs
01593         double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
01594         double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
01595         double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
01596         double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
01597         double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
01598         double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
01599         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
01600         break;
01601       }
01602       else if (len == 13){ // quadratic piramids
01603         double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
01604         double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
01605         double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
01606         double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
01607         double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
01608         double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
01609         double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
01610         double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
01611         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
01612         aVal = Max(aVal,Max(L7,L8));
01613         break;
01614       }
01615       else if (len == 15){ // quadratic pentaidres
01616         double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
01617         double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
01618         double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
01619         double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
01620         double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
01621         double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
01622         double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
01623         double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
01624         double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
01625         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
01626         aVal = Max(aVal,Max(Max(L7,L8),L9));
01627         break;
01628       }
01629       else if (len == 20){ // quadratic hexaider
01630         double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
01631         double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
01632         double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
01633         double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
01634         double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
01635         double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
01636         double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
01637         double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
01638         double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
01639         double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
01640         double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
01641         double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
01642         aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
01643         aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
01644         aVal = Max(aVal,Max(L11,L12));
01645         break;
01646 
01647       }
01648 
01649     default: aVal=-1;
01650     }
01651 
01652     if (aVal <0){
01653       return 0.;
01654     }
01655 
01656     if ( myPrecision >= 0 )
01657     {
01658       double prec = pow( 10., (double)( myPrecision ) );
01659       aVal = floor( aVal * prec + 0.5 ) / prec;
01660     }
01661 
01662     return aVal;
01663 
01664   }
01665   return 0.;
01666 }
01667 
01668 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
01669 {
01670   // meaningless as it is not quality control functor
01671   return Value;
01672 }
01673 
01674 SMDSAbs_ElementType Length2D::GetType() const
01675 {
01676   return SMDSAbs_Face;
01677 }
01678 
01679 Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
01680   myLength(theLength)
01681 {
01682   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
01683   if(thePntId1 > thePntId2){
01684     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
01685   }
01686 }
01687 
01688 bool Length2D::Value::operator<(const Length2D::Value& x) const{
01689   if(myPntId[0] < x.myPntId[0]) return true;
01690   if(myPntId[0] == x.myPntId[0])
01691     if(myPntId[1] < x.myPntId[1]) return true;
01692   return false;
01693 }
01694 
01695 void Length2D::GetValues(TValues& theValues){
01696   TValues aValues;
01697   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
01698   for(; anIter->more(); ){
01699     const SMDS_MeshFace* anElem = anIter->next();
01700 
01701     if(anElem->IsQuadratic()) {
01702       const SMDS_VtkFace* F =
01703         dynamic_cast<const SMDS_VtkFace*>(anElem);
01704       // use special nodes iterator
01705       SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
01706       long aNodeId[4];
01707       gp_Pnt P[4];
01708 
01709       double aLength;
01710       const SMDS_MeshElement* aNode;
01711       if(anIter->more()){
01712         aNode = anIter->next();
01713         const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
01714         P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
01715         aNodeId[0] = aNodeId[1] = aNode->GetID();
01716         aLength = 0;
01717       }
01718       for(; anIter->more(); ){
01719         const SMDS_MeshNode* N1 = static_cast<const SMDS_MeshNode*> (anIter->next());
01720         P[2] = gp_Pnt(N1->X(),N1->Y(),N1->Z());
01721         aNodeId[2] = N1->GetID();
01722         aLength = P[1].Distance(P[2]);
01723         if(!anIter->more()) break;
01724         const SMDS_MeshNode* N2 = static_cast<const SMDS_MeshNode*> (anIter->next());
01725         P[3] = gp_Pnt(N2->X(),N2->Y(),N2->Z());
01726         aNodeId[3] = N2->GetID();
01727         aLength += P[2].Distance(P[3]);
01728         Value aValue1(aLength,aNodeId[1],aNodeId[2]);
01729         Value aValue2(aLength,aNodeId[2],aNodeId[3]);
01730         P[1] = P[3];
01731         aNodeId[1] = aNodeId[3];
01732         theValues.insert(aValue1);
01733         theValues.insert(aValue2);
01734       }
01735       aLength += P[2].Distance(P[0]);
01736       Value aValue1(aLength,aNodeId[1],aNodeId[2]);
01737       Value aValue2(aLength,aNodeId[2],aNodeId[0]);
01738       theValues.insert(aValue1);
01739       theValues.insert(aValue2);
01740     }
01741     else {
01742       SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
01743       long aNodeId[2];
01744       gp_Pnt P[3];
01745 
01746       double aLength;
01747       const SMDS_MeshElement* aNode;
01748       if(aNodesIter->more()){
01749         aNode = aNodesIter->next();
01750         const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
01751         P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
01752         aNodeId[0] = aNodeId[1] = aNode->GetID();
01753         aLength = 0;
01754       }
01755       for(; aNodesIter->more(); ){
01756         aNode = aNodesIter->next();
01757         const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
01758         long anId = aNode->GetID();
01759         
01760         P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
01761         
01762         aLength = P[1].Distance(P[2]);
01763         
01764         Value aValue(aLength,aNodeId[1],anId);
01765         aNodeId[1] = anId;
01766         P[1] = P[2];
01767         theValues.insert(aValue);
01768       }
01769 
01770       aLength = P[0].Distance(P[1]);
01771 
01772       Value aValue(aLength,aNodeId[0],aNodeId[1]);
01773       theValues.insert(aValue);
01774     }
01775   }
01776 }
01777 
01778 /*
01779   Class       : MultiConnection
01780   Description : Functor for calculating number of faces conneted to the edge
01781 */
01782 double MultiConnection::GetValue( const TSequenceOfXYZ& P )
01783 {
01784   return 0;
01785 }
01786 double MultiConnection::GetValue( long theId )
01787 {
01788   return getNbMultiConnection( myMesh, theId );
01789 }
01790 
01791 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
01792 {
01793   // meaningless as it is not quality control functor
01794   return Value;
01795 }
01796 
01797 SMDSAbs_ElementType MultiConnection::GetType() const
01798 {
01799   return SMDSAbs_Edge;
01800 }
01801 
01802 /*
01803   Class       : MultiConnection2D
01804   Description : Functor for calculating number of faces conneted to the edge
01805 */
01806 double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
01807 {
01808   return 0;
01809 }
01810 
01811 double MultiConnection2D::GetValue( long theElementId )
01812 {
01813   int aResult = 0;
01814 
01815   const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId);
01816   SMDSAbs_ElementType aType = aFaceElem->GetType();
01817 
01818   switch (aType) {
01819   case SMDSAbs_Face:
01820     {
01821       int i = 0, len = aFaceElem->NbNodes();
01822       SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
01823       if (!anIter) break;
01824 
01825       const SMDS_MeshNode *aNode, *aNode0;
01826       TColStd_MapOfInteger aMap, aMapPrev;
01827 
01828       for (i = 0; i <= len; i++) {
01829         aMapPrev = aMap;
01830         aMap.Clear();
01831 
01832         int aNb = 0;
01833         if (anIter->more()) {
01834           aNode = (SMDS_MeshNode*)anIter->next();
01835         } else {
01836           if (i == len)
01837             aNode = aNode0;
01838           else
01839             break;
01840         }
01841         if (!aNode) break;
01842         if (i == 0) aNode0 = aNode;
01843 
01844         SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
01845         while (anElemIter->more()) {
01846           const SMDS_MeshElement* anElem = anElemIter->next();
01847           if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
01848             int anId = anElem->GetID();
01849 
01850             aMap.Add(anId);
01851             if (aMapPrev.Contains(anId)) {
01852               aNb++;
01853             }
01854           }
01855         }
01856         aResult = Max(aResult, aNb);
01857       }
01858     }
01859     break;
01860   default:
01861     aResult = 0;
01862   }
01863 
01864   return aResult;
01865 }
01866 
01867 double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
01868 {
01869   // meaningless as it is not quality control functor
01870   return Value;
01871 }
01872 
01873 SMDSAbs_ElementType MultiConnection2D::GetType() const
01874 {
01875   return SMDSAbs_Face;
01876 }
01877 
01878 MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
01879 {
01880   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
01881   if(thePntId1 > thePntId2){
01882     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
01883   }
01884 }
01885 
01886 bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const{
01887   if(myPntId[0] < x.myPntId[0]) return true;
01888   if(myPntId[0] == x.myPntId[0])
01889     if(myPntId[1] < x.myPntId[1]) return true;
01890   return false;
01891 }
01892 
01893 void MultiConnection2D::GetValues(MValues& theValues){
01894   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
01895   for(; anIter->more(); ){
01896     const SMDS_MeshFace* anElem = anIter->next();
01897     SMDS_ElemIteratorPtr aNodesIter;
01898     if ( anElem->IsQuadratic() )
01899       aNodesIter = dynamic_cast<const SMDS_VtkFace*>
01900         (anElem)->interlacedNodesElemIterator();
01901     else
01902       aNodesIter = anElem->nodesIterator();
01903     long aNodeId[3];
01904 
01905     //int aNbConnects=0;
01906     const SMDS_MeshNode* aNode0;
01907     const SMDS_MeshNode* aNode1;
01908     const SMDS_MeshNode* aNode2;
01909     if(aNodesIter->more()){
01910       aNode0 = (SMDS_MeshNode*) aNodesIter->next();
01911       aNode1 = aNode0;
01912       const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
01913       aNodeId[0] = aNodeId[1] = aNodes->GetID();
01914     }
01915     for(; aNodesIter->more(); ) {
01916       aNode2 = (SMDS_MeshNode*) aNodesIter->next();
01917       long anId = aNode2->GetID();
01918       aNodeId[2] = anId;
01919 
01920       Value aValue(aNodeId[1],aNodeId[2]);
01921       MValues::iterator aItr = theValues.find(aValue);
01922       if (aItr != theValues.end()){
01923         aItr->second += 1;
01924         //aNbConnects = nb;
01925       }
01926       else {
01927         theValues[aValue] = 1;
01928         //aNbConnects = 1;
01929       }
01930       //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
01931       aNodeId[1] = aNodeId[2];
01932       aNode1 = aNode2;
01933     }
01934     Value aValue(aNodeId[0],aNodeId[2]);
01935     MValues::iterator aItr = theValues.find(aValue);
01936     if (aItr != theValues.end()) {
01937       aItr->second += 1;
01938       //aNbConnects = nb;
01939     }
01940     else {
01941       theValues[aValue] = 1;
01942       //aNbConnects = 1;
01943     }
01944     //cout << "NodeIds: "<<aNodeId[0]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
01945   }
01946 
01947 }
01948 
01949 /*
01950                             PREDICATES
01951 */
01952 
01953 /*
01954   Class       : BadOrientedVolume
01955   Description : Predicate bad oriented volumes
01956 */
01957 
01958 BadOrientedVolume::BadOrientedVolume()
01959 {
01960   myMesh = 0;
01961 }
01962 
01963 void BadOrientedVolume::SetMesh( const SMDS_Mesh* theMesh )
01964 {
01965   myMesh = theMesh;
01966 }
01967 
01968 bool BadOrientedVolume::IsSatisfy( long theId )
01969 {
01970   if ( myMesh == 0 )
01971     return false;
01972 
01973   SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
01974   return !vTool.IsForward();
01975 }
01976 
01977 SMDSAbs_ElementType BadOrientedVolume::GetType() const
01978 {
01979   return SMDSAbs_Volume;
01980 }
01981 
01982 /*
01983   Class       : BareBorderVolume
01984 */
01985 
01986 bool BareBorderVolume::IsSatisfy(long theElementId )
01987 {
01988   SMDS_VolumeTool  myTool;
01989   if ( myTool.Set( myMesh->FindElement(theElementId)))
01990   {
01991     for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
01992       if ( myTool.IsFreeFace( iF ))
01993       {
01994         const SMDS_MeshNode** n = myTool.GetFaceNodes(iF);
01995         vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF));
01996         if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*Nomedium=*/false))
01997           return true;
01998       }
01999   }
02000   return false;
02001 }
02002 
02003 /*
02004   Class       : BareBorderFace
02005 */
02006 
02007 bool BareBorderFace::IsSatisfy(long theElementId )
02008 {
02009   bool ok = false;
02010   if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
02011   {
02012     if ( face->GetType() == SMDSAbs_Face )
02013     {
02014       int nbN = face->NbCornerNodes();
02015       for ( int i = 0; i < nbN && !ok; ++i )
02016       {
02017         // check if a link is shared by another face
02018         const SMDS_MeshNode* n1 = face->GetNode( i );
02019         const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
02020         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
02021         bool isShared = false;
02022         while ( !isShared && fIt->more() )
02023         {
02024           const SMDS_MeshElement* f = fIt->next();
02025           isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
02026         }
02027         if ( !isShared )
02028         {
02029           const int iQuad = face->IsQuadratic();
02030           myLinkNodes.resize( 2 + iQuad);
02031           myLinkNodes[0] = n1;
02032           myLinkNodes[1] = n2;
02033           if ( iQuad )
02034             myLinkNodes[2] = face->GetNode( i+nbN );
02035           ok = !myMesh->FindElement( myLinkNodes, SMDSAbs_Edge, /*noMedium=*/false);
02036         }
02037       }
02038     }
02039   }
02040   return ok;
02041 }
02042 
02043 /*
02044   Class       : OverConstrainedVolume
02045 */
02046 
02047 bool OverConstrainedVolume::IsSatisfy(long theElementId )
02048 {
02049   // An element is over-constrained if it has N-1 free borders where
02050   // N is the number of edges/faces for a 2D/3D element.
02051   SMDS_VolumeTool  myTool;
02052   if ( myTool.Set( myMesh->FindElement(theElementId)))
02053   {
02054     int nbSharedFaces = 0;
02055     for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
02056       if ( !myTool.IsFreeFace( iF ) && ++nbSharedFaces > 1 )
02057         break;
02058     return ( nbSharedFaces == 1 );
02059   }
02060   return false;
02061 }
02062 
02063 /*
02064   Class       : OverConstrainedFace
02065 */
02066 
02067 bool OverConstrainedFace::IsSatisfy(long theElementId )
02068 {
02069   // An element is over-constrained if it has N-1 free borders where
02070   // N is the number of edges/faces for a 2D/3D element.
02071   if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
02072     if ( face->GetType() == SMDSAbs_Face )
02073     {
02074       int nbSharedBorders = 0;
02075       int nbN = face->NbCornerNodes();
02076       for ( int i = 0; i < nbN; ++i )
02077       {
02078         // check if a link is shared by another face
02079         const SMDS_MeshNode* n1 = face->GetNode( i );
02080         const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
02081         SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
02082         bool isShared = false;
02083         while ( !isShared && fIt->more() )
02084         {
02085           const SMDS_MeshElement* f = fIt->next();
02086           isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
02087         }
02088         if ( isShared && ++nbSharedBorders > 1 )
02089           break;
02090       }
02091       return ( nbSharedBorders == 1 );
02092     }
02093   return false;
02094 }
02095 
02096 /*
02097   Class       : CoincidentNodes
02098   Description : Predicate of Coincident nodes
02099 */
02100 
02101 CoincidentNodes::CoincidentNodes()
02102 {
02103   myToler = 1e-5;
02104 }
02105 
02106 bool CoincidentNodes::IsSatisfy( long theElementId )
02107 {
02108   return myCoincidentIDs.Contains( theElementId );
02109 }
02110 
02111 SMDSAbs_ElementType CoincidentNodes::GetType() const
02112 {
02113   return SMDSAbs_Node;
02114 }
02115 
02116 void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
02117 {
02118   myMeshModifTracer.SetMesh( theMesh );
02119   if ( myMeshModifTracer.IsMeshModified() )
02120   {
02121     TIDSortedNodeSet nodesToCheck;
02122     SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true);
02123     while ( nIt->more() )
02124       nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
02125 
02126     list< list< const SMDS_MeshNode*> > nodeGroups;
02127     SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler );
02128 
02129     myCoincidentIDs.Clear();
02130     list< list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin();
02131     for ( ; groupIt != nodeGroups.end(); ++groupIt )
02132     {
02133       list< const SMDS_MeshNode*>& coincNodes = *groupIt;
02134       list< const SMDS_MeshNode*>::iterator n = coincNodes.begin();
02135       for ( ; n != coincNodes.end(); ++n )
02136         myCoincidentIDs.Add( (*n)->GetID() );
02137     }
02138   }
02139 }
02140 
02141 /*
02142   Class       : CoincidentElements
02143   Description : Predicate of Coincident Elements
02144   Note        : This class is suitable only for visualization of Coincident Elements
02145 */
02146 
02147 CoincidentElements::CoincidentElements()
02148 {
02149   myMesh = 0;
02150 }
02151 
02152 void CoincidentElements::SetMesh( const SMDS_Mesh* theMesh )
02153 {
02154   myMesh = theMesh;
02155 }
02156 
02157 bool CoincidentElements::IsSatisfy( long theElementId )
02158 {
02159   if ( !myMesh ) return false;
02160 
02161   if ( const SMDS_MeshElement* e = myMesh->FindElement( theElementId ))
02162   {
02163     if ( e->GetType() != GetType() ) return false;
02164     set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() );
02165     const int nbNodes = e->NbNodes();
02166     SMDS_ElemIteratorPtr invIt = (*elemNodes.begin())->GetInverseElementIterator( GetType() );
02167     while ( invIt->more() )
02168     {
02169       const SMDS_MeshElement* e2 = invIt->next();
02170       if ( e2 == e || e2->NbNodes() != nbNodes ) continue;
02171 
02172       bool sameNodes = true;
02173       for ( size_t i = 0; i < elemNodes.size() && sameNodes; ++i )
02174         sameNodes = ( elemNodes.count( e2->GetNode( i )));
02175       if ( sameNodes )
02176         return true;
02177     }
02178   }
02179   return false;
02180 }
02181 
02182 SMDSAbs_ElementType CoincidentElements1D::GetType() const
02183 {
02184   return SMDSAbs_Edge;
02185 }
02186 SMDSAbs_ElementType CoincidentElements2D::GetType() const
02187 {
02188   return SMDSAbs_Face;
02189 }
02190 SMDSAbs_ElementType CoincidentElements3D::GetType() const
02191 {
02192   return SMDSAbs_Volume;
02193 }
02194 
02195 
02196 /*
02197   Class       : FreeBorders
02198   Description : Predicate for free borders
02199 */
02200 
02201 FreeBorders::FreeBorders()
02202 {
02203   myMesh = 0;
02204 }
02205 
02206 void FreeBorders::SetMesh( const SMDS_Mesh* theMesh )
02207 {
02208   myMesh = theMesh;
02209 }
02210 
02211 bool FreeBorders::IsSatisfy( long theId )
02212 {
02213   return getNbMultiConnection( myMesh, theId ) == 1;
02214 }
02215 
02216 SMDSAbs_ElementType FreeBorders::GetType() const
02217 {
02218   return SMDSAbs_Edge;
02219 }
02220 
02221 
02222 /*
02223   Class       : FreeEdges
02224   Description : Predicate for free Edges
02225 */
02226 FreeEdges::FreeEdges()
02227 {
02228   myMesh = 0;
02229 }
02230 
02231 void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
02232 {
02233   myMesh = theMesh;
02234 }
02235 
02236 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId  )
02237 {
02238   TColStd_MapOfInteger aMap;
02239   for ( int i = 0; i < 2; i++ )
02240   {
02241     SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator(SMDSAbs_Face);
02242     while( anElemIter->more() )
02243     {
02244       if ( const SMDS_MeshElement* anElem = anElemIter->next())
02245       {
02246         const int anId = anElem->GetID();
02247         if ( anId != theFaceId && !aMap.Add( anId ))
02248           return false;
02249       }
02250     }
02251   }
02252   return true;
02253 }
02254 
02255 bool FreeEdges::IsSatisfy( long theId )
02256 {
02257   if ( myMesh == 0 )
02258     return false;
02259 
02260   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
02261   if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
02262     return false;
02263 
02264   SMDS_ElemIteratorPtr anIter;
02265   if ( aFace->IsQuadratic() ) {
02266     anIter = dynamic_cast<const SMDS_VtkFace*>
02267       (aFace)->interlacedNodesElemIterator();
02268   }
02269   else {
02270     anIter = aFace->nodesIterator();
02271   }
02272   if ( !anIter )
02273     return false;
02274 
02275   int i = 0, nbNodes = aFace->NbNodes();
02276   std::vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
02277   while( anIter->more() )
02278   {
02279     const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
02280     if ( aNode == 0 )
02281       return false;
02282     aNodes[ i++ ] = aNode;
02283   }
02284   aNodes[ nbNodes ] = aNodes[ 0 ];
02285 
02286   for ( i = 0; i < nbNodes; i++ )
02287     if ( IsFreeEdge( &aNodes[ i ], theId ) )
02288       return true;
02289 
02290   return false;
02291 }
02292 
02293 SMDSAbs_ElementType FreeEdges::GetType() const
02294 {
02295   return SMDSAbs_Face;
02296 }
02297 
02298 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
02299   myElemId(theElemId)
02300 {
02301   myPntId[0] = thePntId1;  myPntId[1] = thePntId2;
02302   if(thePntId1 > thePntId2){
02303     myPntId[1] = thePntId1;  myPntId[0] = thePntId2;
02304   }
02305 }
02306 
02307 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
02308   if(myPntId[0] < x.myPntId[0]) return true;
02309   if(myPntId[0] == x.myPntId[0])
02310     if(myPntId[1] < x.myPntId[1]) return true;
02311   return false;
02312 }
02313 
02314 inline void UpdateBorders(const FreeEdges::Border& theBorder,
02315                           FreeEdges::TBorders& theRegistry,
02316                           FreeEdges::TBorders& theContainer)
02317 {
02318   if(theRegistry.find(theBorder) == theRegistry.end()){
02319     theRegistry.insert(theBorder);
02320     theContainer.insert(theBorder);
02321   }else{
02322     theContainer.erase(theBorder);
02323   }
02324 }
02325 
02326 void FreeEdges::GetBoreders(TBorders& theBorders)
02327 {
02328   TBorders aRegistry;
02329   SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
02330   for(; anIter->more(); ){
02331     const SMDS_MeshFace* anElem = anIter->next();
02332     long anElemId = anElem->GetID();
02333     SMDS_ElemIteratorPtr aNodesIter;
02334     if ( anElem->IsQuadratic() )
02335       aNodesIter = static_cast<const SMDS_VtkFace*>(anElem)->
02336         interlacedNodesElemIterator();
02337     else
02338       aNodesIter = anElem->nodesIterator();
02339     long aNodeId[2];
02340     const SMDS_MeshElement* aNode;
02341     if(aNodesIter->more()){
02342       aNode = aNodesIter->next();
02343       aNodeId[0] = aNodeId[1] = aNode->GetID();
02344     }
02345     for(; aNodesIter->more(); ){
02346       aNode = aNodesIter->next();
02347       long anId = aNode->GetID();
02348       Border aBorder(anElemId,aNodeId[1],anId);
02349       aNodeId[1] = anId;
02350       UpdateBorders(aBorder,aRegistry,theBorders);
02351     }
02352     Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
02353     UpdateBorders(aBorder,aRegistry,theBorders);
02354   }
02355 }
02356 
02357 
02358 /*
02359   Class       : FreeNodes
02360   Description : Predicate for free nodes
02361 */
02362 
02363 FreeNodes::FreeNodes()
02364 {
02365   myMesh = 0;
02366 }
02367 
02368 void FreeNodes::SetMesh( const SMDS_Mesh* theMesh )
02369 {
02370   myMesh = theMesh;
02371 }
02372 
02373 bool FreeNodes::IsSatisfy( long theNodeId )
02374 {
02375   const SMDS_MeshNode* aNode = myMesh->FindNode( theNodeId );
02376   if (!aNode)
02377     return false;
02378 
02379   return (aNode->NbInverseElements() < 1);
02380 }
02381 
02382 SMDSAbs_ElementType FreeNodes::GetType() const
02383 {
02384   return SMDSAbs_Node;
02385 }
02386 
02387 
02388 /*
02389   Class       : FreeFaces
02390   Description : Predicate for free faces
02391 */
02392 
02393 FreeFaces::FreeFaces()
02394 {
02395   myMesh = 0;
02396 }
02397 
02398 void FreeFaces::SetMesh( const SMDS_Mesh* theMesh )
02399 {
02400   myMesh = theMesh;
02401 }
02402 
02403 bool FreeFaces::IsSatisfy( long theId )
02404 {
02405   if (!myMesh) return false;
02406   // check that faces nodes refers to less than two common volumes
02407   const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
02408   if ( !aFace || aFace->GetType() != SMDSAbs_Face )
02409     return false;
02410 
02411   int nbNode = aFace->NbNodes();
02412 
02413   // collect volumes check that number of volumss with count equal nbNode not less than 2
02414   typedef map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
02415   typedef map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
02416   TMapOfVolume mapOfVol;
02417 
02418   SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator();
02419   while ( nodeItr->more() ) {
02420     const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
02421     if ( !aNode ) continue;
02422     SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume);
02423     while ( volItr->more() ) {
02424       SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
02425       TItrMapOfVolume itr = mapOfVol.insert(make_pair(aVol, 0)).first;
02426       (*itr).second++;
02427     } 
02428   }
02429   int nbVol = 0;
02430   TItrMapOfVolume volItr = mapOfVol.begin();
02431   TItrMapOfVolume volEnd = mapOfVol.end();
02432   for ( ; volItr != volEnd; ++volItr )
02433     if ( (*volItr).second >= nbNode )
02434        nbVol++;
02435   // face is not free if number of volumes constructed on thier nodes more than one
02436   return (nbVol < 2);
02437 }
02438 
02439 SMDSAbs_ElementType FreeFaces::GetType() const
02440 {
02441   return SMDSAbs_Face;
02442 }
02443 
02444 /*
02445   Class       : LinearOrQuadratic
02446   Description : Predicate to verify whether a mesh element is linear
02447 */
02448 
02449 LinearOrQuadratic::LinearOrQuadratic()
02450 {
02451   myMesh = 0;
02452 }
02453 
02454 void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh )
02455 {
02456   myMesh = theMesh;
02457 }
02458 
02459 bool LinearOrQuadratic::IsSatisfy( long theId )
02460 {
02461   if (!myMesh) return false;
02462   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
02463   if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) )
02464     return false;
02465   return (!anElem->IsQuadratic());
02466 }
02467 
02468 void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType )
02469 {
02470   myType = theType;
02471 }
02472 
02473 SMDSAbs_ElementType LinearOrQuadratic::GetType() const
02474 {
02475   return myType;
02476 }
02477 
02478 /*
02479   Class       : GroupColor
02480   Description : Functor for check color of group to whic mesh element belongs to
02481 */
02482 
02483 GroupColor::GroupColor()
02484 {
02485 }
02486 
02487 bool GroupColor::IsSatisfy( long theId )
02488 {
02489   return (myIDs.find( theId ) != myIDs.end());
02490 }
02491 
02492 void GroupColor::SetType( SMDSAbs_ElementType theType )
02493 {
02494   myType = theType;
02495 }
02496 
02497 SMDSAbs_ElementType GroupColor::GetType() const
02498 {
02499   return myType;
02500 }
02501 
02502 static bool isEqual( const Quantity_Color& theColor1,
02503                      const Quantity_Color& theColor2 )
02504 {
02505   // tolerance to compare colors
02506   const double tol = 5*1e-3;
02507   return ( fabs( theColor1.Red() - theColor2.Red() ) < tol &&
02508            fabs( theColor1.Green() - theColor2.Green() ) < tol &&
02509            fabs( theColor1.Blue() - theColor2.Blue() ) < tol );
02510 }
02511 
02512 
02513 void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
02514 {
02515   myIDs.clear();
02516   
02517   const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
02518   if ( !aMesh )
02519     return;
02520 
02521   int nbGrp = aMesh->GetNbGroups();
02522   if ( !nbGrp )
02523     return;
02524   
02525   // iterates on groups and find necessary elements ids
02526   const std::set<SMESHDS_GroupBase*>& aGroups = aMesh->GetGroups();
02527   set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
02528   for (; GrIt != aGroups.end(); GrIt++) {
02529     SMESHDS_GroupBase* aGrp = (*GrIt);
02530     if ( !aGrp )
02531       continue;
02532     // check type and color of group
02533     if ( !isEqual( myColor, aGrp->GetColor() ) )
02534       continue;
02535     if ( myType != SMDSAbs_All && myType != (SMDSAbs_ElementType)aGrp->GetType() )
02536       continue;
02537 
02538     SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType();
02539     if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) {
02540       // add elements IDS into control
02541       int aSize = aGrp->Extent();
02542       for (int i = 0; i < aSize; i++)
02543         myIDs.insert( aGrp->GetID(i+1) );
02544     }
02545   }
02546 }
02547 
02548 void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
02549 {
02550   TCollection_AsciiString aStr = theStr;
02551   aStr.RemoveAll( ' ' );
02552   aStr.RemoveAll( '\t' );
02553   for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) )
02554     aStr.Remove( aPos, 2 );
02555   Standard_Real clr[3];
02556   clr[0] = clr[1] = clr[2] = 0.;
02557   for ( int i = 0; i < 3; i++ ) {
02558     TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 );
02559     if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() )
02560       clr[i] = tmpStr.RealValue();
02561   }
02562   myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB );
02563 }
02564 
02565 //=======================================================================
02566 // name    : GetRangeStr
02567 // Purpose : Get range as a string.
02568 //           Example: "1,2,3,50-60,63,67,70-"
02569 //=======================================================================
02570 void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
02571 {
02572   theResStr.Clear();
02573   theResStr += TCollection_AsciiString( myColor.Red() );
02574   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() );
02575   theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
02576 }
02577 
02578 /*
02579   Class       : ElemGeomType
02580   Description : Predicate to check element geometry type
02581 */
02582 
02583 ElemGeomType::ElemGeomType()
02584 {
02585   myMesh = 0;
02586   myType = SMDSAbs_All;
02587   myGeomType = SMDSGeom_TRIANGLE;
02588 }
02589 
02590 void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh )
02591 {
02592   myMesh = theMesh;
02593 }
02594 
02595 bool ElemGeomType::IsSatisfy( long theId )
02596 {
02597   if (!myMesh) return false;
02598   const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
02599   if ( !anElem )
02600     return false;
02601   const SMDSAbs_ElementType anElemType = anElem->GetType();
02602   if ( myType != SMDSAbs_All && anElemType != myType )
02603     return false;
02604   const int aNbNode = anElem->NbNodes();
02605   bool isOk = false;
02606   switch( anElemType )
02607   {
02608   case SMDSAbs_Node:
02609     isOk = (myGeomType == SMDSGeom_POINT);
02610     break;
02611 
02612   case SMDSAbs_Edge:
02613     isOk = (myGeomType == SMDSGeom_EDGE);
02614     break;
02615 
02616   case SMDSAbs_Face:
02617     if ( myGeomType == SMDSGeom_TRIANGLE )
02618       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 6 : aNbNode == 3));
02619     else if ( myGeomType == SMDSGeom_QUADRANGLE )
02620       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? ( aNbNode == 8 || aNbNode == 9 ) : aNbNode == 4));
02621     else if ( myGeomType == SMDSGeom_POLYGON )
02622       isOk = anElem->IsPoly();
02623     break;
02624 
02625   case SMDSAbs_Volume:
02626     if ( myGeomType == SMDSGeom_TETRA )
02627       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 10 : aNbNode == 4));
02628     else if ( myGeomType == SMDSGeom_PYRAMID )
02629       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 13 : aNbNode == 5));
02630     else if ( myGeomType == SMDSGeom_PENTA )
02631       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? aNbNode == 15 : aNbNode == 6));
02632     else if ( myGeomType == SMDSGeom_HEXA )
02633       isOk = (!anElem->IsPoly() && (anElem->IsQuadratic() ? ( aNbNode == 20 || aNbNode == 27 ): aNbNode == 8));
02634     else if ( myGeomType == SMDSGeom_HEXAGONAL_PRISM )
02635       isOk = (anElem->GetEntityType() == SMDSEntity_Hexagonal_Prism );
02636      else if ( myGeomType == SMDSGeom_POLYHEDRA )
02637       isOk = anElem->IsPoly();
02638     break;
02639     default: break;
02640   }
02641   return isOk;
02642 }
02643 
02644 void ElemGeomType::SetType( SMDSAbs_ElementType theType )
02645 {
02646   myType = theType;
02647 }
02648 
02649 SMDSAbs_ElementType ElemGeomType::GetType() const
02650 {
02651   return myType;
02652 }
02653 
02654 void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType )
02655 {
02656   myGeomType = theType;
02657 }
02658 
02659 SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
02660 {
02661   return myGeomType;
02662 }
02663 
02664 //================================================================================
02668 //================================================================================
02669 
02670 CoplanarFaces::CoplanarFaces()
02671   : myFaceID(0), myToler(0)
02672 {
02673 }
02674 void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
02675 {
02676   myMeshModifTracer.SetMesh( theMesh );
02677   if ( myMeshModifTracer.IsMeshModified() )
02678   {
02679     // Build a set of coplanar face ids
02680 
02681     myCoplanarIDs.clear();
02682 
02683     if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler )
02684       return;
02685 
02686     const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID );
02687     if ( !face || face->GetType() != SMDSAbs_Face )
02688       return;
02689 
02690     bool normOK;
02691     gp_Vec myNorm = getNormale( static_cast<const SMDS_MeshFace*>(face), &normOK );
02692     if (!normOK)
02693       return;
02694 
02695     const double radianTol = myToler * M_PI / 180.;
02696     typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TFaceIt;
02697     std::set<const SMDS_MeshElement*> checkedFaces, checkedNodes;
02698     std::list<const SMDS_MeshElement*> faceQueue( 1, face );
02699     while ( !faceQueue.empty() )
02700     {
02701       face = faceQueue.front();
02702       if ( checkedFaces.insert( face ).second )
02703       {
02704         gp_Vec norm = getNormale( static_cast<const SMDS_MeshFace*>(face), &normOK );
02705         if (!normOK || myNorm.Angle( norm ) <= radianTol)
02706         {
02707           myCoplanarIDs.insert( face->GetID() );
02708           std::set<const SMDS_MeshElement*> neighborFaces;
02709           for ( int i = 0; i < face->NbCornerNodes(); ++i )
02710           {
02711             const SMDS_MeshNode* n = face->GetNode( i );
02712             if ( checkedNodes.insert( n ).second )
02713               neighborFaces.insert( TFaceIt( n->GetInverseElementIterator(SMDSAbs_Face)),
02714                                     TFaceIt());
02715           }
02716           faceQueue.insert( faceQueue.end(), neighborFaces.begin(), neighborFaces.end() );
02717         }
02718       }
02719       faceQueue.pop_front();
02720     }
02721   }
02722 }
02723 bool CoplanarFaces::IsSatisfy( long theElementId )
02724 {
02725   return myCoplanarIDs.count( theElementId );
02726 }
02727 
02728 /*
02729   *Class       : RangeOfIds
02730   *Description : Predicate for Range of Ids.
02731   *              Range may be specified with two ways.
02732   *              1. Using AddToRange method
02733   *              2. With SetRangeStr method. Parameter of this method is a string
02734   *                 like as "1,2,3,50-60,63,67,70-"
02735 */
02736 
02737 //=======================================================================
02738 // name    : RangeOfIds
02739 // Purpose : Constructor
02740 //=======================================================================
02741 RangeOfIds::RangeOfIds()
02742 {
02743   myMesh = 0;
02744   myType = SMDSAbs_All;
02745 }
02746 
02747 //=======================================================================
02748 // name    : SetMesh
02749 // Purpose : Set mesh
02750 //=======================================================================
02751 void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh )
02752 {
02753   myMesh = theMesh;
02754 }
02755 
02756 //=======================================================================
02757 // name    : AddToRange
02758 // Purpose : Add ID to the range
02759 //=======================================================================
02760 bool RangeOfIds::AddToRange( long theEntityId )
02761 {
02762   myIds.Add( theEntityId );
02763   return true;
02764 }
02765 
02766 //=======================================================================
02767 // name    : GetRangeStr
02768 // Purpose : Get range as a string.
02769 //           Example: "1,2,3,50-60,63,67,70-"
02770 //=======================================================================
02771 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
02772 {
02773   theResStr.Clear();
02774 
02775   TColStd_SequenceOfInteger     anIntSeq;
02776   TColStd_SequenceOfAsciiString aStrSeq;
02777 
02778   TColStd_MapIteratorOfMapOfInteger anIter( myIds );
02779   for ( ; anIter.More(); anIter.Next() )
02780   {
02781     int anId = anIter.Key();
02782     TCollection_AsciiString aStr( anId );
02783     anIntSeq.Append( anId );
02784     aStrSeq.Append( aStr );
02785   }
02786 
02787   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
02788   {
02789     int aMinId = myMin( i );
02790     int aMaxId = myMax( i );
02791 
02792     TCollection_AsciiString aStr;
02793     if ( aMinId != IntegerFirst() )
02794       aStr += aMinId;
02795 
02796     aStr += "-";
02797 
02798     if ( aMaxId != IntegerLast() )
02799       aStr += aMaxId;
02800 
02801     // find position of the string in result sequence and insert string in it
02802     if ( anIntSeq.Length() == 0 )
02803     {
02804       anIntSeq.Append( aMinId );
02805       aStrSeq.Append( aStr );
02806     }
02807     else
02808     {
02809       if ( aMinId < anIntSeq.First() )
02810       {
02811         anIntSeq.Prepend( aMinId );
02812         aStrSeq.Prepend( aStr );
02813       }
02814       else if ( aMinId > anIntSeq.Last() )
02815       {
02816         anIntSeq.Append( aMinId );
02817         aStrSeq.Append( aStr );
02818       }
02819       else
02820         for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
02821           if ( aMinId < anIntSeq( j ) )
02822           {
02823             anIntSeq.InsertBefore( j, aMinId );
02824             aStrSeq.InsertBefore( j, aStr );
02825             break;
02826           }
02827     }
02828   }
02829 
02830   if ( aStrSeq.Length() == 0 )
02831     return;
02832 
02833   theResStr = aStrSeq( 1 );
02834   for ( int j = 2, k = aStrSeq.Length(); j <= k; j++  )
02835   {
02836     theResStr += ",";
02837     theResStr += aStrSeq( j );
02838   }
02839 }
02840 
02841 //=======================================================================
02842 // name    : SetRangeStr
02843 // Purpose : Define range with string
02844 //           Example of entry string: "1,2,3,50-60,63,67,70-"
02845 //=======================================================================
02846 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
02847 {
02848   myMin.Clear();
02849   myMax.Clear();
02850   myIds.Clear();
02851 
02852   TCollection_AsciiString aStr = theStr;
02853   aStr.RemoveAll( ' ' );
02854   aStr.RemoveAll( '\t' );
02855 
02856   for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) )
02857     aStr.Remove( aPos, 2 );
02858 
02859   TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
02860   int i = 1;
02861   while ( tmpStr != "" )
02862   {
02863     tmpStr = aStr.Token( ",", i++ );
02864     int aPos = tmpStr.Search( '-' );
02865 
02866     if ( aPos == -1 )
02867     {
02868       if ( tmpStr.IsIntegerValue() )
02869         myIds.Add( tmpStr.IntegerValue() );
02870       else
02871         return false;
02872     }
02873     else
02874     {
02875       TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
02876       TCollection_AsciiString aMinStr = tmpStr;
02877 
02878       while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
02879       while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
02880 
02881       if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) ||
02882            (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) )
02883         return false;
02884 
02885       myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
02886       myMax.Append( aMaxStr.IsEmpty() ? IntegerLast()  : aMaxStr.IntegerValue() );
02887     }
02888   }
02889 
02890   return true;
02891 }
02892 
02893 //=======================================================================
02894 // name    : GetType
02895 // Purpose : Get type of supported entities
02896 //=======================================================================
02897 SMDSAbs_ElementType RangeOfIds::GetType() const
02898 {
02899   return myType;
02900 }
02901 
02902 //=======================================================================
02903 // name    : SetType
02904 // Purpose : Set type of supported entities
02905 //=======================================================================
02906 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
02907 {
02908   myType = theType;
02909 }
02910 
02911 //=======================================================================
02912 // name    : IsSatisfy
02913 // Purpose : Verify whether entity satisfies to this rpedicate
02914 //=======================================================================
02915 bool RangeOfIds::IsSatisfy( long theId )
02916 {
02917   if ( !myMesh )
02918     return false;
02919 
02920   if ( myType == SMDSAbs_Node )
02921   {
02922     if ( myMesh->FindNode( theId ) == 0 )
02923       return false;
02924   }
02925   else
02926   {
02927     const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
02928     if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All ))
02929       return false;
02930   }
02931 
02932   if ( myIds.Contains( theId ) )
02933     return true;
02934 
02935   for ( int i = 1, n = myMin.Length(); i <= n; i++ )
02936     if ( theId >= myMin( i ) && theId <= myMax( i ) )
02937       return true;
02938 
02939   return false;
02940 }
02941 
02942 /*
02943   Class       : Comparator
02944   Description : Base class for comparators
02945 */
02946 Comparator::Comparator():
02947   myMargin(0)
02948 {}
02949 
02950 Comparator::~Comparator()
02951 {}
02952 
02953 void Comparator::SetMesh( const SMDS_Mesh* theMesh )
02954 {
02955   if ( myFunctor )
02956     myFunctor->SetMesh( theMesh );
02957 }
02958 
02959 void Comparator::SetMargin( double theValue )
02960 {
02961   myMargin = theValue;
02962 }
02963 
02964 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
02965 {
02966   myFunctor = theFunct;
02967 }
02968 
02969 SMDSAbs_ElementType Comparator::GetType() const
02970 {
02971   return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
02972 }
02973 
02974 double Comparator::GetMargin()
02975 {
02976   return myMargin;
02977 }
02978 
02979 
02980 /*
02981   Class       : LessThan
02982   Description : Comparator "<"
02983 */
02984 bool LessThan::IsSatisfy( long theId )
02985 {
02986   return myFunctor && myFunctor->GetValue( theId ) < myMargin;
02987 }
02988 
02989 
02990 /*
02991   Class       : MoreThan
02992   Description : Comparator ">"
02993 */
02994 bool MoreThan::IsSatisfy( long theId )
02995 {
02996   return myFunctor && myFunctor->GetValue( theId ) > myMargin;
02997 }
02998 
02999 
03000 /*
03001   Class       : EqualTo
03002   Description : Comparator "="
03003 */
03004 EqualTo::EqualTo():
03005   myToler(Precision::Confusion())
03006 {}
03007 
03008 bool EqualTo::IsSatisfy( long theId )
03009 {
03010   return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
03011 }
03012 
03013 void EqualTo::SetTolerance( double theToler )
03014 {
03015   myToler = theToler;
03016 }
03017 
03018 double EqualTo::GetTolerance()
03019 {
03020   return myToler;
03021 }
03022 
03023 /*
03024   Class       : LogicalNOT
03025   Description : Logical NOT predicate
03026 */
03027 LogicalNOT::LogicalNOT()
03028 {}
03029 
03030 LogicalNOT::~LogicalNOT()
03031 {}
03032 
03033 bool LogicalNOT::IsSatisfy( long theId )
03034 {
03035   return myPredicate && !myPredicate->IsSatisfy( theId );
03036 }
03037 
03038 void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh )
03039 {
03040   if ( myPredicate )
03041     myPredicate->SetMesh( theMesh );
03042 }
03043 
03044 void LogicalNOT::SetPredicate( PredicatePtr thePred )
03045 {
03046   myPredicate = thePred;
03047 }
03048 
03049 SMDSAbs_ElementType LogicalNOT::GetType() const
03050 {
03051   return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
03052 }
03053 
03054 
03055 /*
03056   Class       : LogicalBinary
03057   Description : Base class for binary logical predicate
03058 */
03059 LogicalBinary::LogicalBinary()
03060 {}
03061 
03062 LogicalBinary::~LogicalBinary()
03063 {}
03064 
03065 void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh )
03066 {
03067   if ( myPredicate1 )
03068     myPredicate1->SetMesh( theMesh );
03069 
03070   if ( myPredicate2 )
03071     myPredicate2->SetMesh( theMesh );
03072 }
03073 
03074 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
03075 {
03076   myPredicate1 = thePredicate;
03077 }
03078 
03079 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
03080 {
03081   myPredicate2 = thePredicate;
03082 }
03083 
03084 SMDSAbs_ElementType LogicalBinary::GetType() const
03085 {
03086   if ( !myPredicate1 || !myPredicate2 )
03087     return SMDSAbs_All;
03088 
03089   SMDSAbs_ElementType aType1 = myPredicate1->GetType();
03090   SMDSAbs_ElementType aType2 = myPredicate2->GetType();
03091 
03092   return aType1 == aType2 ? aType1 : SMDSAbs_All;
03093 }
03094 
03095 
03096 /*
03097   Class       : LogicalAND
03098   Description : Logical AND
03099 */
03100 bool LogicalAND::IsSatisfy( long theId )
03101 {
03102   return
03103     myPredicate1 &&
03104     myPredicate2 &&
03105     myPredicate1->IsSatisfy( theId ) &&
03106     myPredicate2->IsSatisfy( theId );
03107 }
03108 
03109 
03110 /*
03111   Class       : LogicalOR
03112   Description : Logical OR
03113 */
03114 bool LogicalOR::IsSatisfy( long theId )
03115 {
03116   return
03117     myPredicate1 &&
03118     myPredicate2 &&
03119     (myPredicate1->IsSatisfy( theId ) ||
03120     myPredicate2->IsSatisfy( theId ));
03121 }
03122 
03123 
03124 /*
03125                               FILTER
03126 */
03127 
03128 Filter::Filter()
03129 {}
03130 
03131 Filter::~Filter()
03132 {}
03133 
03134 void Filter::SetPredicate( PredicatePtr thePredicate )
03135 {
03136   myPredicate = thePredicate;
03137 }
03138 
03139 template<class TElement, class TIterator, class TPredicate>
03140 inline void FillSequence(const TIterator& theIterator,
03141                          TPredicate& thePredicate,
03142                          Filter::TIdSequence& theSequence)
03143 {
03144   if ( theIterator ) {
03145     while( theIterator->more() ) {
03146       TElement anElem = theIterator->next();
03147       long anId = anElem->GetID();
03148       if ( thePredicate->IsSatisfy( anId ) )
03149         theSequence.push_back( anId );
03150     }
03151   }
03152 }
03153 
03154 void
03155 Filter::
03156 GetElementsId( const SMDS_Mesh* theMesh,
03157                PredicatePtr thePredicate,
03158                TIdSequence& theSequence )
03159 {
03160   theSequence.clear();
03161 
03162   if ( !theMesh || !thePredicate )
03163     return;
03164 
03165   thePredicate->SetMesh( theMesh );
03166 
03167   SMDSAbs_ElementType aType = thePredicate->GetType();
03168   switch(aType){
03169   case SMDSAbs_Node:
03170     FillSequence<const SMDS_MeshNode*>(theMesh->nodesIterator(),thePredicate,theSequence);
03171     break;
03172   case SMDSAbs_Edge:
03173     FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),thePredicate,theSequence);
03174     break;
03175   case SMDSAbs_Face:
03176     FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),thePredicate,theSequence);
03177     break;
03178   case SMDSAbs_Volume:
03179     FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),thePredicate,theSequence);
03180     break;
03181   case SMDSAbs_All:
03182     FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),thePredicate,theSequence);
03183     FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),thePredicate,theSequence);
03184     FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),thePredicate,theSequence);
03185     break;
03186   }
03187 }
03188 
03189 void
03190 Filter::GetElementsId( const SMDS_Mesh* theMesh,
03191                        Filter::TIdSequence& theSequence )
03192 {
03193   GetElementsId(theMesh,myPredicate,theSequence);
03194 }
03195 
03196 /*
03197                               ManifoldPart
03198 */
03199 
03200 typedef std::set<SMDS_MeshFace*>                    TMapOfFacePtr;
03201 
03202 /*
03203    Internal class Link
03204 */
03205 
03206 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
03207                           SMDS_MeshNode* theNode2 )
03208 {
03209   myNode1 = theNode1;
03210   myNode2 = theNode2;
03211 }
03212 
03213 ManifoldPart::Link::~Link()
03214 {
03215   myNode1 = 0;
03216   myNode2 = 0;
03217 }
03218 
03219 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
03220 {
03221   if ( myNode1 == theLink.myNode1 &&
03222        myNode2 == theLink.myNode2 )
03223     return true;
03224   else if ( myNode1 == theLink.myNode2 &&
03225             myNode2 == theLink.myNode1 )
03226     return true;
03227   else
03228     return false;
03229 }
03230 
03231 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
03232 {
03233   if(myNode1 < x.myNode1) return true;
03234   if(myNode1 == x.myNode1)
03235     if(myNode2 < x.myNode2) return true;
03236   return false;
03237 }
03238 
03239 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
03240                             const ManifoldPart::Link& theLink2 )
03241 {
03242   return theLink1.IsEqual( theLink2 );
03243 }
03244 
03245 ManifoldPart::ManifoldPart()
03246 {
03247   myMesh = 0;
03248   myAngToler = Precision::Angular();
03249   myIsOnlyManifold = true;
03250 }
03251 
03252 ManifoldPart::~ManifoldPart()
03253 {
03254   myMesh = 0;
03255 }
03256 
03257 void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh )
03258 {
03259   myMesh = theMesh;
03260   process();
03261 }
03262 
03263 SMDSAbs_ElementType ManifoldPart::GetType() const
03264 { return SMDSAbs_Face; }
03265 
03266 bool ManifoldPart::IsSatisfy( long theElementId )
03267 {
03268   return myMapIds.Contains( theElementId );
03269 }
03270 
03271 void ManifoldPart::SetAngleTolerance( const double theAngToler )
03272 { myAngToler = theAngToler; }
03273 
03274 double ManifoldPart::GetAngleTolerance() const
03275 { return myAngToler; }
03276 
03277 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
03278 { myIsOnlyManifold = theIsOnly; }
03279 
03280 void ManifoldPart::SetStartElem( const long  theStartId )
03281 { myStartElemId = theStartId; }
03282 
03283 bool ManifoldPart::process()
03284 {
03285   myMapIds.Clear();
03286   myMapBadGeomIds.Clear();
03287 
03288   myAllFacePtr.clear();
03289   myAllFacePtrIntDMap.clear();
03290   if ( !myMesh )
03291     return false;
03292 
03293   // collect all faces into own map
03294   SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
03295   for (; anFaceItr->more(); )
03296   {
03297     SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
03298     myAllFacePtr.push_back( aFacePtr );
03299     myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
03300   }
03301 
03302   SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
03303   if ( !aStartFace )
03304     return false;
03305 
03306   // the map of non manifold links and bad geometry
03307   TMapOfLink aMapOfNonManifold;
03308   TColStd_MapOfInteger aMapOfTreated;
03309 
03310   // begin cycle on faces from start index and run on vector till the end
03311   //  and from begin to start index to cover whole vector
03312   const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
03313   bool isStartTreat = false;
03314   for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
03315   {
03316     if ( fi == aStartIndx )
03317       isStartTreat = true;
03318     // as result next time when fi will be equal to aStartIndx
03319 
03320     SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
03321     if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
03322       continue;
03323 
03324     aMapOfTreated.Add( aFacePtr->GetID() );
03325     TColStd_MapOfInteger aResFaces;
03326     if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
03327                          aMapOfNonManifold, aResFaces ) )
03328       continue;
03329     TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
03330     for ( ; anItr.More(); anItr.Next() )
03331     {
03332       int aFaceId = anItr.Key();
03333       aMapOfTreated.Add( aFaceId );
03334       myMapIds.Add( aFaceId );
03335     }
03336 
03337     if ( fi == ( myAllFacePtr.size() - 1 ) )
03338       fi = 0;
03339   } // end run on vector of faces
03340   return !myMapIds.IsEmpty();
03341 }
03342 
03343 static void getLinks( const SMDS_MeshFace* theFace,
03344                       ManifoldPart::TVectorOfLink& theLinks )
03345 {
03346   int aNbNode = theFace->NbNodes();
03347   SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
03348   int i = 1;
03349   SMDS_MeshNode* aNode = 0;
03350   for ( ; aNodeItr->more() && i <= aNbNode; )
03351   {
03352 
03353     SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
03354     if ( i == 1 )
03355       aNode = aN1;
03356     i++;
03357     SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
03358     i++;
03359     ManifoldPart::Link aLink( aN1, aN2 );
03360     theLinks.push_back( aLink );
03361   }
03362 }
03363 
03364 bool ManifoldPart::findConnected
03365                  ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
03366                   SMDS_MeshFace*                           theStartFace,
03367                   ManifoldPart::TMapOfLink&                theNonManifold,
03368                   TColStd_MapOfInteger&                    theResFaces )
03369 {
03370   theResFaces.Clear();
03371   if ( !theAllFacePtrInt.size() )
03372     return false;
03373 
03374   if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
03375   {
03376     myMapBadGeomIds.Add( theStartFace->GetID() );
03377     return false;
03378   }
03379 
03380   ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
03381   ManifoldPart::TVectorOfLink aSeqOfBoundary;
03382   theResFaces.Add( theStartFace->GetID() );
03383   ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
03384 
03385   expandBoundary( aMapOfBoundary, aSeqOfBoundary,
03386                  aDMapLinkFace, theNonManifold, theStartFace );
03387 
03388   bool isDone = false;
03389   while ( !isDone && aMapOfBoundary.size() != 0 )
03390   {
03391     bool isToReset = false;
03392     ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
03393     for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
03394     {
03395       ManifoldPart::Link aLink = *pLink;
03396       if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
03397         continue;
03398       // each link could be treated only once
03399       aMapToSkip.insert( aLink );
03400 
03401       ManifoldPart::TVectorOfFacePtr aFaces;
03402       // find next
03403       if ( myIsOnlyManifold &&
03404            (theNonManifold.find( aLink ) != theNonManifold.end()) )
03405         continue;
03406       else
03407       {
03408         getFacesByLink( aLink, aFaces );
03409         // filter the element to keep only indicated elements
03410         ManifoldPart::TVectorOfFacePtr aFiltered;
03411         ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
03412         for ( ; pFace != aFaces.end(); ++pFace )
03413         {
03414           SMDS_MeshFace* aFace = *pFace;
03415           if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
03416             aFiltered.push_back( aFace );
03417         }
03418         aFaces = aFiltered;
03419         if ( aFaces.size() < 2 )  // no neihgbour faces
03420           continue;
03421         else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
03422         {
03423           theNonManifold.insert( aLink );
03424           continue;
03425         }
03426       }
03427 
03428       // compare normal with normals of neighbor element
03429       SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
03430       ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
03431       for ( ; pFace != aFaces.end(); ++pFace )
03432       {
03433         SMDS_MeshFace* aNextFace = *pFace;
03434         if ( aPrevFace == aNextFace )
03435           continue;
03436         int anNextFaceID = aNextFace->GetID();
03437         if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
03438          // should not be with non manifold restriction. probably bad topology
03439           continue;
03440         // check if face was treated and skipped
03441         if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
03442              !isInPlane( aPrevFace, aNextFace ) )
03443           continue;
03444         // add new element to connected and extend the boundaries.
03445         theResFaces.Add( anNextFaceID );
03446         expandBoundary( aMapOfBoundary, aSeqOfBoundary,
03447                         aDMapLinkFace, theNonManifold, aNextFace );
03448         isToReset = true;
03449       }
03450     }
03451     isDone = !isToReset;
03452   }
03453 
03454   return !theResFaces.IsEmpty();
03455 }
03456 
03457 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
03458                               const SMDS_MeshFace* theFace2 )
03459 {
03460   gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
03461   gp_XYZ aNorm2XYZ = getNormale( theFace2 );
03462   if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
03463   {
03464     myMapBadGeomIds.Add( theFace2->GetID() );
03465     return false;
03466   }
03467   if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
03468     return true;
03469 
03470   return false;
03471 }
03472 
03473 void ManifoldPart::expandBoundary
03474                    ( ManifoldPart::TMapOfLink&            theMapOfBoundary,
03475                      ManifoldPart::TVectorOfLink&         theSeqOfBoundary,
03476                      ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
03477                      ManifoldPart::TMapOfLink&            theNonManifold,
03478                      SMDS_MeshFace*                       theNextFace ) const
03479 {
03480   ManifoldPart::TVectorOfLink aLinks;
03481   getLinks( theNextFace, aLinks );
03482   int aNbLink = (int)aLinks.size();
03483   for ( int i = 0; i < aNbLink; i++ )
03484   {
03485     ManifoldPart::Link aLink = aLinks[ i ];
03486     if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
03487       continue;
03488     if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
03489     {
03490       if ( myIsOnlyManifold )
03491       {
03492         // remove from boundary
03493         theMapOfBoundary.erase( aLink );
03494         ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
03495         for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
03496         {
03497           ManifoldPart::Link aBoundLink = *pLink;
03498           if ( aBoundLink.IsEqual( aLink ) )
03499           {
03500             theSeqOfBoundary.erase( pLink );
03501             break;
03502           }
03503         }
03504       }
03505     }
03506     else
03507     {
03508       theMapOfBoundary.insert( aLink );
03509       theSeqOfBoundary.push_back( aLink );
03510       theDMapLinkFacePtr[ aLink ] = theNextFace;
03511     }
03512   }
03513 }
03514 
03515 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
03516                                    ManifoldPart::TVectorOfFacePtr& theFaces ) const
03517 {
03518   std::set<SMDS_MeshCell *> aSetOfFaces;
03519   // take all faces that shared first node
03520   SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
03521   for ( ; anItr->more(); )
03522   {
03523     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
03524     if ( !aFace )
03525       continue;
03526     aSetOfFaces.insert( aFace );
03527   }
03528   // take all faces that shared second node
03529   anItr = theLink.myNode2->facesIterator();
03530   // find the common part of two sets
03531   for ( ; anItr->more(); )
03532   {
03533     SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
03534     if ( aSetOfFaces.count( aFace ) )
03535       theFaces.push_back( aFace );
03536   }
03537 }
03538 
03539 
03540 /*
03541    ElementsOnSurface
03542 */
03543 
03544 ElementsOnSurface::ElementsOnSurface()
03545 {
03546   myMesh = 0;
03547   myIds.Clear();
03548   myType = SMDSAbs_All;
03549   mySurf.Nullify();
03550   myToler = Precision::Confusion();
03551   myUseBoundaries = false;
03552 }
03553 
03554 ElementsOnSurface::~ElementsOnSurface()
03555 {
03556   myMesh = 0;
03557 }
03558 
03559 void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
03560 {
03561   if ( myMesh == theMesh )
03562     return;
03563   myMesh = theMesh;
03564   process();
03565 }
03566 
03567 bool ElementsOnSurface::IsSatisfy( long theElementId )
03568 {
03569   return myIds.Contains( theElementId );
03570 }
03571 
03572 SMDSAbs_ElementType ElementsOnSurface::GetType() const
03573 { return myType; }
03574 
03575 void ElementsOnSurface::SetTolerance( const double theToler )
03576 {
03577   if ( myToler != theToler )
03578     myIds.Clear();
03579   myToler = theToler;
03580 }
03581 
03582 double ElementsOnSurface::GetTolerance() const
03583 { return myToler; }
03584 
03585 void ElementsOnSurface::SetUseBoundaries( bool theUse )
03586 {
03587   if ( myUseBoundaries != theUse ) {
03588     myUseBoundaries = theUse;
03589     SetSurface( mySurf, myType );
03590   }
03591 }
03592 
03593 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
03594                                     const SMDSAbs_ElementType theType )
03595 {
03596   myIds.Clear();
03597   myType = theType;
03598   mySurf.Nullify();
03599   if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
03600     return;
03601   mySurf = TopoDS::Face( theShape );
03602   BRepAdaptor_Surface SA( mySurf, myUseBoundaries );
03603   Standard_Real
03604     u1 = SA.FirstUParameter(),
03605     u2 = SA.LastUParameter(),
03606     v1 = SA.FirstVParameter(),
03607     v2 = SA.LastVParameter();
03608   Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf );
03609   myProjector.Init( surf, u1,u2, v1,v2 );
03610   process();
03611 }
03612 
03613 void ElementsOnSurface::process()
03614 {
03615   myIds.Clear();
03616   if ( mySurf.IsNull() )
03617     return;
03618 
03619   if ( myMesh == 0 )
03620     return;
03621 
03622   if ( myType == SMDSAbs_Face || myType == SMDSAbs_All )
03623   {
03624     myIds.ReSize( myMesh->NbFaces() );
03625     SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
03626     for(; anIter->more(); )
03627       process( anIter->next() );
03628   }
03629 
03630   if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All )
03631   {
03632     myIds.ReSize( myIds.Extent() + myMesh->NbEdges() );
03633     SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator();
03634     for(; anIter->more(); )
03635       process( anIter->next() );
03636   }
03637 
03638   if ( myType == SMDSAbs_Node )
03639   {
03640     myIds.ReSize( myMesh->NbNodes() );
03641     SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator();
03642     for(; anIter->more(); )
03643       process( anIter->next() );
03644   }
03645 }
03646 
03647 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
03648 {
03649   SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
03650   bool isSatisfy = true;
03651   for ( ; aNodeItr->more(); )
03652   {
03653     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
03654     if ( !isOnSurface( aNode ) )
03655     {
03656       isSatisfy = false;
03657       break;
03658     }
03659   }
03660   if ( isSatisfy )
03661     myIds.Add( theElemPtr->GetID() );
03662 }
03663 
03664 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode )
03665 {
03666   if ( mySurf.IsNull() )
03667     return false;
03668 
03669   gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
03670   //  double aToler2 = myToler * myToler;
03671 //   if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
03672 //   {
03673 //     gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
03674 //     if ( aPln.SquareDistance( aPnt ) > aToler2 )
03675 //       return false;
03676 //   }
03677 //   else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
03678 //   {
03679 //     gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
03680 //     double aRad = aCyl.Radius();
03681 //     gp_Ax3 anAxis = aCyl.Position();
03682 //     gp_XYZ aLoc = aCyl.Location().XYZ();
03683 //     double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
03684 //     double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
03685 //     if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
03686 //       return false;
03687 //   }
03688 //   else
03689 //     return false;
03690   myProjector.Perform( aPnt );
03691   bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler );
03692 
03693   return isOn;
03694 }
03695 
03696 
03697 /*
03698   ElementsOnShape
03699 */
03700 
03701 ElementsOnShape::ElementsOnShape()
03702   : //myMesh(0),
03703     myType(SMDSAbs_All),
03704     myToler(Precision::Confusion()),
03705     myAllNodesFlag(false)
03706 {
03707   myCurShapeType = TopAbs_SHAPE;
03708 }
03709 
03710 ElementsOnShape::~ElementsOnShape()
03711 {
03712 }
03713 
03714 void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
03715 {
03716   myMeshModifTracer.SetMesh( theMesh );
03717   if ( myMeshModifTracer.IsMeshModified())
03718     SetShape(myShape, myType);
03719 }
03720 
03721 bool ElementsOnShape::IsSatisfy (long theElementId)
03722 {
03723   return myIds.Contains(theElementId);
03724 }
03725 
03726 SMDSAbs_ElementType ElementsOnShape::GetType() const
03727 {
03728   return myType;
03729 }
03730 
03731 void ElementsOnShape::SetTolerance (const double theToler)
03732 {
03733   if (myToler != theToler) {
03734     myToler = theToler;
03735     SetShape(myShape, myType);
03736   }
03737 }
03738 
03739 double ElementsOnShape::GetTolerance() const
03740 {
03741   return myToler;
03742 }
03743 
03744 void ElementsOnShape::SetAllNodes (bool theAllNodes)
03745 {
03746   if (myAllNodesFlag != theAllNodes) {
03747     myAllNodesFlag = theAllNodes;
03748     SetShape(myShape, myType);
03749   }
03750 }
03751 
03752 void ElementsOnShape::SetShape (const TopoDS_Shape&       theShape,
03753                                 const SMDSAbs_ElementType theType)
03754 {
03755   myType = theType;
03756   myShape = theShape;
03757   myIds.Clear();
03758 
03759   const SMDS_Mesh* myMesh = myMeshModifTracer.GetMesh();
03760   
03761   if ( !myMesh ) return;
03762 
03763   switch (myType)
03764   {
03765   case SMDSAbs_All:
03766     myIds.ReSize(myMesh->NbEdges() + myMesh->NbFaces() + myMesh->NbVolumes());
03767     break;
03768   case SMDSAbs_Node:
03769     myIds.ReSize(myMesh->NbNodes());
03770     break;
03771   case SMDSAbs_Edge:
03772     myIds.ReSize(myMesh->NbEdges());
03773     break;
03774   case SMDSAbs_Face:
03775     myIds.ReSize(myMesh->NbFaces());
03776     break;
03777   case SMDSAbs_Volume:
03778     myIds.ReSize(myMesh->NbVolumes());
03779     break;
03780   default:
03781     break;
03782   }
03783 
03784   myShapesMap.Clear();
03785   addShape(myShape);
03786 }
03787 
03788 void ElementsOnShape::addShape (const TopoDS_Shape& theShape)
03789 {
03790   if (theShape.IsNull() || myMeshModifTracer.GetMesh() == 0)
03791     return;
03792 
03793   if (!myShapesMap.Add(theShape)) return;
03794 
03795   myCurShapeType = theShape.ShapeType();
03796   switch (myCurShapeType)
03797   {
03798   case TopAbs_COMPOUND:
03799   case TopAbs_COMPSOLID:
03800   case TopAbs_SHELL:
03801   case TopAbs_WIRE:
03802     {
03803       TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
03804       for (; anIt.More(); anIt.Next()) addShape(anIt.Value());
03805     }
03806     break;
03807   case TopAbs_SOLID:
03808     {
03809       myCurSC.Load(theShape);
03810       process();
03811     }
03812     break;
03813   case TopAbs_FACE:
03814     {
03815       TopoDS_Face aFace = TopoDS::Face(theShape);
03816       BRepAdaptor_Surface SA (aFace, true);
03817       Standard_Real
03818         u1 = SA.FirstUParameter(),
03819         u2 = SA.LastUParameter(),
03820         v1 = SA.FirstVParameter(),
03821         v2 = SA.LastVParameter();
03822       Handle(Geom_Surface) surf = BRep_Tool::Surface(aFace);
03823       myCurProjFace.Init(surf, u1,u2, v1,v2);
03824       myCurFace = aFace;
03825       process();
03826     }
03827     break;
03828   case TopAbs_EDGE:
03829     {
03830       TopoDS_Edge anEdge = TopoDS::Edge(theShape);
03831       Standard_Real u1, u2;
03832       Handle(Geom_Curve) curve = BRep_Tool::Curve(anEdge, u1, u2);
03833       myCurProjEdge.Init(curve, u1, u2);
03834       process();
03835     }
03836     break;
03837   case TopAbs_VERTEX:
03838     {
03839       TopoDS_Vertex aV = TopoDS::Vertex(theShape);
03840       myCurPnt = BRep_Tool::Pnt(aV);
03841       process();
03842     }
03843     break;
03844   default:
03845     break;
03846   }
03847 }
03848 
03849 void ElementsOnShape::process()
03850 {
03851   const SMDS_Mesh* myMesh = myMeshModifTracer.GetMesh();
03852   if (myShape.IsNull() || myMesh == 0)
03853     return;
03854 
03855   if (myType == SMDSAbs_Node)
03856   {
03857     SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator();
03858     while (anIter->more())
03859       process(anIter->next());
03860   }
03861   else
03862   {
03863     if (myType == SMDSAbs_Edge || myType == SMDSAbs_All)
03864     {
03865       SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator();
03866       while (anIter->more())
03867         process(anIter->next());
03868     }
03869 
03870     if (myType == SMDSAbs_Face || myType == SMDSAbs_All)
03871     {
03872       SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
03873       while (anIter->more()) {
03874         process(anIter->next());
03875       }
03876     }
03877 
03878     if (myType == SMDSAbs_Volume || myType == SMDSAbs_All)
03879     {
03880       SMDS_VolumeIteratorPtr anIter = myMesh->volumesIterator();
03881       while (anIter->more())
03882         process(anIter->next());
03883     }
03884   }
03885 }
03886 
03887 void ElementsOnShape::process (const SMDS_MeshElement* theElemPtr)
03888 {
03889   if (myShape.IsNull())
03890     return;
03891 
03892   SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
03893   bool isSatisfy = myAllNodesFlag;
03894 
03895   gp_XYZ centerXYZ (0, 0, 0);
03896 
03897   while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
03898   {
03899     SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
03900     gp_Pnt aPnt (aNode->X(), aNode->Y(), aNode->Z());
03901     centerXYZ += aPnt.XYZ();
03902 
03903     switch (myCurShapeType)
03904     {
03905     case TopAbs_SOLID:
03906       {
03907         myCurSC.Perform(aPnt, myToler);
03908         isSatisfy = (myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON);
03909       }
03910       break;
03911     case TopAbs_FACE:
03912       {
03913         myCurProjFace.Perform(aPnt);
03914         isSatisfy = (myCurProjFace.IsDone() && myCurProjFace.LowerDistance() <= myToler);
03915         if (isSatisfy)
03916         {
03917           // check relatively the face
03918           Quantity_Parameter u, v;
03919           myCurProjFace.LowerDistanceParameters(u, v);
03920           gp_Pnt2d aProjPnt (u, v);
03921           BRepClass_FaceClassifier aClsf (myCurFace, aProjPnt, myToler);
03922           isSatisfy = (aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON);
03923         }
03924       }
03925       break;
03926     case TopAbs_EDGE:
03927       {
03928         myCurProjEdge.Perform(aPnt);
03929         isSatisfy = (myCurProjEdge.NbPoints() > 0 && myCurProjEdge.LowerDistance() <= myToler);
03930       }
03931       break;
03932     case TopAbs_VERTEX:
03933       {
03934         isSatisfy = (aPnt.Distance(myCurPnt) <= myToler);
03935       }
03936       break;
03937     default:
03938       {
03939         isSatisfy = false;
03940       }
03941     }
03942   }
03943 
03944   if (isSatisfy && myCurShapeType == TopAbs_SOLID) { // Check the center point for volumes MantisBug 0020168
03945     centerXYZ /= theElemPtr->NbNodes();
03946     gp_Pnt aCenterPnt (centerXYZ);
03947     myCurSC.Perform(aCenterPnt, myToler);
03948     if ( !(myCurSC.State() == TopAbs_IN || myCurSC.State() == TopAbs_ON))
03949       isSatisfy = false;
03950   }
03951 
03952   if (isSatisfy)
03953     myIds.Add(theElemPtr->GetID());
03954 }
03955 
03956 TSequenceOfXYZ::TSequenceOfXYZ()
03957 {}
03958 
03959 TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n)
03960 {}
03961 
03962 TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t)
03963 {}
03964 
03965 TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray)
03966 {}
03967 
03968 template <class InputIterator>
03969 TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd)
03970 {}
03971 
03972 TSequenceOfXYZ::~TSequenceOfXYZ()
03973 {}
03974 
03975 TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
03976 {
03977   myArray = theSequenceOfXYZ.myArray;
03978   return *this;
03979 }
03980 
03981 gp_XYZ& TSequenceOfXYZ::operator()(size_type n)
03982 {
03983   return myArray[n-1];
03984 }
03985 
03986 const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const
03987 {
03988   return myArray[n-1];
03989 }
03990 
03991 void TSequenceOfXYZ::clear()
03992 {
03993   myArray.clear();
03994 }
03995 
03996 void TSequenceOfXYZ::reserve(size_type n)
03997 {
03998   myArray.reserve(n);
03999 }
04000 
04001 void TSequenceOfXYZ::push_back(const gp_XYZ& v)
04002 {
04003   myArray.push_back(v);
04004 }
04005 
04006 TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
04007 {
04008   return myArray.size();
04009 }
04010 
04011 TMeshModifTracer::TMeshModifTracer():
04012   myMeshModifTime(0), myMesh(0)
04013 {
04014 }
04015 void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh )
04016 {
04017   if ( theMesh != myMesh )
04018     myMeshModifTime = 0;
04019   myMesh = theMesh;
04020 }
04021 bool TMeshModifTracer::IsMeshModified()
04022 {
04023   bool modified = false;
04024   if ( myMesh )
04025   {
04026     modified = ( myMeshModifTime != myMesh->GetMTime() );
04027     myMeshModifTime = myMesh->GetMTime();
04028   }
04029   return modified;
04030 }