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
25 #include <Precision.hxx>
26 #include <TColgp_SequenceOfXYZ.hxx>
27 #include <TColStd_MapOfInteger.hxx>
29 #include "SMDS_Mesh.hxx"
30 #include "SMDS_Iterator.hxx"
31 #include "SMDS_MeshElement.hxx"
32 #include "SMDS_MeshNode.hxx"
34 #include "SMESH_Controls.hxx"
40 static inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
42 return gp_Vec( P1 - P2 ).Angle( gp_Vec( P3 - P2 ) );
45 static inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
47 gp_Vec aVec1( P2 - P1 );
48 gp_Vec aVec2( P3 - P1 );
49 return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
52 static inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
54 return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
57 static inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
59 double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
63 static int getNbMultiConnection( SMDS_Mesh* theMesh, const int theId )
68 const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
69 if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge || anEdge->NbNodes() != 2 )
72 TColStd_MapOfInteger aMap;
75 SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
78 while( anIter->more() )
80 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
83 SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
84 while( anElemIter->more() )
86 const SMDS_MeshElement* anElem = anElemIter->next();
87 if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge )
89 int anId = anElem->GetID();
91 if ( anIter->more() ) // i.e. first node
93 else if ( aMap.Contains( anId ) )
104 using namespace SMESH::Controls;
111 Class : NumericalFunctor
112 Description : Base class for numerical functors
114 NumericalFunctor::NumericalFunctor():
118 void NumericalFunctor::SetMesh( SMDS_Mesh* theMesh )
123 bool NumericalFunctor::getPoints( const int theId,
124 TColgp_SequenceOfXYZ& theRes ) const
131 // Get nodes of the face
132 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
133 if ( anElem == 0 || anElem->GetType() != GetType() )
136 int nbNodes = anElem->NbNodes();
138 SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
141 while( anIter->more() )
143 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
145 theRes.Append( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
155 Description : Functor for calculation of minimum angle
157 double MinimumAngle::GetValue( long theId )
159 TColgp_SequenceOfXYZ P;
160 if ( !getPoints( theId, P ) || P.Length() != 3 && P.Length() != 4 )
165 if ( P.Length() == 3 )
167 double A0 = getAngle( P( 3 ), P( 1 ), P( 2 ) );
168 double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
169 double A2 = getAngle( P( 2 ), P( 3 ), P( 1 ) );
171 aMin = Min( A0, Min( A1, A2 ) );
175 double A0 = getAngle( P( 4 ), P( 1 ), P( 2 ) );
176 double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
177 double A2 = getAngle( P( 2 ), P( 3 ), P( 4 ) );
178 double A3 = getAngle( P( 3 ), P( 4 ), P( 1 ) );
180 aMin = Min( Min( A0, A1 ), Min( A2, A3 ) );
183 return aMin * 180 / PI;
186 SMDSAbs_ElementType MinimumAngle::GetType() const
194 Description : Functor for calculating aspect ratio
196 double AspectRatio::GetValue( long theId )
198 TColgp_SequenceOfXYZ P;
199 if ( !getPoints( theId, P ) || P.Length() != 3 && P.Length() != 4 )
202 int nbNodes = P.Length();
204 // Compute lengths of the sides
206 double aLen[ nbNodes ];
207 for ( int i = 0; i < nbNodes - 1; i++ )
208 aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
209 aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
211 // Compute aspect ratio
215 double aMaxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
216 double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
217 static double aCoef = sqrt( 3. ) / 4;
219 return anArea != 0 ? aCoef * aMaxLen * aMaxLen / anArea : 0;
223 double aMaxLen = Max( Max( aLen[ 0 ], aLen[ 1 ] ), Max( aLen[ 2 ], aLen[ 3 ] ) );
224 double aMinLen = Min( Min( aLen[ 0 ], aLen[ 1 ] ), Min( aLen[ 2 ], aLen[ 3 ] ) );
226 return aMinLen != 0 ? aMaxLen / aMinLen : 0;
230 SMDSAbs_ElementType AspectRatio::GetType() const
238 Description : Functor for calculating warping
240 double Warping::GetValue( long theId )
242 TColgp_SequenceOfXYZ P;
243 if ( !getPoints( theId, P ) || P.Length() != 4 )
246 gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4;
248 double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
249 double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
250 double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
251 double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
253 return Max( Max( A1, A2 ), Max( A3, A4 ) );
256 double Warping::ComputeA( const gp_XYZ& thePnt1,
257 const gp_XYZ& thePnt2,
258 const gp_XYZ& thePnt3,
259 const gp_XYZ& theG ) const
261 double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
262 double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
263 double L = Min( aLen1, aLen2 ) * 0.5;
265 gp_XYZ GI = ( thePnt2 - thePnt1 ) / 2. - theG;
266 gp_XYZ GJ = ( thePnt3 - thePnt2 ) / 2. - theG;
267 gp_XYZ N = GI.Crossed( GJ );
270 double H = gp_Vec( thePnt2 - theG ).Dot( gp_Vec( N ) );
271 return asin( fabs( H / L ) ) * 180 / PI;
274 SMDSAbs_ElementType Warping::GetType() const
282 Description : Functor for calculating taper
284 double Taper::GetValue( long theId )
286 TColgp_SequenceOfXYZ P;
287 if ( !getPoints( theId, P ) || P.Length() != 4 )
291 double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2;
292 double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2;
293 double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2;
294 double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2;
296 double JA = 0.25 * ( J1 + J2 + J3 + J4 );
298 double T1 = fabs( ( J1 - JA ) / JA );
299 double T2 = fabs( ( J2 - JA ) / JA );
300 double T3 = fabs( ( J3 - JA ) / JA );
301 double T4 = fabs( ( J4 - JA ) / JA );
303 return Max( Max( T1, T2 ), Max( T3, T4 ) );
306 SMDSAbs_ElementType Taper::GetType() const
314 Description : Functor for calculating skew in degrees
316 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
318 gp_XYZ p12 = ( p2 + p1 ) / 2;
319 gp_XYZ p23 = ( p3 + p2 ) / 2;
320 gp_XYZ p31 = ( p3 + p1 ) / 2;
322 return gp_Vec( p31 - p2 ).Angle( p12 - p23 );
325 double Skew::GetValue( long theId )
327 TColgp_SequenceOfXYZ P;
328 if ( !getPoints( theId, P ) || P.Length() != 3 && P.Length() != 4 )
332 static double PI2 = PI / 2;
333 if ( P.Length() == 3 )
335 double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
336 double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
337 double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
339 return Max( A0, Max( A1, A2 ) ) * 180 / PI;
343 gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2;
344 gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2;
345 gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2;
346 gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2;
348 double A = fabs( PI2 - gp_Vec( p34 - p12 ).Angle( p23 - p41 ) );
354 SMDSAbs_ElementType Skew::GetType() const
362 Description : Functor for calculating area
364 double Area::GetValue( long theId )
366 TColgp_SequenceOfXYZ P;
367 if ( !getPoints( theId, P ) || P.Length() != 3 && P.Length() != 4 )
370 if ( P.Length() == 3 )
371 return getArea( P( 1 ), P( 2 ), P( 3 ) );
373 return getArea( P( 1 ), P( 2 ), P( 3 ) ) + getArea( P( 1 ), P( 3 ), P( 4 ) );
376 SMDSAbs_ElementType Area::GetType() const
384 Description : Functor for calculating length off edge
386 double Length::GetValue( long theId )
388 TColgp_SequenceOfXYZ P;
389 return getPoints( theId, P ) && P.Length() == 2 ? getDistance( P( 1 ), P( 2 ) ) : 0;
392 SMDSAbs_ElementType Length::GetType() const
399 Class : MultiConnection
400 Description : Functor for calculating number of faces conneted to the edge
402 double MultiConnection::GetValue( long theId )
404 return getNbMultiConnection( myMesh, theId );
407 SMDSAbs_ElementType MultiConnection::GetType() const
419 Description : Predicate for free borders
422 FreeBorders::FreeBorders()
427 void FreeBorders::SetMesh( SMDS_Mesh* theMesh )
432 bool FreeBorders::IsSatisfy( long theId )
434 return getNbMultiConnection( myMesh, theId ) == 1;
437 SMDSAbs_ElementType FreeBorders::GetType() const
445 Description : Predicate for free Edges
447 FreeEdges::FreeEdges()
452 void FreeEdges::SetMesh( SMDS_Mesh* theMesh )
457 bool FreeEdges::IsSatisfy( long theId )
459 return getNbMultiConnection( myMesh, theId ) == 1;
462 SMDSAbs_ElementType FreeEdges::GetType() const
467 FreeEdges::Border::Border(long thePntId1, long thePntId2){
468 PntId[0] = thePntId1; PntId[1] = thePntId2;
469 if(thePntId1 > thePntId2){
470 PntId[1] = thePntId1; PntId[0] = thePntId2;
474 //bool operator<(const FreeEdges::Border& x, const FreeEdges::Border& y){
475 // if(x.PntId[0] < y.PntId[0]) return true;
476 // if(x.PntId[0] == y.PntId[0])
477 // if(x.PntId[1] < y.PntId[1]) return true;
483 struct EdgeBorder: public FreeEdges::Border{
485 EdgeBorder(long theElemId, long thePntId1, long thePntId2):
486 FreeEdges::Border(thePntId1,thePntId2),
491 bool operator<(const FreeEdges::Border& x, const FreeEdges::Border& y){
492 if(x.PntId[0] < y.PntId[0]) return true;
493 if(x.PntId[0] == y.PntId[0])
494 if(x.PntId[1] < y.PntId[1]) return true;
498 typedef std::set<EdgeBorder> EdgeBorderS;
500 inline void UpdateBorders(const EdgeBorder& theBorder,
501 EdgeBorderS& theRegistry,
502 EdgeBorderS& theContainer)
504 if(theRegistry.find(theBorder) == theRegistry.end()){
505 theRegistry.insert(theBorder);
506 theContainer.insert(theBorder);
508 theContainer.erase(theBorder);
515 void FreeEdges::GetBoreders(Borders& theBorders)
517 EdgeBorderS aRegistry;
518 EdgeBorderS aContainer;
519 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
520 for(; anIter->more(); ){
521 const SMDS_MeshFace* anElem = anIter->next();
522 long anElemId = anElem->GetID();
523 SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
525 const SMDS_MeshElement* aNode;
526 if(aNodesIter->more()){
527 aNode = aNodesIter->next();
528 aNodeId[0] = aNodeId[1] = aNode->GetID();
530 for(; aNodesIter->more(); ){
531 aNode = aNodesIter->next();
532 EdgeBorder aBorder(anElemId,aNodeId[1],aNode->GetID());
533 aNodeId[1] = aNode->GetID();
534 //std::cout<<aBorder.PntId[0]<<"; "<<aBorder.PntId[1]<<"; "<<aBorder.ElemId<<"\n";
535 UpdateBorders(aBorder,aRegistry,aContainer);
537 EdgeBorder aBorder(anElemId,aNodeId[0],aNodeId[1]);
538 UpdateBorders(aBorder,aRegistry,aContainer);
540 //std::cout<<"aContainer.size() = "<<aContainer.size()<<"\n";
541 if(aContainer.size()){
542 EdgeBorderS::const_iterator anIter = aContainer.begin();
543 for(; anIter != aContainer.end(); anIter++){
544 const EdgeBorder& aBorder = *anIter;
545 //std::cout<<aBorder.PntId[0]<<"; "<<aBorder.PntId[1]<<"; "<<aBorder.ElemId<<"\n";
546 theBorders.insert(Borders::value_type(aBorder.ElemId,aBorder));
554 Description : Base class for comparators
556 Comparator::Comparator():
560 Comparator::~Comparator()
563 void Comparator::SetMesh( SMDS_Mesh* theMesh )
566 myFunctor->SetMesh( theMesh );
569 void Comparator::SetMargin( double theValue )
574 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
576 myFunctor = theFunct;
579 SMDSAbs_ElementType Comparator::GetType() const
581 return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
587 Description : Comparator "<"
589 bool LessThan::IsSatisfy( long theId )
591 return myFunctor && myFunctor->GetValue( theId ) < myMargin;
597 Description : Comparator ">"
599 bool MoreThan::IsSatisfy( long theId )
601 return myFunctor && myFunctor->GetValue( theId ) > myMargin;
607 Description : Comparator "="
610 myToler(Precision::Confusion())
613 bool EqualTo::IsSatisfy( long theId )
615 return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
618 void EqualTo::SetTolerance( double theToler )
626 Description : Logical NOT predicate
628 LogicalNOT::LogicalNOT()
631 LogicalNOT::~LogicalNOT()
634 bool LogicalNOT::IsSatisfy( long theId )
636 return myPredicate && !myPredicate->IsSatisfy( theId );
639 void LogicalNOT::SetMesh( SMDS_Mesh* theMesh )
642 myPredicate->SetMesh( theMesh );
645 void LogicalNOT::SetPredicate( PredicatePtr thePred )
647 myPredicate = thePred;
650 SMDSAbs_ElementType LogicalNOT::GetType() const
652 return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
657 Class : LogicalBinary
658 Description : Base class for binary logical predicate
660 LogicalBinary::LogicalBinary()
663 LogicalBinary::~LogicalBinary()
666 void LogicalBinary::SetMesh( SMDS_Mesh* theMesh )
669 myPredicate1->SetMesh( theMesh );
672 myPredicate2->SetMesh( theMesh );
675 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
677 myPredicate1 = thePredicate;
680 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
682 myPredicate2 = thePredicate;
685 SMDSAbs_ElementType LogicalBinary::GetType() const
687 if ( !myPredicate1 || !myPredicate2 )
690 SMDSAbs_ElementType aType1 = myPredicate1->GetType();
691 SMDSAbs_ElementType aType2 = myPredicate2->GetType();
693 return aType1 == aType2 ? aType1 : SMDSAbs_All;
699 Description : Logical AND
701 bool LogicalAND::IsSatisfy( long theId )
706 myPredicate1->IsSatisfy( theId ) &&
707 myPredicate2->IsSatisfy( theId );
713 Description : Logical OR
715 bool LogicalOR::IsSatisfy( long theId )
720 myPredicate1->IsSatisfy( theId ) ||
721 myPredicate2->IsSatisfy( theId );
735 void Filter::SetPredicate( PredicatePtr thePredicate )
737 myPredicate = thePredicate;
741 Filter::GetElementsId( SMDS_Mesh* theMesh )
743 TIdSequence aSequence;
744 if ( !theMesh || !myPredicate ) return aSequence;
746 myPredicate->SetMesh( theMesh );
748 SMDSAbs_ElementType aType = myPredicate->GetType();
751 SMDS_EdgeIteratorPtr anIter = theMesh->edgesIterator();
753 while( anIter->more() ) {
754 const SMDS_MeshElement* anElem = anIter->next();
755 long anId = anElem->GetID();
756 if ( myPredicate->IsSatisfy( anId ) )
757 aSequence.push_back( anId );
762 SMDS_FaceIteratorPtr anIter = theMesh->facesIterator();
764 while( anIter->more() ) {
765 const SMDS_MeshElement* anElem = anIter->next();
766 long anId = anElem->GetID();
767 if ( myPredicate->IsSatisfy( anId ) )
768 aSequence.push_back( anId );