1 // Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
2 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License.
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 // See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
20 #include "SMESH_Controls.hxx"
24 #include <BRep_Tool.hxx>
26 #include <gp_Cylinder.hxx>
32 #include <Geom_Plane.hxx>
33 #include <Geom_CylindricalSurface.hxx>
34 #include <Precision.hxx>
35 #include <TColgp_Array1OfXYZ.hxx>
36 #include <TColgp_SequenceOfXYZ.hxx>
37 #include <TColStd_MapOfInteger.hxx>
38 #include <TColStd_SequenceOfAsciiString.hxx>
39 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
42 #include <TopoDS_Face.hxx>
43 #include <TopoDS_Shape.hxx>
45 #include "SMDS_Mesh.hxx"
46 #include "SMDS_Iterator.hxx"
47 #include "SMDS_MeshElement.hxx"
48 #include "SMDS_MeshNode.hxx"
56 static inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
58 gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
60 return v1.Magnitude() < gp::Resolution() ||
61 v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
64 static inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
66 gp_Vec aVec1( P2 - P1 );
67 gp_Vec aVec2( P3 - P1 );
68 return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
71 static inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
73 return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
76 static inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
78 double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
82 static int getNbMultiConnection( SMDS_Mesh* theMesh, const int theId )
87 const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
88 if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge || anEdge->NbNodes() != 2 )
91 TColStd_MapOfInteger aMap;
94 SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
97 while( anIter->more() )
99 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
102 SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
103 while( anElemIter->more() )
105 const SMDS_MeshElement* anElem = anElemIter->next();
106 if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge )
108 int anId = anElem->GetID();
110 if ( anIter->more() ) // i.e. first node
112 else if ( aMap.Contains( anId ) )
123 using namespace SMESH::Controls;
130 Class : NumericalFunctor
131 Description : Base class for numerical functors
133 NumericalFunctor::NumericalFunctor():
139 void NumericalFunctor::SetMesh( SMDS_Mesh* theMesh )
144 bool NumericalFunctor::GetPoints(const int theId,
145 TColgp_SequenceOfXYZ& theRes ) const
152 return GetPoints( myMesh->FindElement( theId ), theRes );
155 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
156 TColgp_SequenceOfXYZ& theRes )
163 // Get nodes of the element
164 SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
167 while( anIter->more() )
169 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
171 theRes.Append( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
178 long NumericalFunctor::GetPrecision() const
183 void NumericalFunctor::SetPrecision( const long thePrecision )
185 myPrecision = thePrecision;
188 double NumericalFunctor::GetValue( long theId )
190 TColgp_SequenceOfXYZ P;
191 if ( GetPoints( theId, P ))
193 double aVal = GetValue( P );
194 if ( myPrecision >= 0 )
196 double prec = pow( 10., (double)( myPrecision ) );
197 aVal = floor( aVal * prec + 0.5 ) / prec;
208 Description : Functor for calculation of minimum angle
211 double MinimumAngle::GetValue( const TColgp_SequenceOfXYZ& P )
215 if ( P.Length() == 3 )
217 double A0 = getAngle( P( 3 ), P( 1 ), P( 2 ) );
218 double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
219 double A2 = getAngle( P( 2 ), P( 3 ), P( 1 ) );
221 aMin = Min( A0, Min( A1, A2 ) );
223 else if ( P.Length() == 4 )
225 double A0 = getAngle( P( 4 ), P( 1 ), P( 2 ) );
226 double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
227 double A2 = getAngle( P( 2 ), P( 3 ), P( 4 ) );
228 double A3 = getAngle( P( 3 ), P( 4 ), P( 1 ) );
230 aMin = Min( Min( A0, A1 ), Min( A2, A3 ) );
235 return aMin * 180 / PI;
238 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
240 const double aBestAngle = PI / nbNodes;
241 return ( fabs( aBestAngle - Value ));
244 SMDSAbs_ElementType MinimumAngle::GetType() const
252 Description : Functor for calculating aspect ratio
254 double AspectRatio::GetValue( const TColgp_SequenceOfXYZ& P )
256 int nbNodes = P.Length();
258 if ( nbNodes != 3 && nbNodes != 4 )
261 // Compute lengths of the sides
263 double aLen[ nbNodes ];
264 for ( int i = 0; i < nbNodes - 1; i++ )
265 aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
266 aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
268 // Compute aspect ratio
272 double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
273 if ( anArea <= Precision::Confusion() )
275 double aMaxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
276 static double aCoef = sqrt( 3. ) / 4;
278 return aCoef * aMaxLen * aMaxLen / anArea;
282 double aMinLen = Min( Min( aLen[ 0 ], aLen[ 1 ] ), Min( aLen[ 2 ], aLen[ 3 ] ) );
283 if ( aMinLen <= Precision::Confusion() )
285 double aMaxLen = Max( Max( aLen[ 0 ], aLen[ 1 ] ), Max( aLen[ 2 ], aLen[ 3 ] ) );
287 return aMaxLen / aMinLen;
291 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
293 // the aspect ratio is in the range [1.0,infinity]
296 return Value / 1000.;
299 SMDSAbs_ElementType AspectRatio::GetType() const
306 Class : AspectRatio3D
307 Description : Functor for calculating aspect ratio
310 static inline double getHalfPerimeter(double theTria[3]){
311 return (theTria[0] + theTria[1] + theTria[2])/2.0;
314 static inline double getArea(double theHalfPerim, double theTria[3]){
315 return sqrt(theHalfPerim*
316 (theHalfPerim-theTria[0])*
317 (theHalfPerim-theTria[1])*
318 (theHalfPerim-theTria[2]));
321 static inline double getVolume(double theLen[6]){
322 double a2 = theLen[0]*theLen[0];
323 double b2 = theLen[1]*theLen[1];
324 double c2 = theLen[2]*theLen[2];
325 double d2 = theLen[3]*theLen[3];
326 double e2 = theLen[4]*theLen[4];
327 double f2 = theLen[5]*theLen[5];
328 double P = 4.0*a2*b2*d2;
329 double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
330 double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
331 return sqrt(P-Q+R)/12.0;
334 static inline double getHeight( const gp_Pnt& P1, const gp_Pnt& P2,
335 const gp_Pnt& P3, const gp_Pnt& P4)
337 gp_Vec aVec1( P2.XYZ() - P1.XYZ() );
338 gp_Vec aVec2( P3.XYZ() - P1.XYZ() );
339 gp_Vec aNorm = aVec1 ^ aVec2;
340 aNorm /= aNorm.Magnitude();
341 gp_Vec aVec3( P4.XYZ() - P1.XYZ() );
342 double aDist = aVec1 * aVec2;
343 return fabs( aDist );
346 static inline double getMaxHeight( const TColgp_SequenceOfXYZ& P )
348 double aHeight = getHeight(P(1),P(2),P(3),P(4));
349 aHeight = max(aHeight,getHeight(P(1),P(2),P(4),P(3)));
350 aHeight = max(aHeight,getHeight(P(1),P(3),P(4),P(2)));
351 aHeight = max(aHeight,getHeight(P(2),P(3),P(4),P(1)));
355 double AspectRatio3D::GetValue( const TColgp_SequenceOfXYZ& P )
357 double aQuality = 0.0;
358 int nbNodes = P.Length();
362 getDistance(P(1),P(2)), // a
363 getDistance(P(2),P(3)), // b
364 getDistance(P(3),P(1)), // c
365 getDistance(P(2),P(4)), // d
366 getDistance(P(3),P(4)), // e
367 getDistance(P(1),P(4)) // f
369 double aTria[4][3] = {
370 {aLen[0],aLen[1],aLen[2]}, // abc
371 {aLen[0],aLen[3],aLen[5]}, // adf
372 {aLen[1],aLen[3],aLen[4]}, // bde
373 {aLen[2],aLen[4],aLen[5]} // cef
375 double aHalfPerim = getHalfPerimeter(aTria[0]);
376 double anArea = getArea(aHalfPerim,aTria[0]);
377 aHalfPerim = getHalfPerimeter(aTria[1]);
378 anArea += getArea(aHalfPerim,aTria[1]);
379 aHalfPerim = getHalfPerimeter(aTria[2]);
380 anArea += getArea(aHalfPerim,aTria[2]);
381 double aVolume = getVolume(aLen);
382 double aHeight = getMaxHeight(P);
383 aQuality = 1.0/3.0*aHeight*anArea/aVolume;
390 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
392 // the aspect ratio is in the range [1.0,infinity]
395 return Value / 1000.;
398 SMDSAbs_ElementType AspectRatio3D::GetType() const
400 return SMDSAbs_Volume;
406 Description : Functor for calculating warping
408 double Warping::GetValue( const TColgp_SequenceOfXYZ& P )
410 if ( P.Length() != 4 )
413 gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4;
415 double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
416 double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
417 double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
418 double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
420 return Max( Max( A1, A2 ), Max( A3, A4 ) );
423 double Warping::ComputeA( const gp_XYZ& thePnt1,
424 const gp_XYZ& thePnt2,
425 const gp_XYZ& thePnt3,
426 const gp_XYZ& theG ) const
428 double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
429 double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
430 double L = Min( aLen1, aLen2 ) * 0.5;
431 if ( L < Precision::Confusion())
434 gp_XYZ GI = ( thePnt2 - thePnt1 ) / 2. - theG;
435 gp_XYZ GJ = ( thePnt3 - thePnt2 ) / 2. - theG;
436 gp_XYZ N = GI.Crossed( GJ );
438 if ( N.Modulus() < gp::Resolution() )
443 double H = ( thePnt2 - theG ).Dot( N );
444 return asin( fabs( H / L ) ) * 180 / PI;
447 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
449 // the warp is in the range [0.0,PI/2]
450 // 0.0 = good (no warp)
451 // PI/2 = bad (face pliee)
455 SMDSAbs_ElementType Warping::GetType() const
463 Description : Functor for calculating taper
465 double Taper::GetValue( const TColgp_SequenceOfXYZ& P )
467 if ( P.Length() != 4 )
471 double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2;
472 double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2;
473 double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2;
474 double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2;
476 double JA = 0.25 * ( J1 + J2 + J3 + J4 );
477 if ( JA <= Precision::Confusion() )
480 double T1 = fabs( ( J1 - JA ) / JA );
481 double T2 = fabs( ( J2 - JA ) / JA );
482 double T3 = fabs( ( J3 - JA ) / JA );
483 double T4 = fabs( ( J4 - JA ) / JA );
485 return Max( Max( T1, T2 ), Max( T3, T4 ) );
488 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
490 // the taper is in the range [0.0,1.0]
491 // 0.0 = good (no taper)
492 // 1.0 = bad (les cotes opposes sont allignes)
496 SMDSAbs_ElementType Taper::GetType() const
504 Description : Functor for calculating skew in degrees
506 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
508 gp_XYZ p12 = ( p2 + p1 ) / 2;
509 gp_XYZ p23 = ( p3 + p2 ) / 2;
510 gp_XYZ p31 = ( p3 + p1 ) / 2;
512 gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
514 return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
517 double Skew::GetValue( const TColgp_SequenceOfXYZ& P )
519 if ( P.Length() != 3 && P.Length() != 4 )
523 static double PI2 = PI / 2;
524 if ( P.Length() == 3 )
526 double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
527 double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
528 double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
530 return Max( A0, Max( A1, A2 ) ) * 180 / PI;
534 gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2;
535 gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2;
536 gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2;
537 gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2;
539 gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
540 double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
541 ? 0 : fabs( PI2 - v1.Angle( v2 ) );
547 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
549 // the skew is in the range [0.0,PI/2].
555 SMDSAbs_ElementType Skew::GetType() const
563 Description : Functor for calculating area
565 double Area::GetValue( const TColgp_SequenceOfXYZ& P )
567 if ( P.Length() == 3 )
568 return getArea( P( 1 ), P( 2 ), P( 3 ) );
569 else if ( P.Length() == 4 )
570 return getArea( P( 1 ), P( 2 ), P( 3 ) ) + getArea( P( 1 ), P( 3 ), P( 4 ) );
575 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
580 SMDSAbs_ElementType Area::GetType() const
588 Description : Functor for calculating length off edge
590 double Length::GetValue( const TColgp_SequenceOfXYZ& P )
592 return ( P.Length() == 2 ? getDistance( P( 1 ), P( 2 ) ) : 0 );
595 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
600 SMDSAbs_ElementType Length::GetType() const
607 Class : MultiConnection
608 Description : Functor for calculating number of faces conneted to the edge
610 double MultiConnection::GetValue( const TColgp_SequenceOfXYZ& P )
614 double MultiConnection::GetValue( long theId )
616 return getNbMultiConnection( myMesh, theId );
619 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
624 SMDSAbs_ElementType MultiConnection::GetType() const
636 Description : Predicate for free borders
639 FreeBorders::FreeBorders()
644 void FreeBorders::SetMesh( SMDS_Mesh* theMesh )
649 bool FreeBorders::IsSatisfy( long theId )
651 return getNbMultiConnection( myMesh, theId ) == 1;
654 SMDSAbs_ElementType FreeBorders::GetType() const
662 Description : Predicate for free Edges
664 FreeEdges::FreeEdges()
669 void FreeEdges::SetMesh( SMDS_Mesh* theMesh )
674 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId )
676 TColStd_MapOfInteger aMap;
677 for ( int i = 0; i < 2; i++ )
679 SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator();
680 while( anElemIter->more() )
682 const SMDS_MeshElement* anElem = anElemIter->next();
683 if ( anElem != 0 && anElem->GetType() == SMDSAbs_Face )
685 int anId = anElem->GetID();
689 else if ( aMap.Contains( anId ) && anId != theFaceId )
697 bool FreeEdges::IsSatisfy( long theId )
702 const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
703 if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
706 int nbNodes = aFace->NbNodes();
707 const SMDS_MeshNode* aNodes[ nbNodes ];
709 SMDS_ElemIteratorPtr anIter = aFace->nodesIterator();
712 while( anIter->more() )
714 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
717 aNodes[ i++ ] = aNode;
721 for ( int i = 0; i < nbNodes - 1; i++ )
722 if ( IsFreeEdge( &aNodes[ i ], theId ) )
725 aNodes[ 1 ] = aNodes[ nbNodes - 1 ];
727 return IsFreeEdge( &aNodes[ 0 ], theId );
731 SMDSAbs_ElementType FreeEdges::GetType() const
736 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
739 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
740 if(thePntId1 > thePntId2){
741 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
745 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
746 if(myPntId[0] < x.myPntId[0]) return true;
747 if(myPntId[0] == x.myPntId[0])
748 if(myPntId[1] < x.myPntId[1]) return true;
752 inline void UpdateBorders(const FreeEdges::Border& theBorder,
753 FreeEdges::TBorders& theRegistry,
754 FreeEdges::TBorders& theContainer)
756 if(theRegistry.find(theBorder) == theRegistry.end()){
757 theRegistry.insert(theBorder);
758 theContainer.insert(theBorder);
760 theContainer.erase(theBorder);
764 void FreeEdges::GetBoreders(TBorders& theBorders)
767 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
768 for(; anIter->more(); ){
769 const SMDS_MeshFace* anElem = anIter->next();
770 long anElemId = anElem->GetID();
771 SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
773 const SMDS_MeshElement* aNode;
774 if(aNodesIter->more()){
775 aNode = aNodesIter->next();
776 aNodeId[0] = aNodeId[1] = aNode->GetID();
778 for(; aNodesIter->more(); ){
779 aNode = aNodesIter->next();
780 long anId = aNode->GetID();
781 Border aBorder(anElemId,aNodeId[1],anId);
783 //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
784 UpdateBorders(aBorder,aRegistry,theBorders);
786 Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
787 //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
788 UpdateBorders(aBorder,aRegistry,theBorders);
790 //std::cout<<"theBorders.size() = "<<theBorders.size()<<endl;
795 Description : Predicate for Range of Ids.
796 Range may be specified with two ways.
797 1. Using AddToRange method
798 2. With SetRangeStr method. Parameter of this method is a string
799 like as "1,2,3,50-60,63,67,70-"
802 //=======================================================================
804 // Purpose : Constructor
805 //=======================================================================
806 RangeOfIds::RangeOfIds()
809 myType = SMDSAbs_All;
812 //=======================================================================
814 // Purpose : Set mesh
815 //=======================================================================
816 void RangeOfIds::SetMesh( SMDS_Mesh* theMesh )
821 //=======================================================================
823 // Purpose : Add ID to the range
824 //=======================================================================
825 bool RangeOfIds::AddToRange( long theEntityId )
827 myIds.Add( theEntityId );
831 //=======================================================================
832 // name : GetRangeStr
833 // Purpose : Get range as a string.
834 // Example: "1,2,3,50-60,63,67,70-"
835 //=======================================================================
836 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
840 TColStd_SequenceOfInteger anIntSeq;
841 TColStd_SequenceOfAsciiString aStrSeq;
843 TColStd_MapIteratorOfMapOfInteger anIter( myIds );
844 for ( ; anIter.More(); anIter.Next() )
846 int anId = anIter.Key();
847 TCollection_AsciiString aStr( anId );
848 anIntSeq.Append( anId );
849 aStrSeq.Append( aStr );
852 for ( int i = 1, n = myMin.Length(); i <= n; i++ )
854 int aMinId = myMin( i );
855 int aMaxId = myMax( i );
857 TCollection_AsciiString aStr;
858 if ( aMinId != IntegerFirst() )
863 if ( aMaxId != IntegerLast() )
866 // find position of the string in result sequence and insert string in it
867 if ( anIntSeq.Length() == 0 )
869 anIntSeq.Append( aMinId );
870 aStrSeq.Append( aStr );
874 if ( aMinId < anIntSeq.First() )
876 anIntSeq.Prepend( aMinId );
877 aStrSeq.Prepend( aStr );
879 else if ( aMinId > anIntSeq.Last() )
881 anIntSeq.Append( aMinId );
882 aStrSeq.Append( aStr );
885 for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
886 if ( aMinId < anIntSeq( j ) )
888 anIntSeq.InsertBefore( j, aMinId );
889 aStrSeq.InsertBefore( j, aStr );
895 if ( aStrSeq.Length() == 0 )
898 theResStr = aStrSeq( 1 );
899 for ( int j = 2, k = aStrSeq.Length(); j <= k; j++ )
902 theResStr += aStrSeq( j );
906 //=======================================================================
907 // name : SetRangeStr
908 // Purpose : Define range with string
909 // Example of entry string: "1,2,3,50-60,63,67,70-"
910 //=======================================================================
911 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
917 TCollection_AsciiString aStr = theStr;
918 aStr.RemoveAll( ' ' );
919 aStr.RemoveAll( '\t' );
921 for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) )
922 aStr.Remove( aPos, 2 );
924 TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
926 while ( tmpStr != "" )
928 tmpStr = aStr.Token( ",", i++ );
929 int aPos = tmpStr.Search( '-' );
933 if ( tmpStr.IsIntegerValue() )
934 myIds.Add( tmpStr.IntegerValue() );
940 TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
941 TCollection_AsciiString aMinStr = tmpStr;
943 while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
944 while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
946 if ( !aMinStr.IsEmpty() && !aMinStr.IsIntegerValue() ||
947 !aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue() )
950 myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
951 myMax.Append( aMaxStr.IsEmpty() ? IntegerLast() : aMaxStr.IntegerValue() );
958 //=======================================================================
960 // Purpose : Get type of supported entities
961 //=======================================================================
962 SMDSAbs_ElementType RangeOfIds::GetType() const
967 //=======================================================================
969 // Purpose : Set type of supported entities
970 //=======================================================================
971 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
976 //=======================================================================
978 // Purpose : Verify whether entity satisfies to this rpedicate
979 //=======================================================================
980 bool RangeOfIds::IsSatisfy( long theId )
985 if ( myType == SMDSAbs_Node )
987 if ( myMesh->FindNode( theId ) == 0 )
992 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
993 if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All )
997 if ( myIds.Contains( theId ) )
1000 for ( int i = 1, n = myMin.Length(); i <= n; i++ )
1001 if ( theId >= myMin( i ) && theId <= myMax( i ) )
1009 Description : Base class for comparators
1011 Comparator::Comparator():
1015 Comparator::~Comparator()
1018 void Comparator::SetMesh( SMDS_Mesh* theMesh )
1021 myFunctor->SetMesh( theMesh );
1024 void Comparator::SetMargin( double theValue )
1026 myMargin = theValue;
1029 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
1031 myFunctor = theFunct;
1034 SMDSAbs_ElementType Comparator::GetType() const
1036 return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
1039 double Comparator::GetMargin()
1047 Description : Comparator "<"
1049 bool LessThan::IsSatisfy( long theId )
1051 return myFunctor && myFunctor->GetValue( theId ) < myMargin;
1057 Description : Comparator ">"
1059 bool MoreThan::IsSatisfy( long theId )
1061 return myFunctor && myFunctor->GetValue( theId ) > myMargin;
1067 Description : Comparator "="
1070 myToler(Precision::Confusion())
1073 bool EqualTo::IsSatisfy( long theId )
1075 return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
1078 void EqualTo::SetTolerance( double theToler )
1083 double EqualTo::GetTolerance()
1090 Description : Logical NOT predicate
1092 LogicalNOT::LogicalNOT()
1095 LogicalNOT::~LogicalNOT()
1098 bool LogicalNOT::IsSatisfy( long theId )
1100 return myPredicate && !myPredicate->IsSatisfy( theId );
1103 void LogicalNOT::SetMesh( SMDS_Mesh* theMesh )
1106 myPredicate->SetMesh( theMesh );
1109 void LogicalNOT::SetPredicate( PredicatePtr thePred )
1111 myPredicate = thePred;
1114 SMDSAbs_ElementType LogicalNOT::GetType() const
1116 return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
1121 Class : LogicalBinary
1122 Description : Base class for binary logical predicate
1124 LogicalBinary::LogicalBinary()
1127 LogicalBinary::~LogicalBinary()
1130 void LogicalBinary::SetMesh( SMDS_Mesh* theMesh )
1133 myPredicate1->SetMesh( theMesh );
1136 myPredicate2->SetMesh( theMesh );
1139 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
1141 myPredicate1 = thePredicate;
1144 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
1146 myPredicate2 = thePredicate;
1149 SMDSAbs_ElementType LogicalBinary::GetType() const
1151 if ( !myPredicate1 || !myPredicate2 )
1154 SMDSAbs_ElementType aType1 = myPredicate1->GetType();
1155 SMDSAbs_ElementType aType2 = myPredicate2->GetType();
1157 return aType1 == aType2 ? aType1 : SMDSAbs_All;
1163 Description : Logical AND
1165 bool LogicalAND::IsSatisfy( long theId )
1170 myPredicate1->IsSatisfy( theId ) &&
1171 myPredicate2->IsSatisfy( theId );
1177 Description : Logical OR
1179 bool LogicalOR::IsSatisfy( long theId )
1184 myPredicate1->IsSatisfy( theId ) ||
1185 myPredicate2->IsSatisfy( theId );
1199 void Filter::SetPredicate( PredicatePtr thePredicate )
1201 myPredicate = thePredicate;
1205 template<class TElement, class TIterator, class TPredicate>
1206 void FillSequence(const TIterator& theIterator,
1207 TPredicate& thePredicate,
1208 Filter::TIdSequence& theSequence)
1210 if ( theIterator ) {
1211 while( theIterator->more() ) {
1212 TElement anElem = theIterator->next();
1213 long anId = anElem->GetID();
1214 if ( thePredicate->IsSatisfy( anId ) )
1215 theSequence.push_back( anId );
1221 Filter::GetElementsId( SMDS_Mesh* theMesh )
1223 TIdSequence aSequence;
1224 if ( !theMesh || !myPredicate ) return aSequence;
1226 myPredicate->SetMesh( theMesh );
1228 SMDSAbs_ElementType aType = myPredicate->GetType();
1231 FillSequence<const SMDS_MeshNode*>(theMesh->nodesIterator(),myPredicate,aSequence);
1235 FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
1239 FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
1242 case SMDSAbs_Volume:{
1243 FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
1247 FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
1248 FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
1249 FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
1260 typedef std::set<SMDS_MeshFace*> TMapOfFacePtr;
1266 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
1267 SMDS_MeshNode* theNode2 )
1273 ManifoldPart::Link::~Link()
1279 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
1281 if ( myNode1 == theLink.myNode1 &&
1282 myNode2 == theLink.myNode2 )
1284 else if ( myNode1 == theLink.myNode2 &&
1285 myNode2 == theLink.myNode1 )
1291 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
1293 if(myNode1 < x.myNode1) return true;
1294 if(myNode1 == x.myNode1)
1295 if(myNode2 < x.myNode2) return true;
1299 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
1300 const ManifoldPart::Link& theLink2 )
1302 return theLink1.IsEqual( theLink2 );
1305 ManifoldPart::ManifoldPart()
1308 myAngToler = Precision::Angular();
1309 myIsOnlyManifold = true;
1312 ManifoldPart::~ManifoldPart()
1317 void ManifoldPart::SetMesh( SMDS_Mesh* theMesh )
1323 SMDSAbs_ElementType ManifoldPart::GetType() const
1324 { return SMDSAbs_Face; }
1326 bool ManifoldPart::IsSatisfy( long theElementId )
1328 return myMapIds.Contains( theElementId );
1331 void ManifoldPart::SetAngleTolerance( const double theAngToler )
1332 { myAngToler = theAngToler; }
1334 double ManifoldPart::GetAngleTolerance() const
1335 { return myAngToler; }
1337 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
1338 { myIsOnlyManifold = theIsOnly; }
1340 void ManifoldPart::SetStartElem( const long theStartId )
1341 { myStartElemId = theStartId; }
1343 bool ManifoldPart::process()
1346 myMapBadGeomIds.Clear();
1348 myAllFacePtr.clear();
1349 myAllFacePtrIntDMap.clear();
1353 // collect all faces into own map
1354 SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
1355 for (; anFaceItr->more(); )
1357 SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
1358 myAllFacePtr.push_back( aFacePtr );
1359 myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
1362 SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
1366 // the map of non manifold links and bad geometry
1367 TMapOfLink aMapOfNonManifold;
1368 TColStd_MapOfInteger aMapOfTreated;
1370 // begin cycle on faces from start index and run on vector till the end
1371 // and from begin to start index to cover whole vector
1372 const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
1373 bool isStartTreat = false;
1374 for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
1376 if ( fi == aStartIndx )
1377 isStartTreat = true;
1378 // as result next time when fi will be equal to aStartIndx
1380 SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
1381 if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
1384 aMapOfTreated.Add( aFacePtr->GetID() );
1385 TColStd_MapOfInteger aResFaces;
1386 if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
1387 aMapOfNonManifold, aResFaces ) )
1389 TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
1390 for ( ; anItr.More(); anItr.Next() )
1392 int aFaceId = anItr.Key();
1393 aMapOfTreated.Add( aFaceId );
1394 myMapIds.Add( aFaceId );
1397 if ( fi == ( myAllFacePtr.size() - 1 ) )
1399 } // end run on vector of faces
1400 return !myMapIds.IsEmpty();
1403 static void getLinks( const SMDS_MeshFace* theFace,
1404 ManifoldPart::TVectorOfLink& theLinks )
1406 int aNbNode = theFace->NbNodes();
1407 SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
1409 SMDS_MeshNode* aNode = 0;
1410 for ( ; aNodeItr->more() && i <= aNbNode; )
1413 SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
1417 SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
1419 ManifoldPart::Link aLink( aN1, aN2 );
1420 theLinks.push_back( aLink );
1424 static gp_XYZ getNormale( const SMDS_MeshFace* theFace )
1427 int aNbNode = theFace->NbNodes();
1428 TColgp_Array1OfXYZ anArrOfXYZ(1,4);
1429 gp_XYZ p1, p2, p3, p4;
1430 SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
1432 for ( ; aNodeItr->more() && i <= 4; i++ )
1434 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
1435 anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
1438 gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1);
1439 gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1);
1443 gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1);
1446 double len = n.Modulus();
1453 bool ManifoldPart::findConnected
1454 ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
1455 SMDS_MeshFace* theStartFace,
1456 ManifoldPart::TMapOfLink& theNonManifold,
1457 TColStd_MapOfInteger& theResFaces )
1459 theResFaces.Clear();
1460 if ( !theAllFacePtrInt.size() )
1463 if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
1465 myMapBadGeomIds.Add( theStartFace->GetID() );
1469 ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
1470 ManifoldPart::TVectorOfLink aSeqOfBoundary;
1471 theResFaces.Add( theStartFace->GetID() );
1472 ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
1474 expandBoundary( aMapOfBoundary, aSeqOfBoundary,
1475 aDMapLinkFace, theNonManifold, theStartFace );
1477 bool isDone = false;
1478 while ( !isDone && aMapOfBoundary.size() != 0 )
1480 bool isToReset = false;
1481 ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
1482 for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
1484 ManifoldPart::Link aLink = *pLink;
1485 if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
1487 // each link could be treated only once
1488 aMapToSkip.insert( aLink );
1490 ManifoldPart::TVectorOfFacePtr aFaces;
1492 if ( myIsOnlyManifold &&
1493 (theNonManifold.find( aLink ) != theNonManifold.end()) )
1497 getFacesByLink( aLink, aFaces );
1498 // filter the element to keep only indicated elements
1499 ManifoldPart::TVectorOfFacePtr aFiltered;
1500 ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
1501 for ( ; pFace != aFaces.end(); ++pFace )
1503 SMDS_MeshFace* aFace = *pFace;
1504 if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
1505 aFiltered.push_back( aFace );
1508 if ( aFaces.size() < 2 ) // no neihgbour faces
1510 else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
1512 theNonManifold.insert( aLink );
1517 // compare normal with normals of neighbor element
1518 SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
1519 ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
1520 for ( ; pFace != aFaces.end(); ++pFace )
1522 SMDS_MeshFace* aNextFace = *pFace;
1523 if ( aPrevFace == aNextFace )
1525 int anNextFaceID = aNextFace->GetID();
1526 if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
1527 // should not be with non manifold restriction. probably bad topology
1529 // check if face was treated and skipped
1530 if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
1531 !isInPlane( aPrevFace, aNextFace ) )
1533 // add new element to connected and extend the boundaries.
1534 theResFaces.Add( anNextFaceID );
1535 expandBoundary( aMapOfBoundary, aSeqOfBoundary,
1536 aDMapLinkFace, theNonManifold, aNextFace );
1540 isDone = !isToReset;
1543 return !theResFaces.IsEmpty();
1546 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
1547 const SMDS_MeshFace* theFace2 )
1549 gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
1550 gp_XYZ aNorm2XYZ = getNormale( theFace2 );
1551 if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
1553 myMapBadGeomIds.Add( theFace2->GetID() );
1556 if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
1562 void ManifoldPart::expandBoundary
1563 ( ManifoldPart::TMapOfLink& theMapOfBoundary,
1564 ManifoldPart::TVectorOfLink& theSeqOfBoundary,
1565 ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
1566 ManifoldPart::TMapOfLink& theNonManifold,
1567 SMDS_MeshFace* theNextFace ) const
1569 ManifoldPart::TVectorOfLink aLinks;
1570 getLinks( theNextFace, aLinks );
1571 int aNbLink = aLinks.size();
1572 for ( int i = 0; i < aNbLink; i++ )
1574 ManifoldPart::Link aLink = aLinks[ i ];
1575 if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
1577 if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
1579 if ( myIsOnlyManifold )
1581 // remove from boundary
1582 theMapOfBoundary.erase( aLink );
1583 ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
1584 for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
1586 ManifoldPart::Link aBoundLink = *pLink;
1587 if ( aBoundLink.IsEqual( aLink ) )
1589 theSeqOfBoundary.erase( pLink );
1597 theMapOfBoundary.insert( aLink );
1598 theSeqOfBoundary.push_back( aLink );
1599 theDMapLinkFacePtr[ aLink ] = theNextFace;
1604 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
1605 ManifoldPart::TVectorOfFacePtr& theFaces ) const
1607 SMDS_Mesh::SetOfFaces aSetOfFaces;
1608 // take all faces that shared first node
1609 SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
1610 for ( ; anItr->more(); )
1612 SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
1615 aSetOfFaces.insert( aFace );
1617 // take all faces that shared second node
1618 anItr = theLink.myNode2->facesIterator();
1619 // find the common part of two sets
1620 for ( ; anItr->more(); )
1622 SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
1623 if ( aSetOfFaces.find( aFace ) != aSetOfFaces.end() )
1624 theFaces.push_back( aFace );
1633 ElementsOnSurface::ElementsOnSurface()
1637 myType = SMDSAbs_All;
1639 myToler = Precision::Confusion();
1642 ElementsOnSurface::~ElementsOnSurface()
1647 void ElementsOnSurface::SetMesh( SMDS_Mesh* theMesh )
1649 if ( myMesh == theMesh )
1656 bool ElementsOnSurface::IsSatisfy( long theElementId )
1658 return myIds.Contains( theElementId );
1661 SMDSAbs_ElementType ElementsOnSurface::GetType() const
1664 void ElementsOnSurface::SetTolerance( const double theToler )
1665 { myToler = theToler; }
1667 double ElementsOnSurface::GetTolerance() const
1672 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
1673 const SMDSAbs_ElementType theType )
1677 if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
1682 TopoDS_Face aFace = TopoDS::Face( theShape );
1683 mySurf = BRep_Tool::Surface( aFace );
1686 void ElementsOnSurface::process()
1689 if ( mySurf.IsNull() )
1695 if ( myType == SMDSAbs_Face || myType == SMDSAbs_All )
1697 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1698 for(; anIter->more(); )
1699 process( anIter->next() );
1702 if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All )
1704 SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator();
1705 for(; anIter->more(); )
1706 process( anIter->next() );
1709 if ( myType == SMDSAbs_Node )
1711 SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator();
1712 for(; anIter->more(); )
1713 process( anIter->next() );
1717 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
1719 SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
1720 bool isSatisfy = true;
1721 for ( ; aNodeItr->more(); )
1723 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
1724 if ( !isOnSurface( aNode ) )
1731 myIds.Add( theElemPtr->GetID() );
1734 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) const
1736 if ( mySurf.IsNull() )
1739 gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
1740 double aToler2 = myToler * myToler;
1741 if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
1743 gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
1744 if ( aPln.SquareDistance( aPnt ) > aToler2 )
1747 else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
1749 gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
1750 double aRad = aCyl.Radius();
1751 gp_Ax3 anAxis = aCyl.Position();
1752 gp_XYZ aLoc = aCyl.Location().XYZ();
1753 double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
1754 double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
1755 if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )