1 // Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 #include "SMESH_ControlsDef.hxx"
25 #include "SMDS_BallElement.hxx"
26 #include "SMDS_Iterator.hxx"
27 #include "SMDS_Mesh.hxx"
28 #include "SMDS_MeshElement.hxx"
29 #include "SMDS_MeshNode.hxx"
30 #include "SMDS_QuadraticEdge.hxx"
31 #include "SMDS_QuadraticFaceOfNodes.hxx"
32 #include "SMDS_VolumeTool.hxx"
33 #include "SMESHDS_GroupBase.hxx"
34 #include "SMESHDS_Mesh.hxx"
35 #include "SMESH_OctreeNode.hxx"
36 #include "SMESH_MeshAlgos.hxx"
38 #include <Basics_Utils.hxx>
40 #include <BRepAdaptor_Surface.hxx>
41 #include <BRepClass_FaceClassifier.hxx>
42 #include <BRep_Tool.hxx>
43 #include <Geom_CylindricalSurface.hxx>
44 #include <Geom_Plane.hxx>
45 #include <Geom_Surface.hxx>
46 #include <Precision.hxx>
47 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
48 #include <TColStd_MapOfInteger.hxx>
49 #include <TColStd_SequenceOfAsciiString.hxx>
50 #include <TColgp_Array1OfXYZ.hxx>
54 #include <TopoDS_Edge.hxx>
55 #include <TopoDS_Face.hxx>
56 #include <TopoDS_Iterator.hxx>
57 #include <TopoDS_Shape.hxx>
58 #include <TopoDS_Vertex.hxx>
60 #include <gp_Cylinder.hxx>
67 #include <vtkMeshQuality.h>
71 #include <TopTools_MapOfShape.hxx>
79 const double theEps = 1e-100;
80 const double theInf = 1e+100;
82 inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode )
84 return gp_XYZ(aNode->X(), aNode->Y(), aNode->Z() );
87 inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
89 gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
91 return v1.Magnitude() < gp::Resolution() ||
92 v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
95 inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
97 gp_Vec aVec1( P2 - P1 );
98 gp_Vec aVec2( P3 - P1 );
99 return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
102 inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
104 return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
109 inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
111 double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
115 int getNbMultiConnection( const SMDS_Mesh* theMesh, const int theId )
120 const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
121 if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge/* || anEdge->NbNodes() != 2 */)
124 // for each pair of nodes in anEdge (there are 2 pairs in a quadratic edge)
125 // count elements containing both nodes of the pair.
126 // Note that there may be such cases for a quadratic edge (a horizontal line):
131 // +-----+------+ +-----+------+
134 // result sould be 2 in both cases
136 int aResult0 = 0, aResult1 = 0;
137 // last node, it is a medium one in a quadratic edge
138 const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 );
139 const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 );
140 const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 );
141 if ( aNode1 == aLastNode ) aNode1 = 0;
143 SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator();
144 while( anElemIter->more() ) {
145 const SMDS_MeshElement* anElem = anElemIter->next();
146 if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
147 SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
148 while ( anIter->more() ) {
149 if ( const SMDS_MeshElement* anElemNode = anIter->next() ) {
150 if ( anElemNode == aNode0 ) {
152 if ( !aNode1 ) break; // not a quadratic edge
154 else if ( anElemNode == aNode1 )
160 int aResult = std::max ( aResult0, aResult1 );
162 // TColStd_MapOfInteger aMap;
164 // SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
165 // if ( anIter != 0 ) {
166 // while( anIter->more() ) {
167 // const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
170 // SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
171 // while( anElemIter->more() ) {
172 // const SMDS_MeshElement* anElem = anElemIter->next();
173 // if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
174 // int anId = anElem->GetID();
176 // if ( anIter->more() ) // i.e. first node
178 // else if ( aMap.Contains( anId ) )
188 gp_XYZ getNormale( const SMDS_MeshFace* theFace, bool* ok=0 )
190 int aNbNode = theFace->NbNodes();
192 gp_XYZ q1 = gpXYZ( theFace->GetNode(1)) - gpXYZ( theFace->GetNode(0));
193 gp_XYZ q2 = gpXYZ( theFace->GetNode(2)) - gpXYZ( theFace->GetNode(0));
196 gp_XYZ q3 = gpXYZ( theFace->GetNode(3)) - gpXYZ( theFace->GetNode(0));
199 double len = n.Modulus();
200 bool zeroLen = ( len <= numeric_limits<double>::min());
204 if (ok) *ok = !zeroLen;
212 using namespace SMESH::Controls;
218 //================================================================================
220 Class : NumericalFunctor
221 Description : Base class for numerical functors
223 //================================================================================
225 NumericalFunctor::NumericalFunctor():
231 void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh )
236 bool NumericalFunctor::GetPoints(const int theId,
237 TSequenceOfXYZ& theRes ) const
244 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
245 if ( !anElem || anElem->GetType() != this->GetType() )
248 return GetPoints( anElem, theRes );
251 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
252 TSequenceOfXYZ& theRes )
259 theRes.reserve( anElem->NbNodes() );
261 // Get nodes of the element
262 SMDS_ElemIteratorPtr anIter;
264 if ( anElem->IsQuadratic() ) {
265 switch ( anElem->GetType() ) {
267 anIter = dynamic_cast<const SMDS_VtkEdge*>
268 (anElem)->interlacedNodesElemIterator();
271 anIter = dynamic_cast<const SMDS_VtkFace*>
272 (anElem)->interlacedNodesElemIterator();
275 anIter = anElem->nodesIterator();
280 anIter = anElem->nodesIterator();
284 while( anIter->more() ) {
285 if ( const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( anIter->next() ))
286 theRes.push_back( gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
293 long NumericalFunctor::GetPrecision() const
298 void NumericalFunctor::SetPrecision( const long thePrecision )
300 myPrecision = thePrecision;
301 myPrecisionValue = pow( 10., (double)( myPrecision ) );
304 double NumericalFunctor::GetValue( long theId )
308 myCurrElement = myMesh->FindElement( theId );
311 if ( GetPoints( theId, P )) // elem type is checked here
312 aVal = Round( GetValue( P ));
317 double NumericalFunctor::Round( const double & aVal )
319 return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal;
322 //================================================================================
324 * \brief Return histogram of functor values
325 * \param nbIntervals - number of intervals
326 * \param nbEvents - number of mesh elements having values within i-th interval
327 * \param funValues - boundaries of intervals
328 * \param elements - elements to check vulue of; empty list means "of all"
329 * \param minmax - boundaries of diapason of values to divide into intervals
331 //================================================================================
333 void NumericalFunctor::GetHistogram(int nbIntervals,
334 std::vector<int>& nbEvents,
335 std::vector<double>& funValues,
336 const vector<int>& elements,
337 const double* minmax,
338 const bool isLogarithmic)
340 if ( nbIntervals < 1 ||
342 !myMesh->GetMeshInfo().NbElements( GetType() ))
344 nbEvents.resize( nbIntervals, 0 );
345 funValues.resize( nbIntervals+1 );
347 // get all values sorted
348 std::multiset< double > values;
349 if ( elements.empty() )
351 SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator(GetType());
352 while ( elemIt->more() )
353 values.insert( GetValue( elemIt->next()->GetID() ));
357 vector<int>::const_iterator id = elements.begin();
358 for ( ; id != elements.end(); ++id )
359 values.insert( GetValue( *id ));
364 funValues[0] = minmax[0];
365 funValues[nbIntervals] = minmax[1];
369 funValues[0] = *values.begin();
370 funValues[nbIntervals] = *values.rbegin();
372 // case nbIntervals == 1
373 if ( nbIntervals == 1 )
375 nbEvents[0] = values.size();
379 if (funValues.front() == funValues.back())
381 nbEvents.resize( 1 );
382 nbEvents[0] = values.size();
383 funValues[1] = funValues.back();
384 funValues.resize( 2 );
387 std::multiset< double >::iterator min = values.begin(), max;
388 for ( int i = 0; i < nbIntervals; ++i )
390 // find end value of i-th interval
391 double r = (i+1) / double(nbIntervals);
392 if (isLogarithmic && funValues.front() > 1e-07 && funValues.back() > 1e-07) {
393 double logmin = log10(funValues.front());
394 double lval = logmin + r * (log10(funValues.back()) - logmin);
395 funValues[i+1] = pow(10.0, lval);
398 funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
401 // count values in the i-th interval if there are any
402 if ( min != values.end() && *min <= funValues[i+1] )
404 // find the first value out of the interval
405 max = values.upper_bound( funValues[i+1] ); // max is greater than funValues[i+1], or end()
406 nbEvents[i] = std::distance( min, max );
410 // add values larger than minmax[1]
411 nbEvents.back() += std::distance( min, values.end() );
414 //=======================================================================
417 Description : Functor calculating volume of a 3D element
419 //================================================================================
421 double Volume::GetValue( long theElementId )
423 if ( theElementId && myMesh ) {
424 SMDS_VolumeTool aVolumeTool;
425 if ( aVolumeTool.Set( myMesh->FindElement( theElementId )))
426 return aVolumeTool.GetSize();
431 double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const
436 SMDSAbs_ElementType Volume::GetType() const
438 return SMDSAbs_Volume;
441 //=======================================================================
443 Class : MaxElementLength2D
444 Description : Functor calculating maximum length of 2D element
446 //================================================================================
448 double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P )
454 if( len == 3 ) { // triangles
455 double L1 = getDistance(P( 1 ),P( 2 ));
456 double L2 = getDistance(P( 2 ),P( 3 ));
457 double L3 = getDistance(P( 3 ),P( 1 ));
458 aVal = Max(L1,Max(L2,L3));
460 else if( len == 4 ) { // quadrangles
461 double L1 = getDistance(P( 1 ),P( 2 ));
462 double L2 = getDistance(P( 2 ),P( 3 ));
463 double L3 = getDistance(P( 3 ),P( 4 ));
464 double L4 = getDistance(P( 4 ),P( 1 ));
465 double D1 = getDistance(P( 1 ),P( 3 ));
466 double D2 = getDistance(P( 2 ),P( 4 ));
467 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
469 else if( len == 6 ) { // quadratic triangles
470 double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
471 double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
472 double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
473 aVal = Max(L1,Max(L2,L3));
475 else if( len == 8 || len == 9 ) { // quadratic quadrangles
476 double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
477 double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
478 double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
479 double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
480 double D1 = getDistance(P( 1 ),P( 5 ));
481 double D2 = getDistance(P( 3 ),P( 7 ));
482 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
485 if( myPrecision >= 0 )
487 double prec = pow( 10., (double)myPrecision );
488 aVal = floor( aVal * prec + 0.5 ) / prec;
493 double MaxElementLength2D::GetValue( long theElementId )
496 return GetPoints( theElementId, P ) ? GetValue(P) : 0.0;
499 double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const
504 SMDSAbs_ElementType MaxElementLength2D::GetType() const
509 //=======================================================================
511 Class : MaxElementLength3D
512 Description : Functor calculating maximum length of 3D element
514 //================================================================================
516 double MaxElementLength3D::GetValue( long theElementId )
519 if( GetPoints( theElementId, P ) ) {
521 const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
522 SMDSAbs_ElementType aType = aElem->GetType();
526 if( len == 4 ) { // tetras
527 double L1 = getDistance(P( 1 ),P( 2 ));
528 double L2 = getDistance(P( 2 ),P( 3 ));
529 double L3 = getDistance(P( 3 ),P( 1 ));
530 double L4 = getDistance(P( 1 ),P( 4 ));
531 double L5 = getDistance(P( 2 ),P( 4 ));
532 double L6 = getDistance(P( 3 ),P( 4 ));
533 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
536 else if( len == 5 ) { // pyramids
537 double L1 = getDistance(P( 1 ),P( 2 ));
538 double L2 = getDistance(P( 2 ),P( 3 ));
539 double L3 = getDistance(P( 3 ),P( 4 ));
540 double L4 = getDistance(P( 4 ),P( 1 ));
541 double L5 = getDistance(P( 1 ),P( 5 ));
542 double L6 = getDistance(P( 2 ),P( 5 ));
543 double L7 = getDistance(P( 3 ),P( 5 ));
544 double L8 = getDistance(P( 4 ),P( 5 ));
545 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
546 aVal = Max(aVal,Max(L7,L8));
549 else if( len == 6 ) { // pentas
550 double L1 = getDistance(P( 1 ),P( 2 ));
551 double L2 = getDistance(P( 2 ),P( 3 ));
552 double L3 = getDistance(P( 3 ),P( 1 ));
553 double L4 = getDistance(P( 4 ),P( 5 ));
554 double L5 = getDistance(P( 5 ),P( 6 ));
555 double L6 = getDistance(P( 6 ),P( 4 ));
556 double L7 = getDistance(P( 1 ),P( 4 ));
557 double L8 = getDistance(P( 2 ),P( 5 ));
558 double L9 = getDistance(P( 3 ),P( 6 ));
559 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
560 aVal = Max(aVal,Max(Max(L7,L8),L9));
563 else if( len == 8 ) { // hexas
564 double L1 = getDistance(P( 1 ),P( 2 ));
565 double L2 = getDistance(P( 2 ),P( 3 ));
566 double L3 = getDistance(P( 3 ),P( 4 ));
567 double L4 = getDistance(P( 4 ),P( 1 ));
568 double L5 = getDistance(P( 5 ),P( 6 ));
569 double L6 = getDistance(P( 6 ),P( 7 ));
570 double L7 = getDistance(P( 7 ),P( 8 ));
571 double L8 = getDistance(P( 8 ),P( 5 ));
572 double L9 = getDistance(P( 1 ),P( 5 ));
573 double L10= getDistance(P( 2 ),P( 6 ));
574 double L11= getDistance(P( 3 ),P( 7 ));
575 double L12= getDistance(P( 4 ),P( 8 ));
576 double D1 = getDistance(P( 1 ),P( 7 ));
577 double D2 = getDistance(P( 2 ),P( 8 ));
578 double D3 = getDistance(P( 3 ),P( 5 ));
579 double D4 = getDistance(P( 4 ),P( 6 ));
580 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
581 aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
582 aVal = Max(aVal,Max(L11,L12));
583 aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
586 else if( len == 12 ) { // hexagonal prism
587 for ( int i1 = 1; i1 < 12; ++i1 )
588 for ( int i2 = i1+1; i1 <= 12; ++i1 )
589 aVal = Max( aVal, getDistance(P( i1 ),P( i2 )));
592 else if( len == 10 ) { // quadratic tetras
593 double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
594 double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
595 double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
596 double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
597 double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
598 double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
599 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
602 else if( len == 13 ) { // quadratic pyramids
603 double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
604 double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
605 double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
606 double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
607 double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
608 double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
609 double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
610 double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
611 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
612 aVal = Max(aVal,Max(L7,L8));
615 else if( len == 15 ) { // quadratic pentas
616 double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
617 double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
618 double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
619 double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
620 double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
621 double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
622 double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
623 double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
624 double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
625 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
626 aVal = Max(aVal,Max(Max(L7,L8),L9));
629 else if( len == 20 || len == 27 ) { // quadratic hexas
630 double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
631 double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
632 double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
633 double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
634 double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
635 double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
636 double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
637 double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
638 double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
639 double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
640 double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
641 double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
642 double D1 = getDistance(P( 1 ),P( 7 ));
643 double D2 = getDistance(P( 2 ),P( 8 ));
644 double D3 = getDistance(P( 3 ),P( 5 ));
645 double D4 = getDistance(P( 4 ),P( 6 ));
646 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
647 aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
648 aVal = Max(aVal,Max(L11,L12));
649 aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
652 else if( len > 1 && aElem->IsPoly() ) { // polys
653 // get the maximum distance between all pairs of nodes
654 for( int i = 1; i <= len; i++ ) {
655 for( int j = 1; j <= len; j++ ) {
656 if( j > i ) { // optimization of the loop
657 double D = getDistance( P(i), P(j) );
658 aVal = Max( aVal, D );
665 if( myPrecision >= 0 )
667 double prec = pow( 10., (double)myPrecision );
668 aVal = floor( aVal * prec + 0.5 ) / prec;
675 double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const
680 SMDSAbs_ElementType MaxElementLength3D::GetType() const
682 return SMDSAbs_Volume;
685 //=======================================================================
688 Description : Functor for calculation of minimum angle
690 //================================================================================
692 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
699 aMin = getAngle(P( P.size() ), P( 1 ), P( 2 ));
700 aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 )));
702 for (int i=2; i<P.size();i++){
703 double A0 = getAngle( P( i-1 ), P( i ), P( i+1 ) );
707 return aMin * 180.0 / M_PI;
710 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
712 //const double aBestAngle = PI / nbNodes;
713 const double aBestAngle = 180.0 - ( 360.0 / double(nbNodes) );
714 return ( fabs( aBestAngle - Value ));
717 SMDSAbs_ElementType MinimumAngle::GetType() const
723 //================================================================================
726 Description : Functor for calculating aspect ratio
728 //================================================================================
730 double AspectRatio::GetValue( long theId )
733 myCurrElement = myMesh->FindElement( theId );
734 if ( myCurrElement && myCurrElement->GetVtkType() == VTK_QUAD )
737 vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid();
738 if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() ))
739 aVal = Round( vtkMeshQuality::QuadAspectRatio( avtkCell ));
744 if ( GetPoints( myCurrElement, P ))
745 aVal = Round( GetValue( P ));
750 double AspectRatio::GetValue( const TSequenceOfXYZ& P )
752 // According to "Mesh quality control" by Nadir Bouhamau referring to
753 // Pascal Jean Frey and Paul-Louis George. Maillages, applications aux elements finis.
754 // Hermes Science publications, Paris 1999 ISBN 2-7462-0024-4
757 int nbNodes = P.size();
762 // Compute aspect ratio
764 if ( nbNodes == 3 ) {
765 // Compute lengths of the sides
766 std::vector< double > aLen (nbNodes);
767 for ( int i = 0; i < nbNodes - 1; i++ )
768 aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
769 aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
770 // Q = alfa * h * p / S, where
772 // alfa = sqrt( 3 ) / 6
773 // h - length of the longest edge
774 // p - half perimeter
775 // S - triangle surface
776 const double alfa = sqrt( 3. ) / 6.;
777 double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
778 double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
779 double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
780 if ( anArea <= theEps )
782 return alfa * maxLen * half_perimeter / anArea;
784 else if ( nbNodes == 6 ) { // quadratic triangles
785 // Compute lengths of the sides
786 std::vector< double > aLen (3);
787 aLen[0] = getDistance( P(1), P(3) );
788 aLen[1] = getDistance( P(3), P(5) );
789 aLen[2] = getDistance( P(5), P(1) );
790 // Q = alfa * h * p / S, where
792 // alfa = sqrt( 3 ) / 6
793 // h - length of the longest edge
794 // p - half perimeter
795 // S - triangle surface
796 const double alfa = sqrt( 3. ) / 6.;
797 double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
798 double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
799 double anArea = getArea( P(1), P(3), P(5) );
800 if ( anArea <= theEps )
802 return alfa * maxLen * half_perimeter / anArea;
804 else if( nbNodes == 4 ) { // quadrangle
805 // Compute lengths of the sides
806 std::vector< double > aLen (4);
807 aLen[0] = getDistance( P(1), P(2) );
808 aLen[1] = getDistance( P(2), P(3) );
809 aLen[2] = getDistance( P(3), P(4) );
810 aLen[3] = getDistance( P(4), P(1) );
811 // Compute lengths of the diagonals
812 std::vector< double > aDia (2);
813 aDia[0] = getDistance( P(1), P(3) );
814 aDia[1] = getDistance( P(2), P(4) );
815 // Compute areas of all triangles which can be built
816 // taking three nodes of the quadrangle
817 std::vector< double > anArea (4);
818 anArea[0] = getArea( P(1), P(2), P(3) );
819 anArea[1] = getArea( P(1), P(2), P(4) );
820 anArea[2] = getArea( P(1), P(3), P(4) );
821 anArea[3] = getArea( P(2), P(3), P(4) );
822 // Q = alpha * L * C1 / C2, where
824 // alpha = sqrt( 1/32 )
825 // L = max( L1, L2, L3, L4, D1, D2 )
826 // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
827 // C2 = min( S1, S2, S3, S4 )
828 // Li - lengths of the edges
829 // Di - lengths of the diagonals
830 // Si - areas of the triangles
831 const double alpha = sqrt( 1 / 32. );
832 double L = Max( aLen[ 0 ],
836 Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
837 double C1 = sqrt( ( aLen[0] * aLen[0] +
840 aLen[3] * aLen[3] ) / 4. );
841 double C2 = Min( anArea[ 0 ],
843 Min( anArea[ 2 ], anArea[ 3 ] ) ) );
846 return alpha * L * C1 / C2;
848 else if( nbNodes == 8 || nbNodes == 9 ) { // nbNodes==8 - quadratic quadrangle
849 // Compute lengths of the sides
850 std::vector< double > aLen (4);
851 aLen[0] = getDistance( P(1), P(3) );
852 aLen[1] = getDistance( P(3), P(5) );
853 aLen[2] = getDistance( P(5), P(7) );
854 aLen[3] = getDistance( P(7), P(1) );
855 // Compute lengths of the diagonals
856 std::vector< double > aDia (2);
857 aDia[0] = getDistance( P(1), P(5) );
858 aDia[1] = getDistance( P(3), P(7) );
859 // Compute areas of all triangles which can be built
860 // taking three nodes of the quadrangle
861 std::vector< double > anArea (4);
862 anArea[0] = getArea( P(1), P(3), P(5) );
863 anArea[1] = getArea( P(1), P(3), P(7) );
864 anArea[2] = getArea( P(1), P(5), P(7) );
865 anArea[3] = getArea( P(3), P(5), P(7) );
866 // Q = alpha * L * C1 / C2, where
868 // alpha = sqrt( 1/32 )
869 // L = max( L1, L2, L3, L4, D1, D2 )
870 // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
871 // C2 = min( S1, S2, S3, S4 )
872 // Li - lengths of the edges
873 // Di - lengths of the diagonals
874 // Si - areas of the triangles
875 const double alpha = sqrt( 1 / 32. );
876 double L = Max( aLen[ 0 ],
880 Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
881 double C1 = sqrt( ( aLen[0] * aLen[0] +
884 aLen[3] * aLen[3] ) / 4. );
885 double C2 = Min( anArea[ 0 ],
887 Min( anArea[ 2 ], anArea[ 3 ] ) ) );
890 return alpha * L * C1 / C2;
895 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
897 // the aspect ratio is in the range [1.0,infinity]
898 // < 1.0 = very bad, zero area
901 return ( Value < 0.9 ) ? 1000 : Value / 1000.;
904 SMDSAbs_ElementType AspectRatio::GetType() const
910 //================================================================================
912 Class : AspectRatio3D
913 Description : Functor for calculating aspect ratio
915 //================================================================================
919 inline double getHalfPerimeter(double theTria[3]){
920 return (theTria[0] + theTria[1] + theTria[2])/2.0;
923 inline double getArea(double theHalfPerim, double theTria[3]){
924 return sqrt(theHalfPerim*
925 (theHalfPerim-theTria[0])*
926 (theHalfPerim-theTria[1])*
927 (theHalfPerim-theTria[2]));
930 inline double getVolume(double theLen[6]){
931 double a2 = theLen[0]*theLen[0];
932 double b2 = theLen[1]*theLen[1];
933 double c2 = theLen[2]*theLen[2];
934 double d2 = theLen[3]*theLen[3];
935 double e2 = theLen[4]*theLen[4];
936 double f2 = theLen[5]*theLen[5];
937 double P = 4.0*a2*b2*d2;
938 double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
939 double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
940 return sqrt(P-Q+R)/12.0;
943 inline double getVolume2(double theLen[6]){
944 double a2 = theLen[0]*theLen[0];
945 double b2 = theLen[1]*theLen[1];
946 double c2 = theLen[2]*theLen[2];
947 double d2 = theLen[3]*theLen[3];
948 double e2 = theLen[4]*theLen[4];
949 double f2 = theLen[5]*theLen[5];
951 double P = a2*e2*(b2+c2+d2+f2-a2-e2);
952 double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
953 double R = c2*d2*(a2+b2+e2+f2-c2-d2);
954 double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
956 return sqrt(P+Q+R-S)/12.0;
959 inline double getVolume(const TSequenceOfXYZ& P){
960 gp_Vec aVec1( P( 2 ) - P( 1 ) );
961 gp_Vec aVec2( P( 3 ) - P( 1 ) );
962 gp_Vec aVec3( P( 4 ) - P( 1 ) );
963 gp_Vec anAreaVec( aVec1 ^ aVec2 );
964 return fabs(aVec3 * anAreaVec) / 6.0;
967 inline double getMaxHeight(double theLen[6])
969 double aHeight = std::max(theLen[0],theLen[1]);
970 aHeight = std::max(aHeight,theLen[2]);
971 aHeight = std::max(aHeight,theLen[3]);
972 aHeight = std::max(aHeight,theLen[4]);
973 aHeight = std::max(aHeight,theLen[5]);
979 double AspectRatio3D::GetValue( long theId )
982 myCurrElement = myMesh->FindElement( theId );
983 if ( myCurrElement && myCurrElement->GetVtkType() == VTK_TETRA )
985 // Action from CoTech | ACTION 31.3:
986 // EURIWARE BO: Homogenize the formulas used to calculate the Controls in SMESH to fit with
987 // those of ParaView. The library used by ParaView for those calculations can be reused in SMESH.
988 vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid();
989 if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() ))
990 aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell ));
995 if ( GetPoints( myCurrElement, P ))
996 aVal = Round( GetValue( P ));
1001 double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
1003 double aQuality = 0.0;
1004 if(myCurrElement->IsPoly()) return aQuality;
1006 int nbNodes = P.size();
1008 if(myCurrElement->IsQuadratic()) {
1009 if(nbNodes==10) nbNodes=4; // quadratic tetrahedron
1010 else if(nbNodes==13) nbNodes=5; // quadratic pyramid
1011 else if(nbNodes==15) nbNodes=6; // quadratic pentahedron
1012 else if(nbNodes==20) nbNodes=8; // quadratic hexahedron
1013 else if(nbNodes==27) nbNodes=8; // quadratic hexahedron
1014 else return aQuality;
1020 getDistance(P( 1 ),P( 2 )), // a
1021 getDistance(P( 2 ),P( 3 )), // b
1022 getDistance(P( 3 ),P( 1 )), // c
1023 getDistance(P( 2 ),P( 4 )), // d
1024 getDistance(P( 3 ),P( 4 )), // e
1025 getDistance(P( 1 ),P( 4 )) // f
1027 double aTria[4][3] = {
1028 {aLen[0],aLen[1],aLen[2]}, // abc
1029 {aLen[0],aLen[3],aLen[5]}, // adf
1030 {aLen[1],aLen[3],aLen[4]}, // bde
1031 {aLen[2],aLen[4],aLen[5]} // cef
1033 double aSumArea = 0.0;
1034 double aHalfPerimeter = getHalfPerimeter(aTria[0]);
1035 double anArea = getArea(aHalfPerimeter,aTria[0]);
1037 aHalfPerimeter = getHalfPerimeter(aTria[1]);
1038 anArea = getArea(aHalfPerimeter,aTria[1]);
1040 aHalfPerimeter = getHalfPerimeter(aTria[2]);
1041 anArea = getArea(aHalfPerimeter,aTria[2]);
1043 aHalfPerimeter = getHalfPerimeter(aTria[3]);
1044 anArea = getArea(aHalfPerimeter,aTria[3]);
1046 double aVolume = getVolume(P);
1047 //double aVolume = getVolume(aLen);
1048 double aHeight = getMaxHeight(aLen);
1049 static double aCoeff = sqrt(2.0)/12.0;
1050 if ( aVolume > DBL_MIN )
1051 aQuality = aCoeff*aHeight*aSumArea/aVolume;
1056 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
1057 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1060 gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
1061 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1064 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
1065 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1068 gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
1069 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1075 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
1076 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1079 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
1080 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1083 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
1084 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1087 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
1088 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1091 gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
1092 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1095 gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
1096 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1102 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
1103 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1106 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
1107 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1110 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
1111 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1114 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
1115 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1118 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
1119 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1122 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
1123 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1126 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
1127 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1130 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
1131 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1134 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
1135 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1138 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
1139 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1142 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
1143 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1146 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
1147 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1150 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
1151 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1154 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
1155 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1158 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
1159 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1162 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
1163 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1166 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
1167 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1170 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
1171 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1174 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
1175 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1178 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
1179 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1182 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
1183 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1186 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1187 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1190 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
1191 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1194 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
1195 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1198 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1199 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1202 gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
1203 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1206 gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
1207 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1210 gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
1211 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1214 gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
1215 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1218 gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
1219 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1222 gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
1223 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1226 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
1227 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1230 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
1231 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1237 gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )};
1238 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1241 gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )};
1242 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1245 gp_XYZ aXYZ[8] = {P( 3 ),P( 4 ),P( 6 ),P( 1 ),P( 9 ),P( 10 ),P( 12 ),P( 7 )};
1246 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1249 } // switch(nbNodes)
1251 if ( nbNodes > 4 ) {
1252 // avaluate aspect ratio of quadranle faces
1253 AspectRatio aspect2D;
1254 SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes );
1255 int nbFaces = SMDS_VolumeTool::NbFaces( type );
1256 TSequenceOfXYZ points(4);
1257 for ( int i = 0; i < nbFaces; ++i ) { // loop on faces of a volume
1258 if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 )
1260 const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true );
1261 for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face
1262 points( p + 1 ) = P( pInd[ p ] + 1 );
1263 aQuality = std::max( aQuality, aspect2D.GetValue( points ));
1269 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
1271 // the aspect ratio is in the range [1.0,infinity]
1274 return Value / 1000.;
1277 SMDSAbs_ElementType AspectRatio3D::GetType() const
1279 return SMDSAbs_Volume;
1283 //================================================================================
1286 Description : Functor for calculating warping
1288 //================================================================================
1290 double Warping::GetValue( const TSequenceOfXYZ& P )
1292 if ( P.size() != 4 )
1295 gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4.;
1297 double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
1298 double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
1299 double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
1300 double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
1302 double val = Max( Max( A1, A2 ), Max( A3, A4 ) );
1304 const double eps = 0.1; // val is in degrees
1306 return val < eps ? 0. : val;
1309 double Warping::ComputeA( const gp_XYZ& thePnt1,
1310 const gp_XYZ& thePnt2,
1311 const gp_XYZ& thePnt3,
1312 const gp_XYZ& theG ) const
1314 double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
1315 double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
1316 double L = Min( aLen1, aLen2 ) * 0.5;
1320 gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG;
1321 gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG;
1322 gp_XYZ N = GI.Crossed( GJ );
1324 if ( N.Modulus() < gp::Resolution() )
1329 double H = ( thePnt2 - theG ).Dot( N );
1330 return asin( fabs( H / L ) ) * 180. / M_PI;
1333 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
1335 // the warp is in the range [0.0,PI/2]
1336 // 0.0 = good (no warp)
1337 // PI/2 = bad (face pliee)
1341 SMDSAbs_ElementType Warping::GetType() const
1343 return SMDSAbs_Face;
1347 //================================================================================
1350 Description : Functor for calculating taper
1352 //================================================================================
1354 double Taper::GetValue( const TSequenceOfXYZ& P )
1356 if ( P.size() != 4 )
1360 double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) );
1361 double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) );
1362 double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) );
1363 double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) );
1365 double JA = 0.25 * ( J1 + J2 + J3 + J4 );
1369 double T1 = fabs( ( J1 - JA ) / JA );
1370 double T2 = fabs( ( J2 - JA ) / JA );
1371 double T3 = fabs( ( J3 - JA ) / JA );
1372 double T4 = fabs( ( J4 - JA ) / JA );
1374 double val = Max( Max( T1, T2 ), Max( T3, T4 ) );
1376 const double eps = 0.01;
1378 return val < eps ? 0. : val;
1381 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
1383 // the taper is in the range [0.0,1.0]
1384 // 0.0 = good (no taper)
1385 // 1.0 = bad (les cotes opposes sont allignes)
1389 SMDSAbs_ElementType Taper::GetType() const
1391 return SMDSAbs_Face;
1394 //================================================================================
1397 Description : Functor for calculating skew in degrees
1399 //================================================================================
1401 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
1403 gp_XYZ p12 = ( p2 + p1 ) / 2.;
1404 gp_XYZ p23 = ( p3 + p2 ) / 2.;
1405 gp_XYZ p31 = ( p3 + p1 ) / 2.;
1407 gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
1409 return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 );
1412 double Skew::GetValue( const TSequenceOfXYZ& P )
1414 if ( P.size() != 3 && P.size() != 4 )
1418 const double PI2 = M_PI / 2.;
1419 if ( P.size() == 3 )
1421 double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
1422 double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
1423 double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
1425 return Max( A0, Max( A1, A2 ) ) * 180. / M_PI;
1429 gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2.;
1430 gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2.;
1431 gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2.;
1432 gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2.;
1434 gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
1435 double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
1436 ? 0. : fabs( PI2 - v1.Angle( v2 ) );
1438 double val = A * 180. / M_PI;
1440 const double eps = 0.1; // val is in degrees
1442 return val < eps ? 0. : val;
1446 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
1448 // the skew is in the range [0.0,PI/2].
1454 SMDSAbs_ElementType Skew::GetType() const
1456 return SMDSAbs_Face;
1460 //================================================================================
1463 Description : Functor for calculating area
1465 //================================================================================
1467 double Area::GetValue( const TSequenceOfXYZ& P )
1470 if ( P.size() > 2 ) {
1471 gp_Vec aVec1( P(2) - P(1) );
1472 gp_Vec aVec2( P(3) - P(1) );
1473 gp_Vec SumVec = aVec1 ^ aVec2;
1474 for (int i=4; i<=P.size(); i++) {
1475 gp_Vec aVec1( P(i-1) - P(1) );
1476 gp_Vec aVec2( P(i) - P(1) );
1477 gp_Vec tmp = aVec1 ^ aVec2;
1480 val = SumVec.Magnitude() * 0.5;
1485 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
1487 // meaningless as it is not a quality control functor
1491 SMDSAbs_ElementType Area::GetType() const
1493 return SMDSAbs_Face;
1496 //================================================================================
1499 Description : Functor for calculating length of edge
1501 //================================================================================
1503 double Length::GetValue( const TSequenceOfXYZ& P )
1505 switch ( P.size() ) {
1506 case 2: return getDistance( P( 1 ), P( 2 ) );
1507 case 3: return getDistance( P( 1 ), P( 2 ) ) + getDistance( P( 2 ), P( 3 ) );
1512 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
1514 // meaningless as it is not quality control functor
1518 SMDSAbs_ElementType Length::GetType() const
1520 return SMDSAbs_Edge;
1523 //================================================================================
1526 Description : Functor for calculating length of edge
1528 //================================================================================
1530 double Length2D::GetValue( long theElementId )
1534 //cout<<"Length2D::GetValue"<<endl;
1535 if (GetPoints(theElementId,P)){
1536 //for(int jj=1; jj<=P.size(); jj++)
1537 // cout<<"jj="<<jj<<" P("<<P(jj).X()<<","<<P(jj).Y()<<","<<P(jj).Z()<<")"<<endl;
1539 double aVal;// = GetValue( P );
1540 const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
1541 SMDSAbs_ElementType aType = aElem->GetType();
1550 aVal = getDistance( P( 1 ), P( 2 ) );
1553 else if (len == 3){ // quadratic edge
1554 aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
1558 if (len == 3){ // triangles
1559 double L1 = getDistance(P( 1 ),P( 2 ));
1560 double L2 = getDistance(P( 2 ),P( 3 ));
1561 double L3 = getDistance(P( 3 ),P( 1 ));
1562 aVal = Min(L1,Min(L2,L3));
1565 else if (len == 4){ // quadrangles
1566 double L1 = getDistance(P( 1 ),P( 2 ));
1567 double L2 = getDistance(P( 2 ),P( 3 ));
1568 double L3 = getDistance(P( 3 ),P( 4 ));
1569 double L4 = getDistance(P( 4 ),P( 1 ));
1570 aVal = Min(Min(L1,L2),Min(L3,L4));
1573 if (len == 6){ // quadratic triangles
1574 double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1575 double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1576 double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
1577 aVal = Min(L1,Min(L2,L3));
1578 //cout<<"L1="<<L1<<" L2="<<L2<<"L3="<<L3<<" aVal="<<aVal<<endl;
1581 else if (len == 8){ // quadratic quadrangles
1582 double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1583 double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1584 double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
1585 double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
1586 aVal = Min(Min(L1,L2),Min(L3,L4));
1589 case SMDSAbs_Volume:
1590 if (len == 4){ // tetraidrs
1591 double L1 = getDistance(P( 1 ),P( 2 ));
1592 double L2 = getDistance(P( 2 ),P( 3 ));
1593 double L3 = getDistance(P( 3 ),P( 1 ));
1594 double L4 = getDistance(P( 1 ),P( 4 ));
1595 double L5 = getDistance(P( 2 ),P( 4 ));
1596 double L6 = getDistance(P( 3 ),P( 4 ));
1597 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1600 else if (len == 5){ // piramids
1601 double L1 = getDistance(P( 1 ),P( 2 ));
1602 double L2 = getDistance(P( 2 ),P( 3 ));
1603 double L3 = getDistance(P( 3 ),P( 4 ));
1604 double L4 = getDistance(P( 4 ),P( 1 ));
1605 double L5 = getDistance(P( 1 ),P( 5 ));
1606 double L6 = getDistance(P( 2 ),P( 5 ));
1607 double L7 = getDistance(P( 3 ),P( 5 ));
1608 double L8 = getDistance(P( 4 ),P( 5 ));
1610 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1611 aVal = Min(aVal,Min(L7,L8));
1614 else if (len == 6){ // pentaidres
1615 double L1 = getDistance(P( 1 ),P( 2 ));
1616 double L2 = getDistance(P( 2 ),P( 3 ));
1617 double L3 = getDistance(P( 3 ),P( 1 ));
1618 double L4 = getDistance(P( 4 ),P( 5 ));
1619 double L5 = getDistance(P( 5 ),P( 6 ));
1620 double L6 = getDistance(P( 6 ),P( 4 ));
1621 double L7 = getDistance(P( 1 ),P( 4 ));
1622 double L8 = getDistance(P( 2 ),P( 5 ));
1623 double L9 = getDistance(P( 3 ),P( 6 ));
1625 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1626 aVal = Min(aVal,Min(Min(L7,L8),L9));
1629 else if (len == 8){ // hexaider
1630 double L1 = getDistance(P( 1 ),P( 2 ));
1631 double L2 = getDistance(P( 2 ),P( 3 ));
1632 double L3 = getDistance(P( 3 ),P( 4 ));
1633 double L4 = getDistance(P( 4 ),P( 1 ));
1634 double L5 = getDistance(P( 5 ),P( 6 ));
1635 double L6 = getDistance(P( 6 ),P( 7 ));
1636 double L7 = getDistance(P( 7 ),P( 8 ));
1637 double L8 = getDistance(P( 8 ),P( 5 ));
1638 double L9 = getDistance(P( 1 ),P( 5 ));
1639 double L10= getDistance(P( 2 ),P( 6 ));
1640 double L11= getDistance(P( 3 ),P( 7 ));
1641 double L12= getDistance(P( 4 ),P( 8 ));
1643 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1644 aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
1645 aVal = Min(aVal,Min(L11,L12));
1650 if (len == 10){ // quadratic tetraidrs
1651 double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
1652 double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
1653 double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
1654 double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1655 double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
1656 double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
1657 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1660 else if (len == 13){ // quadratic piramids
1661 double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
1662 double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
1663 double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1664 double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1665 double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1666 double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
1667 double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
1668 double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
1669 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1670 aVal = Min(aVal,Min(L7,L8));
1673 else if (len == 15){ // quadratic pentaidres
1674 double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
1675 double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
1676 double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1677 double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1678 double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
1679 double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
1680 double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
1681 double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
1682 double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
1683 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1684 aVal = Min(aVal,Min(Min(L7,L8),L9));
1687 else if (len == 20){ // quadratic hexaider
1688 double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
1689 double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
1690 double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
1691 double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
1692 double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
1693 double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
1694 double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
1695 double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
1696 double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
1697 double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
1698 double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
1699 double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
1700 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1701 aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
1702 aVal = Min(aVal,Min(L11,L12));
1714 if ( myPrecision >= 0 )
1716 double prec = pow( 10., (double)( myPrecision ) );
1717 aVal = floor( aVal * prec + 0.5 ) / prec;
1726 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1728 // meaningless as it is not a quality control functor
1732 SMDSAbs_ElementType Length2D::GetType() const
1734 return SMDSAbs_Face;
1737 Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
1740 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
1741 if(thePntId1 > thePntId2){
1742 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
1746 bool Length2D::Value::operator<(const Length2D::Value& x) const{
1747 if(myPntId[0] < x.myPntId[0]) return true;
1748 if(myPntId[0] == x.myPntId[0])
1749 if(myPntId[1] < x.myPntId[1]) return true;
1753 void Length2D::GetValues(TValues& theValues){
1755 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1756 for(; anIter->more(); ){
1757 const SMDS_MeshFace* anElem = anIter->next();
1759 if(anElem->IsQuadratic()) {
1760 const SMDS_VtkFace* F =
1761 dynamic_cast<const SMDS_VtkFace*>(anElem);
1762 // use special nodes iterator
1763 SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
1768 const SMDS_MeshElement* aNode;
1770 aNode = anIter->next();
1771 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1772 P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1773 aNodeId[0] = aNodeId[1] = aNode->GetID();
1776 for(; anIter->more(); ){
1777 const SMDS_MeshNode* N1 = static_cast<const SMDS_MeshNode*> (anIter->next());
1778 P[2] = gp_Pnt(N1->X(),N1->Y(),N1->Z());
1779 aNodeId[2] = N1->GetID();
1780 aLength = P[1].Distance(P[2]);
1781 if(!anIter->more()) break;
1782 const SMDS_MeshNode* N2 = static_cast<const SMDS_MeshNode*> (anIter->next());
1783 P[3] = gp_Pnt(N2->X(),N2->Y(),N2->Z());
1784 aNodeId[3] = N2->GetID();
1785 aLength += P[2].Distance(P[3]);
1786 Value aValue1(aLength,aNodeId[1],aNodeId[2]);
1787 Value aValue2(aLength,aNodeId[2],aNodeId[3]);
1789 aNodeId[1] = aNodeId[3];
1790 theValues.insert(aValue1);
1791 theValues.insert(aValue2);
1793 aLength += P[2].Distance(P[0]);
1794 Value aValue1(aLength,aNodeId[1],aNodeId[2]);
1795 Value aValue2(aLength,aNodeId[2],aNodeId[0]);
1796 theValues.insert(aValue1);
1797 theValues.insert(aValue2);
1800 SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
1805 const SMDS_MeshElement* aNode;
1806 if(aNodesIter->more()){
1807 aNode = aNodesIter->next();
1808 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1809 P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1810 aNodeId[0] = aNodeId[1] = aNode->GetID();
1813 for(; aNodesIter->more(); ){
1814 aNode = aNodesIter->next();
1815 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1816 long anId = aNode->GetID();
1818 P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1820 aLength = P[1].Distance(P[2]);
1822 Value aValue(aLength,aNodeId[1],anId);
1825 theValues.insert(aValue);
1828 aLength = P[0].Distance(P[1]);
1830 Value aValue(aLength,aNodeId[0],aNodeId[1]);
1831 theValues.insert(aValue);
1836 //================================================================================
1838 Class : MultiConnection
1839 Description : Functor for calculating number of faces conneted to the edge
1841 //================================================================================
1843 double MultiConnection::GetValue( const TSequenceOfXYZ& P )
1847 double MultiConnection::GetValue( long theId )
1849 return getNbMultiConnection( myMesh, theId );
1852 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
1854 // meaningless as it is not quality control functor
1858 SMDSAbs_ElementType MultiConnection::GetType() const
1860 return SMDSAbs_Edge;
1863 //================================================================================
1865 Class : MultiConnection2D
1866 Description : Functor for calculating number of faces conneted to the edge
1868 //================================================================================
1870 double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
1875 double MultiConnection2D::GetValue( long theElementId )
1879 const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId);
1880 SMDSAbs_ElementType aType = aFaceElem->GetType();
1885 int i = 0, len = aFaceElem->NbNodes();
1886 SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
1889 const SMDS_MeshNode *aNode, *aNode0;
1890 TColStd_MapOfInteger aMap, aMapPrev;
1892 for (i = 0; i <= len; i++) {
1897 if (anIter->more()) {
1898 aNode = (SMDS_MeshNode*)anIter->next();
1906 if (i == 0) aNode0 = aNode;
1908 SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
1909 while (anElemIter->more()) {
1910 const SMDS_MeshElement* anElem = anElemIter->next();
1911 if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
1912 int anId = anElem->GetID();
1915 if (aMapPrev.Contains(anId)) {
1920 aResult = Max(aResult, aNb);
1931 double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1933 // meaningless as it is not quality control functor
1937 SMDSAbs_ElementType MultiConnection2D::GetType() const
1939 return SMDSAbs_Face;
1942 MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
1944 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
1945 if(thePntId1 > thePntId2){
1946 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
1950 bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const{
1951 if(myPntId[0] < x.myPntId[0]) return true;
1952 if(myPntId[0] == x.myPntId[0])
1953 if(myPntId[1] < x.myPntId[1]) return true;
1957 void MultiConnection2D::GetValues(MValues& theValues){
1958 if ( !myMesh ) return;
1959 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1960 for(; anIter->more(); ){
1961 const SMDS_MeshFace* anElem = anIter->next();
1962 SMDS_ElemIteratorPtr aNodesIter;
1963 if ( anElem->IsQuadratic() )
1964 aNodesIter = dynamic_cast<const SMDS_VtkFace*>
1965 (anElem)->interlacedNodesElemIterator();
1967 aNodesIter = anElem->nodesIterator();
1970 //int aNbConnects=0;
1971 const SMDS_MeshNode* aNode0;
1972 const SMDS_MeshNode* aNode1;
1973 const SMDS_MeshNode* aNode2;
1974 if(aNodesIter->more()){
1975 aNode0 = (SMDS_MeshNode*) aNodesIter->next();
1977 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
1978 aNodeId[0] = aNodeId[1] = aNodes->GetID();
1980 for(; aNodesIter->more(); ) {
1981 aNode2 = (SMDS_MeshNode*) aNodesIter->next();
1982 long anId = aNode2->GetID();
1985 Value aValue(aNodeId[1],aNodeId[2]);
1986 MValues::iterator aItr = theValues.find(aValue);
1987 if (aItr != theValues.end()){
1992 theValues[aValue] = 1;
1995 //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
1996 aNodeId[1] = aNodeId[2];
1999 Value aValue(aNodeId[0],aNodeId[2]);
2000 MValues::iterator aItr = theValues.find(aValue);
2001 if (aItr != theValues.end()) {
2006 theValues[aValue] = 1;
2009 //cout << "NodeIds: "<<aNodeId[0]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
2014 //================================================================================
2016 Class : BallDiameter
2017 Description : Functor returning diameter of a ball element
2019 //================================================================================
2021 double BallDiameter::GetValue( long theId )
2023 double diameter = 0;
2025 if ( const SMDS_BallElement* ball =
2026 dynamic_cast<const SMDS_BallElement*>( myMesh->FindElement( theId )))
2028 diameter = ball->GetDiameter();
2033 double BallDiameter::GetBadRate( double Value, int /*nbNodes*/ ) const
2035 // meaningless as it is not a quality control functor
2039 SMDSAbs_ElementType BallDiameter::GetType() const
2041 return SMDSAbs_Ball;
2049 //================================================================================
2051 Class : BadOrientedVolume
2052 Description : Predicate bad oriented volumes
2054 //================================================================================
2056 BadOrientedVolume::BadOrientedVolume()
2061 void BadOrientedVolume::SetMesh( const SMDS_Mesh* theMesh )
2066 bool BadOrientedVolume::IsSatisfy( long theId )
2071 SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
2072 return !vTool.IsForward();
2075 SMDSAbs_ElementType BadOrientedVolume::GetType() const
2077 return SMDSAbs_Volume;
2081 Class : BareBorderVolume
2084 bool BareBorderVolume::IsSatisfy(long theElementId )
2086 SMDS_VolumeTool myTool;
2087 if ( myTool.Set( myMesh->FindElement(theElementId)))
2089 for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
2090 if ( myTool.IsFreeFace( iF ))
2092 const SMDS_MeshNode** n = myTool.GetFaceNodes(iF);
2093 vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF));
2094 if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*Nomedium=*/false))
2101 //================================================================================
2103 Class : BareBorderFace
2105 //================================================================================
2107 bool BareBorderFace::IsSatisfy(long theElementId )
2110 if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
2112 if ( face->GetType() == SMDSAbs_Face )
2114 int nbN = face->NbCornerNodes();
2115 for ( int i = 0; i < nbN && !ok; ++i )
2117 // check if a link is shared by another face
2118 const SMDS_MeshNode* n1 = face->GetNode( i );
2119 const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
2120 SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
2121 bool isShared = false;
2122 while ( !isShared && fIt->more() )
2124 const SMDS_MeshElement* f = fIt->next();
2125 isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
2129 const int iQuad = face->IsQuadratic();
2130 myLinkNodes.resize( 2 + iQuad);
2131 myLinkNodes[0] = n1;
2132 myLinkNodes[1] = n2;
2134 myLinkNodes[2] = face->GetNode( i+nbN );
2135 ok = !myMesh->FindElement( myLinkNodes, SMDSAbs_Edge, /*noMedium=*/false);
2143 //================================================================================
2145 Class : OverConstrainedVolume
2147 //================================================================================
2149 bool OverConstrainedVolume::IsSatisfy(long theElementId )
2151 // An element is over-constrained if it has N-1 free borders where
2152 // N is the number of edges/faces for a 2D/3D element.
2153 SMDS_VolumeTool myTool;
2154 if ( myTool.Set( myMesh->FindElement(theElementId)))
2156 int nbSharedFaces = 0;
2157 for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
2158 if ( !myTool.IsFreeFace( iF ) && ++nbSharedFaces > 1 )
2160 return ( nbSharedFaces == 1 );
2165 //================================================================================
2167 Class : OverConstrainedFace
2169 //================================================================================
2171 bool OverConstrainedFace::IsSatisfy(long theElementId )
2173 // An element is over-constrained if it has N-1 free borders where
2174 // N is the number of edges/faces for a 2D/3D element.
2175 if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
2176 if ( face->GetType() == SMDSAbs_Face )
2178 int nbSharedBorders = 0;
2179 int nbN = face->NbCornerNodes();
2180 for ( int i = 0; i < nbN; ++i )
2182 // check if a link is shared by another face
2183 const SMDS_MeshNode* n1 = face->GetNode( i );
2184 const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
2185 SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
2186 bool isShared = false;
2187 while ( !isShared && fIt->more() )
2189 const SMDS_MeshElement* f = fIt->next();
2190 isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
2192 if ( isShared && ++nbSharedBorders > 1 )
2195 return ( nbSharedBorders == 1 );
2200 //================================================================================
2202 Class : CoincidentNodes
2203 Description : Predicate of Coincident nodes
2205 //================================================================================
2207 CoincidentNodes::CoincidentNodes()
2212 bool CoincidentNodes::IsSatisfy( long theElementId )
2214 return myCoincidentIDs.Contains( theElementId );
2217 SMDSAbs_ElementType CoincidentNodes::GetType() const
2219 return SMDSAbs_Node;
2222 void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
2224 myMeshModifTracer.SetMesh( theMesh );
2225 if ( myMeshModifTracer.IsMeshModified() )
2227 TIDSortedNodeSet nodesToCheck;
2228 SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true);
2229 while ( nIt->more() )
2230 nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
2232 list< list< const SMDS_MeshNode*> > nodeGroups;
2233 SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler );
2235 myCoincidentIDs.Clear();
2236 list< list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin();
2237 for ( ; groupIt != nodeGroups.end(); ++groupIt )
2239 list< const SMDS_MeshNode*>& coincNodes = *groupIt;
2240 list< const SMDS_MeshNode*>::iterator n = coincNodes.begin();
2241 for ( ; n != coincNodes.end(); ++n )
2242 myCoincidentIDs.Add( (*n)->GetID() );
2247 //================================================================================
2249 Class : CoincidentElements
2250 Description : Predicate of Coincident Elements
2251 Note : This class is suitable only for visualization of Coincident Elements
2253 //================================================================================
2255 CoincidentElements::CoincidentElements()
2260 void CoincidentElements::SetMesh( const SMDS_Mesh* theMesh )
2265 bool CoincidentElements::IsSatisfy( long theElementId )
2267 if ( !myMesh ) return false;
2269 if ( const SMDS_MeshElement* e = myMesh->FindElement( theElementId ))
2271 if ( e->GetType() != GetType() ) return false;
2272 set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() );
2273 const int nbNodes = e->NbNodes();
2274 SMDS_ElemIteratorPtr invIt = (*elemNodes.begin())->GetInverseElementIterator( GetType() );
2275 while ( invIt->more() )
2277 const SMDS_MeshElement* e2 = invIt->next();
2278 if ( e2 == e || e2->NbNodes() != nbNodes ) continue;
2280 bool sameNodes = true;
2281 for ( size_t i = 0; i < elemNodes.size() && sameNodes; ++i )
2282 sameNodes = ( elemNodes.count( e2->GetNode( i )));
2290 SMDSAbs_ElementType CoincidentElements1D::GetType() const
2292 return SMDSAbs_Edge;
2294 SMDSAbs_ElementType CoincidentElements2D::GetType() const
2296 return SMDSAbs_Face;
2298 SMDSAbs_ElementType CoincidentElements3D::GetType() const
2300 return SMDSAbs_Volume;
2304 //================================================================================
2307 Description : Predicate for free borders
2309 //================================================================================
2311 FreeBorders::FreeBorders()
2316 void FreeBorders::SetMesh( const SMDS_Mesh* theMesh )
2321 bool FreeBorders::IsSatisfy( long theId )
2323 return getNbMultiConnection( myMesh, theId ) == 1;
2326 SMDSAbs_ElementType FreeBorders::GetType() const
2328 return SMDSAbs_Edge;
2332 //================================================================================
2335 Description : Predicate for free Edges
2337 //================================================================================
2339 FreeEdges::FreeEdges()
2344 void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
2349 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId )
2351 TColStd_MapOfInteger aMap;
2352 for ( int i = 0; i < 2; i++ )
2354 SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator(SMDSAbs_Face);
2355 while( anElemIter->more() )
2357 if ( const SMDS_MeshElement* anElem = anElemIter->next())
2359 const int anId = anElem->GetID();
2360 if ( anId != theFaceId && !aMap.Add( anId ))
2368 bool FreeEdges::IsSatisfy( long theId )
2373 const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2374 if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
2377 SMDS_ElemIteratorPtr anIter;
2378 if ( aFace->IsQuadratic() ) {
2379 anIter = dynamic_cast<const SMDS_VtkFace*>
2380 (aFace)->interlacedNodesElemIterator();
2383 anIter = aFace->nodesIterator();
2388 int i = 0, nbNodes = aFace->NbNodes();
2389 std::vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
2390 while( anIter->more() )
2392 const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
2395 aNodes[ i++ ] = aNode;
2397 aNodes[ nbNodes ] = aNodes[ 0 ];
2399 for ( i = 0; i < nbNodes; i++ )
2400 if ( IsFreeEdge( &aNodes[ i ], theId ) )
2406 SMDSAbs_ElementType FreeEdges::GetType() const
2408 return SMDSAbs_Face;
2411 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
2414 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
2415 if(thePntId1 > thePntId2){
2416 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
2420 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
2421 if(myPntId[0] < x.myPntId[0]) return true;
2422 if(myPntId[0] == x.myPntId[0])
2423 if(myPntId[1] < x.myPntId[1]) return true;
2427 inline void UpdateBorders(const FreeEdges::Border& theBorder,
2428 FreeEdges::TBorders& theRegistry,
2429 FreeEdges::TBorders& theContainer)
2431 if(theRegistry.find(theBorder) == theRegistry.end()){
2432 theRegistry.insert(theBorder);
2433 theContainer.insert(theBorder);
2435 theContainer.erase(theBorder);
2439 void FreeEdges::GetBoreders(TBorders& theBorders)
2442 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
2443 for(; anIter->more(); ){
2444 const SMDS_MeshFace* anElem = anIter->next();
2445 long anElemId = anElem->GetID();
2446 SMDS_ElemIteratorPtr aNodesIter;
2447 if ( anElem->IsQuadratic() )
2448 aNodesIter = static_cast<const SMDS_VtkFace*>(anElem)->
2449 interlacedNodesElemIterator();
2451 aNodesIter = anElem->nodesIterator();
2453 const SMDS_MeshElement* aNode;
2454 if(aNodesIter->more()){
2455 aNode = aNodesIter->next();
2456 aNodeId[0] = aNodeId[1] = aNode->GetID();
2458 for(; aNodesIter->more(); ){
2459 aNode = aNodesIter->next();
2460 long anId = aNode->GetID();
2461 Border aBorder(anElemId,aNodeId[1],anId);
2463 UpdateBorders(aBorder,aRegistry,theBorders);
2465 Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
2466 UpdateBorders(aBorder,aRegistry,theBorders);
2470 //================================================================================
2473 Description : Predicate for free nodes
2475 //================================================================================
2477 FreeNodes::FreeNodes()
2482 void FreeNodes::SetMesh( const SMDS_Mesh* theMesh )
2487 bool FreeNodes::IsSatisfy( long theNodeId )
2489 const SMDS_MeshNode* aNode = myMesh->FindNode( theNodeId );
2493 return (aNode->NbInverseElements() < 1);
2496 SMDSAbs_ElementType FreeNodes::GetType() const
2498 return SMDSAbs_Node;
2502 //================================================================================
2505 Description : Predicate for free faces
2507 //================================================================================
2509 FreeFaces::FreeFaces()
2514 void FreeFaces::SetMesh( const SMDS_Mesh* theMesh )
2519 bool FreeFaces::IsSatisfy( long theId )
2521 if (!myMesh) return false;
2522 // check that faces nodes refers to less than two common volumes
2523 const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2524 if ( !aFace || aFace->GetType() != SMDSAbs_Face )
2527 int nbNode = aFace->NbNodes();
2529 // collect volumes check that number of volumss with count equal nbNode not less than 2
2530 typedef map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
2531 typedef map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
2532 TMapOfVolume mapOfVol;
2534 SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator();
2535 while ( nodeItr->more() ) {
2536 const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
2537 if ( !aNode ) continue;
2538 SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume);
2539 while ( volItr->more() ) {
2540 SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
2541 TItrMapOfVolume itr = mapOfVol.insert(make_pair(aVol, 0)).first;
2546 TItrMapOfVolume volItr = mapOfVol.begin();
2547 TItrMapOfVolume volEnd = mapOfVol.end();
2548 for ( ; volItr != volEnd; ++volItr )
2549 if ( (*volItr).second >= nbNode )
2551 // face is not free if number of volumes constructed on thier nodes more than one
2555 SMDSAbs_ElementType FreeFaces::GetType() const
2557 return SMDSAbs_Face;
2560 //================================================================================
2562 Class : LinearOrQuadratic
2563 Description : Predicate to verify whether a mesh element is linear
2565 //================================================================================
2567 LinearOrQuadratic::LinearOrQuadratic()
2572 void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh )
2577 bool LinearOrQuadratic::IsSatisfy( long theId )
2579 if (!myMesh) return false;
2580 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2581 if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) )
2583 return (!anElem->IsQuadratic());
2586 void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType )
2591 SMDSAbs_ElementType LinearOrQuadratic::GetType() const
2596 //================================================================================
2599 Description : Functor for check color of group to whic mesh element belongs to
2601 //================================================================================
2603 GroupColor::GroupColor()
2607 bool GroupColor::IsSatisfy( long theId )
2609 return (myIDs.find( theId ) != myIDs.end());
2612 void GroupColor::SetType( SMDSAbs_ElementType theType )
2617 SMDSAbs_ElementType GroupColor::GetType() const
2622 static bool isEqual( const Quantity_Color& theColor1,
2623 const Quantity_Color& theColor2 )
2625 // tolerance to compare colors
2626 const double tol = 5*1e-3;
2627 return ( fabs( theColor1.Red() - theColor2.Red() ) < tol &&
2628 fabs( theColor1.Green() - theColor2.Green() ) < tol &&
2629 fabs( theColor1.Blue() - theColor2.Blue() ) < tol );
2633 void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
2637 const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
2641 int nbGrp = aMesh->GetNbGroups();
2645 // iterates on groups and find necessary elements ids
2646 const std::set<SMESHDS_GroupBase*>& aGroups = aMesh->GetGroups();
2647 set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
2648 for (; GrIt != aGroups.end(); GrIt++) {
2649 SMESHDS_GroupBase* aGrp = (*GrIt);
2652 // check type and color of group
2653 if ( !isEqual( myColor, aGrp->GetColor() ) )
2655 if ( myType != SMDSAbs_All && myType != (SMDSAbs_ElementType)aGrp->GetType() )
2658 SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType();
2659 if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) {
2660 // add elements IDS into control
2661 int aSize = aGrp->Extent();
2662 for (int i = 0; i < aSize; i++)
2663 myIDs.insert( aGrp->GetID(i+1) );
2668 void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
2670 Kernel_Utils::Localizer loc;
2671 TCollection_AsciiString aStr = theStr;
2672 aStr.RemoveAll( ' ' );
2673 aStr.RemoveAll( '\t' );
2674 for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) )
2675 aStr.Remove( aPos, 2 );
2676 Standard_Real clr[3];
2677 clr[0] = clr[1] = clr[2] = 0.;
2678 for ( int i = 0; i < 3; i++ ) {
2679 TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 );
2680 if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() )
2681 clr[i] = tmpStr.RealValue();
2683 myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB );
2686 //=======================================================================
2687 // name : GetRangeStr
2688 // Purpose : Get range as a string.
2689 // Example: "1,2,3,50-60,63,67,70-"
2690 //=======================================================================
2692 void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
2695 theResStr += TCollection_AsciiString( myColor.Red() );
2696 theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() );
2697 theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
2700 //================================================================================
2702 Class : ElemGeomType
2703 Description : Predicate to check element geometry type
2705 //================================================================================
2707 ElemGeomType::ElemGeomType()
2710 myType = SMDSAbs_All;
2711 myGeomType = SMDSGeom_TRIANGLE;
2714 void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh )
2719 bool ElemGeomType::IsSatisfy( long theId )
2721 if (!myMesh) return false;
2722 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2725 const SMDSAbs_ElementType anElemType = anElem->GetType();
2726 if ( myType != SMDSAbs_All && anElemType != myType )
2728 bool isOk = ( anElem->GetGeomType() == myGeomType );
2732 void ElemGeomType::SetType( SMDSAbs_ElementType theType )
2737 SMDSAbs_ElementType ElemGeomType::GetType() const
2742 void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType )
2744 myGeomType = theType;
2747 SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
2752 //================================================================================
2754 Class : ElemEntityType
2755 Description : Predicate to check element entity type
2757 //================================================================================
2759 ElemEntityType::ElemEntityType():
2761 myType( SMDSAbs_All ),
2762 myEntityType( SMDSEntity_0D )
2766 void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh )
2771 bool ElemEntityType::IsSatisfy( long theId )
2773 if ( !myMesh ) return false;
2774 if ( myType == SMDSAbs_Node )
2775 return myMesh->FindNode( theId );
2776 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2778 myEntityType == anElem->GetEntityType() );
2781 void ElemEntityType::SetType( SMDSAbs_ElementType theType )
2786 SMDSAbs_ElementType ElemEntityType::GetType() const
2791 void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType )
2793 myEntityType = theEntityType;
2796 SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const
2798 return myEntityType;
2801 //================================================================================
2803 * \brief Class ConnectedElements
2805 //================================================================================
2807 ConnectedElements::ConnectedElements():
2808 myNodeID(0), myType( SMDSAbs_All ), myOkIDsReady( false ) {}
2810 SMDSAbs_ElementType ConnectedElements::GetType() const
2813 int ConnectedElements::GetNode() const
2814 { return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ
2816 std::vector<double> ConnectedElements::GetPoint() const
2819 void ConnectedElements::clearOkIDs()
2820 { myOkIDsReady = false; myOkIDs.clear(); }
2822 void ConnectedElements::SetType( SMDSAbs_ElementType theType )
2824 if ( myType != theType || myMeshModifTracer.IsMeshModified() )
2829 void ConnectedElements::SetMesh( const SMDS_Mesh* theMesh )
2831 myMeshModifTracer.SetMesh( theMesh );
2832 if ( myMeshModifTracer.IsMeshModified() )
2835 if ( !myXYZ.empty() )
2836 SetPoint( myXYZ[0], myXYZ[1], myXYZ[2] ); // find a node near myXYZ it in a new mesh
2840 void ConnectedElements::SetNode( int nodeID )
2845 bool isSameDomain = false;
2846 if ( myOkIDsReady && myMeshModifTracer.GetMesh() && !myMeshModifTracer.IsMeshModified() )
2847 if ( const SMDS_MeshNode* n = myMeshModifTracer.GetMesh()->FindNode( myNodeID ))
2849 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( myType );
2850 while ( !isSameDomain && eIt->more() )
2851 isSameDomain = IsSatisfy( eIt->next()->GetID() );
2853 if ( !isSameDomain )
2857 void ConnectedElements::SetPoint( double x, double y, double z )
2865 bool isSameDomain = false;
2867 // find myNodeID by myXYZ if possible
2868 if ( myMeshModifTracer.GetMesh() )
2870 auto_ptr<SMESH_ElementSearcher> searcher
2871 ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() ));
2873 vector< const SMDS_MeshElement* > foundElems;
2874 searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems );
2876 if ( !foundElems.empty() )
2878 myNodeID = foundElems[0]->GetNode(0)->GetID();
2879 if ( myOkIDsReady && !myMeshModifTracer.IsMeshModified() )
2880 isSameDomain = IsSatisfy( foundElems[0]->GetID() );
2883 if ( !isSameDomain )
2887 bool ConnectedElements::IsSatisfy( long theElementId )
2889 // Here we do NOT check if the mesh has changed, we do it in Set...() only!!!
2891 if ( !myOkIDsReady )
2893 if ( !myMeshModifTracer.GetMesh() )
2895 const SMDS_MeshNode* node0 = myMeshModifTracer.GetMesh()->FindNode( myNodeID );
2899 list< const SMDS_MeshNode* > nodeQueue( 1, node0 );
2900 std::set< int > checkedNodeIDs;
2902 // foreach node in nodeQueue:
2903 // foreach element sharing a node:
2904 // add ID of an element of myType to myOkIDs;
2905 // push all element nodes absent from checkedNodeIDs to nodeQueue;
2906 while ( !nodeQueue.empty() )
2908 const SMDS_MeshNode* node = nodeQueue.front();
2909 nodeQueue.pop_front();
2911 // loop on elements sharing the node
2912 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
2913 while ( eIt->more() )
2915 // keep elements of myType
2916 const SMDS_MeshElement* element = eIt->next();
2917 if ( element->GetType() == myType )
2918 myOkIDs.insert( myOkIDs.end(), element->GetID() );
2920 // enqueue nodes of the element
2921 SMDS_ElemIteratorPtr nIt = element->nodesIterator();
2922 while ( nIt->more() )
2924 const SMDS_MeshNode* n = static_cast< const SMDS_MeshNode* >( nIt->next() );
2925 if ( checkedNodeIDs.insert( n->GetID() ).second )
2926 nodeQueue.push_back( n );
2930 if ( myType == SMDSAbs_Node )
2931 std::swap( myOkIDs, checkedNodeIDs );
2933 size_t totalNbElems = myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType );
2934 if ( myOkIDs.size() == totalNbElems )
2937 myOkIDsReady = true;
2940 return myOkIDs.empty() ? true : myOkIDs.count( theElementId );
2943 //================================================================================
2945 * \brief Class CoplanarFaces
2947 //================================================================================
2949 CoplanarFaces::CoplanarFaces()
2950 : myFaceID(0), myToler(0)
2953 void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
2955 myMeshModifTracer.SetMesh( theMesh );
2956 if ( myMeshModifTracer.IsMeshModified() )
2958 // Build a set of coplanar face ids
2960 myCoplanarIDs.clear();
2962 if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler )
2965 const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID );
2966 if ( !face || face->GetType() != SMDSAbs_Face )
2970 gp_Vec myNorm = getNormale( static_cast<const SMDS_MeshFace*>(face), &normOK );
2974 const double radianTol = myToler * M_PI / 180.;
2975 std::set< SMESH_TLink > checkedLinks;
2977 std::list< pair< const SMDS_MeshElement*, gp_Vec > > faceQueue;
2978 faceQueue.push_back( make_pair( face, myNorm ));
2979 while ( !faceQueue.empty() )
2981 face = faceQueue.front().first;
2982 myNorm = faceQueue.front().second;
2983 faceQueue.pop_front();
2985 for ( int i = 0, nbN = face->NbCornerNodes(); i < nbN; ++i )
2987 const SMDS_MeshNode* n1 = face->GetNode( i );
2988 const SMDS_MeshNode* n2 = face->GetNode(( i+1 )%nbN);
2989 if ( !checkedLinks.insert( SMESH_TLink( n1, n2 )).second )
2991 SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator(SMDSAbs_Face);
2992 while ( fIt->more() )
2994 const SMDS_MeshElement* f = fIt->next();
2995 if ( f->GetNodeIndex( n2 ) > -1 )
2997 gp_Vec norm = getNormale( static_cast<const SMDS_MeshFace*>(f), &normOK );
2998 if (!normOK || myNorm.Angle( norm ) <= radianTol)
3000 myCoplanarIDs.insert( f->GetID() );
3001 faceQueue.push_back( make_pair( f, norm ));
3009 bool CoplanarFaces::IsSatisfy( long theElementId )
3011 return myCoplanarIDs.count( theElementId );
3016 *Description : Predicate for Range of Ids.
3017 * Range may be specified with two ways.
3018 * 1. Using AddToRange method
3019 * 2. With SetRangeStr method. Parameter of this method is a string
3020 * like as "1,2,3,50-60,63,67,70-"
3023 //=======================================================================
3024 // name : RangeOfIds
3025 // Purpose : Constructor
3026 //=======================================================================
3027 RangeOfIds::RangeOfIds()
3030 myType = SMDSAbs_All;
3033 //=======================================================================
3035 // Purpose : Set mesh
3036 //=======================================================================
3037 void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh )
3042 //=======================================================================
3043 // name : AddToRange
3044 // Purpose : Add ID to the range
3045 //=======================================================================
3046 bool RangeOfIds::AddToRange( long theEntityId )
3048 myIds.Add( theEntityId );
3052 //=======================================================================
3053 // name : GetRangeStr
3054 // Purpose : Get range as a string.
3055 // Example: "1,2,3,50-60,63,67,70-"
3056 //=======================================================================
3057 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
3061 TColStd_SequenceOfInteger anIntSeq;
3062 TColStd_SequenceOfAsciiString aStrSeq;
3064 TColStd_MapIteratorOfMapOfInteger anIter( myIds );
3065 for ( ; anIter.More(); anIter.Next() )
3067 int anId = anIter.Key();
3068 TCollection_AsciiString aStr( anId );
3069 anIntSeq.Append( anId );
3070 aStrSeq.Append( aStr );
3073 for ( int i = 1, n = myMin.Length(); i <= n; i++ )
3075 int aMinId = myMin( i );
3076 int aMaxId = myMax( i );
3078 TCollection_AsciiString aStr;
3079 if ( aMinId != IntegerFirst() )
3084 if ( aMaxId != IntegerLast() )
3087 // find position of the string in result sequence and insert string in it
3088 if ( anIntSeq.Length() == 0 )
3090 anIntSeq.Append( aMinId );
3091 aStrSeq.Append( aStr );
3095 if ( aMinId < anIntSeq.First() )
3097 anIntSeq.Prepend( aMinId );
3098 aStrSeq.Prepend( aStr );
3100 else if ( aMinId > anIntSeq.Last() )
3102 anIntSeq.Append( aMinId );
3103 aStrSeq.Append( aStr );
3106 for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
3107 if ( aMinId < anIntSeq( j ) )
3109 anIntSeq.InsertBefore( j, aMinId );
3110 aStrSeq.InsertBefore( j, aStr );
3116 if ( aStrSeq.Length() == 0 )
3119 theResStr = aStrSeq( 1 );
3120 for ( int j = 2, k = aStrSeq.Length(); j <= k; j++ )
3123 theResStr += aStrSeq( j );
3127 //=======================================================================
3128 // name : SetRangeStr
3129 // Purpose : Define range with string
3130 // Example of entry string: "1,2,3,50-60,63,67,70-"
3131 //=======================================================================
3132 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
3138 TCollection_AsciiString aStr = theStr;
3139 //aStr.RemoveAll( ' ' );
3140 //aStr.RemoveAll( '\t' );
3141 for ( int i = 1; i <= aStr.Length(); ++i )
3142 if ( isspace( aStr.Value( i )))
3143 aStr.SetValue( i, ',');
3145 for ( int aPos = aStr.Search( ",," ); aPos != -1; aPos = aStr.Search( ",," ) )
3146 aStr.Remove( aPos, 1 );
3148 TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
3150 while ( tmpStr != "" )
3152 tmpStr = aStr.Token( ",", i++ );
3153 int aPos = tmpStr.Search( '-' );
3157 if ( tmpStr.IsIntegerValue() )
3158 myIds.Add( tmpStr.IntegerValue() );
3164 TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
3165 TCollection_AsciiString aMinStr = tmpStr;
3167 while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
3168 while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
3170 if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) ||
3171 (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) )
3174 myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
3175 myMax.Append( aMaxStr.IsEmpty() ? IntegerLast() : aMaxStr.IntegerValue() );
3182 //=======================================================================
3184 // Purpose : Get type of supported entities
3185 //=======================================================================
3186 SMDSAbs_ElementType RangeOfIds::GetType() const
3191 //=======================================================================
3193 // Purpose : Set type of supported entities
3194 //=======================================================================
3195 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
3200 //=======================================================================
3202 // Purpose : Verify whether entity satisfies to this rpedicate
3203 //=======================================================================
3204 bool RangeOfIds::IsSatisfy( long theId )
3209 if ( myType == SMDSAbs_Node )
3211 if ( myMesh->FindNode( theId ) == 0 )
3216 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3217 if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All ))
3221 if ( myIds.Contains( theId ) )
3224 for ( int i = 1, n = myMin.Length(); i <= n; i++ )
3225 if ( theId >= myMin( i ) && theId <= myMax( i ) )
3233 Description : Base class for comparators
3235 Comparator::Comparator():
3239 Comparator::~Comparator()
3242 void Comparator::SetMesh( const SMDS_Mesh* theMesh )
3245 myFunctor->SetMesh( theMesh );
3248 void Comparator::SetMargin( double theValue )
3250 myMargin = theValue;
3253 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
3255 myFunctor = theFunct;
3258 SMDSAbs_ElementType Comparator::GetType() const
3260 return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
3263 double Comparator::GetMargin()
3271 Description : Comparator "<"
3273 bool LessThan::IsSatisfy( long theId )
3275 return myFunctor && myFunctor->GetValue( theId ) < myMargin;
3281 Description : Comparator ">"
3283 bool MoreThan::IsSatisfy( long theId )
3285 return myFunctor && myFunctor->GetValue( theId ) > myMargin;
3291 Description : Comparator "="
3294 myToler(Precision::Confusion())
3297 bool EqualTo::IsSatisfy( long theId )
3299 return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
3302 void EqualTo::SetTolerance( double theToler )
3307 double EqualTo::GetTolerance()
3314 Description : Logical NOT predicate
3316 LogicalNOT::LogicalNOT()
3319 LogicalNOT::~LogicalNOT()
3322 bool LogicalNOT::IsSatisfy( long theId )
3324 return myPredicate && !myPredicate->IsSatisfy( theId );
3327 void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh )
3330 myPredicate->SetMesh( theMesh );
3333 void LogicalNOT::SetPredicate( PredicatePtr thePred )
3335 myPredicate = thePred;
3338 SMDSAbs_ElementType LogicalNOT::GetType() const
3340 return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
3345 Class : LogicalBinary
3346 Description : Base class for binary logical predicate
3348 LogicalBinary::LogicalBinary()
3351 LogicalBinary::~LogicalBinary()
3354 void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh )
3357 myPredicate1->SetMesh( theMesh );
3360 myPredicate2->SetMesh( theMesh );
3363 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
3365 myPredicate1 = thePredicate;
3368 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
3370 myPredicate2 = thePredicate;
3373 SMDSAbs_ElementType LogicalBinary::GetType() const
3375 if ( !myPredicate1 || !myPredicate2 )
3378 SMDSAbs_ElementType aType1 = myPredicate1->GetType();
3379 SMDSAbs_ElementType aType2 = myPredicate2->GetType();
3381 return aType1 == aType2 ? aType1 : SMDSAbs_All;
3387 Description : Logical AND
3389 bool LogicalAND::IsSatisfy( long theId )
3394 myPredicate1->IsSatisfy( theId ) &&
3395 myPredicate2->IsSatisfy( theId );
3401 Description : Logical OR
3403 bool LogicalOR::IsSatisfy( long theId )
3408 (myPredicate1->IsSatisfy( theId ) ||
3409 myPredicate2->IsSatisfy( theId ));
3418 // #include <tbb/parallel_for.h>
3419 // #include <tbb/enumerable_thread_specific.h>
3421 // namespace Parallel
3423 // typedef tbb::enumerable_thread_specific< TIdSequence > TIdSeq;
3427 // const SMDS_Mesh* myMesh;
3428 // PredicatePtr myPredicate;
3429 // TIdSeq & myOKIds;
3430 // Predicate( const SMDS_Mesh* m, PredicatePtr p, TIdSeq & ids ):
3431 // myMesh(m), myPredicate(p->Duplicate()), myOKIds(ids) {}
3432 // void operator() ( const tbb::blocked_range<size_t>& r ) const
3434 // for ( size_t i = r.begin(); i != r.end(); ++i )
3435 // if ( myPredicate->IsSatisfy( i ))
3436 // myOKIds.local().push_back();
3448 void Filter::SetPredicate( PredicatePtr thePredicate )
3450 myPredicate = thePredicate;
3453 void Filter::GetElementsId( const SMDS_Mesh* theMesh,
3454 PredicatePtr thePredicate,
3455 TIdSequence& theSequence )
3457 theSequence.clear();
3459 if ( !theMesh || !thePredicate )
3462 thePredicate->SetMesh( theMesh );
3464 SMDS_ElemIteratorPtr elemIt = theMesh->elementsIterator( thePredicate->GetType() );
3466 while ( elemIt->more() ) {
3467 const SMDS_MeshElement* anElem = elemIt->next();
3468 long anId = anElem->GetID();
3469 if ( thePredicate->IsSatisfy( anId ) )
3470 theSequence.push_back( anId );
3475 void Filter::GetElementsId( const SMDS_Mesh* theMesh,
3476 Filter::TIdSequence& theSequence )
3478 GetElementsId(theMesh,myPredicate,theSequence);
3485 typedef std::set<SMDS_MeshFace*> TMapOfFacePtr;
3491 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
3492 SMDS_MeshNode* theNode2 )
3498 ManifoldPart::Link::~Link()
3504 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
3506 if ( myNode1 == theLink.myNode1 &&
3507 myNode2 == theLink.myNode2 )
3509 else if ( myNode1 == theLink.myNode2 &&
3510 myNode2 == theLink.myNode1 )
3516 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
3518 if(myNode1 < x.myNode1) return true;
3519 if(myNode1 == x.myNode1)
3520 if(myNode2 < x.myNode2) return true;
3524 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
3525 const ManifoldPart::Link& theLink2 )
3527 return theLink1.IsEqual( theLink2 );
3530 ManifoldPart::ManifoldPart()
3533 myAngToler = Precision::Angular();
3534 myIsOnlyManifold = true;
3537 ManifoldPart::~ManifoldPart()
3542 void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh )
3548 SMDSAbs_ElementType ManifoldPart::GetType() const
3549 { return SMDSAbs_Face; }
3551 bool ManifoldPart::IsSatisfy( long theElementId )
3553 return myMapIds.Contains( theElementId );
3556 void ManifoldPart::SetAngleTolerance( const double theAngToler )
3557 { myAngToler = theAngToler; }
3559 double ManifoldPart::GetAngleTolerance() const
3560 { return myAngToler; }
3562 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
3563 { myIsOnlyManifold = theIsOnly; }
3565 void ManifoldPart::SetStartElem( const long theStartId )
3566 { myStartElemId = theStartId; }
3568 bool ManifoldPart::process()
3571 myMapBadGeomIds.Clear();
3573 myAllFacePtr.clear();
3574 myAllFacePtrIntDMap.clear();
3578 // collect all faces into own map
3579 SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
3580 for (; anFaceItr->more(); )
3582 SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
3583 myAllFacePtr.push_back( aFacePtr );
3584 myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
3587 SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
3591 // the map of non manifold links and bad geometry
3592 TMapOfLink aMapOfNonManifold;
3593 TColStd_MapOfInteger aMapOfTreated;
3595 // begin cycle on faces from start index and run on vector till the end
3596 // and from begin to start index to cover whole vector
3597 const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
3598 bool isStartTreat = false;
3599 for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
3601 if ( fi == aStartIndx )
3602 isStartTreat = true;
3603 // as result next time when fi will be equal to aStartIndx
3605 SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
3606 if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
3609 aMapOfTreated.Add( aFacePtr->GetID() );
3610 TColStd_MapOfInteger aResFaces;
3611 if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
3612 aMapOfNonManifold, aResFaces ) )
3614 TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
3615 for ( ; anItr.More(); anItr.Next() )
3617 int aFaceId = anItr.Key();
3618 aMapOfTreated.Add( aFaceId );
3619 myMapIds.Add( aFaceId );
3622 if ( fi == ( myAllFacePtr.size() - 1 ) )
3624 } // end run on vector of faces
3625 return !myMapIds.IsEmpty();
3628 static void getLinks( const SMDS_MeshFace* theFace,
3629 ManifoldPart::TVectorOfLink& theLinks )
3631 int aNbNode = theFace->NbNodes();
3632 SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
3634 SMDS_MeshNode* aNode = 0;
3635 for ( ; aNodeItr->more() && i <= aNbNode; )
3638 SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
3642 SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
3644 ManifoldPart::Link aLink( aN1, aN2 );
3645 theLinks.push_back( aLink );
3649 bool ManifoldPart::findConnected
3650 ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
3651 SMDS_MeshFace* theStartFace,
3652 ManifoldPart::TMapOfLink& theNonManifold,
3653 TColStd_MapOfInteger& theResFaces )
3655 theResFaces.Clear();
3656 if ( !theAllFacePtrInt.size() )
3659 if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
3661 myMapBadGeomIds.Add( theStartFace->GetID() );
3665 ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
3666 ManifoldPart::TVectorOfLink aSeqOfBoundary;
3667 theResFaces.Add( theStartFace->GetID() );
3668 ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
3670 expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3671 aDMapLinkFace, theNonManifold, theStartFace );
3673 bool isDone = false;
3674 while ( !isDone && aMapOfBoundary.size() != 0 )
3676 bool isToReset = false;
3677 ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
3678 for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
3680 ManifoldPart::Link aLink = *pLink;
3681 if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
3683 // each link could be treated only once
3684 aMapToSkip.insert( aLink );
3686 ManifoldPart::TVectorOfFacePtr aFaces;
3688 if ( myIsOnlyManifold &&
3689 (theNonManifold.find( aLink ) != theNonManifold.end()) )
3693 getFacesByLink( aLink, aFaces );
3694 // filter the element to keep only indicated elements
3695 ManifoldPart::TVectorOfFacePtr aFiltered;
3696 ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3697 for ( ; pFace != aFaces.end(); ++pFace )
3699 SMDS_MeshFace* aFace = *pFace;
3700 if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
3701 aFiltered.push_back( aFace );
3704 if ( aFaces.size() < 2 ) // no neihgbour faces
3706 else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
3708 theNonManifold.insert( aLink );
3713 // compare normal with normals of neighbor element
3714 SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
3715 ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3716 for ( ; pFace != aFaces.end(); ++pFace )
3718 SMDS_MeshFace* aNextFace = *pFace;
3719 if ( aPrevFace == aNextFace )
3721 int anNextFaceID = aNextFace->GetID();
3722 if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
3723 // should not be with non manifold restriction. probably bad topology
3725 // check if face was treated and skipped
3726 if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
3727 !isInPlane( aPrevFace, aNextFace ) )
3729 // add new element to connected and extend the boundaries.
3730 theResFaces.Add( anNextFaceID );
3731 expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3732 aDMapLinkFace, theNonManifold, aNextFace );
3736 isDone = !isToReset;
3739 return !theResFaces.IsEmpty();
3742 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
3743 const SMDS_MeshFace* theFace2 )
3745 gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
3746 gp_XYZ aNorm2XYZ = getNormale( theFace2 );
3747 if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
3749 myMapBadGeomIds.Add( theFace2->GetID() );
3752 if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
3758 void ManifoldPart::expandBoundary
3759 ( ManifoldPart::TMapOfLink& theMapOfBoundary,
3760 ManifoldPart::TVectorOfLink& theSeqOfBoundary,
3761 ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
3762 ManifoldPart::TMapOfLink& theNonManifold,
3763 SMDS_MeshFace* theNextFace ) const
3765 ManifoldPart::TVectorOfLink aLinks;
3766 getLinks( theNextFace, aLinks );
3767 int aNbLink = (int)aLinks.size();
3768 for ( int i = 0; i < aNbLink; i++ )
3770 ManifoldPart::Link aLink = aLinks[ i ];
3771 if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
3773 if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
3775 if ( myIsOnlyManifold )
3777 // remove from boundary
3778 theMapOfBoundary.erase( aLink );
3779 ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
3780 for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
3782 ManifoldPart::Link aBoundLink = *pLink;
3783 if ( aBoundLink.IsEqual( aLink ) )
3785 theSeqOfBoundary.erase( pLink );
3793 theMapOfBoundary.insert( aLink );
3794 theSeqOfBoundary.push_back( aLink );
3795 theDMapLinkFacePtr[ aLink ] = theNextFace;
3800 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
3801 ManifoldPart::TVectorOfFacePtr& theFaces ) const
3803 std::set<SMDS_MeshCell *> aSetOfFaces;
3804 // take all faces that shared first node
3805 SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
3806 for ( ; anItr->more(); )
3808 SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
3811 aSetOfFaces.insert( aFace );
3813 // take all faces that shared second node
3814 anItr = theLink.myNode2->facesIterator();
3815 // find the common part of two sets
3816 for ( ; anItr->more(); )
3818 SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
3819 if ( aSetOfFaces.count( aFace ) )
3820 theFaces.push_back( aFace );
3825 Class : BelongToMeshGroup
3826 Description : Verify whether a mesh element is included into a mesh group
3828 BelongToMeshGroup::BelongToMeshGroup(): myGroup( 0 )
3832 void BelongToMeshGroup::SetGroup( SMESHDS_GroupBase* g )
3837 void BelongToMeshGroup::SetStoreName( const std::string& sn )
3842 void BelongToMeshGroup::SetMesh( const SMDS_Mesh* theMesh )
3844 if ( myGroup && myGroup->GetMesh() != theMesh )
3848 if ( !myGroup && !myStoreName.empty() )
3850 if ( const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh))
3852 const std::set<SMESHDS_GroupBase*>& grps = aMesh->GetGroups();
3853 std::set<SMESHDS_GroupBase*>::const_iterator g = grps.begin();
3854 for ( ; g != grps.end() && !myGroup; ++g )
3855 if ( *g && myStoreName == (*g)->GetStoreName() )
3861 myGroup->IsEmpty(); // make GroupOnFilter update its predicate
3865 bool BelongToMeshGroup::IsSatisfy( long theElementId )
3867 return myGroup ? myGroup->Contains( theElementId ) : false;
3870 SMDSAbs_ElementType BelongToMeshGroup::GetType() const
3872 return myGroup ? myGroup->GetType() : SMDSAbs_All;
3879 ElementsOnSurface::ElementsOnSurface()
3882 myType = SMDSAbs_All;
3884 myToler = Precision::Confusion();
3885 myUseBoundaries = false;
3888 ElementsOnSurface::~ElementsOnSurface()
3892 void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
3894 myMeshModifTracer.SetMesh( theMesh );
3895 if ( myMeshModifTracer.IsMeshModified())
3899 bool ElementsOnSurface::IsSatisfy( long theElementId )
3901 return myIds.Contains( theElementId );
3904 SMDSAbs_ElementType ElementsOnSurface::GetType() const
3907 void ElementsOnSurface::SetTolerance( const double theToler )
3909 if ( myToler != theToler )
3914 double ElementsOnSurface::GetTolerance() const
3917 void ElementsOnSurface::SetUseBoundaries( bool theUse )
3919 if ( myUseBoundaries != theUse ) {
3920 myUseBoundaries = theUse;
3921 SetSurface( mySurf, myType );
3925 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
3926 const SMDSAbs_ElementType theType )
3931 if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
3933 mySurf = TopoDS::Face( theShape );
3934 BRepAdaptor_Surface SA( mySurf, myUseBoundaries );
3936 u1 = SA.FirstUParameter(),
3937 u2 = SA.LastUParameter(),
3938 v1 = SA.FirstVParameter(),
3939 v2 = SA.LastVParameter();
3940 Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf );
3941 myProjector.Init( surf, u1,u2, v1,v2 );
3945 void ElementsOnSurface::process()
3948 if ( mySurf.IsNull() )
3951 if ( !myMeshModifTracer.GetMesh() )
3954 myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
3956 SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType );
3957 for(; anIter->more(); )
3958 process( anIter->next() );
3961 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
3963 SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
3964 bool isSatisfy = true;
3965 for ( ; aNodeItr->more(); )
3967 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
3968 if ( !isOnSurface( aNode ) )
3975 myIds.Add( theElemPtr->GetID() );
3978 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode )
3980 if ( mySurf.IsNull() )
3983 gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
3984 // double aToler2 = myToler * myToler;
3985 // if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
3987 // gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
3988 // if ( aPln.SquareDistance( aPnt ) > aToler2 )
3991 // else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
3993 // gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
3994 // double aRad = aCyl.Radius();
3995 // gp_Ax3 anAxis = aCyl.Position();
3996 // gp_XYZ aLoc = aCyl.Location().XYZ();
3997 // double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
3998 // double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
3999 // if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
4004 myProjector.Perform( aPnt );
4005 bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler );
4015 ElementsOnShape::ElementsOnShape()
4017 myType(SMDSAbs_All),
4018 myToler(Precision::Confusion()),
4019 myAllNodesFlag(false)
4023 ElementsOnShape::~ElementsOnShape()
4028 SMDSAbs_ElementType ElementsOnShape::GetType() const
4033 void ElementsOnShape::SetTolerance (const double theToler)
4035 if (myToler != theToler) {
4037 SetShape(myShape, myType);
4041 double ElementsOnShape::GetTolerance() const
4046 void ElementsOnShape::SetAllNodes (bool theAllNodes)
4048 myAllNodesFlag = theAllNodes;
4051 void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
4053 myMeshModifTracer.SetMesh( theMesh );
4054 if ( myMeshModifTracer.IsMeshModified())
4056 size_t nbNodes = theMesh ? theMesh->NbNodes() : 0;
4057 if ( myNodeIsChecked.size() == nbNodes )
4059 std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4063 SMESHUtils::FreeVector( myNodeIsChecked );
4064 SMESHUtils::FreeVector( myNodeIsOut );
4065 myNodeIsChecked.resize( nbNodes, false );
4066 myNodeIsOut.resize( nbNodes );
4071 bool ElementsOnShape::getNodeIsOut( const SMDS_MeshNode* n, bool& isOut )
4073 if ( n->GetID() >= (int) myNodeIsChecked.size() ||
4074 !myNodeIsChecked[ n->GetID() ])
4077 isOut = myNodeIsOut[ n->GetID() ];
4081 void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool isOut )
4083 if ( n->GetID() < (int) myNodeIsChecked.size() )
4085 myNodeIsChecked[ n->GetID() ] = true;
4086 myNodeIsOut [ n->GetID() ] = isOut;
4090 void ElementsOnShape::SetShape (const TopoDS_Shape& theShape,
4091 const SMDSAbs_ElementType theType)
4095 if ( myShape.IsNull() ) return;
4097 TopTools_IndexedMapOfShape shapesMap;
4098 TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX };
4099 TopExp_Explorer sub;
4100 for ( int i = 0; i < 4; ++i )
4102 if ( shapesMap.IsEmpty() )
4103 for ( sub.Init( myShape, shapeTypes[i] ); sub.More(); sub.Next() )
4104 shapesMap.Add( sub.Current() );
4106 for ( sub.Init( myShape, shapeTypes[i], shapeTypes[i-1] ); sub.More(); sub.Next() )
4107 shapesMap.Add( sub.Current() );
4111 myClassifiers.resize( shapesMap.Extent() );
4112 for ( int i = 0; i < shapesMap.Extent(); ++i )
4113 myClassifiers[ i ] = new TClassifier( shapesMap( i+1 ), myToler );
4115 if ( theType == SMDSAbs_Node )
4117 SMESHUtils::FreeVector( myNodeIsChecked );
4118 SMESHUtils::FreeVector( myNodeIsOut );
4122 std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4126 void ElementsOnShape::clearClassifiers()
4128 for ( size_t i = 0; i < myClassifiers.size(); ++i )
4129 delete myClassifiers[ i ];
4130 myClassifiers.clear();
4133 bool ElementsOnShape::IsSatisfy (long elemId)
4135 const SMDS_Mesh* mesh = myMeshModifTracer.GetMesh();
4136 const SMDS_MeshElement* elem =
4137 ( myType == SMDSAbs_Node ? mesh->FindNode( elemId ) : mesh->FindElement( elemId ));
4138 if ( !elem || myClassifiers.empty() )
4141 bool isSatisfy = myAllNodesFlag, isNodeOut;
4143 gp_XYZ centerXYZ (0, 0, 0);
4145 SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator();
4146 while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
4148 SMESH_TNodeXYZ aPnt( aNodeItr->next() );
4152 if ( !getNodeIsOut( aPnt._node, isNodeOut ))
4154 for ( size_t i = 0; i < myClassifiers.size() && isNodeOut; ++i )
4155 isNodeOut = myClassifiers[i]->IsOut( aPnt );
4157 setNodeIsOut( aPnt._node, isNodeOut );
4159 isSatisfy = !isNodeOut;
4162 // Check the center point for volumes MantisBug 0020168
4165 myClassifiers[0]->ShapeType() == TopAbs_SOLID)
4167 centerXYZ /= elem->NbNodes();
4169 for ( size_t i = 0; i < myClassifiers.size() && !isSatisfy; ++i )
4170 isSatisfy = ! myClassifiers[i]->IsOut( centerXYZ );
4176 TopAbs_ShapeEnum ElementsOnShape::TClassifier::ShapeType() const
4178 return myShape.ShapeType();
4181 bool ElementsOnShape::TClassifier::IsOut(const gp_Pnt& p)
4183 return (this->*myIsOutFun)( p );
4186 void ElementsOnShape::TClassifier::Init (const TopoDS_Shape& theShape, double theTol)
4190 switch ( myShape.ShapeType() )
4192 case TopAbs_SOLID: {
4193 if ( isBox( theShape ))
4195 myIsOutFun = & ElementsOnShape::TClassifier::isOutOfBox;
4199 mySolidClfr.Load(theShape);
4200 myIsOutFun = & ElementsOnShape::TClassifier::isOutOfSolid;
4205 Standard_Real u1,u2,v1,v2;
4206 Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape ));
4207 surf->Bounds( u1,u2,v1,v2 );
4208 myProjFace.Init(surf, u1,u2, v1,v2, myTol );
4209 myIsOutFun = & ElementsOnShape::TClassifier::isOutOfFace;
4213 Standard_Real u1, u2;
4214 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge(theShape), u1, u2);
4215 myProjEdge.Init(curve, u1, u2);
4216 myIsOutFun = & ElementsOnShape::TClassifier::isOutOfEdge;
4219 case TopAbs_VERTEX:{
4220 myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) );
4221 myIsOutFun = & ElementsOnShape::TClassifier::isOutOfVertex;
4225 throw SALOME_Exception("Programmer error in usage of ElementsOnShape::TClassifier");
4229 bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p)
4231 mySolidClfr.Perform( p, myTol );
4232 return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON );
4235 bool ElementsOnShape::TClassifier::isOutOfBox (const gp_Pnt& p)
4237 return myBox.IsOut( p.XYZ() );
4240 bool ElementsOnShape::TClassifier::isOutOfFace (const gp_Pnt& p)
4242 myProjFace.Perform( p );
4243 if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol )
4245 // check relatively to the face
4246 Quantity_Parameter u, v;
4247 myProjFace.LowerDistanceParameters(u, v);
4248 gp_Pnt2d aProjPnt (u, v);
4249 BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol );
4250 if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON )
4256 bool ElementsOnShape::TClassifier::isOutOfEdge (const gp_Pnt& p)
4258 myProjEdge.Perform( p );
4259 return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
4262 bool ElementsOnShape::TClassifier::isOutOfVertex(const gp_Pnt& p)
4264 return ( myVertexXYZ.Distance( p ) > myTol );
4267 bool ElementsOnShape::TClassifier::isBox (const TopoDS_Shape& theShape)
4269 TopTools_IndexedMapOfShape vMap;
4270 TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
4271 if ( vMap.Extent() != 8 )
4275 for ( int i = 1; i <= 8; ++i )
4276 myBox.Add( BRep_Tool::Pnt( TopoDS::Vertex( vMap( i ))).XYZ() );
4278 gp_XYZ pMin = myBox.CornerMin(), pMax = myBox.CornerMax();
4279 for ( int i = 1; i <= 8; ++i )
4281 gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vMap( i )));
4282 for ( int iC = 1; iC <= 3; ++ iC )
4284 double d1 = Abs( pMin.Coord( iC ) - p.Coord( iC ));
4285 double d2 = Abs( pMax.Coord( iC ) - p.Coord( iC ));
4286 if ( Min( d1, d2 ) > myTol )
4290 myBox.Enlarge( myTol );
4296 Class : BelongToGeom
4297 Description : Predicate for verifying whether entity belongs to
4298 specified geometrical support
4301 BelongToGeom::BelongToGeom()
4303 myType(SMDSAbs_All),
4304 myIsSubshape(false),
4305 myTolerance(Precision::Confusion())
4308 void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh )
4310 myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
4314 void BelongToGeom::SetGeom( const TopoDS_Shape& theShape )
4320 static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap,
4321 const TopoDS_Shape& theShape)
4323 if (theMap.Contains(theShape)) return true;
4325 if (theShape.ShapeType() == TopAbs_COMPOUND ||
4326 theShape.ShapeType() == TopAbs_COMPSOLID)
4328 TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
4329 for (; anIt.More(); anIt.Next())
4331 if (!IsSubShape(theMap, anIt.Value())) {
4341 void BelongToGeom::init()
4343 if (!myMeshDS || myShape.IsNull()) return;
4345 // is sub-shape of main shape?
4346 TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
4347 if (aMainShape.IsNull()) {
4348 myIsSubshape = false;
4351 TopTools_IndexedMapOfShape aMap;
4352 TopExp::MapShapes(aMainShape, aMap);
4353 myIsSubshape = IsSubShape(aMap, myShape);
4356 //if (!myIsSubshape) // to be always ready to check an element not bound to geometry
4358 myElementsOnShapePtr.reset(new ElementsOnShape());
4359 myElementsOnShapePtr->SetTolerance(myTolerance);
4360 myElementsOnShapePtr->SetAllNodes(true); // "belong", while false means "lays on"
4361 myElementsOnShapePtr->SetMesh(myMeshDS);
4362 myElementsOnShapePtr->SetShape(myShape, myType);
4366 static bool IsContains( const SMESHDS_Mesh* theMeshDS,
4367 const TopoDS_Shape& theShape,
4368 const SMDS_MeshElement* theElem,
4369 TopAbs_ShapeEnum theFindShapeEnum,
4370 TopAbs_ShapeEnum theAvoidShapeEnum = TopAbs_SHAPE )
4372 TopExp_Explorer anExp( theShape,theFindShapeEnum,theAvoidShapeEnum );
4374 while( anExp.More() )
4376 const TopoDS_Shape& aShape = anExp.Current();
4377 if( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape ) ){
4378 if( aSubMesh->Contains( theElem ) )
4386 bool BelongToGeom::IsSatisfy (long theId)
4388 if (myMeshDS == 0 || myShape.IsNull())
4393 return myElementsOnShapePtr->IsSatisfy(theId);
4397 if (myType == SMDSAbs_Node)
4399 if( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ) )
4401 if ( aNode->getshapeId() < 1 )
4402 return myElementsOnShapePtr->IsSatisfy(theId);
4404 const SMDS_PositionPtr& aPosition = aNode->GetPosition();
4405 SMDS_TypeOfPosition aTypeOfPosition = aPosition->GetTypeOfPosition();
4406 switch( aTypeOfPosition )
4408 case SMDS_TOP_VERTEX : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_VERTEX ));
4409 case SMDS_TOP_EDGE : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_EDGE ));
4410 case SMDS_TOP_FACE : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_FACE ));
4411 case SMDS_TOP_3DSPACE: return ( IsContains( myMeshDS,myShape,aNode,TopAbs_SOLID ) ||
4412 IsContains( myMeshDS,myShape,aNode,TopAbs_SHELL ));
4418 if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ))
4420 if ( anElem->getshapeId() < 1 )
4421 return myElementsOnShapePtr->IsSatisfy(theId);
4423 if( myType == SMDSAbs_All )
4425 return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE ) ||
4426 IsContains( myMeshDS,myShape,anElem,TopAbs_FACE ) ||
4427 IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )||
4428 IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL ));
4430 else if( myType == anElem->GetType() )
4434 case SMDSAbs_Edge : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE ));
4435 case SMDSAbs_Face : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_FACE ));
4436 case SMDSAbs_Volume: return ( IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )||
4437 IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL ));
4446 void BelongToGeom::SetType (SMDSAbs_ElementType theType)
4452 SMDSAbs_ElementType BelongToGeom::GetType() const
4457 TopoDS_Shape BelongToGeom::GetShape()
4462 const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const
4467 void BelongToGeom::SetTolerance (double theTolerance)
4469 myTolerance = theTolerance;
4474 double BelongToGeom::GetTolerance()
4481 Description : Predicate for verifying whether entiy lying or partially lying on
4482 specified geometrical support
4485 LyingOnGeom::LyingOnGeom()
4487 myType(SMDSAbs_All),
4488 myIsSubshape(false),
4489 myTolerance(Precision::Confusion())
4492 void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh )
4494 myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
4498 void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape )
4504 void LyingOnGeom::init()
4506 if (!myMeshDS || myShape.IsNull()) return;
4508 // is sub-shape of main shape?
4509 TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
4510 if (aMainShape.IsNull()) {
4511 myIsSubshape = false;
4514 myIsSubshape = myMeshDS->IsGroupOfSubShapes( myShape );
4519 TopTools_IndexedMapOfShape shapes;
4520 TopExp::MapShapes( myShape, shapes );
4521 mySubShapesIDs.Clear();
4522 for ( int i = 1; i <= shapes.Extent(); ++i )
4524 int subID = myMeshDS->ShapeToIndex( shapes( i ));
4526 mySubShapesIDs.Add( subID );
4531 myElementsOnShapePtr.reset(new ElementsOnShape());
4532 myElementsOnShapePtr->SetTolerance(myTolerance);
4533 myElementsOnShapePtr->SetAllNodes(false); // lays on, while true means "belong"
4534 myElementsOnShapePtr->SetMesh(myMeshDS);
4535 myElementsOnShapePtr->SetShape(myShape, myType);
4539 bool LyingOnGeom::IsSatisfy( long theId )
4541 if ( myMeshDS == 0 || myShape.IsNull() )
4546 return myElementsOnShapePtr->IsSatisfy(theId);
4551 const SMDS_MeshElement* elem =
4552 ( myType == SMDSAbs_Node ) ? myMeshDS->FindNode( theId ) : myMeshDS->FindElement( theId );
4554 if ( mySubShapesIDs.Contains( elem->getshapeId() ))
4557 if ( elem->GetType() != SMDSAbs_Node )
4559 SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator();
4560 while ( nodeItr->more() )
4562 const SMDS_MeshElement* aNode = nodeItr->next();
4563 if ( mySubShapesIDs.Contains( aNode->getshapeId() ))
4571 void LyingOnGeom::SetType( SMDSAbs_ElementType theType )
4577 SMDSAbs_ElementType LyingOnGeom::GetType() const
4582 TopoDS_Shape LyingOnGeom::GetShape()
4587 const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const
4592 void LyingOnGeom::SetTolerance (double theTolerance)
4594 myTolerance = theTolerance;
4599 double LyingOnGeom::GetTolerance()
4604 bool LyingOnGeom::Contains( const SMESHDS_Mesh* theMeshDS,
4605 const TopoDS_Shape& theShape,
4606 const SMDS_MeshElement* theElem,
4607 TopAbs_ShapeEnum theFindShapeEnum,
4608 TopAbs_ShapeEnum theAvoidShapeEnum )
4610 // if (IsContains(theMeshDS, theShape, theElem, theFindShapeEnum, theAvoidShapeEnum))
4613 // TopTools_MapOfShape aSubShapes;
4614 // TopExp_Explorer exp( theShape, theFindShapeEnum, theAvoidShapeEnum );
4615 // for ( ; exp.More(); exp.Next() )
4617 // const TopoDS_Shape& aShape = exp.Current();
4618 // if ( !aSubShapes.Add( aShape )) continue;
4620 // if ( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape ))
4622 // if ( aSubMesh->Contains( theElem ))
4625 // SMDS_ElemIteratorPtr nodeItr = theElem->nodesIterator();
4626 // while ( nodeItr->more() )
4628 // const SMDS_MeshElement* aNode = nodeItr->next();
4629 // if ( aSubMesh->Contains( aNode ))
4637 TSequenceOfXYZ::TSequenceOfXYZ()
4640 TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n)
4643 TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t)
4646 TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray)
4649 template <class InputIterator>
4650 TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd)
4653 TSequenceOfXYZ::~TSequenceOfXYZ()
4656 TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
4658 myArray = theSequenceOfXYZ.myArray;
4662 gp_XYZ& TSequenceOfXYZ::operator()(size_type n)
4664 return myArray[n-1];
4667 const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const
4669 return myArray[n-1];
4672 void TSequenceOfXYZ::clear()
4677 void TSequenceOfXYZ::reserve(size_type n)
4682 void TSequenceOfXYZ::push_back(const gp_XYZ& v)
4684 myArray.push_back(v);
4687 TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
4689 return myArray.size();
4692 TMeshModifTracer::TMeshModifTracer():
4693 myMeshModifTime(0), myMesh(0)
4696 void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh )
4698 if ( theMesh != myMesh )
4699 myMeshModifTime = 0;
4702 bool TMeshModifTracer::IsMeshModified()
4704 bool modified = false;
4707 modified = ( myMeshModifTime != myMesh->GetMTime() );
4708 myMeshModifTime = myMesh->GetMTime();