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_ControlsDef.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 <TColStd_MapOfInteger.hxx>
37 #include <TColStd_SequenceOfAsciiString.hxx>
38 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
41 #include <TopoDS_Face.hxx>
42 #include <TopoDS_Shape.hxx>
44 #include "SMDS_Mesh.hxx"
45 #include "SMDS_Iterator.hxx"
46 #include "SMDS_MeshElement.hxx"
47 #include "SMDS_MeshNode.hxx"
55 inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
57 gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
59 return v1.Magnitude() < gp::Resolution() ||
60 v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
63 inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
65 gp_Vec aVec1( P2 - P1 );
66 gp_Vec aVec2( P3 - P1 );
67 return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
70 inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
72 return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
77 inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
79 double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
83 int getNbMultiConnection( SMDS_Mesh* theMesh, const int theId )
88 const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
89 if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge || anEdge->NbNodes() != 2 )
92 TColStd_MapOfInteger aMap;
95 SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
97 while( anIter->more() ) {
98 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
101 SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
102 while( anElemIter->more() ) {
103 const SMDS_MeshElement* anElem = anElemIter->next();
104 if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
105 int anId = anElem->GetID();
107 if ( anIter->more() ) // i.e. first node
109 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 TSequenceOfXYZ& theRes ) const
152 return GetPoints( myMesh->FindElement( theId ), theRes );
155 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
156 TSequenceOfXYZ& 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.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
179 long NumericalFunctor::GetPrecision() const
184 void NumericalFunctor::SetPrecision( const long thePrecision )
186 myPrecision = thePrecision;
189 double NumericalFunctor::GetValue( long theId )
192 if ( GetPoints( theId, P ))
194 double aVal = GetValue( P );
195 if ( myPrecision >= 0 )
197 double prec = pow( 10., (double)( myPrecision ) );
198 aVal = floor( aVal * prec + 0.5 ) / prec;
208 Description : Functor for calculation of minimum angle
211 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
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.size() == 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 TSequenceOfXYZ& P )
256 int nbNodes = P.size();
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
311 inline double getHalfPerimeter(double theTria[3]){
312 return (theTria[0] + theTria[1] + theTria[2])/2.0;
315 inline double getArea(double theHalfPerim, double theTria[3]){
316 return sqrt(theHalfPerim*
317 (theHalfPerim-theTria[0])*
318 (theHalfPerim-theTria[1])*
319 (theHalfPerim-theTria[2]));
322 inline double getVolume(double theLen[6]){
323 double a2 = theLen[0]*theLen[0];
324 double b2 = theLen[1]*theLen[1];
325 double c2 = theLen[2]*theLen[2];
326 double d2 = theLen[3]*theLen[3];
327 double e2 = theLen[4]*theLen[4];
328 double f2 = theLen[5]*theLen[5];
329 double P = 4.0*a2*b2*d2;
330 double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
331 double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
332 return sqrt(P-Q+R)/12.0;
335 inline double getVolume2(double theLen[6]){
336 double a2 = theLen[0]*theLen[0];
337 double b2 = theLen[1]*theLen[1];
338 double c2 = theLen[2]*theLen[2];
339 double d2 = theLen[3]*theLen[3];
340 double e2 = theLen[4]*theLen[4];
341 double f2 = theLen[5]*theLen[5];
343 double P = a2*e2*(b2+c2+d2+f2-a2-e2);
344 double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
345 double R = c2*d2*(a2+b2+e2+f2-c2-d2);
346 double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
348 return sqrt(P+Q+R-S)/12.0;
351 inline double getVolume(const TSequenceOfXYZ& P){
352 gp_Vec aVec1( P( 2 ) - P( 1 ) );
353 gp_Vec aVec2( P( 3 ) - P( 1 ) );
354 gp_Vec aVec3( P( 4 ) - P( 1 ) );
355 gp_Vec anAreaVec( aVec1 ^ aVec2 );
356 return abs(aVec3 * anAreaVec) / 6.0;
359 inline double getMaxHeight(double theLen[6])
361 double aHeight = max(theLen[0],theLen[1]);
362 aHeight = max(aHeight,theLen[2]);
363 aHeight = max(aHeight,theLen[3]);
364 aHeight = max(aHeight,theLen[4]);
365 aHeight = max(aHeight,theLen[5]);
371 double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
373 double aQuality = 0.0;
374 int nbNodes = P.size();
378 getDistance(P( 1 ),P( 2 )), // a
379 getDistance(P( 2 ),P( 3 )), // b
380 getDistance(P( 3 ),P( 1 )), // c
381 getDistance(P( 2 ),P( 4 )), // d
382 getDistance(P( 3 ),P( 4 )), // e
383 getDistance(P( 1 ),P( 4 )) // f
385 double aTria[4][3] = {
386 {aLen[0],aLen[1],aLen[2]}, // abc
387 {aLen[0],aLen[3],aLen[5]}, // adf
388 {aLen[1],aLen[3],aLen[4]}, // bde
389 {aLen[2],aLen[4],aLen[5]} // cef
391 double aSumArea = 0.0;
392 double aHalfPerimeter = getHalfPerimeter(aTria[0]);
393 double anArea = getArea(aHalfPerimeter,aTria[0]);
395 aHalfPerimeter = getHalfPerimeter(aTria[1]);
396 anArea = getArea(aHalfPerimeter,aTria[1]);
398 aHalfPerimeter = getHalfPerimeter(aTria[2]);
399 anArea = getArea(aHalfPerimeter,aTria[2]);
401 aHalfPerimeter = getHalfPerimeter(aTria[3]);
402 anArea = getArea(aHalfPerimeter,aTria[3]);
404 double aVolume = getVolume(P);
405 //double aVolume = getVolume(aLen);
406 double aHeight = getMaxHeight(aLen);
407 static double aCoeff = sqrt(6.0)/36.0;
408 aQuality = aCoeff*aHeight*aSumArea/aVolume;
413 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
414 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
417 gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
418 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
421 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
422 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
425 gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
426 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
432 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
433 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
436 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
437 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
440 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
441 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
444 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
445 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
448 gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
449 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
452 gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
453 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
459 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
460 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
463 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
464 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
467 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
468 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
471 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
472 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
475 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
476 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
479 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
480 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
483 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
484 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
487 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
488 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
491 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
492 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
495 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
496 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
499 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
500 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
503 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
504 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
507 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
508 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
511 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
512 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
515 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
516 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
519 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
520 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
523 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
524 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
527 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
528 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
531 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
532 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
535 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
536 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
539 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
540 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
543 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
544 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
547 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
548 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
551 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
552 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
555 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
556 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
559 gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
560 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
563 gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
564 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
567 gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
568 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
571 gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
572 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
575 gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
576 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
579 gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
580 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
583 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
584 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
587 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
588 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
596 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
598 // the aspect ratio is in the range [1.0,infinity]
601 return Value / 1000.;
604 SMDSAbs_ElementType AspectRatio3D::GetType() const
606 return SMDSAbs_Volume;
612 Description : Functor for calculating warping
614 double Warping::GetValue( const TSequenceOfXYZ& P )
619 gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4;
621 double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
622 double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
623 double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
624 double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
626 return Max( Max( A1, A2 ), Max( A3, A4 ) );
629 double Warping::ComputeA( const gp_XYZ& thePnt1,
630 const gp_XYZ& thePnt2,
631 const gp_XYZ& thePnt3,
632 const gp_XYZ& theG ) const
634 double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
635 double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
636 double L = Min( aLen1, aLen2 ) * 0.5;
637 if ( L < Precision::Confusion())
640 gp_XYZ GI = ( thePnt2 - thePnt1 ) / 2. - theG;
641 gp_XYZ GJ = ( thePnt3 - thePnt2 ) / 2. - theG;
642 gp_XYZ N = GI.Crossed( GJ );
644 if ( N.Modulus() < gp::Resolution() )
649 double H = ( thePnt2 - theG ).Dot( N );
650 return asin( fabs( H / L ) ) * 180 / PI;
653 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
655 // the warp is in the range [0.0,PI/2]
656 // 0.0 = good (no warp)
657 // PI/2 = bad (face pliee)
661 SMDSAbs_ElementType Warping::GetType() const
669 Description : Functor for calculating taper
671 double Taper::GetValue( const TSequenceOfXYZ& P )
677 double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2;
678 double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2;
679 double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2;
680 double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2;
682 double JA = 0.25 * ( J1 + J2 + J3 + J4 );
683 if ( JA <= Precision::Confusion() )
686 double T1 = fabs( ( J1 - JA ) / JA );
687 double T2 = fabs( ( J2 - JA ) / JA );
688 double T3 = fabs( ( J3 - JA ) / JA );
689 double T4 = fabs( ( J4 - JA ) / JA );
691 return Max( Max( T1, T2 ), Max( T3, T4 ) );
694 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
696 // the taper is in the range [0.0,1.0]
697 // 0.0 = good (no taper)
698 // 1.0 = bad (les cotes opposes sont allignes)
702 SMDSAbs_ElementType Taper::GetType() const
710 Description : Functor for calculating skew in degrees
712 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
714 gp_XYZ p12 = ( p2 + p1 ) / 2;
715 gp_XYZ p23 = ( p3 + p2 ) / 2;
716 gp_XYZ p31 = ( p3 + p1 ) / 2;
718 gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
720 return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
723 double Skew::GetValue( const TSequenceOfXYZ& P )
725 if ( P.size() != 3 && P.size() != 4 )
729 static double PI2 = PI / 2;
732 double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
733 double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
734 double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
736 return Max( A0, Max( A1, A2 ) ) * 180 / PI;
740 gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2;
741 gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2;
742 gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2;
743 gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2;
745 gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
746 double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
747 ? 0 : fabs( PI2 - v1.Angle( v2 ) );
753 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
755 // the skew is in the range [0.0,PI/2].
761 SMDSAbs_ElementType Skew::GetType() const
769 Description : Functor for calculating area
771 double Area::GetValue( const TSequenceOfXYZ& P )
774 return getArea( P( 1 ), P( 2 ), P( 3 ) );
775 else if ( P.size() == 4 )
776 return getArea( P( 1 ), P( 2 ), P( 3 ) ) + getArea( P( 1 ), P( 3 ), P( 4 ) );
781 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
786 SMDSAbs_ElementType Area::GetType() const
794 Description : Functor for calculating length off edge
796 double Length::GetValue( const TSequenceOfXYZ& P )
798 return ( P.size() == 2 ? getDistance( P( 1 ), P( 2 ) ) : 0 );
801 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
806 SMDSAbs_ElementType Length::GetType() const
813 Description : Functor for calculating length of edge
816 double Length2D::GetValue( long theElementId)
820 if (GetPoints(theElementId,P)){
822 double aVal;// = GetValue( P );
823 const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
824 SMDSAbs_ElementType aType = aElem->GetType();
833 aVal = getDistance( P( 1 ), P( 2 ) );
837 if (len == 3){ // triangles
838 double L1 = getDistance(P( 1 ),P( 2 ));
839 double L2 = getDistance(P( 2 ),P( 3 ));
840 double L3 = getDistance(P( 3 ),P( 1 ));
841 aVal = Max(L1,Max(L2,L3));
844 else if (len == 4){ // quadrangles
845 double L1 = getDistance(P( 1 ),P( 2 ));
846 double L2 = getDistance(P( 2 ),P( 3 ));
847 double L3 = getDistance(P( 3 ),P( 4 ));
848 double L4 = getDistance(P( 4 ),P( 1 ));
849 aVal = Max(Max(L1,L2),Max(L3,L4));
853 if (len == 4){ // tetraidrs
854 double L1 = getDistance(P( 1 ),P( 2 ));
855 double L2 = getDistance(P( 2 ),P( 3 ));
856 double L3 = getDistance(P( 3 ),P( 1 ));
857 double L4 = getDistance(P( 1 ),P( 4 ));
858 double L5 = getDistance(P( 2 ),P( 4 ));
859 double L6 = getDistance(P( 3 ),P( 4 ));
860 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
863 else if (len == 5){ // piramids
864 double L1 = getDistance(P( 1 ),P( 2 ));
865 double L2 = getDistance(P( 2 ),P( 3 ));
866 double L3 = getDistance(P( 3 ),P( 1 ));
867 double L4 = getDistance(P( 4 ),P( 1 ));
868 double L5 = getDistance(P( 1 ),P( 5 ));
869 double L6 = getDistance(P( 2 ),P( 5 ));
870 double L7 = getDistance(P( 3 ),P( 5 ));
871 double L8 = getDistance(P( 4 ),P( 5 ));
873 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
874 aVal = Max(aVal,Max(L7,L8));
877 else if (len == 6){ // pentaidres
878 double L1 = getDistance(P( 1 ),P( 2 ));
879 double L2 = getDistance(P( 2 ),P( 3 ));
880 double L3 = getDistance(P( 3 ),P( 1 ));
881 double L4 = getDistance(P( 4 ),P( 5 ));
882 double L5 = getDistance(P( 5 ),P( 6 ));
883 double L6 = getDistance(P( 6 ),P( 4 ));
884 double L7 = getDistance(P( 1 ),P( 4 ));
885 double L8 = getDistance(P( 2 ),P( 5 ));
886 double L9 = getDistance(P( 3 ),P( 6 ));
888 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
889 aVal = Max(aVal,Max(Max(L7,L8),L9));
892 else if (len == 8){ // hexaider
893 double L1 = getDistance(P( 1 ),P( 2 ));
894 double L2 = getDistance(P( 2 ),P( 3 ));
895 double L3 = getDistance(P( 3 ),P( 4 ));
896 double L4 = getDistance(P( 4 ),P( 1 ));
897 double L5 = getDistance(P( 5 ),P( 6 ));
898 double L6 = getDistance(P( 6 ),P( 7 ));
899 double L7 = getDistance(P( 7 ),P( 8 ));
900 double L8 = getDistance(P( 8 ),P( 5 ));
901 double L9 = getDistance(P( 1 ),P( 5 ));
902 double L10= getDistance(P( 2 ),P( 6 ));
903 double L11= getDistance(P( 3 ),P( 7 ));
904 double L12= getDistance(P( 4 ),P( 8 ));
906 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
907 aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
908 aVal = Max(aVal,Max(L11,L12));
920 if ( myPrecision >= 0 )
922 double prec = pow( 10., (double)( myPrecision ) );
923 aVal = floor( aVal * prec + 0.5 ) / prec;
932 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
937 SMDSAbs_ElementType Length2D::GetType() const
942 Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
945 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
946 if(thePntId1 > thePntId2){
947 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
951 bool Length2D::Value::operator<(const Length2D::Value& x) const{
952 if(myPntId[0] < x.myPntId[0]) return true;
953 if(myPntId[0] == x.myPntId[0])
954 if(myPntId[1] < x.myPntId[1]) return true;
958 void Length2D::GetValues(TValues& theValues){
960 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
962 for(; anIter->more(); ){
963 const SMDS_MeshFace* anElem = anIter->next();
964 long anElemId = anElem->GetID();
965 SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
970 const SMDS_MeshElement* aNode;
971 if(aNodesIter->more()){
972 aNode = aNodesIter->next();
973 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
974 P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
975 aNodeId[0] = aNodeId[1] = aNode->GetID();
978 for(; aNodesIter->more(); ){
979 aNode = aNodesIter->next();
980 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
981 long anId = aNode->GetID();
983 P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
985 aLength = P[1].Distance(P[2]);
987 Value aValue(aLength,aNodeId[1],anId);
990 theValues.insert(aValue);
993 aLength = P[0].Distance(P[1]);
995 Value aValue(aLength,aNodeId[0],aNodeId[1]);
996 theValues.insert(aValue);
1001 Class : MultiConnection
1002 Description : Functor for calculating number of faces conneted to the edge
1004 double MultiConnection::GetValue( const TSequenceOfXYZ& P )
1008 double MultiConnection::GetValue( long theId )
1010 return getNbMultiConnection( myMesh, theId );
1013 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
1018 SMDSAbs_ElementType MultiConnection::GetType() const
1020 return SMDSAbs_Edge;
1024 Class : MultiConnection2D
1025 Description : Functor for calculating number of faces conneted to the edge
1027 double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
1032 double MultiConnection2D::GetValue( long theElementId )
1037 if (GetPoints(theElementId,P)){
1039 const SMDS_MeshElement* anFaceElem = myMesh->FindElement( theElementId );
1040 SMDSAbs_ElementType aType = anFaceElem->GetType();
1044 TColStd_MapOfInteger aMap;
1052 if (len == 3){ // triangles
1053 int Nb[3] = {0,0,0};
1056 SMDS_ElemIteratorPtr anIter = anFaceElem->nodesIterator();
1057 if ( anIter != 0 ) {
1058 while( anIter->more() ) {
1059 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
1063 SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
1064 while( anElemIter->more() ) {
1065 const SMDS_MeshElement* anElem = anElemIter->next();
1066 if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
1067 int anId = anElem->GetID();
1069 if ( anIter->more() ) // i.e. first node
1071 else if ( aMap.Contains( anId ) ){
1075 else if ( anElem != 0 && anElem->GetType() == SMDSAbs_Edge ) i++;
1080 aResult = Max(Max(Nb[0],Nb[1]),Nb[2]);
1083 case SMDSAbs_Volume:
1088 return aResult;//getNbMultiConnection( myMesh, theId );
1091 double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1096 SMDSAbs_ElementType MultiConnection2D::GetType() const
1098 return SMDSAbs_Face;
1101 MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
1103 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
1104 if(thePntId1 > thePntId2){
1105 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
1109 bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const{
1110 if(myPntId[0] < x.myPntId[0]) return true;
1111 if(myPntId[0] == x.myPntId[0])
1112 if(myPntId[1] < x.myPntId[1]) return true;
1116 void MultiConnection2D::GetValues(MValues& theValues){
1117 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1118 for(; anIter->more(); ){
1119 const SMDS_MeshFace* anElem = anIter->next();
1120 long anElemId = anElem->GetID();
1121 SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
1124 //int aNbConnects=0;
1125 const SMDS_MeshNode* aNode0;
1126 const SMDS_MeshNode* aNode1;
1127 const SMDS_MeshNode* aNode2;
1128 if(aNodesIter->more()){
1129 aNode0 = (SMDS_MeshNode*) aNodesIter->next();
1131 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
1132 aNodeId[0] = aNodeId[1] = aNodes->GetID();
1134 for(; aNodesIter->more(); ){
1135 aNode2 = (SMDS_MeshNode*) aNodesIter->next();
1136 long anId = aNode2->GetID();
1139 Value aValue(aNodeId[1],aNodeId[2]);
1140 MValues::iterator aItr = theValues.find(aValue);
1141 if (aItr != theValues.end()){
1145 theValues[aValue] = 1;
1148 //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
1149 aNodeId[1] = aNodeId[2];
1152 Value aValue(aNodeId[0],aNodeId[2]);
1153 MValues::iterator aItr = theValues.find(aValue);
1154 if (aItr != theValues.end()){
1158 theValues[aValue] = 1;
1161 //cout << "NodeIds: "<<aNodeId[0]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
1172 Description : Predicate for free borders
1175 FreeBorders::FreeBorders()
1180 void FreeBorders::SetMesh( SMDS_Mesh* theMesh )
1185 bool FreeBorders::IsSatisfy( long theId )
1187 return getNbMultiConnection( myMesh, theId ) == 1;
1190 SMDSAbs_ElementType FreeBorders::GetType() const
1192 return SMDSAbs_Edge;
1198 Description : Predicate for free Edges
1200 FreeEdges::FreeEdges()
1205 void FreeEdges::SetMesh( SMDS_Mesh* theMesh )
1210 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId )
1212 TColStd_MapOfInteger aMap;
1213 for ( int i = 0; i < 2; i++ )
1215 SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator();
1216 while( anElemIter->more() )
1218 const SMDS_MeshElement* anElem = anElemIter->next();
1219 if ( anElem != 0 && anElem->GetType() == SMDSAbs_Face )
1221 int anId = anElem->GetID();
1225 else if ( aMap.Contains( anId ) && anId != theFaceId )
1233 bool FreeEdges::IsSatisfy( long theId )
1238 const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
1239 if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
1242 int nbNodes = aFace->NbNodes();
1243 const SMDS_MeshNode* aNodes[ nbNodes ];
1245 SMDS_ElemIteratorPtr anIter = aFace->nodesIterator();
1248 while( anIter->more() )
1250 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
1253 aNodes[ i++ ] = aNode;
1257 for ( int i = 0; i < nbNodes - 1; i++ )
1258 if ( IsFreeEdge( &aNodes[ i ], theId ) )
1261 aNodes[ 1 ] = aNodes[ nbNodes - 1 ];
1263 return IsFreeEdge( &aNodes[ 0 ], theId );
1267 SMDSAbs_ElementType FreeEdges::GetType() const
1269 return SMDSAbs_Face;
1272 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
1275 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
1276 if(thePntId1 > thePntId2){
1277 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
1281 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
1282 if(myPntId[0] < x.myPntId[0]) return true;
1283 if(myPntId[0] == x.myPntId[0])
1284 if(myPntId[1] < x.myPntId[1]) return true;
1288 inline void UpdateBorders(const FreeEdges::Border& theBorder,
1289 FreeEdges::TBorders& theRegistry,
1290 FreeEdges::TBorders& theContainer)
1292 if(theRegistry.find(theBorder) == theRegistry.end()){
1293 theRegistry.insert(theBorder);
1294 theContainer.insert(theBorder);
1296 theContainer.erase(theBorder);
1300 void FreeEdges::GetBoreders(TBorders& theBorders)
1303 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1304 for(; anIter->more(); ){
1305 const SMDS_MeshFace* anElem = anIter->next();
1306 long anElemId = anElem->GetID();
1307 SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
1309 const SMDS_MeshElement* aNode;
1310 if(aNodesIter->more()){
1311 aNode = aNodesIter->next();
1312 aNodeId[0] = aNodeId[1] = aNode->GetID();
1314 for(; aNodesIter->more(); ){
1315 aNode = aNodesIter->next();
1316 long anId = aNode->GetID();
1317 Border aBorder(anElemId,aNodeId[1],anId);
1319 //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
1320 UpdateBorders(aBorder,aRegistry,theBorders);
1322 Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
1323 //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
1324 UpdateBorders(aBorder,aRegistry,theBorders);
1326 //std::cout<<"theBorders.size() = "<<theBorders.size()<<endl;
1331 Description : Predicate for Range of Ids.
1332 Range may be specified with two ways.
1333 1. Using AddToRange method
1334 2. With SetRangeStr method. Parameter of this method is a string
1335 like as "1,2,3,50-60,63,67,70-"
1338 //=======================================================================
1339 // name : RangeOfIds
1340 // Purpose : Constructor
1341 //=======================================================================
1342 RangeOfIds::RangeOfIds()
1345 myType = SMDSAbs_All;
1348 //=======================================================================
1350 // Purpose : Set mesh
1351 //=======================================================================
1352 void RangeOfIds::SetMesh( SMDS_Mesh* theMesh )
1357 //=======================================================================
1358 // name : AddToRange
1359 // Purpose : Add ID to the range
1360 //=======================================================================
1361 bool RangeOfIds::AddToRange( long theEntityId )
1363 myIds.Add( theEntityId );
1367 //=======================================================================
1368 // name : GetRangeStr
1369 // Purpose : Get range as a string.
1370 // Example: "1,2,3,50-60,63,67,70-"
1371 //=======================================================================
1372 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
1376 TColStd_SequenceOfInteger anIntSeq;
1377 TColStd_SequenceOfAsciiString aStrSeq;
1379 TColStd_MapIteratorOfMapOfInteger anIter( myIds );
1380 for ( ; anIter.More(); anIter.Next() )
1382 int anId = anIter.Key();
1383 TCollection_AsciiString aStr( anId );
1384 anIntSeq.Append( anId );
1385 aStrSeq.Append( aStr );
1388 for ( int i = 1, n = myMin.Length(); i <= n; i++ )
1390 int aMinId = myMin( i );
1391 int aMaxId = myMax( i );
1393 TCollection_AsciiString aStr;
1394 if ( aMinId != IntegerFirst() )
1399 if ( aMaxId != IntegerLast() )
1402 // find position of the string in result sequence and insert string in it
1403 if ( anIntSeq.Length() == 0 )
1405 anIntSeq.Append( aMinId );
1406 aStrSeq.Append( aStr );
1410 if ( aMinId < anIntSeq.First() )
1412 anIntSeq.Prepend( aMinId );
1413 aStrSeq.Prepend( aStr );
1415 else if ( aMinId > anIntSeq.Last() )
1417 anIntSeq.Append( aMinId );
1418 aStrSeq.Append( aStr );
1421 for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
1422 if ( aMinId < anIntSeq( j ) )
1424 anIntSeq.InsertBefore( j, aMinId );
1425 aStrSeq.InsertBefore( j, aStr );
1431 if ( aStrSeq.Length() == 0 )
1434 theResStr = aStrSeq( 1 );
1435 for ( int j = 2, k = aStrSeq.Length(); j <= k; j++ )
1438 theResStr += aStrSeq( j );
1442 //=======================================================================
1443 // name : SetRangeStr
1444 // Purpose : Define range with string
1445 // Example of entry string: "1,2,3,50-60,63,67,70-"
1446 //=======================================================================
1447 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
1453 TCollection_AsciiString aStr = theStr;
1454 aStr.RemoveAll( ' ' );
1455 aStr.RemoveAll( '\t' );
1457 for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) )
1458 aStr.Remove( aPos, 2 );
1460 TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
1462 while ( tmpStr != "" )
1464 tmpStr = aStr.Token( ",", i++ );
1465 int aPos = tmpStr.Search( '-' );
1469 if ( tmpStr.IsIntegerValue() )
1470 myIds.Add( tmpStr.IntegerValue() );
1476 TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
1477 TCollection_AsciiString aMinStr = tmpStr;
1479 while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
1480 while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
1482 if ( !aMinStr.IsEmpty() && !aMinStr.IsIntegerValue() ||
1483 !aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue() )
1486 myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
1487 myMax.Append( aMaxStr.IsEmpty() ? IntegerLast() : aMaxStr.IntegerValue() );
1494 //=======================================================================
1496 // Purpose : Get type of supported entities
1497 //=======================================================================
1498 SMDSAbs_ElementType RangeOfIds::GetType() const
1503 //=======================================================================
1505 // Purpose : Set type of supported entities
1506 //=======================================================================
1507 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
1512 //=======================================================================
1514 // Purpose : Verify whether entity satisfies to this rpedicate
1515 //=======================================================================
1516 bool RangeOfIds::IsSatisfy( long theId )
1521 if ( myType == SMDSAbs_Node )
1523 if ( myMesh->FindNode( theId ) == 0 )
1528 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
1529 if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All )
1533 if ( myIds.Contains( theId ) )
1536 for ( int i = 1, n = myMin.Length(); i <= n; i++ )
1537 if ( theId >= myMin( i ) && theId <= myMax( i ) )
1545 Description : Base class for comparators
1547 Comparator::Comparator():
1551 Comparator::~Comparator()
1554 void Comparator::SetMesh( SMDS_Mesh* theMesh )
1557 myFunctor->SetMesh( theMesh );
1560 void Comparator::SetMargin( double theValue )
1562 myMargin = theValue;
1565 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
1567 myFunctor = theFunct;
1570 SMDSAbs_ElementType Comparator::GetType() const
1572 return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
1575 double Comparator::GetMargin()
1583 Description : Comparator "<"
1585 bool LessThan::IsSatisfy( long theId )
1587 return myFunctor && myFunctor->GetValue( theId ) < myMargin;
1593 Description : Comparator ">"
1595 bool MoreThan::IsSatisfy( long theId )
1597 return myFunctor && myFunctor->GetValue( theId ) > myMargin;
1603 Description : Comparator "="
1606 myToler(Precision::Confusion())
1609 bool EqualTo::IsSatisfy( long theId )
1611 return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
1614 void EqualTo::SetTolerance( double theToler )
1619 double EqualTo::GetTolerance()
1626 Description : Logical NOT predicate
1628 LogicalNOT::LogicalNOT()
1631 LogicalNOT::~LogicalNOT()
1634 bool LogicalNOT::IsSatisfy( long theId )
1636 return myPredicate && !myPredicate->IsSatisfy( theId );
1639 void LogicalNOT::SetMesh( SMDS_Mesh* theMesh )
1642 myPredicate->SetMesh( theMesh );
1645 void LogicalNOT::SetPredicate( PredicatePtr thePred )
1647 myPredicate = thePred;
1650 SMDSAbs_ElementType LogicalNOT::GetType() const
1652 return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
1657 Class : LogicalBinary
1658 Description : Base class for binary logical predicate
1660 LogicalBinary::LogicalBinary()
1663 LogicalBinary::~LogicalBinary()
1666 void LogicalBinary::SetMesh( SMDS_Mesh* theMesh )
1669 myPredicate1->SetMesh( theMesh );
1672 myPredicate2->SetMesh( theMesh );
1675 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
1677 myPredicate1 = thePredicate;
1680 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
1682 myPredicate2 = thePredicate;
1685 SMDSAbs_ElementType LogicalBinary::GetType() const
1687 if ( !myPredicate1 || !myPredicate2 )
1690 SMDSAbs_ElementType aType1 = myPredicate1->GetType();
1691 SMDSAbs_ElementType aType2 = myPredicate2->GetType();
1693 return aType1 == aType2 ? aType1 : SMDSAbs_All;
1699 Description : Logical AND
1701 bool LogicalAND::IsSatisfy( long theId )
1706 myPredicate1->IsSatisfy( theId ) &&
1707 myPredicate2->IsSatisfy( theId );
1713 Description : Logical OR
1715 bool LogicalOR::IsSatisfy( long theId )
1720 myPredicate1->IsSatisfy( theId ) ||
1721 myPredicate2->IsSatisfy( theId );
1735 void Filter::SetPredicate( PredicatePtr thePredicate )
1737 myPredicate = thePredicate;
1741 template<class TElement, class TIterator, class TPredicate>
1742 void FillSequence(const TIterator& theIterator,
1743 TPredicate& thePredicate,
1744 Filter::TIdSequence& theSequence)
1746 if ( theIterator ) {
1747 while( theIterator->more() ) {
1748 TElement anElem = theIterator->next();
1749 long anId = anElem->GetID();
1750 if ( thePredicate->IsSatisfy( anId ) )
1751 theSequence.push_back( anId );
1757 Filter::GetElementsId( SMDS_Mesh* theMesh )
1759 TIdSequence aSequence;
1760 if ( !theMesh || !myPredicate ) return aSequence;
1762 myPredicate->SetMesh( theMesh );
1764 SMDSAbs_ElementType aType = myPredicate->GetType();
1767 FillSequence<const SMDS_MeshNode*>(theMesh->nodesIterator(),myPredicate,aSequence);
1771 FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
1775 FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
1778 case SMDSAbs_Volume:{
1779 FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
1783 FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
1784 FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
1785 FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
1796 typedef std::set<SMDS_MeshFace*> TMapOfFacePtr;
1802 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
1803 SMDS_MeshNode* theNode2 )
1809 ManifoldPart::Link::~Link()
1815 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
1817 if ( myNode1 == theLink.myNode1 &&
1818 myNode2 == theLink.myNode2 )
1820 else if ( myNode1 == theLink.myNode2 &&
1821 myNode2 == theLink.myNode1 )
1827 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
1829 if(myNode1 < x.myNode1) return true;
1830 if(myNode1 == x.myNode1)
1831 if(myNode2 < x.myNode2) return true;
1835 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
1836 const ManifoldPart::Link& theLink2 )
1838 return theLink1.IsEqual( theLink2 );
1841 ManifoldPart::ManifoldPart()
1844 myAngToler = Precision::Angular();
1845 myIsOnlyManifold = true;
1848 ManifoldPart::~ManifoldPart()
1853 void ManifoldPart::SetMesh( SMDS_Mesh* theMesh )
1859 SMDSAbs_ElementType ManifoldPart::GetType() const
1860 { return SMDSAbs_Face; }
1862 bool ManifoldPart::IsSatisfy( long theElementId )
1864 return myMapIds.Contains( theElementId );
1867 void ManifoldPart::SetAngleTolerance( const double theAngToler )
1868 { myAngToler = theAngToler; }
1870 double ManifoldPart::GetAngleTolerance() const
1871 { return myAngToler; }
1873 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
1874 { myIsOnlyManifold = theIsOnly; }
1876 void ManifoldPart::SetStartElem( const long theStartId )
1877 { myStartElemId = theStartId; }
1879 bool ManifoldPart::process()
1882 myMapBadGeomIds.Clear();
1884 myAllFacePtr.clear();
1885 myAllFacePtrIntDMap.clear();
1889 // collect all faces into own map
1890 SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
1891 for (; anFaceItr->more(); )
1893 SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
1894 myAllFacePtr.push_back( aFacePtr );
1895 myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
1898 SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
1902 // the map of non manifold links and bad geometry
1903 TMapOfLink aMapOfNonManifold;
1904 TColStd_MapOfInteger aMapOfTreated;
1906 // begin cycle on faces from start index and run on vector till the end
1907 // and from begin to start index to cover whole vector
1908 const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
1909 bool isStartTreat = false;
1910 for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
1912 if ( fi == aStartIndx )
1913 isStartTreat = true;
1914 // as result next time when fi will be equal to aStartIndx
1916 SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
1917 if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
1920 aMapOfTreated.Add( aFacePtr->GetID() );
1921 TColStd_MapOfInteger aResFaces;
1922 if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
1923 aMapOfNonManifold, aResFaces ) )
1925 TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
1926 for ( ; anItr.More(); anItr.Next() )
1928 int aFaceId = anItr.Key();
1929 aMapOfTreated.Add( aFaceId );
1930 myMapIds.Add( aFaceId );
1933 if ( fi == ( myAllFacePtr.size() - 1 ) )
1935 } // end run on vector of faces
1936 return !myMapIds.IsEmpty();
1939 static void getLinks( const SMDS_MeshFace* theFace,
1940 ManifoldPart::TVectorOfLink& theLinks )
1942 int aNbNode = theFace->NbNodes();
1943 SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
1945 SMDS_MeshNode* aNode = 0;
1946 for ( ; aNodeItr->more() && i <= aNbNode; )
1949 SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
1953 SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
1955 ManifoldPart::Link aLink( aN1, aN2 );
1956 theLinks.push_back( aLink );
1960 static gp_XYZ getNormale( const SMDS_MeshFace* theFace )
1963 int aNbNode = theFace->NbNodes();
1964 TColgp_Array1OfXYZ anArrOfXYZ(1,4);
1965 SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
1967 for ( ; aNodeItr->more() && i <= 4; i++ )
1969 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
1970 anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
1973 gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1);
1974 gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1);
1978 gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1);
1981 double len = n.Modulus();
1988 bool ManifoldPart::findConnected
1989 ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
1990 SMDS_MeshFace* theStartFace,
1991 ManifoldPart::TMapOfLink& theNonManifold,
1992 TColStd_MapOfInteger& theResFaces )
1994 theResFaces.Clear();
1995 if ( !theAllFacePtrInt.size() )
1998 if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
2000 myMapBadGeomIds.Add( theStartFace->GetID() );
2004 ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
2005 ManifoldPart::TVectorOfLink aSeqOfBoundary;
2006 theResFaces.Add( theStartFace->GetID() );
2007 ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
2009 expandBoundary( aMapOfBoundary, aSeqOfBoundary,
2010 aDMapLinkFace, theNonManifold, theStartFace );
2012 bool isDone = false;
2013 while ( !isDone && aMapOfBoundary.size() != 0 )
2015 bool isToReset = false;
2016 ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
2017 for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
2019 ManifoldPart::Link aLink = *pLink;
2020 if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
2022 // each link could be treated only once
2023 aMapToSkip.insert( aLink );
2025 ManifoldPart::TVectorOfFacePtr aFaces;
2027 if ( myIsOnlyManifold &&
2028 (theNonManifold.find( aLink ) != theNonManifold.end()) )
2032 getFacesByLink( aLink, aFaces );
2033 // filter the element to keep only indicated elements
2034 ManifoldPart::TVectorOfFacePtr aFiltered;
2035 ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
2036 for ( ; pFace != aFaces.end(); ++pFace )
2038 SMDS_MeshFace* aFace = *pFace;
2039 if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
2040 aFiltered.push_back( aFace );
2043 if ( aFaces.size() < 2 ) // no neihgbour faces
2045 else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
2047 theNonManifold.insert( aLink );
2052 // compare normal with normals of neighbor element
2053 SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
2054 ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
2055 for ( ; pFace != aFaces.end(); ++pFace )
2057 SMDS_MeshFace* aNextFace = *pFace;
2058 if ( aPrevFace == aNextFace )
2060 int anNextFaceID = aNextFace->GetID();
2061 if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
2062 // should not be with non manifold restriction. probably bad topology
2064 // check if face was treated and skipped
2065 if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
2066 !isInPlane( aPrevFace, aNextFace ) )
2068 // add new element to connected and extend the boundaries.
2069 theResFaces.Add( anNextFaceID );
2070 expandBoundary( aMapOfBoundary, aSeqOfBoundary,
2071 aDMapLinkFace, theNonManifold, aNextFace );
2075 isDone = !isToReset;
2078 return !theResFaces.IsEmpty();
2081 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
2082 const SMDS_MeshFace* theFace2 )
2084 gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
2085 gp_XYZ aNorm2XYZ = getNormale( theFace2 );
2086 if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
2088 myMapBadGeomIds.Add( theFace2->GetID() );
2091 if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
2097 void ManifoldPart::expandBoundary
2098 ( ManifoldPart::TMapOfLink& theMapOfBoundary,
2099 ManifoldPart::TVectorOfLink& theSeqOfBoundary,
2100 ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
2101 ManifoldPart::TMapOfLink& theNonManifold,
2102 SMDS_MeshFace* theNextFace ) const
2104 ManifoldPart::TVectorOfLink aLinks;
2105 getLinks( theNextFace, aLinks );
2106 int aNbLink = aLinks.size();
2107 for ( int i = 0; i < aNbLink; i++ )
2109 ManifoldPart::Link aLink = aLinks[ i ];
2110 if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
2112 if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
2114 if ( myIsOnlyManifold )
2116 // remove from boundary
2117 theMapOfBoundary.erase( aLink );
2118 ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
2119 for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
2121 ManifoldPart::Link aBoundLink = *pLink;
2122 if ( aBoundLink.IsEqual( aLink ) )
2124 theSeqOfBoundary.erase( pLink );
2132 theMapOfBoundary.insert( aLink );
2133 theSeqOfBoundary.push_back( aLink );
2134 theDMapLinkFacePtr[ aLink ] = theNextFace;
2139 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
2140 ManifoldPart::TVectorOfFacePtr& theFaces ) const
2142 SMDS_Mesh::SetOfFaces aSetOfFaces;
2143 // take all faces that shared first node
2144 SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
2145 for ( ; anItr->more(); )
2147 SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
2150 aSetOfFaces.insert( aFace );
2152 // take all faces that shared second node
2153 anItr = theLink.myNode2->facesIterator();
2154 // find the common part of two sets
2155 for ( ; anItr->more(); )
2157 SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
2158 if ( aSetOfFaces.find( aFace ) != aSetOfFaces.end() )
2159 theFaces.push_back( aFace );
2168 ElementsOnSurface::ElementsOnSurface()
2172 myType = SMDSAbs_All;
2174 myToler = Precision::Confusion();
2177 ElementsOnSurface::~ElementsOnSurface()
2182 void ElementsOnSurface::SetMesh( SMDS_Mesh* theMesh )
2184 if ( myMesh == theMesh )
2191 bool ElementsOnSurface::IsSatisfy( long theElementId )
2193 return myIds.Contains( theElementId );
2196 SMDSAbs_ElementType ElementsOnSurface::GetType() const
2199 void ElementsOnSurface::SetTolerance( const double theToler )
2200 { myToler = theToler; }
2202 double ElementsOnSurface::GetTolerance() const
2207 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
2208 const SMDSAbs_ElementType theType )
2212 if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
2217 TopoDS_Face aFace = TopoDS::Face( theShape );
2218 mySurf = BRep_Tool::Surface( aFace );
2221 void ElementsOnSurface::process()
2224 if ( mySurf.IsNull() )
2230 if ( myType == SMDSAbs_Face || myType == SMDSAbs_All )
2232 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
2233 for(; anIter->more(); )
2234 process( anIter->next() );
2237 if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All )
2239 SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator();
2240 for(; anIter->more(); )
2241 process( anIter->next() );
2244 if ( myType == SMDSAbs_Node )
2246 SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator();
2247 for(; anIter->more(); )
2248 process( anIter->next() );
2252 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
2254 SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
2255 bool isSatisfy = true;
2256 for ( ; aNodeItr->more(); )
2258 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
2259 if ( !isOnSurface( aNode ) )
2266 myIds.Add( theElemPtr->GetID() );
2269 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) const
2271 if ( mySurf.IsNull() )
2274 gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
2275 double aToler2 = myToler * myToler;
2276 if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
2278 gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
2279 if ( aPln.SquareDistance( aPnt ) > aToler2 )
2282 else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
2284 gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
2285 double aRad = aCyl.Radius();
2286 gp_Ax3 anAxis = aCyl.Position();
2287 gp_XYZ aLoc = aCyl.Location().XYZ();
2288 double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
2289 double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
2290 if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )