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"
48 #include "SMDS_VolumeTool.hxx"
56 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 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 inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
73 return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
78 inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
80 double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
84 int getNbMultiConnection( SMDS_Mesh* theMesh, const int theId )
89 const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
90 if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge || anEdge->NbNodes() != 2 )
93 TColStd_MapOfInteger aMap;
96 SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
98 while( anIter->more() ) {
99 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
102 SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
103 while( anElemIter->more() ) {
104 const SMDS_MeshElement* anElem = anElemIter->next();
105 if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
106 int anId = anElem->GetID();
108 if ( anIter->more() ) // i.e. first node
110 else if ( aMap.Contains( anId ) )
124 using namespace SMESH::Controls;
131 Class : NumericalFunctor
132 Description : Base class for numerical functors
134 NumericalFunctor::NumericalFunctor():
140 void NumericalFunctor::SetMesh( SMDS_Mesh* theMesh )
145 bool NumericalFunctor::GetPoints(const int theId,
146 TSequenceOfXYZ& theRes ) const
153 return GetPoints( myMesh->FindElement( theId ), theRes );
156 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
157 TSequenceOfXYZ& theRes )
164 // Get nodes of the element
165 SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
168 while( anIter->more() )
170 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
172 theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
180 long NumericalFunctor::GetPrecision() const
185 void NumericalFunctor::SetPrecision( const long thePrecision )
187 myPrecision = thePrecision;
190 double NumericalFunctor::GetValue( long theId )
193 if ( GetPoints( theId, P ))
195 double aVal = GetValue( P );
196 if ( myPrecision >= 0 )
198 double prec = pow( 10., (double)( myPrecision ) );
199 aVal = floor( aVal * prec + 0.5 ) / prec;
209 Description : Functor for calculation of minimum angle
212 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
218 double A0 = getAngle( P( 3 ), P( 1 ), P( 2 ) );
219 double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
220 double A2 = getAngle( P( 2 ), P( 3 ), P( 1 ) );
222 aMin = Min( A0, Min( A1, A2 ) );
224 else if ( P.size() == 4 )
226 double A0 = getAngle( P( 4 ), P( 1 ), P( 2 ) );
227 double A1 = getAngle( P( 1 ), P( 2 ), P( 3 ) );
228 double A2 = getAngle( P( 2 ), P( 3 ), P( 4 ) );
229 double A3 = getAngle( P( 3 ), P( 4 ), P( 1 ) );
231 aMin = Min( Min( A0, A1 ), Min( A2, A3 ) );
236 return aMin * 180 / PI;
239 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
241 const double aBestAngle = PI / nbNodes;
242 return ( fabs( aBestAngle - Value ));
245 SMDSAbs_ElementType MinimumAngle::GetType() const
253 Description : Functor for calculating aspect ratio
255 double AspectRatio::GetValue( const TSequenceOfXYZ& P )
257 int nbNodes = P.size();
259 if ( nbNodes != 3 && nbNodes != 4 )
262 // Compute lengths of the sides
264 double aLen[ nbNodes ];
265 for ( int i = 0; i < nbNodes - 1; i++ )
266 aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
267 aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
269 // Compute aspect ratio
273 double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
274 if ( anArea <= Precision::Confusion() )
276 double aMaxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
277 static double aCoef = sqrt( 3. ) / 4;
279 return aCoef * aMaxLen * aMaxLen / anArea;
283 double aMinLen = Min( Min( aLen[ 0 ], aLen[ 1 ] ), Min( aLen[ 2 ], aLen[ 3 ] ) );
284 if ( aMinLen <= Precision::Confusion() )
286 double aMaxLen = Max( Max( aLen[ 0 ], aLen[ 1 ] ), Max( aLen[ 2 ], aLen[ 3 ] ) );
288 return aMaxLen / aMinLen;
292 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
294 // the aspect ratio is in the range [1.0,infinity]
297 return Value / 1000.;
300 SMDSAbs_ElementType AspectRatio::GetType() const
307 Class : AspectRatio3D
308 Description : Functor for calculating aspect ratio
312 inline double getHalfPerimeter(double theTria[3]){
313 return (theTria[0] + theTria[1] + theTria[2])/2.0;
316 inline double getArea(double theHalfPerim, double theTria[3]){
317 return sqrt(theHalfPerim*
318 (theHalfPerim-theTria[0])*
319 (theHalfPerim-theTria[1])*
320 (theHalfPerim-theTria[2]));
323 inline double getVolume(double theLen[6]){
324 double a2 = theLen[0]*theLen[0];
325 double b2 = theLen[1]*theLen[1];
326 double c2 = theLen[2]*theLen[2];
327 double d2 = theLen[3]*theLen[3];
328 double e2 = theLen[4]*theLen[4];
329 double f2 = theLen[5]*theLen[5];
330 double P = 4.0*a2*b2*d2;
331 double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
332 double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
333 return sqrt(P-Q+R)/12.0;
336 inline double getVolume2(double theLen[6]){
337 double a2 = theLen[0]*theLen[0];
338 double b2 = theLen[1]*theLen[1];
339 double c2 = theLen[2]*theLen[2];
340 double d2 = theLen[3]*theLen[3];
341 double e2 = theLen[4]*theLen[4];
342 double f2 = theLen[5]*theLen[5];
344 double P = a2*e2*(b2+c2+d2+f2-a2-e2);
345 double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
346 double R = c2*d2*(a2+b2+e2+f2-c2-d2);
347 double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
349 return sqrt(P+Q+R-S)/12.0;
352 inline double getVolume(const TSequenceOfXYZ& P){
353 gp_Vec aVec1( P( 2 ) - P( 1 ) );
354 gp_Vec aVec2( P( 3 ) - P( 1 ) );
355 gp_Vec aVec3( P( 4 ) - P( 1 ) );
356 gp_Vec anAreaVec( aVec1 ^ aVec2 );
357 return Abs(aVec3 * anAreaVec) / 6.0;
360 inline double getMaxHeight(double theLen[6])
362 double aHeight = max(theLen[0],theLen[1]);
363 aHeight = max(aHeight,theLen[2]);
364 aHeight = max(aHeight,theLen[3]);
365 aHeight = max(aHeight,theLen[4]);
366 aHeight = max(aHeight,theLen[5]);
372 double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
374 double aQuality = 0.0;
375 int nbNodes = P.size();
379 getDistance(P( 1 ),P( 2 )), // a
380 getDistance(P( 2 ),P( 3 )), // b
381 getDistance(P( 3 ),P( 1 )), // c
382 getDistance(P( 2 ),P( 4 )), // d
383 getDistance(P( 3 ),P( 4 )), // e
384 getDistance(P( 1 ),P( 4 )) // f
386 double aTria[4][3] = {
387 {aLen[0],aLen[1],aLen[2]}, // abc
388 {aLen[0],aLen[3],aLen[5]}, // adf
389 {aLen[1],aLen[3],aLen[4]}, // bde
390 {aLen[2],aLen[4],aLen[5]} // cef
392 double aSumArea = 0.0;
393 double aHalfPerimeter = getHalfPerimeter(aTria[0]);
394 double anArea = getArea(aHalfPerimeter,aTria[0]);
396 aHalfPerimeter = getHalfPerimeter(aTria[1]);
397 anArea = getArea(aHalfPerimeter,aTria[1]);
399 aHalfPerimeter = getHalfPerimeter(aTria[2]);
400 anArea = getArea(aHalfPerimeter,aTria[2]);
402 aHalfPerimeter = getHalfPerimeter(aTria[3]);
403 anArea = getArea(aHalfPerimeter,aTria[3]);
405 double aVolume = getVolume(P);
406 //double aVolume = getVolume(aLen);
407 double aHeight = getMaxHeight(aLen);
408 static double aCoeff = sqrt(6.0)/36.0;
409 aQuality = aCoeff*aHeight*aSumArea/aVolume;
414 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
415 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
418 gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
419 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
422 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
423 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
426 gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
427 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
433 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
434 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
437 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
438 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
441 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
442 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
445 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
446 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
449 gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
450 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
453 gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
454 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
460 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
461 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
464 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
465 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
468 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
469 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
472 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
473 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
476 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
477 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
480 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
481 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
484 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
485 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
488 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
489 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
492 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
493 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
496 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
497 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
500 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
501 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
504 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
505 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
508 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
509 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
512 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
513 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
516 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
517 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
520 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
521 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
524 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
525 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
528 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
529 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
532 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
533 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
536 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
537 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
540 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
541 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
544 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
545 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
548 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
549 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
552 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
553 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
556 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
557 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
560 gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
561 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
564 gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
565 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
568 gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
569 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
572 gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
573 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
576 gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
577 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
580 gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
581 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
584 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
585 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
588 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
589 aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
597 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
599 // the aspect ratio is in the range [1.0,infinity]
602 return Value / 1000.;
605 SMDSAbs_ElementType AspectRatio3D::GetType() const
607 return SMDSAbs_Volume;
613 Description : Functor for calculating warping
615 double Warping::GetValue( const TSequenceOfXYZ& P )
620 gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4;
622 double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
623 double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
624 double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
625 double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
627 return Max( Max( A1, A2 ), Max( A3, A4 ) );
630 double Warping::ComputeA( const gp_XYZ& thePnt1,
631 const gp_XYZ& thePnt2,
632 const gp_XYZ& thePnt3,
633 const gp_XYZ& theG ) const
635 double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
636 double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
637 double L = Min( aLen1, aLen2 ) * 0.5;
638 if ( L < Precision::Confusion())
641 gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG;
642 gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG;
643 gp_XYZ N = GI.Crossed( GJ );
645 if ( N.Modulus() < gp::Resolution() )
650 double H = ( thePnt2 - theG ).Dot( N );
651 return asin( fabs( H / L ) ) * 180 / PI;
654 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
656 // the warp is in the range [0.0,PI/2]
657 // 0.0 = good (no warp)
658 // PI/2 = bad (face pliee)
662 SMDSAbs_ElementType Warping::GetType() const
670 Description : Functor for calculating taper
672 double Taper::GetValue( const TSequenceOfXYZ& P )
678 double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2;
679 double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2;
680 double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2;
681 double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2;
683 double JA = 0.25 * ( J1 + J2 + J3 + J4 );
684 if ( JA <= Precision::Confusion() )
687 double T1 = fabs( ( J1 - JA ) / JA );
688 double T2 = fabs( ( J2 - JA ) / JA );
689 double T3 = fabs( ( J3 - JA ) / JA );
690 double T4 = fabs( ( J4 - JA ) / JA );
692 return Max( Max( T1, T2 ), Max( T3, T4 ) );
695 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
697 // the taper is in the range [0.0,1.0]
698 // 0.0 = good (no taper)
699 // 1.0 = bad (les cotes opposes sont allignes)
703 SMDSAbs_ElementType Taper::GetType() const
711 Description : Functor for calculating skew in degrees
713 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
715 gp_XYZ p12 = ( p2 + p1 ) / 2;
716 gp_XYZ p23 = ( p3 + p2 ) / 2;
717 gp_XYZ p31 = ( p3 + p1 ) / 2;
719 gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
721 return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
724 double Skew::GetValue( const TSequenceOfXYZ& P )
726 if ( P.size() != 3 && P.size() != 4 )
730 static double PI2 = PI / 2;
733 double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
734 double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
735 double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
737 return Max( A0, Max( A1, A2 ) ) * 180 / PI;
741 gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2;
742 gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2;
743 gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2;
744 gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2;
746 gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
747 double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
748 ? 0 : fabs( PI2 - v1.Angle( v2 ) );
754 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
756 // the skew is in the range [0.0,PI/2].
762 SMDSAbs_ElementType Skew::GetType() const
770 Description : Functor for calculating area
772 double Area::GetValue( const TSequenceOfXYZ& P )
775 return getArea( P( 1 ), P( 2 ), P( 3 ) );
776 else if ( P.size() == 4 )
777 return getArea( P( 1 ), P( 2 ), P( 3 ) ) + getArea( P( 1 ), P( 3 ), P( 4 ) );
782 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
787 SMDSAbs_ElementType Area::GetType() const
795 Description : Functor for calculating length off edge
797 double Length::GetValue( const TSequenceOfXYZ& P )
799 return ( P.size() == 2 ? getDistance( P( 1 ), P( 2 ) ) : 0 );
802 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
807 SMDSAbs_ElementType Length::GetType() const
814 Description : Functor for calculating length of edge
817 double Length2D::GetValue( long theElementId)
821 if (GetPoints(theElementId,P)){
823 double aVal;// = GetValue( P );
824 const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
825 SMDSAbs_ElementType aType = aElem->GetType();
834 aVal = getDistance( P( 1 ), P( 2 ) );
838 if (len == 3){ // triangles
839 double L1 = getDistance(P( 1 ),P( 2 ));
840 double L2 = getDistance(P( 2 ),P( 3 ));
841 double L3 = getDistance(P( 3 ),P( 1 ));
842 aVal = Max(L1,Max(L2,L3));
845 else if (len == 4){ // quadrangles
846 double L1 = getDistance(P( 1 ),P( 2 ));
847 double L2 = getDistance(P( 2 ),P( 3 ));
848 double L3 = getDistance(P( 3 ),P( 4 ));
849 double L4 = getDistance(P( 4 ),P( 1 ));
850 aVal = Max(Max(L1,L2),Max(L3,L4));
854 if (len == 4){ // tetraidrs
855 double L1 = getDistance(P( 1 ),P( 2 ));
856 double L2 = getDistance(P( 2 ),P( 3 ));
857 double L3 = getDistance(P( 3 ),P( 1 ));
858 double L4 = getDistance(P( 1 ),P( 4 ));
859 double L5 = getDistance(P( 2 ),P( 4 ));
860 double L6 = getDistance(P( 3 ),P( 4 ));
861 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
864 else if (len == 5){ // piramids
865 double L1 = getDistance(P( 1 ),P( 2 ));
866 double L2 = getDistance(P( 2 ),P( 3 ));
867 double L3 = getDistance(P( 3 ),P( 1 ));
868 double L4 = getDistance(P( 4 ),P( 1 ));
869 double L5 = getDistance(P( 1 ),P( 5 ));
870 double L6 = getDistance(P( 2 ),P( 5 ));
871 double L7 = getDistance(P( 3 ),P( 5 ));
872 double L8 = getDistance(P( 4 ),P( 5 ));
874 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
875 aVal = Max(aVal,Max(L7,L8));
878 else if (len == 6){ // pentaidres
879 double L1 = getDistance(P( 1 ),P( 2 ));
880 double L2 = getDistance(P( 2 ),P( 3 ));
881 double L3 = getDistance(P( 3 ),P( 1 ));
882 double L4 = getDistance(P( 4 ),P( 5 ));
883 double L5 = getDistance(P( 5 ),P( 6 ));
884 double L6 = getDistance(P( 6 ),P( 4 ));
885 double L7 = getDistance(P( 1 ),P( 4 ));
886 double L8 = getDistance(P( 2 ),P( 5 ));
887 double L9 = getDistance(P( 3 ),P( 6 ));
889 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
890 aVal = Max(aVal,Max(Max(L7,L8),L9));
893 else if (len == 8){ // hexaider
894 double L1 = getDistance(P( 1 ),P( 2 ));
895 double L2 = getDistance(P( 2 ),P( 3 ));
896 double L3 = getDistance(P( 3 ),P( 4 ));
897 double L4 = getDistance(P( 4 ),P( 1 ));
898 double L5 = getDistance(P( 5 ),P( 6 ));
899 double L6 = getDistance(P( 6 ),P( 7 ));
900 double L7 = getDistance(P( 7 ),P( 8 ));
901 double L8 = getDistance(P( 8 ),P( 5 ));
902 double L9 = getDistance(P( 1 ),P( 5 ));
903 double L10= getDistance(P( 2 ),P( 6 ));
904 double L11= getDistance(P( 3 ),P( 7 ));
905 double L12= getDistance(P( 4 ),P( 8 ));
907 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
908 aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
909 aVal = Max(aVal,Max(L11,L12));
921 if ( myPrecision >= 0 )
923 double prec = pow( 10., (double)( myPrecision ) );
924 aVal = floor( aVal * prec + 0.5 ) / prec;
933 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
938 SMDSAbs_ElementType Length2D::GetType() const
943 Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
946 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
947 if(thePntId1 > thePntId2){
948 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
952 bool Length2D::Value::operator<(const Length2D::Value& x) const{
953 if(myPntId[0] < x.myPntId[0]) return true;
954 if(myPntId[0] == x.myPntId[0])
955 if(myPntId[1] < x.myPntId[1]) return true;
959 void Length2D::GetValues(TValues& theValues){
961 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
962 for(; anIter->more(); ){
963 const SMDS_MeshFace* anElem = anIter->next();
964 SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
969 const SMDS_MeshElement* aNode;
970 if(aNodesIter->more()){
971 aNode = aNodesIter->next();
972 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
973 P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
974 aNodeId[0] = aNodeId[1] = aNode->GetID();
977 for(; aNodesIter->more(); ){
978 aNode = aNodesIter->next();
979 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
980 long anId = aNode->GetID();
982 P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
984 aLength = P[1].Distance(P[2]);
986 Value aValue(aLength,aNodeId[1],anId);
989 theValues.insert(aValue);
992 aLength = P[0].Distance(P[1]);
994 Value aValue(aLength,aNodeId[0],aNodeId[1]);
995 theValues.insert(aValue);
1000 Class : MultiConnection
1001 Description : Functor for calculating number of faces conneted to the edge
1003 double MultiConnection::GetValue( const TSequenceOfXYZ& P )
1007 double MultiConnection::GetValue( long theId )
1009 return getNbMultiConnection( myMesh, theId );
1012 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
1017 SMDSAbs_ElementType MultiConnection::GetType() const
1019 return SMDSAbs_Edge;
1023 Class : MultiConnection2D
1024 Description : Functor for calculating number of faces conneted to the edge
1026 double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
1031 double MultiConnection2D::GetValue( long theElementId )
1036 if (GetPoints(theElementId,P)){
1038 const SMDS_MeshElement* anFaceElem = myMesh->FindElement( theElementId );
1039 SMDSAbs_ElementType aType = anFaceElem->GetType();
1043 TColStd_MapOfInteger aMap;
1051 if (len == 3){ // triangles
1052 int Nb[3] = {0,0,0};
1055 SMDS_ElemIteratorPtr anIter = anFaceElem->nodesIterator();
1056 if ( anIter != 0 ) {
1057 while( anIter->more() ) {
1058 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
1062 SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
1063 while( anElemIter->more() ) {
1064 const SMDS_MeshElement* anElem = anElemIter->next();
1065 if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
1066 int anId = anElem->GetID();
1068 if ( anIter->more() ) // i.e. first node
1070 else if ( aMap.Contains( anId ) ){
1074 else if ( anElem != 0 && anElem->GetType() == SMDSAbs_Edge ) i++;
1079 aResult = Max(Max(Nb[0],Nb[1]),Nb[2]);
1082 case SMDSAbs_Volume:
1087 return aResult;//getNbMultiConnection( myMesh, theId );
1090 double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1095 SMDSAbs_ElementType MultiConnection2D::GetType() const
1097 return SMDSAbs_Face;
1100 MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
1102 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
1103 if(thePntId1 > thePntId2){
1104 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
1108 bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const{
1109 if(myPntId[0] < x.myPntId[0]) return true;
1110 if(myPntId[0] == x.myPntId[0])
1111 if(myPntId[1] < x.myPntId[1]) return true;
1115 void MultiConnection2D::GetValues(MValues& theValues){
1116 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1117 for(; anIter->more(); ){
1118 const SMDS_MeshFace* anElem = anIter->next();
1119 long anElemId = anElem->GetID();
1120 SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
1123 //int aNbConnects=0;
1124 const SMDS_MeshNode* aNode0;
1125 const SMDS_MeshNode* aNode1;
1126 const SMDS_MeshNode* aNode2;
1127 if(aNodesIter->more()){
1128 aNode0 = (SMDS_MeshNode*) aNodesIter->next();
1130 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
1131 aNodeId[0] = aNodeId[1] = aNodes->GetID();
1133 for(; aNodesIter->more(); ){
1134 aNode2 = (SMDS_MeshNode*) aNodesIter->next();
1135 long anId = aNode2->GetID();
1138 Value aValue(aNodeId[1],aNodeId[2]);
1139 MValues::iterator aItr = theValues.find(aValue);
1140 if (aItr != theValues.end()){
1144 theValues[aValue] = 1;
1147 //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
1148 aNodeId[1] = aNodeId[2];
1151 Value aValue(aNodeId[0],aNodeId[2]);
1152 MValues::iterator aItr = theValues.find(aValue);
1153 if (aItr != theValues.end()){
1157 theValues[aValue] = 1;
1160 //cout << "NodeIds: "<<aNodeId[0]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
1170 Class : BadOrientedVolume
1171 Description : Predicate bad oriented volumes
1174 BadOrientedVolume::BadOrientedVolume()
1179 void BadOrientedVolume::SetMesh( SMDS_Mesh* theMesh )
1184 bool BadOrientedVolume::IsSatisfy( long theId )
1189 SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
1190 return !vTool.IsForward();
1193 SMDSAbs_ElementType BadOrientedVolume::GetType() const
1195 return SMDSAbs_Volume;
1202 Description : Predicate for free borders
1205 FreeBorders::FreeBorders()
1210 void FreeBorders::SetMesh( SMDS_Mesh* theMesh )
1215 bool FreeBorders::IsSatisfy( long theId )
1217 return getNbMultiConnection( myMesh, theId ) == 1;
1220 SMDSAbs_ElementType FreeBorders::GetType() const
1222 return SMDSAbs_Edge;
1228 Description : Predicate for free Edges
1230 FreeEdges::FreeEdges()
1235 void FreeEdges::SetMesh( SMDS_Mesh* theMesh )
1240 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId )
1242 TColStd_MapOfInteger aMap;
1243 for ( int i = 0; i < 2; i++ )
1245 SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator();
1246 while( anElemIter->more() )
1248 const SMDS_MeshElement* anElem = anElemIter->next();
1249 if ( anElem != 0 && anElem->GetType() == SMDSAbs_Face )
1251 int anId = anElem->GetID();
1255 else if ( aMap.Contains( anId ) && anId != theFaceId )
1263 bool FreeEdges::IsSatisfy( long theId )
1268 const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
1269 if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
1272 int nbNodes = aFace->NbNodes();
1273 const SMDS_MeshNode* aNodes[ nbNodes ];
1275 SMDS_ElemIteratorPtr anIter = aFace->nodesIterator();
1278 while( anIter->more() )
1280 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
1283 aNodes[ i++ ] = aNode;
1287 for ( int i = 0; i < nbNodes - 1; i++ )
1288 if ( IsFreeEdge( &aNodes[ i ], theId ) )
1291 aNodes[ 1 ] = aNodes[ nbNodes - 1 ];
1293 return IsFreeEdge( &aNodes[ 0 ], theId );
1297 SMDSAbs_ElementType FreeEdges::GetType() const
1299 return SMDSAbs_Face;
1302 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
1305 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
1306 if(thePntId1 > thePntId2){
1307 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
1311 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
1312 if(myPntId[0] < x.myPntId[0]) return true;
1313 if(myPntId[0] == x.myPntId[0])
1314 if(myPntId[1] < x.myPntId[1]) return true;
1318 inline void UpdateBorders(const FreeEdges::Border& theBorder,
1319 FreeEdges::TBorders& theRegistry,
1320 FreeEdges::TBorders& theContainer)
1322 if(theRegistry.find(theBorder) == theRegistry.end()){
1323 theRegistry.insert(theBorder);
1324 theContainer.insert(theBorder);
1326 theContainer.erase(theBorder);
1330 void FreeEdges::GetBoreders(TBorders& theBorders)
1333 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1334 for(; anIter->more(); ){
1335 const SMDS_MeshFace* anElem = anIter->next();
1336 long anElemId = anElem->GetID();
1337 SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
1339 const SMDS_MeshElement* aNode;
1340 if(aNodesIter->more()){
1341 aNode = aNodesIter->next();
1342 aNodeId[0] = aNodeId[1] = aNode->GetID();
1344 for(; aNodesIter->more(); ){
1345 aNode = aNodesIter->next();
1346 long anId = aNode->GetID();
1347 Border aBorder(anElemId,aNodeId[1],anId);
1349 //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
1350 UpdateBorders(aBorder,aRegistry,theBorders);
1352 Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
1353 //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
1354 UpdateBorders(aBorder,aRegistry,theBorders);
1356 //std::cout<<"theBorders.size() = "<<theBorders.size()<<endl;
1361 Description : Predicate for Range of Ids.
1362 Range may be specified with two ways.
1363 1. Using AddToRange method
1364 2. With SetRangeStr method. Parameter of this method is a string
1365 like as "1,2,3,50-60,63,67,70-"
1368 //=======================================================================
1369 // name : RangeOfIds
1370 // Purpose : Constructor
1371 //=======================================================================
1372 RangeOfIds::RangeOfIds()
1375 myType = SMDSAbs_All;
1378 //=======================================================================
1380 // Purpose : Set mesh
1381 //=======================================================================
1382 void RangeOfIds::SetMesh( SMDS_Mesh* theMesh )
1387 //=======================================================================
1388 // name : AddToRange
1389 // Purpose : Add ID to the range
1390 //=======================================================================
1391 bool RangeOfIds::AddToRange( long theEntityId )
1393 myIds.Add( theEntityId );
1397 //=======================================================================
1398 // name : GetRangeStr
1399 // Purpose : Get range as a string.
1400 // Example: "1,2,3,50-60,63,67,70-"
1401 //=======================================================================
1402 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
1406 TColStd_SequenceOfInteger anIntSeq;
1407 TColStd_SequenceOfAsciiString aStrSeq;
1409 TColStd_MapIteratorOfMapOfInteger anIter( myIds );
1410 for ( ; anIter.More(); anIter.Next() )
1412 int anId = anIter.Key();
1413 TCollection_AsciiString aStr( anId );
1414 anIntSeq.Append( anId );
1415 aStrSeq.Append( aStr );
1418 for ( int i = 1, n = myMin.Length(); i <= n; i++ )
1420 int aMinId = myMin( i );
1421 int aMaxId = myMax( i );
1423 TCollection_AsciiString aStr;
1424 if ( aMinId != IntegerFirst() )
1429 if ( aMaxId != IntegerLast() )
1432 // find position of the string in result sequence and insert string in it
1433 if ( anIntSeq.Length() == 0 )
1435 anIntSeq.Append( aMinId );
1436 aStrSeq.Append( aStr );
1440 if ( aMinId < anIntSeq.First() )
1442 anIntSeq.Prepend( aMinId );
1443 aStrSeq.Prepend( aStr );
1445 else if ( aMinId > anIntSeq.Last() )
1447 anIntSeq.Append( aMinId );
1448 aStrSeq.Append( aStr );
1451 for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
1452 if ( aMinId < anIntSeq( j ) )
1454 anIntSeq.InsertBefore( j, aMinId );
1455 aStrSeq.InsertBefore( j, aStr );
1461 if ( aStrSeq.Length() == 0 )
1464 theResStr = aStrSeq( 1 );
1465 for ( int j = 2, k = aStrSeq.Length(); j <= k; j++ )
1468 theResStr += aStrSeq( j );
1472 //=======================================================================
1473 // name : SetRangeStr
1474 // Purpose : Define range with string
1475 // Example of entry string: "1,2,3,50-60,63,67,70-"
1476 //=======================================================================
1477 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
1483 TCollection_AsciiString aStr = theStr;
1484 aStr.RemoveAll( ' ' );
1485 aStr.RemoveAll( '\t' );
1487 for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) )
1488 aStr.Remove( aPos, 2 );
1490 TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
1492 while ( tmpStr != "" )
1494 tmpStr = aStr.Token( ",", i++ );
1495 int aPos = tmpStr.Search( '-' );
1499 if ( tmpStr.IsIntegerValue() )
1500 myIds.Add( tmpStr.IntegerValue() );
1506 TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
1507 TCollection_AsciiString aMinStr = tmpStr;
1509 while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
1510 while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
1512 if ( !aMinStr.IsEmpty() && !aMinStr.IsIntegerValue() ||
1513 !aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue() )
1516 myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
1517 myMax.Append( aMaxStr.IsEmpty() ? IntegerLast() : aMaxStr.IntegerValue() );
1524 //=======================================================================
1526 // Purpose : Get type of supported entities
1527 //=======================================================================
1528 SMDSAbs_ElementType RangeOfIds::GetType() const
1533 //=======================================================================
1535 // Purpose : Set type of supported entities
1536 //=======================================================================
1537 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
1542 //=======================================================================
1544 // Purpose : Verify whether entity satisfies to this rpedicate
1545 //=======================================================================
1546 bool RangeOfIds::IsSatisfy( long theId )
1551 if ( myType == SMDSAbs_Node )
1553 if ( myMesh->FindNode( theId ) == 0 )
1558 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
1559 if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All )
1563 if ( myIds.Contains( theId ) )
1566 for ( int i = 1, n = myMin.Length(); i <= n; i++ )
1567 if ( theId >= myMin( i ) && theId <= myMax( i ) )
1575 Description : Base class for comparators
1577 Comparator::Comparator():
1581 Comparator::~Comparator()
1584 void Comparator::SetMesh( SMDS_Mesh* theMesh )
1587 myFunctor->SetMesh( theMesh );
1590 void Comparator::SetMargin( double theValue )
1592 myMargin = theValue;
1595 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
1597 myFunctor = theFunct;
1600 SMDSAbs_ElementType Comparator::GetType() const
1602 return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
1605 double Comparator::GetMargin()
1613 Description : Comparator "<"
1615 bool LessThan::IsSatisfy( long theId )
1617 return myFunctor && myFunctor->GetValue( theId ) < myMargin;
1623 Description : Comparator ">"
1625 bool MoreThan::IsSatisfy( long theId )
1627 return myFunctor && myFunctor->GetValue( theId ) > myMargin;
1633 Description : Comparator "="
1636 myToler(Precision::Confusion())
1639 bool EqualTo::IsSatisfy( long theId )
1641 return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
1644 void EqualTo::SetTolerance( double theToler )
1649 double EqualTo::GetTolerance()
1656 Description : Logical NOT predicate
1658 LogicalNOT::LogicalNOT()
1661 LogicalNOT::~LogicalNOT()
1664 bool LogicalNOT::IsSatisfy( long theId )
1666 return myPredicate && !myPredicate->IsSatisfy( theId );
1669 void LogicalNOT::SetMesh( SMDS_Mesh* theMesh )
1672 myPredicate->SetMesh( theMesh );
1675 void LogicalNOT::SetPredicate( PredicatePtr thePred )
1677 myPredicate = thePred;
1680 SMDSAbs_ElementType LogicalNOT::GetType() const
1682 return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
1687 Class : LogicalBinary
1688 Description : Base class for binary logical predicate
1690 LogicalBinary::LogicalBinary()
1693 LogicalBinary::~LogicalBinary()
1696 void LogicalBinary::SetMesh( SMDS_Mesh* theMesh )
1699 myPredicate1->SetMesh( theMesh );
1702 myPredicate2->SetMesh( theMesh );
1705 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
1707 myPredicate1 = thePredicate;
1710 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
1712 myPredicate2 = thePredicate;
1715 SMDSAbs_ElementType LogicalBinary::GetType() const
1717 if ( !myPredicate1 || !myPredicate2 )
1720 SMDSAbs_ElementType aType1 = myPredicate1->GetType();
1721 SMDSAbs_ElementType aType2 = myPredicate2->GetType();
1723 return aType1 == aType2 ? aType1 : SMDSAbs_All;
1729 Description : Logical AND
1731 bool LogicalAND::IsSatisfy( long theId )
1736 myPredicate1->IsSatisfy( theId ) &&
1737 myPredicate2->IsSatisfy( theId );
1743 Description : Logical OR
1745 bool LogicalOR::IsSatisfy( long theId )
1750 myPredicate1->IsSatisfy( theId ) ||
1751 myPredicate2->IsSatisfy( theId );
1765 void Filter::SetPredicate( PredicatePtr thePredicate )
1767 myPredicate = thePredicate;
1771 template<class TElement, class TIterator, class TPredicate>
1772 void FillSequence(const TIterator& theIterator,
1773 TPredicate& thePredicate,
1774 Filter::TIdSequence& theSequence)
1776 if ( theIterator ) {
1777 while( theIterator->more() ) {
1778 TElement anElem = theIterator->next();
1779 long anId = anElem->GetID();
1780 if ( thePredicate->IsSatisfy( anId ) )
1781 theSequence.push_back( anId );
1787 Filter::GetElementsId( SMDS_Mesh* theMesh )
1789 TIdSequence aSequence;
1790 if ( !theMesh || !myPredicate ) return aSequence;
1792 myPredicate->SetMesh( theMesh );
1794 SMDSAbs_ElementType aType = myPredicate->GetType();
1797 FillSequence<const SMDS_MeshNode*>(theMesh->nodesIterator(),myPredicate,aSequence);
1801 FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
1805 FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
1808 case SMDSAbs_Volume:{
1809 FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
1813 FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),myPredicate,aSequence);
1814 FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),myPredicate,aSequence);
1815 FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),myPredicate,aSequence);
1826 typedef std::set<SMDS_MeshFace*> TMapOfFacePtr;
1832 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
1833 SMDS_MeshNode* theNode2 )
1839 ManifoldPart::Link::~Link()
1845 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
1847 if ( myNode1 == theLink.myNode1 &&
1848 myNode2 == theLink.myNode2 )
1850 else if ( myNode1 == theLink.myNode2 &&
1851 myNode2 == theLink.myNode1 )
1857 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
1859 if(myNode1 < x.myNode1) return true;
1860 if(myNode1 == x.myNode1)
1861 if(myNode2 < x.myNode2) return true;
1865 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
1866 const ManifoldPart::Link& theLink2 )
1868 return theLink1.IsEqual( theLink2 );
1871 ManifoldPart::ManifoldPart()
1874 myAngToler = Precision::Angular();
1875 myIsOnlyManifold = true;
1878 ManifoldPart::~ManifoldPart()
1883 void ManifoldPart::SetMesh( SMDS_Mesh* theMesh )
1889 SMDSAbs_ElementType ManifoldPart::GetType() const
1890 { return SMDSAbs_Face; }
1892 bool ManifoldPart::IsSatisfy( long theElementId )
1894 return myMapIds.Contains( theElementId );
1897 void ManifoldPart::SetAngleTolerance( const double theAngToler )
1898 { myAngToler = theAngToler; }
1900 double ManifoldPart::GetAngleTolerance() const
1901 { return myAngToler; }
1903 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
1904 { myIsOnlyManifold = theIsOnly; }
1906 void ManifoldPart::SetStartElem( const long theStartId )
1907 { myStartElemId = theStartId; }
1909 bool ManifoldPart::process()
1912 myMapBadGeomIds.Clear();
1914 myAllFacePtr.clear();
1915 myAllFacePtrIntDMap.clear();
1919 // collect all faces into own map
1920 SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
1921 for (; anFaceItr->more(); )
1923 SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
1924 myAllFacePtr.push_back( aFacePtr );
1925 myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
1928 SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
1932 // the map of non manifold links and bad geometry
1933 TMapOfLink aMapOfNonManifold;
1934 TColStd_MapOfInteger aMapOfTreated;
1936 // begin cycle on faces from start index and run on vector till the end
1937 // and from begin to start index to cover whole vector
1938 const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
1939 bool isStartTreat = false;
1940 for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
1942 if ( fi == aStartIndx )
1943 isStartTreat = true;
1944 // as result next time when fi will be equal to aStartIndx
1946 SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
1947 if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
1950 aMapOfTreated.Add( aFacePtr->GetID() );
1951 TColStd_MapOfInteger aResFaces;
1952 if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
1953 aMapOfNonManifold, aResFaces ) )
1955 TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
1956 for ( ; anItr.More(); anItr.Next() )
1958 int aFaceId = anItr.Key();
1959 aMapOfTreated.Add( aFaceId );
1960 myMapIds.Add( aFaceId );
1963 if ( fi == ( myAllFacePtr.size() - 1 ) )
1965 } // end run on vector of faces
1966 return !myMapIds.IsEmpty();
1969 static void getLinks( const SMDS_MeshFace* theFace,
1970 ManifoldPart::TVectorOfLink& theLinks )
1972 int aNbNode = theFace->NbNodes();
1973 SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
1975 SMDS_MeshNode* aNode = 0;
1976 for ( ; aNodeItr->more() && i <= aNbNode; )
1979 SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
1983 SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
1985 ManifoldPart::Link aLink( aN1, aN2 );
1986 theLinks.push_back( aLink );
1990 static gp_XYZ getNormale( const SMDS_MeshFace* theFace )
1993 int aNbNode = theFace->NbNodes();
1994 TColgp_Array1OfXYZ anArrOfXYZ(1,4);
1995 SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
1997 for ( ; aNodeItr->more() && i <= 4; i++ )
1999 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
2000 anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
2003 gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1);
2004 gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1);
2008 gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1);
2011 double len = n.Modulus();
2018 bool ManifoldPart::findConnected
2019 ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
2020 SMDS_MeshFace* theStartFace,
2021 ManifoldPart::TMapOfLink& theNonManifold,
2022 TColStd_MapOfInteger& theResFaces )
2024 theResFaces.Clear();
2025 if ( !theAllFacePtrInt.size() )
2028 if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
2030 myMapBadGeomIds.Add( theStartFace->GetID() );
2034 ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
2035 ManifoldPart::TVectorOfLink aSeqOfBoundary;
2036 theResFaces.Add( theStartFace->GetID() );
2037 ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
2039 expandBoundary( aMapOfBoundary, aSeqOfBoundary,
2040 aDMapLinkFace, theNonManifold, theStartFace );
2042 bool isDone = false;
2043 while ( !isDone && aMapOfBoundary.size() != 0 )
2045 bool isToReset = false;
2046 ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
2047 for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
2049 ManifoldPart::Link aLink = *pLink;
2050 if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
2052 // each link could be treated only once
2053 aMapToSkip.insert( aLink );
2055 ManifoldPart::TVectorOfFacePtr aFaces;
2057 if ( myIsOnlyManifold &&
2058 (theNonManifold.find( aLink ) != theNonManifold.end()) )
2062 getFacesByLink( aLink, aFaces );
2063 // filter the element to keep only indicated elements
2064 ManifoldPart::TVectorOfFacePtr aFiltered;
2065 ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
2066 for ( ; pFace != aFaces.end(); ++pFace )
2068 SMDS_MeshFace* aFace = *pFace;
2069 if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
2070 aFiltered.push_back( aFace );
2073 if ( aFaces.size() < 2 ) // no neihgbour faces
2075 else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
2077 theNonManifold.insert( aLink );
2082 // compare normal with normals of neighbor element
2083 SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
2084 ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
2085 for ( ; pFace != aFaces.end(); ++pFace )
2087 SMDS_MeshFace* aNextFace = *pFace;
2088 if ( aPrevFace == aNextFace )
2090 int anNextFaceID = aNextFace->GetID();
2091 if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
2092 // should not be with non manifold restriction. probably bad topology
2094 // check if face was treated and skipped
2095 if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
2096 !isInPlane( aPrevFace, aNextFace ) )
2098 // add new element to connected and extend the boundaries.
2099 theResFaces.Add( anNextFaceID );
2100 expandBoundary( aMapOfBoundary, aSeqOfBoundary,
2101 aDMapLinkFace, theNonManifold, aNextFace );
2105 isDone = !isToReset;
2108 return !theResFaces.IsEmpty();
2111 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
2112 const SMDS_MeshFace* theFace2 )
2114 gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
2115 gp_XYZ aNorm2XYZ = getNormale( theFace2 );
2116 if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
2118 myMapBadGeomIds.Add( theFace2->GetID() );
2121 if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
2127 void ManifoldPart::expandBoundary
2128 ( ManifoldPart::TMapOfLink& theMapOfBoundary,
2129 ManifoldPart::TVectorOfLink& theSeqOfBoundary,
2130 ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
2131 ManifoldPart::TMapOfLink& theNonManifold,
2132 SMDS_MeshFace* theNextFace ) const
2134 ManifoldPart::TVectorOfLink aLinks;
2135 getLinks( theNextFace, aLinks );
2136 int aNbLink = aLinks.size();
2137 for ( int i = 0; i < aNbLink; i++ )
2139 ManifoldPart::Link aLink = aLinks[ i ];
2140 if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
2142 if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
2144 if ( myIsOnlyManifold )
2146 // remove from boundary
2147 theMapOfBoundary.erase( aLink );
2148 ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
2149 for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
2151 ManifoldPart::Link aBoundLink = *pLink;
2152 if ( aBoundLink.IsEqual( aLink ) )
2154 theSeqOfBoundary.erase( pLink );
2162 theMapOfBoundary.insert( aLink );
2163 theSeqOfBoundary.push_back( aLink );
2164 theDMapLinkFacePtr[ aLink ] = theNextFace;
2169 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
2170 ManifoldPart::TVectorOfFacePtr& theFaces ) const
2172 SMDS_Mesh::SetOfFaces aSetOfFaces;
2173 // take all faces that shared first node
2174 SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
2175 for ( ; anItr->more(); )
2177 SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
2180 aSetOfFaces.Add( aFace );
2182 // take all faces that shared second node
2183 anItr = theLink.myNode2->facesIterator();
2184 // find the common part of two sets
2185 for ( ; anItr->more(); )
2187 SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
2188 if ( aSetOfFaces.Contains( aFace ) )
2189 theFaces.push_back( aFace );
2198 ElementsOnSurface::ElementsOnSurface()
2202 myType = SMDSAbs_All;
2204 myToler = Precision::Confusion();
2207 ElementsOnSurface::~ElementsOnSurface()
2212 void ElementsOnSurface::SetMesh( SMDS_Mesh* theMesh )
2214 if ( myMesh == theMesh )
2221 bool ElementsOnSurface::IsSatisfy( long theElementId )
2223 return myIds.Contains( theElementId );
2226 SMDSAbs_ElementType ElementsOnSurface::GetType() const
2229 void ElementsOnSurface::SetTolerance( const double theToler )
2230 { myToler = theToler; }
2232 double ElementsOnSurface::GetTolerance() const
2237 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
2238 const SMDSAbs_ElementType theType )
2242 if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
2247 TopoDS_Face aFace = TopoDS::Face( theShape );
2248 mySurf = BRep_Tool::Surface( aFace );
2251 void ElementsOnSurface::process()
2254 if ( mySurf.IsNull() )
2260 if ( myType == SMDSAbs_Face || myType == SMDSAbs_All )
2262 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
2263 for(; anIter->more(); )
2264 process( anIter->next() );
2267 if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All )
2269 SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator();
2270 for(; anIter->more(); )
2271 process( anIter->next() );
2274 if ( myType == SMDSAbs_Node )
2276 SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator();
2277 for(; anIter->more(); )
2278 process( anIter->next() );
2282 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
2284 SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
2285 bool isSatisfy = true;
2286 for ( ; aNodeItr->more(); )
2288 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
2289 if ( !isOnSurface( aNode ) )
2296 myIds.Add( theElemPtr->GetID() );
2299 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) const
2301 if ( mySurf.IsNull() )
2304 gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
2305 double aToler2 = myToler * myToler;
2306 if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
2308 gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
2309 if ( aPln.SquareDistance( aPnt ) > aToler2 )
2312 else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
2314 gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
2315 double aRad = aCyl.Radius();
2316 gp_Ax3 anAxis = aCyl.Position();
2317 gp_XYZ aLoc = aCyl.Location().XYZ();
2318 double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
2319 double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
2320 if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )