1 // Copyright (C) 2007-2016 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_GroupOnFilter.hxx"
35 #include "SMESHDS_Mesh.hxx"
36 #include "SMESH_MeshAlgos.hxx"
37 #include "SMESH_OctreeNode.hxx"
39 #include <Basics_Utils.hxx>
41 #include <BRepAdaptor_Surface.hxx>
42 #include <BRepBndLib.hxx>
43 #include <BRepClass_FaceClassifier.hxx>
44 #include <BRep_Tool.hxx>
45 #include <Geom_CylindricalSurface.hxx>
46 #include <Geom_Plane.hxx>
47 #include <Geom_Surface.hxx>
48 #include <NCollection_Map.hxx>
49 #include <Precision.hxx>
50 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
51 #include <TColStd_MapOfInteger.hxx>
52 #include <TColStd_SequenceOfAsciiString.hxx>
53 #include <TColgp_Array1OfXYZ.hxx>
57 #include <TopoDS_Edge.hxx>
58 #include <TopoDS_Face.hxx>
59 #include <TopoDS_Iterator.hxx>
60 #include <TopoDS_Shape.hxx>
61 #include <TopoDS_Vertex.hxx>
63 #include <gp_Cylinder.hxx>
70 #include <vtkMeshQuality.h>
81 const double theEps = 1e-100;
82 const double theInf = 1e+100;
84 inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode )
86 return gp_XYZ(aNode->X(), aNode->Y(), aNode->Z() );
89 inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
91 gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
93 return v1.Magnitude() < gp::Resolution() ||
94 v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
97 inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
99 gp_Vec aVec1( P2 - P1 );
100 gp_Vec aVec2( P3 - P1 );
101 return ( aVec1 ^ aVec2 ).Magnitude() * 0.5;
104 inline double getArea( const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3 )
106 return getArea( P1.XYZ(), P2.XYZ(), P3.XYZ() );
111 inline double getDistance( const gp_XYZ& P1, const gp_XYZ& P2 )
113 double aDist = gp_Pnt( P1 ).Distance( gp_Pnt( P2 ) );
117 int getNbMultiConnection( const SMDS_Mesh* theMesh, const int theId )
122 const SMDS_MeshElement* anEdge = theMesh->FindElement( theId );
123 if ( anEdge == 0 || anEdge->GetType() != SMDSAbs_Edge/* || anEdge->NbNodes() != 2 */)
126 // for each pair of nodes in anEdge (there are 2 pairs in a quadratic edge)
127 // count elements containing both nodes of the pair.
128 // Note that there may be such cases for a quadratic edge (a horizontal line):
133 // +-----+------+ +-----+------+
136 // result sould be 2 in both cases
138 int aResult0 = 0, aResult1 = 0;
139 // last node, it is a medium one in a quadratic edge
140 const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 );
141 const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 );
142 const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 );
143 if ( aNode1 == aLastNode ) aNode1 = 0;
145 SMDS_ElemIteratorPtr anElemIter = aLastNode->GetInverseElementIterator();
146 while( anElemIter->more() ) {
147 const SMDS_MeshElement* anElem = anElemIter->next();
148 if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
149 SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
150 while ( anIter->more() ) {
151 if ( const SMDS_MeshElement* anElemNode = anIter->next() ) {
152 if ( anElemNode == aNode0 ) {
154 if ( !aNode1 ) break; // not a quadratic edge
156 else if ( anElemNode == aNode1 )
162 int aResult = std::max ( aResult0, aResult1 );
167 gp_XYZ getNormale( const SMDS_MeshFace* theFace, bool* ok=0 )
169 int aNbNode = theFace->NbNodes();
171 gp_XYZ q1 = gpXYZ( theFace->GetNode(1)) - gpXYZ( theFace->GetNode(0));
172 gp_XYZ q2 = gpXYZ( theFace->GetNode(2)) - gpXYZ( theFace->GetNode(0));
175 gp_XYZ q3 = gpXYZ( theFace->GetNode(3)) - gpXYZ( theFace->GetNode(0));
178 double len = n.Modulus();
179 bool zeroLen = ( len <= std::numeric_limits<double>::min());
183 if (ok) *ok = !zeroLen;
191 using namespace SMESH::Controls;
197 //================================================================================
199 Class : NumericalFunctor
200 Description : Base class for numerical functors
202 //================================================================================
204 NumericalFunctor::NumericalFunctor():
210 void NumericalFunctor::SetMesh( const SMDS_Mesh* theMesh )
215 bool NumericalFunctor::GetPoints(const int theId,
216 TSequenceOfXYZ& theRes ) const
223 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
224 if ( !anElem || anElem->GetType() != this->GetType() )
227 return GetPoints( anElem, theRes );
230 bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
231 TSequenceOfXYZ& theRes )
238 theRes.reserve( anElem->NbNodes() );
239 theRes.setElement( anElem );
241 // Get nodes of the element
242 SMDS_ElemIteratorPtr anIter;
244 if ( anElem->IsQuadratic() ) {
245 switch ( anElem->GetType() ) {
247 anIter = dynamic_cast<const SMDS_VtkEdge*>
248 (anElem)->interlacedNodesElemIterator();
251 anIter = dynamic_cast<const SMDS_VtkFace*>
252 (anElem)->interlacedNodesElemIterator();
255 anIter = anElem->nodesIterator();
259 anIter = anElem->nodesIterator();
264 while( anIter->more() ) {
265 if ( const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( anIter->next() ))
267 aNode->GetXYZ( xyz );
268 theRes.push_back( gp_XYZ( xyz[0], xyz[1], xyz[2] ));
276 long NumericalFunctor::GetPrecision() const
281 void NumericalFunctor::SetPrecision( const long thePrecision )
283 myPrecision = thePrecision;
284 myPrecisionValue = pow( 10., (double)( myPrecision ) );
287 double NumericalFunctor::GetValue( long theId )
291 myCurrElement = myMesh->FindElement( theId );
294 if ( GetPoints( theId, P )) // elem type is checked here
295 aVal = Round( GetValue( P ));
300 double NumericalFunctor::Round( const double & aVal )
302 return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal;
305 //================================================================================
307 * \brief Return histogram of functor values
308 * \param nbIntervals - number of intervals
309 * \param nbEvents - number of mesh elements having values within i-th interval
310 * \param funValues - boundaries of intervals
311 * \param elements - elements to check vulue of; empty list means "of all"
312 * \param minmax - boundaries of diapason of values to divide into intervals
314 //================================================================================
316 void NumericalFunctor::GetHistogram(int nbIntervals,
317 std::vector<int>& nbEvents,
318 std::vector<double>& funValues,
319 const std::vector<int>& elements,
320 const double* minmax,
321 const bool isLogarithmic)
323 if ( nbIntervals < 1 ||
325 !myMesh->GetMeshInfo().NbElements( GetType() ))
327 nbEvents.resize( nbIntervals, 0 );
328 funValues.resize( nbIntervals+1 );
330 // get all values sorted
331 std::multiset< double > values;
332 if ( elements.empty() )
334 SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator( GetType() );
335 while ( elemIt->more() )
336 values.insert( GetValue( elemIt->next()->GetID() ));
340 std::vector<int>::const_iterator id = elements.begin();
341 for ( ; id != elements.end(); ++id )
342 values.insert( GetValue( *id ));
347 funValues[0] = minmax[0];
348 funValues[nbIntervals] = minmax[1];
352 funValues[0] = *values.begin();
353 funValues[nbIntervals] = *values.rbegin();
355 // case nbIntervals == 1
356 if ( nbIntervals == 1 )
358 nbEvents[0] = values.size();
362 if (funValues.front() == funValues.back())
364 nbEvents.resize( 1 );
365 nbEvents[0] = values.size();
366 funValues[1] = funValues.back();
367 funValues.resize( 2 );
370 std::multiset< double >::iterator min = values.begin(), max;
371 for ( int i = 0; i < nbIntervals; ++i )
373 // find end value of i-th interval
374 double r = (i+1) / double(nbIntervals);
375 if (isLogarithmic && funValues.front() > 1e-07 && funValues.back() > 1e-07) {
376 double logmin = log10(funValues.front());
377 double lval = logmin + r * (log10(funValues.back()) - logmin);
378 funValues[i+1] = pow(10.0, lval);
381 funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
384 // count values in the i-th interval if there are any
385 if ( min != values.end() && *min <= funValues[i+1] )
387 // find the first value out of the interval
388 max = values.upper_bound( funValues[i+1] ); // max is greater than funValues[i+1], or end()
389 nbEvents[i] = std::distance( min, max );
393 // add values larger than minmax[1]
394 nbEvents.back() += std::distance( min, values.end() );
397 //=======================================================================
400 Description : Functor calculating volume of a 3D element
402 //================================================================================
404 double Volume::GetValue( long theElementId )
406 if ( theElementId && myMesh ) {
407 SMDS_VolumeTool aVolumeTool;
408 if ( aVolumeTool.Set( myMesh->FindElement( theElementId )))
409 return aVolumeTool.GetSize();
414 double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const
419 SMDSAbs_ElementType Volume::GetType() const
421 return SMDSAbs_Volume;
424 //=======================================================================
426 Class : MaxElementLength2D
427 Description : Functor calculating maximum length of 2D element
429 //================================================================================
431 double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P )
437 if( len == 3 ) { // triangles
438 double L1 = getDistance(P( 1 ),P( 2 ));
439 double L2 = getDistance(P( 2 ),P( 3 ));
440 double L3 = getDistance(P( 3 ),P( 1 ));
441 aVal = Max(L1,Max(L2,L3));
443 else if( len == 4 ) { // quadrangles
444 double L1 = getDistance(P( 1 ),P( 2 ));
445 double L2 = getDistance(P( 2 ),P( 3 ));
446 double L3 = getDistance(P( 3 ),P( 4 ));
447 double L4 = getDistance(P( 4 ),P( 1 ));
448 double D1 = getDistance(P( 1 ),P( 3 ));
449 double D2 = getDistance(P( 2 ),P( 4 ));
450 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
452 else if( len == 6 ) { // quadratic triangles
453 double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
454 double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
455 double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
456 aVal = Max(L1,Max(L2,L3));
458 else if( len == 8 || len == 9 ) { // quadratic quadrangles
459 double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
460 double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
461 double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
462 double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
463 double D1 = getDistance(P( 1 ),P( 5 ));
464 double D2 = getDistance(P( 3 ),P( 7 ));
465 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
467 // Diagonals are undefined for concave polygons
468 // else if ( P.getElementEntity() == SMDSEntity_Quad_Polygon && P.size() > 2 ) // quad polygon
471 // aVal = getDistance( P( 1 ), P( P.size() )) + getDistance( P( P.size() ), P( P.size()-1 ));
472 // for ( size_t i = 1; i < P.size()-1; i += 2 )
474 // double L = getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 ));
475 // aVal = Max( aVal, L );
478 // for ( int i = P.size()-5; i > 0; i -= 2 )
479 // for ( int j = i + 4; j < P.size() + i - 2; i += 2 )
481 // double D = getDistance( P( i ), P( j ));
482 // aVal = Max( aVal, D );
489 if( myPrecision >= 0 )
491 double prec = pow( 10., (double)myPrecision );
492 aVal = floor( aVal * prec + 0.5 ) / prec;
497 double MaxElementLength2D::GetValue( long theElementId )
500 return GetPoints( theElementId, P ) ? GetValue(P) : 0.0;
503 double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const
508 SMDSAbs_ElementType MaxElementLength2D::GetType() const
513 //=======================================================================
515 Class : MaxElementLength3D
516 Description : Functor calculating maximum length of 3D element
518 //================================================================================
520 double MaxElementLength3D::GetValue( long theElementId )
523 if( GetPoints( theElementId, P ) ) {
525 const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
526 SMDSAbs_EntityType aType = aElem->GetEntityType();
529 case SMDSEntity_Tetra: { // tetras
530 double L1 = getDistance(P( 1 ),P( 2 ));
531 double L2 = getDistance(P( 2 ),P( 3 ));
532 double L3 = getDistance(P( 3 ),P( 1 ));
533 double L4 = getDistance(P( 1 ),P( 4 ));
534 double L5 = getDistance(P( 2 ),P( 4 ));
535 double L6 = getDistance(P( 3 ),P( 4 ));
536 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
539 case SMDSEntity_Pyramid: { // pyramids
540 double L1 = getDistance(P( 1 ),P( 2 ));
541 double L2 = getDistance(P( 2 ),P( 3 ));
542 double L3 = getDistance(P( 3 ),P( 4 ));
543 double L4 = getDistance(P( 4 ),P( 1 ));
544 double L5 = getDistance(P( 1 ),P( 5 ));
545 double L6 = getDistance(P( 2 ),P( 5 ));
546 double L7 = getDistance(P( 3 ),P( 5 ));
547 double L8 = getDistance(P( 4 ),P( 5 ));
548 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
549 aVal = Max(aVal,Max(L7,L8));
552 case SMDSEntity_Penta: { // pentas
553 double L1 = getDistance(P( 1 ),P( 2 ));
554 double L2 = getDistance(P( 2 ),P( 3 ));
555 double L3 = getDistance(P( 3 ),P( 1 ));
556 double L4 = getDistance(P( 4 ),P( 5 ));
557 double L5 = getDistance(P( 5 ),P( 6 ));
558 double L6 = getDistance(P( 6 ),P( 4 ));
559 double L7 = getDistance(P( 1 ),P( 4 ));
560 double L8 = getDistance(P( 2 ),P( 5 ));
561 double L9 = getDistance(P( 3 ),P( 6 ));
562 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
563 aVal = Max(aVal,Max(Max(L7,L8),L9));
566 case SMDSEntity_Hexa: { // hexas
567 double L1 = getDistance(P( 1 ),P( 2 ));
568 double L2 = getDistance(P( 2 ),P( 3 ));
569 double L3 = getDistance(P( 3 ),P( 4 ));
570 double L4 = getDistance(P( 4 ),P( 1 ));
571 double L5 = getDistance(P( 5 ),P( 6 ));
572 double L6 = getDistance(P( 6 ),P( 7 ));
573 double L7 = getDistance(P( 7 ),P( 8 ));
574 double L8 = getDistance(P( 8 ),P( 5 ));
575 double L9 = getDistance(P( 1 ),P( 5 ));
576 double L10= getDistance(P( 2 ),P( 6 ));
577 double L11= getDistance(P( 3 ),P( 7 ));
578 double L12= getDistance(P( 4 ),P( 8 ));
579 double D1 = getDistance(P( 1 ),P( 7 ));
580 double D2 = getDistance(P( 2 ),P( 8 ));
581 double D3 = getDistance(P( 3 ),P( 5 ));
582 double D4 = getDistance(P( 4 ),P( 6 ));
583 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
584 aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
585 aVal = Max(aVal,Max(L11,L12));
586 aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
589 case SMDSEntity_Hexagonal_Prism: { // hexagonal prism
590 for ( int i1 = 1; i1 < 12; ++i1 )
591 for ( int i2 = i1+1; i1 <= 12; ++i1 )
592 aVal = Max( aVal, getDistance(P( i1 ),P( i2 )));
595 case SMDSEntity_Quad_Tetra: { // quadratic tetras
596 double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
597 double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
598 double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
599 double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
600 double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
601 double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
602 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
605 case SMDSEntity_Quad_Pyramid: { // quadratic pyramids
606 double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
607 double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
608 double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
609 double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
610 double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
611 double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
612 double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
613 double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
614 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
615 aVal = Max(aVal,Max(L7,L8));
618 case SMDSEntity_Quad_Penta: { // quadratic pentas
619 double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
620 double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
621 double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
622 double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
623 double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
624 double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
625 double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
626 double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
627 double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
628 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
629 aVal = Max(aVal,Max(Max(L7,L8),L9));
632 case SMDSEntity_Quad_Hexa:
633 case SMDSEntity_TriQuad_Hexa: { // quadratic hexas
634 double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
635 double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
636 double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
637 double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
638 double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
639 double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
640 double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
641 double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
642 double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
643 double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
644 double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
645 double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
646 double D1 = getDistance(P( 1 ),P( 7 ));
647 double D2 = getDistance(P( 2 ),P( 8 ));
648 double D3 = getDistance(P( 3 ),P( 5 ));
649 double D4 = getDistance(P( 4 ),P( 6 ));
650 aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
651 aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
652 aVal = Max(aVal,Max(L11,L12));
653 aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
656 case SMDSEntity_Quad_Polyhedra:
657 case SMDSEntity_Polyhedra: { // polys
658 // get the maximum distance between all pairs of nodes
659 for( int i = 1; i <= len; i++ ) {
660 for( int j = 1; j <= len; j++ ) {
661 if( j > i ) { // optimization of the loop
662 double D = getDistance( P(i), P(j) );
663 aVal = Max( aVal, D );
669 case SMDSEntity_Node:
671 case SMDSEntity_Edge:
672 case SMDSEntity_Quad_Edge:
673 case SMDSEntity_Triangle:
674 case SMDSEntity_Quad_Triangle:
675 case SMDSEntity_BiQuad_Triangle:
676 case SMDSEntity_Quadrangle:
677 case SMDSEntity_Quad_Quadrangle:
678 case SMDSEntity_BiQuad_Quadrangle:
679 case SMDSEntity_Polygon:
680 case SMDSEntity_Quad_Polygon:
681 case SMDSEntity_Ball:
682 case SMDSEntity_Last: return 0;
683 } // switch ( aType )
685 if( myPrecision >= 0 )
687 double prec = pow( 10., (double)myPrecision );
688 aVal = floor( aVal * prec + 0.5 ) / prec;
695 double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const
700 SMDSAbs_ElementType MaxElementLength3D::GetType() const
702 return SMDSAbs_Volume;
705 //=======================================================================
708 Description : Functor for calculation of minimum angle
710 //================================================================================
712 double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
719 aMin = getAngle(P( P.size() ), P( 1 ), P( 2 ));
720 aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 )));
722 for ( size_t i = 2; i < P.size(); i++ )
724 double A0 = getAngle( P( i-1 ), P( i ), P( i+1 ) );
728 return aMin * 180.0 / M_PI;
731 double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
733 //const double aBestAngle = PI / nbNodes;
734 const double aBestAngle = 180.0 - ( 360.0 / double(nbNodes) );
735 return ( fabs( aBestAngle - Value ));
738 SMDSAbs_ElementType MinimumAngle::GetType() const
744 //================================================================================
747 Description : Functor for calculating aspect ratio
749 //================================================================================
751 double AspectRatio::GetValue( long theId )
754 myCurrElement = myMesh->FindElement( theId );
755 if ( myCurrElement && myCurrElement->GetVtkType() == VTK_QUAD )
758 vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid();
759 if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() ))
760 aVal = Round( vtkMeshQuality::QuadAspectRatio( avtkCell ));
765 if ( GetPoints( myCurrElement, P ))
766 aVal = Round( GetValue( P ));
771 double AspectRatio::GetValue( const TSequenceOfXYZ& P )
773 // According to "Mesh quality control" by Nadir Bouhamau referring to
774 // Pascal Jean Frey and Paul-Louis George. Maillages, applications aux elements finis.
775 // Hermes Science publications, Paris 1999 ISBN 2-7462-0024-4
778 int nbNodes = P.size();
783 // Compute aspect ratio
785 if ( nbNodes == 3 ) {
786 // Compute lengths of the sides
787 std::vector< double > aLen (nbNodes);
788 for ( int i = 0; i < nbNodes - 1; i++ )
789 aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
790 aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
791 // Q = alfa * h * p / S, where
793 // alfa = sqrt( 3 ) / 6
794 // h - length of the longest edge
795 // p - half perimeter
796 // S - triangle surface
797 const double alfa = sqrt( 3. ) / 6.;
798 double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
799 double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
800 double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
801 if ( anArea <= theEps )
803 return alfa * maxLen * half_perimeter / anArea;
805 else if ( nbNodes == 6 ) { // quadratic triangles
806 // Compute lengths of the sides
807 std::vector< double > aLen (3);
808 aLen[0] = getDistance( P(1), P(3) );
809 aLen[1] = getDistance( P(3), P(5) );
810 aLen[2] = getDistance( P(5), P(1) );
811 // Q = alfa * h * p / S, where
813 // alfa = sqrt( 3 ) / 6
814 // h - length of the longest edge
815 // p - half perimeter
816 // S - triangle surface
817 const double alfa = sqrt( 3. ) / 6.;
818 double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
819 double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
820 double anArea = getArea( P(1), P(3), P(5) );
821 if ( anArea <= theEps )
823 return alfa * maxLen * half_perimeter / anArea;
825 else if( nbNodes == 4 ) { // quadrangle
826 // Compute lengths of the sides
827 std::vector< double > aLen (4);
828 aLen[0] = getDistance( P(1), P(2) );
829 aLen[1] = getDistance( P(2), P(3) );
830 aLen[2] = getDistance( P(3), P(4) );
831 aLen[3] = getDistance( P(4), P(1) );
832 // Compute lengths of the diagonals
833 std::vector< double > aDia (2);
834 aDia[0] = getDistance( P(1), P(3) );
835 aDia[1] = getDistance( P(2), P(4) );
836 // Compute areas of all triangles which can be built
837 // taking three nodes of the quadrangle
838 std::vector< double > anArea (4);
839 anArea[0] = getArea( P(1), P(2), P(3) );
840 anArea[1] = getArea( P(1), P(2), P(4) );
841 anArea[2] = getArea( P(1), P(3), P(4) );
842 anArea[3] = getArea( P(2), P(3), P(4) );
843 // Q = alpha * L * C1 / C2, where
845 // alpha = sqrt( 1/32 )
846 // L = max( L1, L2, L3, L4, D1, D2 )
847 // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
848 // C2 = min( S1, S2, S3, S4 )
849 // Li - lengths of the edges
850 // Di - lengths of the diagonals
851 // Si - areas of the triangles
852 const double alpha = sqrt( 1 / 32. );
853 double L = Max( aLen[ 0 ],
857 Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
858 double C1 = sqrt( ( aLen[0] * aLen[0] +
861 aLen[3] * aLen[3] ) / 4. );
862 double C2 = Min( anArea[ 0 ],
864 Min( anArea[ 2 ], anArea[ 3 ] ) ) );
867 return alpha * L * C1 / C2;
869 else if( nbNodes == 8 || nbNodes == 9 ) { // nbNodes==8 - quadratic quadrangle
870 // Compute lengths of the sides
871 std::vector< double > aLen (4);
872 aLen[0] = getDistance( P(1), P(3) );
873 aLen[1] = getDistance( P(3), P(5) );
874 aLen[2] = getDistance( P(5), P(7) );
875 aLen[3] = getDistance( P(7), P(1) );
876 // Compute lengths of the diagonals
877 std::vector< double > aDia (2);
878 aDia[0] = getDistance( P(1), P(5) );
879 aDia[1] = getDistance( P(3), P(7) );
880 // Compute areas of all triangles which can be built
881 // taking three nodes of the quadrangle
882 std::vector< double > anArea (4);
883 anArea[0] = getArea( P(1), P(3), P(5) );
884 anArea[1] = getArea( P(1), P(3), P(7) );
885 anArea[2] = getArea( P(1), P(5), P(7) );
886 anArea[3] = getArea( P(3), P(5), P(7) );
887 // Q = alpha * L * C1 / C2, where
889 // alpha = sqrt( 1/32 )
890 // L = max( L1, L2, L3, L4, D1, D2 )
891 // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
892 // C2 = min( S1, S2, S3, S4 )
893 // Li - lengths of the edges
894 // Di - lengths of the diagonals
895 // Si - areas of the triangles
896 const double alpha = sqrt( 1 / 32. );
897 double L = Max( aLen[ 0 ],
901 Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
902 double C1 = sqrt( ( aLen[0] * aLen[0] +
905 aLen[3] * aLen[3] ) / 4. );
906 double C2 = Min( anArea[ 0 ],
908 Min( anArea[ 2 ], anArea[ 3 ] ) ) );
911 return alpha * L * C1 / C2;
916 double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
918 // the aspect ratio is in the range [1.0,infinity]
919 // < 1.0 = very bad, zero area
922 return ( Value < 0.9 ) ? 1000 : Value / 1000.;
925 SMDSAbs_ElementType AspectRatio::GetType() const
931 //================================================================================
933 Class : AspectRatio3D
934 Description : Functor for calculating aspect ratio
936 //================================================================================
940 inline double getHalfPerimeter(double theTria[3]){
941 return (theTria[0] + theTria[1] + theTria[2])/2.0;
944 inline double getArea(double theHalfPerim, double theTria[3]){
945 return sqrt(theHalfPerim*
946 (theHalfPerim-theTria[0])*
947 (theHalfPerim-theTria[1])*
948 (theHalfPerim-theTria[2]));
951 inline double getVolume(double theLen[6]){
952 double a2 = theLen[0]*theLen[0];
953 double b2 = theLen[1]*theLen[1];
954 double c2 = theLen[2]*theLen[2];
955 double d2 = theLen[3]*theLen[3];
956 double e2 = theLen[4]*theLen[4];
957 double f2 = theLen[5]*theLen[5];
958 double P = 4.0*a2*b2*d2;
959 double Q = a2*(b2+d2-e2)-b2*(a2+d2-f2)-d2*(a2+b2-c2);
960 double R = (b2+d2-e2)*(a2+d2-f2)*(a2+d2-f2);
961 return sqrt(P-Q+R)/12.0;
964 inline double getVolume2(double theLen[6]){
965 double a2 = theLen[0]*theLen[0];
966 double b2 = theLen[1]*theLen[1];
967 double c2 = theLen[2]*theLen[2];
968 double d2 = theLen[3]*theLen[3];
969 double e2 = theLen[4]*theLen[4];
970 double f2 = theLen[5]*theLen[5];
972 double P = a2*e2*(b2+c2+d2+f2-a2-e2);
973 double Q = b2*f2*(a2+c2+d2+e2-b2-f2);
974 double R = c2*d2*(a2+b2+e2+f2-c2-d2);
975 double S = a2*b2*d2+b2*c2*e2+a2*c2*f2+d2*e2*f2;
977 return sqrt(P+Q+R-S)/12.0;
980 inline double getVolume(const TSequenceOfXYZ& P){
981 gp_Vec aVec1( P( 2 ) - P( 1 ) );
982 gp_Vec aVec2( P( 3 ) - P( 1 ) );
983 gp_Vec aVec3( P( 4 ) - P( 1 ) );
984 gp_Vec anAreaVec( aVec1 ^ aVec2 );
985 return fabs(aVec3 * anAreaVec) / 6.0;
988 inline double getMaxHeight(double theLen[6])
990 double aHeight = std::max(theLen[0],theLen[1]);
991 aHeight = std::max(aHeight,theLen[2]);
992 aHeight = std::max(aHeight,theLen[3]);
993 aHeight = std::max(aHeight,theLen[4]);
994 aHeight = std::max(aHeight,theLen[5]);
1000 double AspectRatio3D::GetValue( long theId )
1003 myCurrElement = myMesh->FindElement( theId );
1004 if ( myCurrElement && myCurrElement->GetVtkType() == VTK_TETRA )
1006 // Action from CoTech | ACTION 31.3:
1007 // EURIWARE BO: Homogenize the formulas used to calculate the Controls in SMESH to fit with
1008 // those of ParaView. The library used by ParaView for those calculations can be reused in SMESH.
1009 vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid();
1010 if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() ))
1011 aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell ));
1016 if ( GetPoints( myCurrElement, P ))
1017 aVal = Round( GetValue( P ));
1022 double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
1024 double aQuality = 0.0;
1025 if(myCurrElement->IsPoly()) return aQuality;
1027 int nbNodes = P.size();
1029 if(myCurrElement->IsQuadratic()) {
1030 if(nbNodes==10) nbNodes=4; // quadratic tetrahedron
1031 else if(nbNodes==13) nbNodes=5; // quadratic pyramid
1032 else if(nbNodes==15) nbNodes=6; // quadratic pentahedron
1033 else if(nbNodes==20) nbNodes=8; // quadratic hexahedron
1034 else if(nbNodes==27) nbNodes=8; // quadratic hexahedron
1035 else return aQuality;
1041 getDistance(P( 1 ),P( 2 )), // a
1042 getDistance(P( 2 ),P( 3 )), // b
1043 getDistance(P( 3 ),P( 1 )), // c
1044 getDistance(P( 2 ),P( 4 )), // d
1045 getDistance(P( 3 ),P( 4 )), // e
1046 getDistance(P( 1 ),P( 4 )) // f
1048 double aTria[4][3] = {
1049 {aLen[0],aLen[1],aLen[2]}, // abc
1050 {aLen[0],aLen[3],aLen[5]}, // adf
1051 {aLen[1],aLen[3],aLen[4]}, // bde
1052 {aLen[2],aLen[4],aLen[5]} // cef
1054 double aSumArea = 0.0;
1055 double aHalfPerimeter = getHalfPerimeter(aTria[0]);
1056 double anArea = getArea(aHalfPerimeter,aTria[0]);
1058 aHalfPerimeter = getHalfPerimeter(aTria[1]);
1059 anArea = getArea(aHalfPerimeter,aTria[1]);
1061 aHalfPerimeter = getHalfPerimeter(aTria[2]);
1062 anArea = getArea(aHalfPerimeter,aTria[2]);
1064 aHalfPerimeter = getHalfPerimeter(aTria[3]);
1065 anArea = getArea(aHalfPerimeter,aTria[3]);
1067 double aVolume = getVolume(P);
1068 //double aVolume = getVolume(aLen);
1069 double aHeight = getMaxHeight(aLen);
1070 static double aCoeff = sqrt(2.0)/12.0;
1071 if ( aVolume > DBL_MIN )
1072 aQuality = aCoeff*aHeight*aSumArea/aVolume;
1077 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
1078 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1081 gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
1082 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1085 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
1086 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1089 gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
1090 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1096 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
1097 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1100 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
1101 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1104 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
1105 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1108 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
1109 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1112 gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
1113 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1116 gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
1117 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1123 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
1124 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1127 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
1128 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1131 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
1132 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1135 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
1136 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1139 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
1140 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1143 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
1144 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1147 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
1148 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1151 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
1152 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1155 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
1156 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1159 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
1160 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1163 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
1164 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1167 gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
1168 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1171 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
1172 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1175 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
1176 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1179 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
1180 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1183 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
1184 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1187 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
1188 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1191 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
1192 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1195 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
1196 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1199 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
1200 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1203 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
1204 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1207 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1208 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1211 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
1212 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1215 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
1216 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1219 gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
1220 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1223 gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
1224 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1227 gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
1228 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1231 gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
1232 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1235 gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
1236 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1239 gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
1240 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1243 gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
1244 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1247 gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
1248 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1251 gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
1252 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
1258 gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )};
1259 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1262 gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )};
1263 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1266 gp_XYZ aXYZ[8] = {P( 3 ),P( 4 ),P( 6 ),P( 1 ),P( 9 ),P( 10 ),P( 12 ),P( 7 )};
1267 aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
1270 } // switch(nbNodes)
1272 if ( nbNodes > 4 ) {
1273 // avaluate aspect ratio of quadranle faces
1274 AspectRatio aspect2D;
1275 SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes );
1276 int nbFaces = SMDS_VolumeTool::NbFaces( type );
1277 TSequenceOfXYZ points(4);
1278 for ( int i = 0; i < nbFaces; ++i ) { // loop on faces of a volume
1279 if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 )
1281 const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true );
1282 for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face
1283 points( p + 1 ) = P( pInd[ p ] + 1 );
1284 aQuality = std::max( aQuality, aspect2D.GetValue( points ));
1290 double AspectRatio3D::GetBadRate( double Value, int /*nbNodes*/ ) const
1292 // the aspect ratio is in the range [1.0,infinity]
1295 return Value / 1000.;
1298 SMDSAbs_ElementType AspectRatio3D::GetType() const
1300 return SMDSAbs_Volume;
1304 //================================================================================
1307 Description : Functor for calculating warping
1309 //================================================================================
1311 double Warping::GetValue( const TSequenceOfXYZ& P )
1313 if ( P.size() != 4 )
1316 gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4.;
1318 double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
1319 double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
1320 double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
1321 double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
1323 double val = Max( Max( A1, A2 ), Max( A3, A4 ) );
1325 const double eps = 0.1; // val is in degrees
1327 return val < eps ? 0. : val;
1330 double Warping::ComputeA( const gp_XYZ& thePnt1,
1331 const gp_XYZ& thePnt2,
1332 const gp_XYZ& thePnt3,
1333 const gp_XYZ& theG ) const
1335 double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
1336 double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
1337 double L = Min( aLen1, aLen2 ) * 0.5;
1341 gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG;
1342 gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG;
1343 gp_XYZ N = GI.Crossed( GJ );
1345 if ( N.Modulus() < gp::Resolution() )
1350 double H = ( thePnt2 - theG ).Dot( N );
1351 return asin( fabs( H / L ) ) * 180. / M_PI;
1354 double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
1356 // the warp is in the range [0.0,PI/2]
1357 // 0.0 = good (no warp)
1358 // PI/2 = bad (face pliee)
1362 SMDSAbs_ElementType Warping::GetType() const
1364 return SMDSAbs_Face;
1368 //================================================================================
1371 Description : Functor for calculating taper
1373 //================================================================================
1375 double Taper::GetValue( const TSequenceOfXYZ& P )
1377 if ( P.size() != 4 )
1381 double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) );
1382 double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) );
1383 double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) );
1384 double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) );
1386 double JA = 0.25 * ( J1 + J2 + J3 + J4 );
1390 double T1 = fabs( ( J1 - JA ) / JA );
1391 double T2 = fabs( ( J2 - JA ) / JA );
1392 double T3 = fabs( ( J3 - JA ) / JA );
1393 double T4 = fabs( ( J4 - JA ) / JA );
1395 double val = Max( Max( T1, T2 ), Max( T3, T4 ) );
1397 const double eps = 0.01;
1399 return val < eps ? 0. : val;
1402 double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
1404 // the taper is in the range [0.0,1.0]
1405 // 0.0 = good (no taper)
1406 // 1.0 = bad (les cotes opposes sont allignes)
1410 SMDSAbs_ElementType Taper::GetType() const
1412 return SMDSAbs_Face;
1415 //================================================================================
1418 Description : Functor for calculating skew in degrees
1420 //================================================================================
1422 static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
1424 gp_XYZ p12 = ( p2 + p1 ) / 2.;
1425 gp_XYZ p23 = ( p3 + p2 ) / 2.;
1426 gp_XYZ p31 = ( p3 + p1 ) / 2.;
1428 gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
1430 return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 );
1433 double Skew::GetValue( const TSequenceOfXYZ& P )
1435 if ( P.size() != 3 && P.size() != 4 )
1439 const double PI2 = M_PI / 2.;
1440 if ( P.size() == 3 )
1442 double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
1443 double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
1444 double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
1446 return Max( A0, Max( A1, A2 ) ) * 180. / M_PI;
1450 gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2.;
1451 gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2.;
1452 gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2.;
1453 gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2.;
1455 gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
1456 double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
1457 ? 0. : fabs( PI2 - v1.Angle( v2 ) );
1459 double val = A * 180. / M_PI;
1461 const double eps = 0.1; // val is in degrees
1463 return val < eps ? 0. : val;
1467 double Skew::GetBadRate( double Value, int /*nbNodes*/ ) const
1469 // the skew is in the range [0.0,PI/2].
1475 SMDSAbs_ElementType Skew::GetType() const
1477 return SMDSAbs_Face;
1481 //================================================================================
1484 Description : Functor for calculating area
1486 //================================================================================
1488 double Area::GetValue( const TSequenceOfXYZ& P )
1493 gp_Vec aVec1( P(2) - P(1) );
1494 gp_Vec aVec2( P(3) - P(1) );
1495 gp_Vec SumVec = aVec1 ^ aVec2;
1497 for (size_t i=4; i<=P.size(); i++)
1499 gp_Vec aVec1( P(i-1) - P(1) );
1500 gp_Vec aVec2( P(i ) - P(1) );
1501 gp_Vec tmp = aVec1 ^ aVec2;
1504 val = SumVec.Magnitude() * 0.5;
1509 double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
1511 // meaningless as it is not a quality control functor
1515 SMDSAbs_ElementType Area::GetType() const
1517 return SMDSAbs_Face;
1520 //================================================================================
1523 Description : Functor for calculating length of edge
1525 //================================================================================
1527 double Length::GetValue( const TSequenceOfXYZ& P )
1529 switch ( P.size() ) {
1530 case 2: return getDistance( P( 1 ), P( 2 ) );
1531 case 3: return getDistance( P( 1 ), P( 2 ) ) + getDistance( P( 2 ), P( 3 ) );
1536 double Length::GetBadRate( double Value, int /*nbNodes*/ ) const
1538 // meaningless as it is not quality control functor
1542 SMDSAbs_ElementType Length::GetType() const
1544 return SMDSAbs_Edge;
1547 //================================================================================
1550 Description : Functor for calculating minimal length of edge
1552 //================================================================================
1554 double Length2D::GetValue( long theElementId )
1558 if ( GetPoints( theElementId, P ))
1562 SMDSAbs_EntityType aType = P.getElementEntity();
1565 case SMDSEntity_Edge:
1567 aVal = getDistance( P( 1 ), P( 2 ) );
1569 case SMDSEntity_Quad_Edge:
1570 if (len == 3) // quadratic edge
1571 aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
1573 case SMDSEntity_Triangle:
1574 if (len == 3){ // triangles
1575 double L1 = getDistance(P( 1 ),P( 2 ));
1576 double L2 = getDistance(P( 2 ),P( 3 ));
1577 double L3 = getDistance(P( 3 ),P( 1 ));
1578 aVal = Min(L1,Min(L2,L3));
1581 case SMDSEntity_Quadrangle:
1582 if (len == 4){ // quadrangles
1583 double L1 = getDistance(P( 1 ),P( 2 ));
1584 double L2 = getDistance(P( 2 ),P( 3 ));
1585 double L3 = getDistance(P( 3 ),P( 4 ));
1586 double L4 = getDistance(P( 4 ),P( 1 ));
1587 aVal = Min(Min(L1,L2),Min(L3,L4));
1590 case SMDSEntity_Quad_Triangle:
1591 case SMDSEntity_BiQuad_Triangle:
1592 if (len >= 6){ // quadratic triangles
1593 double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1594 double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1595 double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
1596 aVal = Min(L1,Min(L2,L3));
1599 case SMDSEntity_Quad_Quadrangle:
1600 case SMDSEntity_BiQuad_Quadrangle:
1601 if (len >= 8){ // quadratic quadrangles
1602 double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
1603 double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
1604 double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
1605 double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
1606 aVal = Min(Min(L1,L2),Min(L3,L4));
1609 case SMDSEntity_Tetra:
1610 if (len == 4){ // tetrahedra
1611 double L1 = getDistance(P( 1 ),P( 2 ));
1612 double L2 = getDistance(P( 2 ),P( 3 ));
1613 double L3 = getDistance(P( 3 ),P( 1 ));
1614 double L4 = getDistance(P( 1 ),P( 4 ));
1615 double L5 = getDistance(P( 2 ),P( 4 ));
1616 double L6 = getDistance(P( 3 ),P( 4 ));
1617 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1620 case SMDSEntity_Pyramid:
1621 if (len == 5){ // piramids
1622 double L1 = getDistance(P( 1 ),P( 2 ));
1623 double L2 = getDistance(P( 2 ),P( 3 ));
1624 double L3 = getDistance(P( 3 ),P( 4 ));
1625 double L4 = getDistance(P( 4 ),P( 1 ));
1626 double L5 = getDistance(P( 1 ),P( 5 ));
1627 double L6 = getDistance(P( 2 ),P( 5 ));
1628 double L7 = getDistance(P( 3 ),P( 5 ));
1629 double L8 = getDistance(P( 4 ),P( 5 ));
1631 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1632 aVal = Min(aVal,Min(L7,L8));
1635 case SMDSEntity_Penta:
1636 if (len == 6) { // pentaidres
1637 double L1 = getDistance(P( 1 ),P( 2 ));
1638 double L2 = getDistance(P( 2 ),P( 3 ));
1639 double L3 = getDistance(P( 3 ),P( 1 ));
1640 double L4 = getDistance(P( 4 ),P( 5 ));
1641 double L5 = getDistance(P( 5 ),P( 6 ));
1642 double L6 = getDistance(P( 6 ),P( 4 ));
1643 double L7 = getDistance(P( 1 ),P( 4 ));
1644 double L8 = getDistance(P( 2 ),P( 5 ));
1645 double L9 = getDistance(P( 3 ),P( 6 ));
1647 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1648 aVal = Min(aVal,Min(Min(L7,L8),L9));
1651 case SMDSEntity_Hexa:
1652 if (len == 8){ // hexahedron
1653 double L1 = getDistance(P( 1 ),P( 2 ));
1654 double L2 = getDistance(P( 2 ),P( 3 ));
1655 double L3 = getDistance(P( 3 ),P( 4 ));
1656 double L4 = getDistance(P( 4 ),P( 1 ));
1657 double L5 = getDistance(P( 5 ),P( 6 ));
1658 double L6 = getDistance(P( 6 ),P( 7 ));
1659 double L7 = getDistance(P( 7 ),P( 8 ));
1660 double L8 = getDistance(P( 8 ),P( 5 ));
1661 double L9 = getDistance(P( 1 ),P( 5 ));
1662 double L10= getDistance(P( 2 ),P( 6 ));
1663 double L11= getDistance(P( 3 ),P( 7 ));
1664 double L12= getDistance(P( 4 ),P( 8 ));
1666 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1667 aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
1668 aVal = Min(aVal,Min(L11,L12));
1671 case SMDSEntity_Quad_Tetra:
1672 if (len == 10){ // quadratic tetraidrs
1673 double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
1674 double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
1675 double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
1676 double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1677 double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
1678 double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
1679 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1682 case SMDSEntity_Quad_Pyramid:
1683 if (len == 13){ // quadratic piramids
1684 double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
1685 double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
1686 double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
1687 double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1688 double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1689 double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
1690 double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
1691 double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
1692 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1693 aVal = Min(aVal,Min(L7,L8));
1696 case SMDSEntity_Quad_Penta:
1697 if (len == 15){ // quadratic pentaidres
1698 double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
1699 double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
1700 double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
1701 double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
1702 double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
1703 double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
1704 double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
1705 double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
1706 double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
1707 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1708 aVal = Min(aVal,Min(Min(L7,L8),L9));
1711 case SMDSEntity_Quad_Hexa:
1712 case SMDSEntity_TriQuad_Hexa:
1713 if (len >= 20) { // quadratic hexaider
1714 double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
1715 double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
1716 double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
1717 double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
1718 double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
1719 double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
1720 double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
1721 double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
1722 double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
1723 double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
1724 double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
1725 double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
1726 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1727 aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
1728 aVal = Min(aVal,Min(L11,L12));
1731 case SMDSEntity_Polygon:
1733 aVal = getDistance( P(1), P( P.size() ));
1734 for ( size_t i = 1; i < P.size(); ++i )
1735 aVal = Min( aVal, getDistance( P( i ), P( i+1 )));
1738 case SMDSEntity_Quad_Polygon:
1740 aVal = getDistance( P(1), P( P.size() )) + getDistance( P(P.size()), P( P.size()-1 ));
1741 for ( size_t i = 1; i < P.size()-1; i += 2 )
1742 aVal = Min( aVal, getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 )));
1745 case SMDSEntity_Hexagonal_Prism:
1746 if (len == 12) { // hexagonal prism
1747 double L1 = getDistance(P( 1 ),P( 2 ));
1748 double L2 = getDistance(P( 2 ),P( 3 ));
1749 double L3 = getDistance(P( 3 ),P( 4 ));
1750 double L4 = getDistance(P( 4 ),P( 5 ));
1751 double L5 = getDistance(P( 5 ),P( 6 ));
1752 double L6 = getDistance(P( 6 ),P( 1 ));
1754 double L7 = getDistance(P( 7 ), P( 8 ));
1755 double L8 = getDistance(P( 8 ), P( 9 ));
1756 double L9 = getDistance(P( 9 ), P( 10 ));
1757 double L10= getDistance(P( 10 ),P( 11 ));
1758 double L11= getDistance(P( 11 ),P( 12 ));
1759 double L12= getDistance(P( 12 ),P( 7 ));
1761 double L13 = getDistance(P( 1 ),P( 7 ));
1762 double L14 = getDistance(P( 2 ),P( 8 ));
1763 double L15 = getDistance(P( 3 ),P( 9 ));
1764 double L16 = getDistance(P( 4 ),P( 10 ));
1765 double L17 = getDistance(P( 5 ),P( 11 ));
1766 double L18 = getDistance(P( 6 ),P( 12 ));
1767 aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
1768 aVal = Min(aVal, Min(Min(Min(L7,L8),Min(L9,L10)),Min(L11,L12)));
1769 aVal = Min(aVal, Min(Min(Min(L13,L14),Min(L15,L16)),Min(L17,L18)));
1772 case SMDSEntity_Polyhedra:
1784 if ( myPrecision >= 0 )
1786 double prec = pow( 10., (double)( myPrecision ) );
1787 aVal = floor( aVal * prec + 0.5 ) / prec;
1796 double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
1798 // meaningless as it is not a quality control functor
1802 SMDSAbs_ElementType Length2D::GetType() const
1804 return SMDSAbs_Face;
1807 Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
1810 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
1811 if(thePntId1 > thePntId2){
1812 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
1816 bool Length2D::Value::operator<(const Length2D::Value& x) const
1818 if(myPntId[0] < x.myPntId[0]) return true;
1819 if(myPntId[0] == x.myPntId[0])
1820 if(myPntId[1] < x.myPntId[1]) return true;
1824 void Length2D::GetValues(TValues& theValues)
1827 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
1828 for(; anIter->more(); ){
1829 const SMDS_MeshFace* anElem = anIter->next();
1831 if(anElem->IsQuadratic()) {
1832 const SMDS_VtkFace* F =
1833 dynamic_cast<const SMDS_VtkFace*>(anElem);
1834 // use special nodes iterator
1835 SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
1836 long aNodeId[4] = { 0,0,0,0 };
1840 const SMDS_MeshElement* aNode;
1842 aNode = anIter->next();
1843 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1844 P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1845 aNodeId[0] = aNodeId[1] = aNode->GetID();
1848 for(; anIter->more(); ){
1849 const SMDS_MeshNode* N1 = static_cast<const SMDS_MeshNode*> (anIter->next());
1850 P[2] = gp_Pnt(N1->X(),N1->Y(),N1->Z());
1851 aNodeId[2] = N1->GetID();
1852 aLength = P[1].Distance(P[2]);
1853 if(!anIter->more()) break;
1854 const SMDS_MeshNode* N2 = static_cast<const SMDS_MeshNode*> (anIter->next());
1855 P[3] = gp_Pnt(N2->X(),N2->Y(),N2->Z());
1856 aNodeId[3] = N2->GetID();
1857 aLength += P[2].Distance(P[3]);
1858 Value aValue1(aLength,aNodeId[1],aNodeId[2]);
1859 Value aValue2(aLength,aNodeId[2],aNodeId[3]);
1861 aNodeId[1] = aNodeId[3];
1862 theValues.insert(aValue1);
1863 theValues.insert(aValue2);
1865 aLength += P[2].Distance(P[0]);
1866 Value aValue1(aLength,aNodeId[1],aNodeId[2]);
1867 Value aValue2(aLength,aNodeId[2],aNodeId[0]);
1868 theValues.insert(aValue1);
1869 theValues.insert(aValue2);
1872 SMDS_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
1873 long aNodeId[2] = {0,0};
1877 const SMDS_MeshElement* aNode;
1878 if(aNodesIter->more()){
1879 aNode = aNodesIter->next();
1880 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1881 P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1882 aNodeId[0] = aNodeId[1] = aNode->GetID();
1885 for(; aNodesIter->more(); ){
1886 aNode = aNodesIter->next();
1887 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
1888 long anId = aNode->GetID();
1890 P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
1892 aLength = P[1].Distance(P[2]);
1894 Value aValue(aLength,aNodeId[1],anId);
1897 theValues.insert(aValue);
1900 aLength = P[0].Distance(P[1]);
1902 Value aValue(aLength,aNodeId[0],aNodeId[1]);
1903 theValues.insert(aValue);
1908 //================================================================================
1910 Class : MultiConnection
1911 Description : Functor for calculating number of faces conneted to the edge
1913 //================================================================================
1915 double MultiConnection::GetValue( const TSequenceOfXYZ& P )
1919 double MultiConnection::GetValue( long theId )
1921 return getNbMultiConnection( myMesh, theId );
1924 double MultiConnection::GetBadRate( double Value, int /*nbNodes*/ ) const
1926 // meaningless as it is not quality control functor
1930 SMDSAbs_ElementType MultiConnection::GetType() const
1932 return SMDSAbs_Edge;
1935 //================================================================================
1937 Class : MultiConnection2D
1938 Description : Functor for calculating number of faces conneted to the edge
1940 //================================================================================
1942 double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
1947 double MultiConnection2D::GetValue( long theElementId )
1951 const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId);
1952 SMDSAbs_ElementType aType = aFaceElem->GetType();
1957 int i = 0, len = aFaceElem->NbNodes();
1958 SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
1961 const SMDS_MeshNode *aNode, *aNode0 = 0;
1962 TColStd_MapOfInteger aMap, aMapPrev;
1964 for (i = 0; i <= len; i++) {
1969 if (anIter->more()) {
1970 aNode = (SMDS_MeshNode*)anIter->next();
1978 if (i == 0) aNode0 = aNode;
1980 SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
1981 while (anElemIter->more()) {
1982 const SMDS_MeshElement* anElem = anElemIter->next();
1983 if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
1984 int anId = anElem->GetID();
1987 if (aMapPrev.Contains(anId)) {
1992 aResult = Max(aResult, aNb);
2003 double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
2005 // meaningless as it is not quality control functor
2009 SMDSAbs_ElementType MultiConnection2D::GetType() const
2011 return SMDSAbs_Face;
2014 MultiConnection2D::Value::Value(long thePntId1, long thePntId2)
2016 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
2017 if(thePntId1 > thePntId2){
2018 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
2022 bool MultiConnection2D::Value::operator<(const MultiConnection2D::Value& x) const
2024 if(myPntId[0] < x.myPntId[0]) return true;
2025 if(myPntId[0] == x.myPntId[0])
2026 if(myPntId[1] < x.myPntId[1]) return true;
2030 void MultiConnection2D::GetValues(MValues& theValues)
2032 if ( !myMesh ) return;
2033 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
2034 for(; anIter->more(); ){
2035 const SMDS_MeshFace* anElem = anIter->next();
2036 SMDS_ElemIteratorPtr aNodesIter;
2037 if ( anElem->IsQuadratic() )
2038 aNodesIter = dynamic_cast<const SMDS_VtkFace*>
2039 (anElem)->interlacedNodesElemIterator();
2041 aNodesIter = anElem->nodesIterator();
2042 long aNodeId[3] = {0,0,0};
2044 //int aNbConnects=0;
2045 const SMDS_MeshNode* aNode0;
2046 const SMDS_MeshNode* aNode1;
2047 const SMDS_MeshNode* aNode2;
2048 if(aNodesIter->more()){
2049 aNode0 = (SMDS_MeshNode*) aNodesIter->next();
2051 const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
2052 aNodeId[0] = aNodeId[1] = aNodes->GetID();
2054 for(; aNodesIter->more(); ) {
2055 aNode2 = (SMDS_MeshNode*) aNodesIter->next();
2056 long anId = aNode2->GetID();
2059 Value aValue(aNodeId[1],aNodeId[2]);
2060 MValues::iterator aItr = theValues.find(aValue);
2061 if (aItr != theValues.end()){
2066 theValues[aValue] = 1;
2069 //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
2070 aNodeId[1] = aNodeId[2];
2073 Value aValue(aNodeId[0],aNodeId[2]);
2074 MValues::iterator aItr = theValues.find(aValue);
2075 if (aItr != theValues.end()) {
2080 theValues[aValue] = 1;
2083 //cout << "NodeIds: "<<aNodeId[0]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
2088 //================================================================================
2090 Class : BallDiameter
2091 Description : Functor returning diameter of a ball element
2093 //================================================================================
2095 double BallDiameter::GetValue( long theId )
2097 double diameter = 0;
2099 if ( const SMDS_BallElement* ball =
2100 dynamic_cast<const SMDS_BallElement*>( myMesh->FindElement( theId )))
2102 diameter = ball->GetDiameter();
2107 double BallDiameter::GetBadRate( double Value, int /*nbNodes*/ ) const
2109 // meaningless as it is not a quality control functor
2113 SMDSAbs_ElementType BallDiameter::GetType() const
2115 return SMDSAbs_Ball;
2123 //================================================================================
2125 Class : BadOrientedVolume
2126 Description : Predicate bad oriented volumes
2128 //================================================================================
2130 BadOrientedVolume::BadOrientedVolume()
2135 void BadOrientedVolume::SetMesh( const SMDS_Mesh* theMesh )
2140 bool BadOrientedVolume::IsSatisfy( long theId )
2145 SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
2146 return !vTool.IsForward();
2149 SMDSAbs_ElementType BadOrientedVolume::GetType() const
2151 return SMDSAbs_Volume;
2155 Class : BareBorderVolume
2158 bool BareBorderVolume::IsSatisfy(long theElementId )
2160 SMDS_VolumeTool myTool;
2161 if ( myTool.Set( myMesh->FindElement(theElementId)))
2163 for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
2164 if ( myTool.IsFreeFace( iF ))
2166 const SMDS_MeshNode** n = myTool.GetFaceNodes(iF);
2167 std::vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF));
2168 if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*Nomedium=*/false))
2175 //================================================================================
2177 Class : BareBorderFace
2179 //================================================================================
2181 bool BareBorderFace::IsSatisfy(long theElementId )
2184 if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
2186 if ( face->GetType() == SMDSAbs_Face )
2188 int nbN = face->NbCornerNodes();
2189 for ( int i = 0; i < nbN && !ok; ++i )
2191 // check if a link is shared by another face
2192 const SMDS_MeshNode* n1 = face->GetNode( i );
2193 const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
2194 SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
2195 bool isShared = false;
2196 while ( !isShared && fIt->more() )
2198 const SMDS_MeshElement* f = fIt->next();
2199 isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
2203 const int iQuad = face->IsQuadratic();
2204 myLinkNodes.resize( 2 + iQuad);
2205 myLinkNodes[0] = n1;
2206 myLinkNodes[1] = n2;
2208 myLinkNodes[2] = face->GetNode( i+nbN );
2209 ok = !myMesh->FindElement( myLinkNodes, SMDSAbs_Edge, /*noMedium=*/false);
2217 //================================================================================
2219 Class : OverConstrainedVolume
2221 //================================================================================
2223 bool OverConstrainedVolume::IsSatisfy(long theElementId )
2225 // An element is over-constrained if it has N-1 free borders where
2226 // N is the number of edges/faces for a 2D/3D element.
2227 SMDS_VolumeTool myTool;
2228 if ( myTool.Set( myMesh->FindElement(theElementId)))
2230 int nbSharedFaces = 0;
2231 for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
2232 if ( !myTool.IsFreeFace( iF ) && ++nbSharedFaces > 1 )
2234 return ( nbSharedFaces == 1 );
2239 //================================================================================
2241 Class : OverConstrainedFace
2243 //================================================================================
2245 bool OverConstrainedFace::IsSatisfy(long theElementId )
2247 // An element is over-constrained if it has N-1 free borders where
2248 // N is the number of edges/faces for a 2D/3D element.
2249 if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
2250 if ( face->GetType() == SMDSAbs_Face )
2252 int nbSharedBorders = 0;
2253 int nbN = face->NbCornerNodes();
2254 for ( int i = 0; i < nbN; ++i )
2256 // check if a link is shared by another face
2257 const SMDS_MeshNode* n1 = face->GetNode( i );
2258 const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
2259 SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
2260 bool isShared = false;
2261 while ( !isShared && fIt->more() )
2263 const SMDS_MeshElement* f = fIt->next();
2264 isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
2266 if ( isShared && ++nbSharedBorders > 1 )
2269 return ( nbSharedBorders == 1 );
2274 //================================================================================
2276 Class : CoincidentNodes
2277 Description : Predicate of Coincident nodes
2279 //================================================================================
2281 CoincidentNodes::CoincidentNodes()
2286 bool CoincidentNodes::IsSatisfy( long theElementId )
2288 return myCoincidentIDs.Contains( theElementId );
2291 SMDSAbs_ElementType CoincidentNodes::GetType() const
2293 return SMDSAbs_Node;
2296 void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
2298 myMeshModifTracer.SetMesh( theMesh );
2299 if ( myMeshModifTracer.IsMeshModified() )
2301 TIDSortedNodeSet nodesToCheck;
2302 SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true);
2303 while ( nIt->more() )
2304 nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
2306 std::list< std::list< const SMDS_MeshNode*> > nodeGroups;
2307 SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler );
2309 myCoincidentIDs.Clear();
2310 std::list< std::list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin();
2311 for ( ; groupIt != nodeGroups.end(); ++groupIt )
2313 std::list< const SMDS_MeshNode*>& coincNodes = *groupIt;
2314 std::list< const SMDS_MeshNode*>::iterator n = coincNodes.begin();
2315 for ( ; n != coincNodes.end(); ++n )
2316 myCoincidentIDs.Add( (*n)->GetID() );
2321 //================================================================================
2323 Class : CoincidentElements
2324 Description : Predicate of Coincident Elements
2325 Note : This class is suitable only for visualization of Coincident Elements
2327 //================================================================================
2329 CoincidentElements::CoincidentElements()
2334 void CoincidentElements::SetMesh( const SMDS_Mesh* theMesh )
2339 bool CoincidentElements::IsSatisfy( long theElementId )
2341 if ( !myMesh ) return false;
2343 if ( const SMDS_MeshElement* e = myMesh->FindElement( theElementId ))
2345 if ( e->GetType() != GetType() ) return false;
2346 std::set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() );
2347 const int nbNodes = e->NbNodes();
2348 SMDS_ElemIteratorPtr invIt = (*elemNodes.begin())->GetInverseElementIterator( GetType() );
2349 while ( invIt->more() )
2351 const SMDS_MeshElement* e2 = invIt->next();
2352 if ( e2 == e || e2->NbNodes() != nbNodes ) continue;
2354 bool sameNodes = true;
2355 for ( size_t i = 0; i < elemNodes.size() && sameNodes; ++i )
2356 sameNodes = ( elemNodes.count( e2->GetNode( i )));
2364 SMDSAbs_ElementType CoincidentElements1D::GetType() const
2366 return SMDSAbs_Edge;
2368 SMDSAbs_ElementType CoincidentElements2D::GetType() const
2370 return SMDSAbs_Face;
2372 SMDSAbs_ElementType CoincidentElements3D::GetType() const
2374 return SMDSAbs_Volume;
2378 //================================================================================
2381 Description : Predicate for free borders
2383 //================================================================================
2385 FreeBorders::FreeBorders()
2390 void FreeBorders::SetMesh( const SMDS_Mesh* theMesh )
2395 bool FreeBorders::IsSatisfy( long theId )
2397 return getNbMultiConnection( myMesh, theId ) == 1;
2400 SMDSAbs_ElementType FreeBorders::GetType() const
2402 return SMDSAbs_Edge;
2406 //================================================================================
2409 Description : Predicate for free Edges
2411 //================================================================================
2413 FreeEdges::FreeEdges()
2418 void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
2423 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId )
2425 TColStd_MapOfInteger aMap;
2426 for ( int i = 0; i < 2; i++ )
2428 SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator(SMDSAbs_Face);
2429 while( anElemIter->more() )
2431 if ( const SMDS_MeshElement* anElem = anElemIter->next())
2433 const int anId = anElem->GetID();
2434 if ( anId != theFaceId && !aMap.Add( anId ))
2442 bool FreeEdges::IsSatisfy( long theId )
2447 const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2448 if ( aFace == 0 || aFace->GetType() != SMDSAbs_Face || aFace->NbNodes() < 3 )
2451 SMDS_NodeIteratorPtr anIter = aFace->interlacedNodesIterator();
2455 int i = 0, nbNodes = aFace->NbNodes();
2456 std::vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
2457 while( anIter->more() )
2458 if ( ! ( aNodes[ i++ ] = anIter->next() ))
2460 aNodes[ nbNodes ] = aNodes[ 0 ];
2462 for ( i = 0; i < nbNodes; i++ )
2463 if ( IsFreeEdge( &aNodes[ i ], theId ) )
2469 SMDSAbs_ElementType FreeEdges::GetType() const
2471 return SMDSAbs_Face;
2474 FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
2477 myPntId[0] = thePntId1; myPntId[1] = thePntId2;
2478 if(thePntId1 > thePntId2){
2479 myPntId[1] = thePntId1; myPntId[0] = thePntId2;
2483 bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
2484 if(myPntId[0] < x.myPntId[0]) return true;
2485 if(myPntId[0] == x.myPntId[0])
2486 if(myPntId[1] < x.myPntId[1]) return true;
2490 inline void UpdateBorders(const FreeEdges::Border& theBorder,
2491 FreeEdges::TBorders& theRegistry,
2492 FreeEdges::TBorders& theContainer)
2494 if(theRegistry.find(theBorder) == theRegistry.end()){
2495 theRegistry.insert(theBorder);
2496 theContainer.insert(theBorder);
2498 theContainer.erase(theBorder);
2502 void FreeEdges::GetBoreders(TBorders& theBorders)
2505 SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
2506 for(; anIter->more(); ){
2507 const SMDS_MeshFace* anElem = anIter->next();
2508 long anElemId = anElem->GetID();
2509 SMDS_ElemIteratorPtr aNodesIter;
2510 if ( anElem->IsQuadratic() )
2511 aNodesIter = static_cast<const SMDS_VtkFace*>(anElem)->
2512 interlacedNodesElemIterator();
2514 aNodesIter = anElem->nodesIterator();
2515 long aNodeId[2] = {0,0};
2516 const SMDS_MeshElement* aNode;
2517 if(aNodesIter->more()){
2518 aNode = aNodesIter->next();
2519 aNodeId[0] = aNodeId[1] = aNode->GetID();
2521 for(; aNodesIter->more(); ){
2522 aNode = aNodesIter->next();
2523 long anId = aNode->GetID();
2524 Border aBorder(anElemId,aNodeId[1],anId);
2526 UpdateBorders(aBorder,aRegistry,theBorders);
2528 Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
2529 UpdateBorders(aBorder,aRegistry,theBorders);
2533 //================================================================================
2536 Description : Predicate for free nodes
2538 //================================================================================
2540 FreeNodes::FreeNodes()
2545 void FreeNodes::SetMesh( const SMDS_Mesh* theMesh )
2550 bool FreeNodes::IsSatisfy( long theNodeId )
2552 const SMDS_MeshNode* aNode = myMesh->FindNode( theNodeId );
2556 return (aNode->NbInverseElements() < 1);
2559 SMDSAbs_ElementType FreeNodes::GetType() const
2561 return SMDSAbs_Node;
2565 //================================================================================
2568 Description : Predicate for free faces
2570 //================================================================================
2572 FreeFaces::FreeFaces()
2577 void FreeFaces::SetMesh( const SMDS_Mesh* theMesh )
2582 bool FreeFaces::IsSatisfy( long theId )
2584 if (!myMesh) return false;
2585 // check that faces nodes refers to less than two common volumes
2586 const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
2587 if ( !aFace || aFace->GetType() != SMDSAbs_Face )
2590 int nbNode = aFace->NbNodes();
2592 // collect volumes to check that number of volumes with count equal nbNode not less than 2
2593 typedef std::map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
2594 typedef std::map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
2595 TMapOfVolume mapOfVol;
2597 SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator();
2598 while ( nodeItr->more() )
2600 const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
2601 if ( !aNode ) continue;
2602 SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume);
2603 while ( volItr->more() )
2605 SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
2606 TItrMapOfVolume itr = mapOfVol.insert( std::make_pair( aVol, 0 )).first;
2611 TItrMapOfVolume volItr = mapOfVol.begin();
2612 TItrMapOfVolume volEnd = mapOfVol.end();
2613 for ( ; volItr != volEnd; ++volItr )
2614 if ( (*volItr).second >= nbNode )
2616 // face is not free if number of volumes constructed on thier nodes more than one
2620 SMDSAbs_ElementType FreeFaces::GetType() const
2622 return SMDSAbs_Face;
2625 //================================================================================
2627 Class : LinearOrQuadratic
2628 Description : Predicate to verify whether a mesh element is linear
2630 //================================================================================
2632 LinearOrQuadratic::LinearOrQuadratic()
2637 void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh )
2642 bool LinearOrQuadratic::IsSatisfy( long theId )
2644 if (!myMesh) return false;
2645 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2646 if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) )
2648 return (!anElem->IsQuadratic());
2651 void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType )
2656 SMDSAbs_ElementType LinearOrQuadratic::GetType() const
2661 //================================================================================
2664 Description : Functor for check color of group to whic mesh element belongs to
2666 //================================================================================
2668 GroupColor::GroupColor()
2672 bool GroupColor::IsSatisfy( long theId )
2674 return myIDs.count( theId );
2677 void GroupColor::SetType( SMDSAbs_ElementType theType )
2682 SMDSAbs_ElementType GroupColor::GetType() const
2687 static bool isEqual( const Quantity_Color& theColor1,
2688 const Quantity_Color& theColor2 )
2690 // tolerance to compare colors
2691 const double tol = 5*1e-3;
2692 return ( fabs( theColor1.Red() - theColor2.Red() ) < tol &&
2693 fabs( theColor1.Green() - theColor2.Green() ) < tol &&
2694 fabs( theColor1.Blue() - theColor2.Blue() ) < tol );
2697 void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
2701 const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
2705 int nbGrp = aMesh->GetNbGroups();
2709 // iterates on groups and find necessary elements ids
2710 const std::set<SMESHDS_GroupBase*>& aGroups = aMesh->GetGroups();
2711 std::set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
2712 for (; GrIt != aGroups.end(); GrIt++)
2714 SMESHDS_GroupBase* aGrp = (*GrIt);
2717 // check type and color of group
2718 if ( !isEqual( myColor, aGrp->GetColor() ))
2721 // IPAL52867 (prevent infinite recursion via GroupOnFilter)
2722 if ( SMESHDS_GroupOnFilter * gof = dynamic_cast< SMESHDS_GroupOnFilter* >( aGrp ))
2723 if ( gof->GetPredicate().get() == this )
2726 SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType();
2727 if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) {
2728 // add elements IDS into control
2729 int aSize = aGrp->Extent();
2730 for (int i = 0; i < aSize; i++)
2731 myIDs.insert( aGrp->GetID(i+1) );
2736 void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
2738 Kernel_Utils::Localizer loc;
2739 TCollection_AsciiString aStr = theStr;
2740 aStr.RemoveAll( ' ' );
2741 aStr.RemoveAll( '\t' );
2742 for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) )
2743 aStr.Remove( aPos, 2 );
2744 Standard_Real clr[3];
2745 clr[0] = clr[1] = clr[2] = 0.;
2746 for ( int i = 0; i < 3; i++ ) {
2747 TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 );
2748 if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() )
2749 clr[i] = tmpStr.RealValue();
2751 myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB );
2754 //=======================================================================
2755 // name : GetRangeStr
2756 // Purpose : Get range as a string.
2757 // Example: "1,2,3,50-60,63,67,70-"
2758 //=======================================================================
2760 void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
2763 theResStr += TCollection_AsciiString( myColor.Red() );
2764 theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() );
2765 theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
2768 //================================================================================
2770 Class : ElemGeomType
2771 Description : Predicate to check element geometry type
2773 //================================================================================
2775 ElemGeomType::ElemGeomType()
2778 myType = SMDSAbs_All;
2779 myGeomType = SMDSGeom_TRIANGLE;
2782 void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh )
2787 bool ElemGeomType::IsSatisfy( long theId )
2789 if (!myMesh) return false;
2790 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2793 const SMDSAbs_ElementType anElemType = anElem->GetType();
2794 if ( myType != SMDSAbs_All && anElemType != myType )
2796 bool isOk = ( anElem->GetGeomType() == myGeomType );
2800 void ElemGeomType::SetType( SMDSAbs_ElementType theType )
2805 SMDSAbs_ElementType ElemGeomType::GetType() const
2810 void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType )
2812 myGeomType = theType;
2815 SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
2820 //================================================================================
2822 Class : ElemEntityType
2823 Description : Predicate to check element entity type
2825 //================================================================================
2827 ElemEntityType::ElemEntityType():
2829 myType( SMDSAbs_All ),
2830 myEntityType( SMDSEntity_0D )
2834 void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh )
2839 bool ElemEntityType::IsSatisfy( long theId )
2841 if ( !myMesh ) return false;
2842 if ( myType == SMDSAbs_Node )
2843 return myMesh->FindNode( theId );
2844 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
2846 myEntityType == anElem->GetEntityType() );
2849 void ElemEntityType::SetType( SMDSAbs_ElementType theType )
2854 SMDSAbs_ElementType ElemEntityType::GetType() const
2859 void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType )
2861 myEntityType = theEntityType;
2864 SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const
2866 return myEntityType;
2869 //================================================================================
2871 * \brief Class ConnectedElements
2873 //================================================================================
2875 ConnectedElements::ConnectedElements():
2876 myNodeID(0), myType( SMDSAbs_All ), myOkIDsReady( false ) {}
2878 SMDSAbs_ElementType ConnectedElements::GetType() const
2881 int ConnectedElements::GetNode() const
2882 { return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ
2884 std::vector<double> ConnectedElements::GetPoint() const
2887 void ConnectedElements::clearOkIDs()
2888 { myOkIDsReady = false; myOkIDs.clear(); }
2890 void ConnectedElements::SetType( SMDSAbs_ElementType theType )
2892 if ( myType != theType || myMeshModifTracer.IsMeshModified() )
2897 void ConnectedElements::SetMesh( const SMDS_Mesh* theMesh )
2899 myMeshModifTracer.SetMesh( theMesh );
2900 if ( myMeshModifTracer.IsMeshModified() )
2903 if ( !myXYZ.empty() )
2904 SetPoint( myXYZ[0], myXYZ[1], myXYZ[2] ); // find a node near myXYZ it in a new mesh
2908 void ConnectedElements::SetNode( int nodeID )
2913 bool isSameDomain = false;
2914 if ( myOkIDsReady && myMeshModifTracer.GetMesh() && !myMeshModifTracer.IsMeshModified() )
2915 if ( const SMDS_MeshNode* n = myMeshModifTracer.GetMesh()->FindNode( myNodeID ))
2917 SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( myType );
2918 while ( !isSameDomain && eIt->more() )
2919 isSameDomain = IsSatisfy( eIt->next()->GetID() );
2921 if ( !isSameDomain )
2925 void ConnectedElements::SetPoint( double x, double y, double z )
2933 bool isSameDomain = false;
2935 // find myNodeID by myXYZ if possible
2936 if ( myMeshModifTracer.GetMesh() )
2938 SMESHUtils::Deleter<SMESH_ElementSearcher> searcher
2939 ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() ));
2941 std::vector< const SMDS_MeshElement* > foundElems;
2942 searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems );
2944 if ( !foundElems.empty() )
2946 myNodeID = foundElems[0]->GetNode(0)->GetID();
2947 if ( myOkIDsReady && !myMeshModifTracer.IsMeshModified() )
2948 isSameDomain = IsSatisfy( foundElems[0]->GetID() );
2951 if ( !isSameDomain )
2955 bool ConnectedElements::IsSatisfy( long theElementId )
2957 // Here we do NOT check if the mesh has changed, we do it in Set...() only!!!
2959 if ( !myOkIDsReady )
2961 if ( !myMeshModifTracer.GetMesh() )
2963 const SMDS_MeshNode* node0 = myMeshModifTracer.GetMesh()->FindNode( myNodeID );
2967 std::list< const SMDS_MeshNode* > nodeQueue( 1, node0 );
2968 std::set< int > checkedNodeIDs;
2970 // foreach node in nodeQueue:
2971 // foreach element sharing a node:
2972 // add ID of an element of myType to myOkIDs;
2973 // push all element nodes absent from checkedNodeIDs to nodeQueue;
2974 while ( !nodeQueue.empty() )
2976 const SMDS_MeshNode* node = nodeQueue.front();
2977 nodeQueue.pop_front();
2979 // loop on elements sharing the node
2980 SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
2981 while ( eIt->more() )
2983 // keep elements of myType
2984 const SMDS_MeshElement* element = eIt->next();
2985 if ( element->GetType() == myType )
2986 myOkIDs.insert( myOkIDs.end(), element->GetID() );
2988 // enqueue nodes of the element
2989 SMDS_ElemIteratorPtr nIt = element->nodesIterator();
2990 while ( nIt->more() )
2992 const SMDS_MeshNode* n = static_cast< const SMDS_MeshNode* >( nIt->next() );
2993 if ( checkedNodeIDs.insert( n->GetID() ).second )
2994 nodeQueue.push_back( n );
2998 if ( myType == SMDSAbs_Node )
2999 std::swap( myOkIDs, checkedNodeIDs );
3001 size_t totalNbElems = myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType );
3002 if ( myOkIDs.size() == totalNbElems )
3005 myOkIDsReady = true;
3008 return myOkIDs.empty() ? true : myOkIDs.count( theElementId );
3011 //================================================================================
3013 * \brief Class CoplanarFaces
3015 //================================================================================
3019 inline bool isLessAngle( const gp_Vec& v1, const gp_Vec& v2, const double cos )
3021 double dot = v1 * v2; // cos * |v1| * |v2|
3022 double l1 = v1.SquareMagnitude();
3023 double l2 = v2.SquareMagnitude();
3024 return (( dot * cos >= 0 ) &&
3025 ( dot * dot ) / l1 / l2 >= ( cos * cos ));
3028 CoplanarFaces::CoplanarFaces()
3029 : myFaceID(0), myToler(0)
3032 void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
3034 myMeshModifTracer.SetMesh( theMesh );
3035 if ( myMeshModifTracer.IsMeshModified() )
3037 // Build a set of coplanar face ids
3039 myCoplanarIDs.Clear();
3041 if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler )
3044 const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID );
3045 if ( !face || face->GetType() != SMDSAbs_Face )
3049 gp_Vec myNorm = getNormale( static_cast<const SMDS_MeshFace*>(face), &normOK );
3053 const double cosTol = Cos( myToler * M_PI / 180. );
3054 NCollection_Map< SMESH_TLink, SMESH_TLink > checkedLinks;
3056 std::list< std::pair< const SMDS_MeshElement*, gp_Vec > > faceQueue;
3057 faceQueue.push_back( std::make_pair( face, myNorm ));
3058 while ( !faceQueue.empty() )
3060 face = faceQueue.front().first;
3061 myNorm = faceQueue.front().second;
3062 faceQueue.pop_front();
3064 for ( int i = 0, nbN = face->NbCornerNodes(); i < nbN; ++i )
3066 const SMDS_MeshNode* n1 = face->GetNode( i );
3067 const SMDS_MeshNode* n2 = face->GetNode(( i+1 )%nbN);
3068 if ( !checkedLinks.Add( SMESH_TLink( n1, n2 )))
3070 SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator(SMDSAbs_Face);
3071 while ( fIt->more() )
3073 const SMDS_MeshElement* f = fIt->next();
3074 if ( f->GetNodeIndex( n2 ) > -1 )
3076 gp_Vec norm = getNormale( static_cast<const SMDS_MeshFace*>(f), &normOK );
3077 if (!normOK || isLessAngle( myNorm, norm, cosTol))
3079 myCoplanarIDs.Add( f->GetID() );
3080 faceQueue.push_back( std::make_pair( f, norm ));
3088 bool CoplanarFaces::IsSatisfy( long theElementId )
3090 return myCoplanarIDs.Contains( theElementId );
3095 *Description : Predicate for Range of Ids.
3096 * Range may be specified with two ways.
3097 * 1. Using AddToRange method
3098 * 2. With SetRangeStr method. Parameter of this method is a string
3099 * like as "1,2,3,50-60,63,67,70-"
3102 //=======================================================================
3103 // name : RangeOfIds
3104 // Purpose : Constructor
3105 //=======================================================================
3106 RangeOfIds::RangeOfIds()
3109 myType = SMDSAbs_All;
3112 //=======================================================================
3114 // Purpose : Set mesh
3115 //=======================================================================
3116 void RangeOfIds::SetMesh( const SMDS_Mesh* theMesh )
3121 //=======================================================================
3122 // name : AddToRange
3123 // Purpose : Add ID to the range
3124 //=======================================================================
3125 bool RangeOfIds::AddToRange( long theEntityId )
3127 myIds.Add( theEntityId );
3131 //=======================================================================
3132 // name : GetRangeStr
3133 // Purpose : Get range as a string.
3134 // Example: "1,2,3,50-60,63,67,70-"
3135 //=======================================================================
3136 void RangeOfIds::GetRangeStr( TCollection_AsciiString& theResStr )
3140 TColStd_SequenceOfInteger anIntSeq;
3141 TColStd_SequenceOfAsciiString aStrSeq;
3143 TColStd_MapIteratorOfMapOfInteger anIter( myIds );
3144 for ( ; anIter.More(); anIter.Next() )
3146 int anId = anIter.Key();
3147 TCollection_AsciiString aStr( anId );
3148 anIntSeq.Append( anId );
3149 aStrSeq.Append( aStr );
3152 for ( int i = 1, n = myMin.Length(); i <= n; i++ )
3154 int aMinId = myMin( i );
3155 int aMaxId = myMax( i );
3157 TCollection_AsciiString aStr;
3158 if ( aMinId != IntegerFirst() )
3163 if ( aMaxId != IntegerLast() )
3166 // find position of the string in result sequence and insert string in it
3167 if ( anIntSeq.Length() == 0 )
3169 anIntSeq.Append( aMinId );
3170 aStrSeq.Append( aStr );
3174 if ( aMinId < anIntSeq.First() )
3176 anIntSeq.Prepend( aMinId );
3177 aStrSeq.Prepend( aStr );
3179 else if ( aMinId > anIntSeq.Last() )
3181 anIntSeq.Append( aMinId );
3182 aStrSeq.Append( aStr );
3185 for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
3186 if ( aMinId < anIntSeq( j ) )
3188 anIntSeq.InsertBefore( j, aMinId );
3189 aStrSeq.InsertBefore( j, aStr );
3195 if ( aStrSeq.Length() == 0 )
3198 theResStr = aStrSeq( 1 );
3199 for ( int j = 2, k = aStrSeq.Length(); j <= k; j++ )
3202 theResStr += aStrSeq( j );
3206 //=======================================================================
3207 // name : SetRangeStr
3208 // Purpose : Define range with string
3209 // Example of entry string: "1,2,3,50-60,63,67,70-"
3210 //=======================================================================
3211 bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
3217 TCollection_AsciiString aStr = theStr;
3218 for ( int i = 1; i <= aStr.Length(); ++i )
3220 char c = aStr.Value( i );
3221 if ( !isdigit( c ) && c != ',' && c != '-' )
3222 aStr.SetValue( i, ' ');
3224 aStr.RemoveAll( ' ' );
3226 TCollection_AsciiString tmpStr = aStr.Token( ",", 1 );
3228 while ( tmpStr != "" )
3230 tmpStr = aStr.Token( ",", i++ );
3231 int aPos = tmpStr.Search( '-' );
3235 if ( tmpStr.IsIntegerValue() )
3236 myIds.Add( tmpStr.IntegerValue() );
3242 TCollection_AsciiString aMaxStr = tmpStr.Split( aPos );
3243 TCollection_AsciiString aMinStr = tmpStr;
3245 while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
3246 while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
3248 if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) ||
3249 (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) )
3252 myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
3253 myMax.Append( aMaxStr.IsEmpty() ? IntegerLast() : aMaxStr.IntegerValue() );
3260 //=======================================================================
3262 // Purpose : Get type of supported entities
3263 //=======================================================================
3264 SMDSAbs_ElementType RangeOfIds::GetType() const
3269 //=======================================================================
3271 // Purpose : Set type of supported entities
3272 //=======================================================================
3273 void RangeOfIds::SetType( SMDSAbs_ElementType theType )
3278 //=======================================================================
3280 // Purpose : Verify whether entity satisfies to this rpedicate
3281 //=======================================================================
3282 bool RangeOfIds::IsSatisfy( long theId )
3287 if ( myType == SMDSAbs_Node )
3289 if ( myMesh->FindNode( theId ) == 0 )
3294 const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
3295 if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All ))
3299 if ( myIds.Contains( theId ) )
3302 for ( int i = 1, n = myMin.Length(); i <= n; i++ )
3303 if ( theId >= myMin( i ) && theId <= myMax( i ) )
3311 Description : Base class for comparators
3313 Comparator::Comparator():
3317 Comparator::~Comparator()
3320 void Comparator::SetMesh( const SMDS_Mesh* theMesh )
3323 myFunctor->SetMesh( theMesh );
3326 void Comparator::SetMargin( double theValue )
3328 myMargin = theValue;
3331 void Comparator::SetNumFunctor( NumericalFunctorPtr theFunct )
3333 myFunctor = theFunct;
3336 SMDSAbs_ElementType Comparator::GetType() const
3338 return myFunctor ? myFunctor->GetType() : SMDSAbs_All;
3341 double Comparator::GetMargin()
3349 Description : Comparator "<"
3351 bool LessThan::IsSatisfy( long theId )
3353 return myFunctor && myFunctor->GetValue( theId ) < myMargin;
3359 Description : Comparator ">"
3361 bool MoreThan::IsSatisfy( long theId )
3363 return myFunctor && myFunctor->GetValue( theId ) > myMargin;
3369 Description : Comparator "="
3372 myToler(Precision::Confusion())
3375 bool EqualTo::IsSatisfy( long theId )
3377 return myFunctor && fabs( myFunctor->GetValue( theId ) - myMargin ) < myToler;
3380 void EqualTo::SetTolerance( double theToler )
3385 double EqualTo::GetTolerance()
3392 Description : Logical NOT predicate
3394 LogicalNOT::LogicalNOT()
3397 LogicalNOT::~LogicalNOT()
3400 bool LogicalNOT::IsSatisfy( long theId )
3402 return myPredicate && !myPredicate->IsSatisfy( theId );
3405 void LogicalNOT::SetMesh( const SMDS_Mesh* theMesh )
3408 myPredicate->SetMesh( theMesh );
3411 void LogicalNOT::SetPredicate( PredicatePtr thePred )
3413 myPredicate = thePred;
3416 SMDSAbs_ElementType LogicalNOT::GetType() const
3418 return myPredicate ? myPredicate->GetType() : SMDSAbs_All;
3423 Class : LogicalBinary
3424 Description : Base class for binary logical predicate
3426 LogicalBinary::LogicalBinary()
3429 LogicalBinary::~LogicalBinary()
3432 void LogicalBinary::SetMesh( const SMDS_Mesh* theMesh )
3435 myPredicate1->SetMesh( theMesh );
3438 myPredicate2->SetMesh( theMesh );
3441 void LogicalBinary::SetPredicate1( PredicatePtr thePredicate )
3443 myPredicate1 = thePredicate;
3446 void LogicalBinary::SetPredicate2( PredicatePtr thePredicate )
3448 myPredicate2 = thePredicate;
3451 SMDSAbs_ElementType LogicalBinary::GetType() const
3453 if ( !myPredicate1 || !myPredicate2 )
3456 SMDSAbs_ElementType aType1 = myPredicate1->GetType();
3457 SMDSAbs_ElementType aType2 = myPredicate2->GetType();
3459 return aType1 == aType2 ? aType1 : SMDSAbs_All;
3465 Description : Logical AND
3467 bool LogicalAND::IsSatisfy( long theId )
3472 myPredicate1->IsSatisfy( theId ) &&
3473 myPredicate2->IsSatisfy( theId );
3479 Description : Logical OR
3481 bool LogicalOR::IsSatisfy( long theId )
3486 (myPredicate1->IsSatisfy( theId ) ||
3487 myPredicate2->IsSatisfy( theId ));
3496 // #include <tbb/parallel_for.h>
3497 // #include <tbb/enumerable_thread_specific.h>
3499 // namespace Parallel
3501 // typedef tbb::enumerable_thread_specific< TIdSequence > TIdSeq;
3505 // const SMDS_Mesh* myMesh;
3506 // PredicatePtr myPredicate;
3507 // TIdSeq & myOKIds;
3508 // Predicate( const SMDS_Mesh* m, PredicatePtr p, TIdSeq & ids ):
3509 // myMesh(m), myPredicate(p->Duplicate()), myOKIds(ids) {}
3510 // void operator() ( const tbb::blocked_range<size_t>& r ) const
3512 // for ( size_t i = r.begin(); i != r.end(); ++i )
3513 // if ( myPredicate->IsSatisfy( i ))
3514 // myOKIds.local().push_back();
3526 void Filter::SetPredicate( PredicatePtr thePredicate )
3528 myPredicate = thePredicate;
3531 void Filter::GetElementsId( const SMDS_Mesh* theMesh,
3532 PredicatePtr thePredicate,
3533 TIdSequence& theSequence )
3535 theSequence.clear();
3537 if ( !theMesh || !thePredicate )
3540 thePredicate->SetMesh( theMesh );
3542 SMDS_ElemIteratorPtr elemIt = theMesh->elementsIterator( thePredicate->GetType() );
3544 while ( elemIt->more() ) {
3545 const SMDS_MeshElement* anElem = elemIt->next();
3546 long anId = anElem->GetID();
3547 if ( thePredicate->IsSatisfy( anId ) )
3548 theSequence.push_back( anId );
3553 void Filter::GetElementsId( const SMDS_Mesh* theMesh,
3554 Filter::TIdSequence& theSequence )
3556 GetElementsId(theMesh,myPredicate,theSequence);
3563 typedef std::set<SMDS_MeshFace*> TMapOfFacePtr;
3569 ManifoldPart::Link::Link( SMDS_MeshNode* theNode1,
3570 SMDS_MeshNode* theNode2 )
3576 ManifoldPart::Link::~Link()
3582 bool ManifoldPart::Link::IsEqual( const ManifoldPart::Link& theLink ) const
3584 if ( myNode1 == theLink.myNode1 &&
3585 myNode2 == theLink.myNode2 )
3587 else if ( myNode1 == theLink.myNode2 &&
3588 myNode2 == theLink.myNode1 )
3594 bool ManifoldPart::Link::operator<( const ManifoldPart::Link& x ) const
3596 if(myNode1 < x.myNode1) return true;
3597 if(myNode1 == x.myNode1)
3598 if(myNode2 < x.myNode2) return true;
3602 bool ManifoldPart::IsEqual( const ManifoldPart::Link& theLink1,
3603 const ManifoldPart::Link& theLink2 )
3605 return theLink1.IsEqual( theLink2 );
3608 ManifoldPart::ManifoldPart()
3611 myAngToler = Precision::Angular();
3612 myIsOnlyManifold = true;
3615 ManifoldPart::~ManifoldPart()
3620 void ManifoldPart::SetMesh( const SMDS_Mesh* theMesh )
3626 SMDSAbs_ElementType ManifoldPart::GetType() const
3627 { return SMDSAbs_Face; }
3629 bool ManifoldPart::IsSatisfy( long theElementId )
3631 return myMapIds.Contains( theElementId );
3634 void ManifoldPart::SetAngleTolerance( const double theAngToler )
3635 { myAngToler = theAngToler; }
3637 double ManifoldPart::GetAngleTolerance() const
3638 { return myAngToler; }
3640 void ManifoldPart::SetIsOnlyManifold( const bool theIsOnly )
3641 { myIsOnlyManifold = theIsOnly; }
3643 void ManifoldPart::SetStartElem( const long theStartId )
3644 { myStartElemId = theStartId; }
3646 bool ManifoldPart::process()
3649 myMapBadGeomIds.Clear();
3651 myAllFacePtr.clear();
3652 myAllFacePtrIntDMap.clear();
3656 // collect all faces into own map
3657 SMDS_FaceIteratorPtr anFaceItr = myMesh->facesIterator();
3658 for (; anFaceItr->more(); )
3660 SMDS_MeshFace* aFacePtr = (SMDS_MeshFace*)anFaceItr->next();
3661 myAllFacePtr.push_back( aFacePtr );
3662 myAllFacePtrIntDMap[aFacePtr] = myAllFacePtr.size()-1;
3665 SMDS_MeshFace* aStartFace = (SMDS_MeshFace*)myMesh->FindElement( myStartElemId );
3669 // the map of non manifold links and bad geometry
3670 TMapOfLink aMapOfNonManifold;
3671 TColStd_MapOfInteger aMapOfTreated;
3673 // begin cycle on faces from start index and run on vector till the end
3674 // and from begin to start index to cover whole vector
3675 const int aStartIndx = myAllFacePtrIntDMap[aStartFace];
3676 bool isStartTreat = false;
3677 for ( int fi = aStartIndx; !isStartTreat || fi != aStartIndx ; fi++ )
3679 if ( fi == aStartIndx )
3680 isStartTreat = true;
3681 // as result next time when fi will be equal to aStartIndx
3683 SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
3684 if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
3687 aMapOfTreated.Add( aFacePtr->GetID() );
3688 TColStd_MapOfInteger aResFaces;
3689 if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
3690 aMapOfNonManifold, aResFaces ) )
3692 TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
3693 for ( ; anItr.More(); anItr.Next() )
3695 int aFaceId = anItr.Key();
3696 aMapOfTreated.Add( aFaceId );
3697 myMapIds.Add( aFaceId );
3700 if ( fi == int( myAllFacePtr.size() - 1 ))
3702 } // end run on vector of faces
3703 return !myMapIds.IsEmpty();
3706 static void getLinks( const SMDS_MeshFace* theFace,
3707 ManifoldPart::TVectorOfLink& theLinks )
3709 int aNbNode = theFace->NbNodes();
3710 SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
3712 SMDS_MeshNode* aNode = 0;
3713 for ( ; aNodeItr->more() && i <= aNbNode; )
3716 SMDS_MeshNode* aN1 = (SMDS_MeshNode*)aNodeItr->next();
3720 SMDS_MeshNode* aN2 = ( i >= aNbNode ) ? aNode : (SMDS_MeshNode*)aNodeItr->next();
3722 ManifoldPart::Link aLink( aN1, aN2 );
3723 theLinks.push_back( aLink );
3727 bool ManifoldPart::findConnected
3728 ( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
3729 SMDS_MeshFace* theStartFace,
3730 ManifoldPart::TMapOfLink& theNonManifold,
3731 TColStd_MapOfInteger& theResFaces )
3733 theResFaces.Clear();
3734 if ( !theAllFacePtrInt.size() )
3737 if ( getNormale( theStartFace ).SquareModulus() <= gp::Resolution() )
3739 myMapBadGeomIds.Add( theStartFace->GetID() );
3743 ManifoldPart::TMapOfLink aMapOfBoundary, aMapToSkip;
3744 ManifoldPart::TVectorOfLink aSeqOfBoundary;
3745 theResFaces.Add( theStartFace->GetID() );
3746 ManifoldPart::TDataMapOfLinkFacePtr aDMapLinkFace;
3748 expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3749 aDMapLinkFace, theNonManifold, theStartFace );
3751 bool isDone = false;
3752 while ( !isDone && aMapOfBoundary.size() != 0 )
3754 bool isToReset = false;
3755 ManifoldPart::TVectorOfLink::iterator pLink = aSeqOfBoundary.begin();
3756 for ( ; !isToReset && pLink != aSeqOfBoundary.end(); ++pLink )
3758 ManifoldPart::Link aLink = *pLink;
3759 if ( aMapToSkip.find( aLink ) != aMapToSkip.end() )
3761 // each link could be treated only once
3762 aMapToSkip.insert( aLink );
3764 ManifoldPart::TVectorOfFacePtr aFaces;
3766 if ( myIsOnlyManifold &&
3767 (theNonManifold.find( aLink ) != theNonManifold.end()) )
3771 getFacesByLink( aLink, aFaces );
3772 // filter the element to keep only indicated elements
3773 ManifoldPart::TVectorOfFacePtr aFiltered;
3774 ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3775 for ( ; pFace != aFaces.end(); ++pFace )
3777 SMDS_MeshFace* aFace = *pFace;
3778 if ( myAllFacePtrIntDMap.find( aFace ) != myAllFacePtrIntDMap.end() )
3779 aFiltered.push_back( aFace );
3782 if ( aFaces.size() < 2 ) // no neihgbour faces
3784 else if ( myIsOnlyManifold && aFaces.size() > 2 ) // non manifold case
3786 theNonManifold.insert( aLink );
3791 // compare normal with normals of neighbor element
3792 SMDS_MeshFace* aPrevFace = aDMapLinkFace[ aLink ];
3793 ManifoldPart::TVectorOfFacePtr::iterator pFace = aFaces.begin();
3794 for ( ; pFace != aFaces.end(); ++pFace )
3796 SMDS_MeshFace* aNextFace = *pFace;
3797 if ( aPrevFace == aNextFace )
3799 int anNextFaceID = aNextFace->GetID();
3800 if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
3801 // should not be with non manifold restriction. probably bad topology
3803 // check if face was treated and skipped
3804 if ( myMapBadGeomIds.Contains( anNextFaceID ) ||
3805 !isInPlane( aPrevFace, aNextFace ) )
3807 // add new element to connected and extend the boundaries.
3808 theResFaces.Add( anNextFaceID );
3809 expandBoundary( aMapOfBoundary, aSeqOfBoundary,
3810 aDMapLinkFace, theNonManifold, aNextFace );
3814 isDone = !isToReset;
3817 return !theResFaces.IsEmpty();
3820 bool ManifoldPart::isInPlane( const SMDS_MeshFace* theFace1,
3821 const SMDS_MeshFace* theFace2 )
3823 gp_Dir aNorm1 = gp_Dir( getNormale( theFace1 ) );
3824 gp_XYZ aNorm2XYZ = getNormale( theFace2 );
3825 if ( aNorm2XYZ.SquareModulus() <= gp::Resolution() )
3827 myMapBadGeomIds.Add( theFace2->GetID() );
3830 if ( aNorm1.IsParallel( gp_Dir( aNorm2XYZ ), myAngToler ) )
3836 void ManifoldPart::expandBoundary
3837 ( ManifoldPart::TMapOfLink& theMapOfBoundary,
3838 ManifoldPart::TVectorOfLink& theSeqOfBoundary,
3839 ManifoldPart::TDataMapOfLinkFacePtr& theDMapLinkFacePtr,
3840 ManifoldPart::TMapOfLink& theNonManifold,
3841 SMDS_MeshFace* theNextFace ) const
3843 ManifoldPart::TVectorOfLink aLinks;
3844 getLinks( theNextFace, aLinks );
3845 int aNbLink = (int)aLinks.size();
3846 for ( int i = 0; i < aNbLink; i++ )
3848 ManifoldPart::Link aLink = aLinks[ i ];
3849 if ( myIsOnlyManifold && (theNonManifold.find( aLink ) != theNonManifold.end()) )
3851 if ( theMapOfBoundary.find( aLink ) != theMapOfBoundary.end() )
3853 if ( myIsOnlyManifold )
3855 // remove from boundary
3856 theMapOfBoundary.erase( aLink );
3857 ManifoldPart::TVectorOfLink::iterator pLink = theSeqOfBoundary.begin();
3858 for ( ; pLink != theSeqOfBoundary.end(); ++pLink )
3860 ManifoldPart::Link aBoundLink = *pLink;
3861 if ( aBoundLink.IsEqual( aLink ) )
3863 theSeqOfBoundary.erase( pLink );
3871 theMapOfBoundary.insert( aLink );
3872 theSeqOfBoundary.push_back( aLink );
3873 theDMapLinkFacePtr[ aLink ] = theNextFace;
3878 void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
3879 ManifoldPart::TVectorOfFacePtr& theFaces ) const
3881 std::set<SMDS_MeshCell *> aSetOfFaces;
3882 // take all faces that shared first node
3883 SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
3884 for ( ; anItr->more(); )
3886 SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
3889 aSetOfFaces.insert( aFace );
3891 // take all faces that shared second node
3892 anItr = theLink.myNode2->facesIterator();
3893 // find the common part of two sets
3894 for ( ; anItr->more(); )
3896 SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
3897 if ( aSetOfFaces.count( aFace ) )
3898 theFaces.push_back( aFace );
3903 Class : BelongToMeshGroup
3904 Description : Verify whether a mesh element is included into a mesh group
3906 BelongToMeshGroup::BelongToMeshGroup(): myGroup( 0 )
3910 void BelongToMeshGroup::SetGroup( SMESHDS_GroupBase* g )
3915 void BelongToMeshGroup::SetStoreName( const std::string& sn )
3920 void BelongToMeshGroup::SetMesh( const SMDS_Mesh* theMesh )
3922 if ( myGroup && myGroup->GetMesh() != theMesh )
3926 if ( !myGroup && !myStoreName.empty() )
3928 if ( const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh))
3930 const std::set<SMESHDS_GroupBase*>& grps = aMesh->GetGroups();
3931 std::set<SMESHDS_GroupBase*>::const_iterator g = grps.begin();
3932 for ( ; g != grps.end() && !myGroup; ++g )
3933 if ( *g && myStoreName == (*g)->GetStoreName() )
3939 myGroup->IsEmpty(); // make GroupOnFilter update its predicate
3943 bool BelongToMeshGroup::IsSatisfy( long theElementId )
3945 return myGroup ? myGroup->Contains( theElementId ) : false;
3948 SMDSAbs_ElementType BelongToMeshGroup::GetType() const
3950 return myGroup ? myGroup->GetType() : SMDSAbs_All;
3957 ElementsOnSurface::ElementsOnSurface()
3960 myType = SMDSAbs_All;
3962 myToler = Precision::Confusion();
3963 myUseBoundaries = false;
3966 ElementsOnSurface::~ElementsOnSurface()
3970 void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
3972 myMeshModifTracer.SetMesh( theMesh );
3973 if ( myMeshModifTracer.IsMeshModified())
3977 bool ElementsOnSurface::IsSatisfy( long theElementId )
3979 return myIds.Contains( theElementId );
3982 SMDSAbs_ElementType ElementsOnSurface::GetType() const
3985 void ElementsOnSurface::SetTolerance( const double theToler )
3987 if ( myToler != theToler )
3992 double ElementsOnSurface::GetTolerance() const
3995 void ElementsOnSurface::SetUseBoundaries( bool theUse )
3997 if ( myUseBoundaries != theUse ) {
3998 myUseBoundaries = theUse;
3999 SetSurface( mySurf, myType );
4003 void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
4004 const SMDSAbs_ElementType theType )
4009 if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
4011 mySurf = TopoDS::Face( theShape );
4012 BRepAdaptor_Surface SA( mySurf, myUseBoundaries );
4014 u1 = SA.FirstUParameter(),
4015 u2 = SA.LastUParameter(),
4016 v1 = SA.FirstVParameter(),
4017 v2 = SA.LastVParameter();
4018 Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf );
4019 myProjector.Init( surf, u1,u2, v1,v2 );
4023 void ElementsOnSurface::process()
4026 if ( mySurf.IsNull() )
4029 if ( !myMeshModifTracer.GetMesh() )
4032 myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
4034 SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType );
4035 for(; anIter->more(); )
4036 process( anIter->next() );
4039 void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
4041 SMDS_ElemIteratorPtr aNodeItr = theElemPtr->nodesIterator();
4042 bool isSatisfy = true;
4043 for ( ; aNodeItr->more(); )
4045 SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
4046 if ( !isOnSurface( aNode ) )
4053 myIds.Add( theElemPtr->GetID() );
4056 bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode )
4058 if ( mySurf.IsNull() )
4061 gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
4062 // double aToler2 = myToler * myToler;
4063 // if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
4065 // gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
4066 // if ( aPln.SquareDistance( aPnt ) > aToler2 )
4069 // else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
4071 // gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
4072 // double aRad = aCyl.Radius();
4073 // gp_Ax3 anAxis = aCyl.Position();
4074 // gp_XYZ aLoc = aCyl.Location().XYZ();
4075 // double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
4076 // double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
4077 // if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
4082 myProjector.Perform( aPnt );
4083 bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler );
4093 struct ElementsOnShape::Classifier
4095 Classifier(const TopoDS_Shape& s, double tol) { Init(s,tol); }
4096 void Init(const TopoDS_Shape& s, double tol);
4097 bool IsOut(const gp_Pnt& p) { return myIsChecked = true, (this->*myIsOutFun)( p ); }
4098 TopAbs_ShapeEnum ShapeType() const { return myShape.ShapeType(); }
4099 Bnd_B3d* GetBndBox() { return & myBox; }
4100 bool& IsChecked() { return myIsChecked; }
4102 bool isOutOfSolid (const gp_Pnt& p);
4103 bool isOutOfBox (const gp_Pnt& p);
4104 bool isOutOfFace (const gp_Pnt& p);
4105 bool isOutOfEdge (const gp_Pnt& p);
4106 bool isOutOfVertex(const gp_Pnt& p);
4107 bool isBox (const TopoDS_Shape& s);
4109 bool (Classifier::* myIsOutFun)(const gp_Pnt& p);
4110 BRepClass3d_SolidClassifier mySolidClfr;
4112 GeomAPI_ProjectPointOnSurf myProjFace;
4113 GeomAPI_ProjectPointOnCurve myProjEdge;
4115 TopoDS_Shape myShape;
4120 struct ElementsOnShape::OctreeClassifier : public SMESH_Octree
4122 OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers );
4123 void GetClassifiersAtPoint( const gp_XYZ& p,
4124 std::vector< ElementsOnShape::Classifier* >& classifiers );
4126 OctreeClassifier() {}
4127 SMESH_Octree* newChild() const { return new OctreeClassifier; }
4128 void buildChildrenData();
4129 Bnd_B3d* buildRootBox();
4131 std::vector< ElementsOnShape::Classifier* > myClassifiers;
4135 ElementsOnShape::ElementsOnShape():
4137 myType(SMDSAbs_All),
4138 myToler(Precision::Confusion()),
4139 myAllNodesFlag(false)
4143 ElementsOnShape::~ElementsOnShape()
4148 SMDSAbs_ElementType ElementsOnShape::GetType() const
4153 void ElementsOnShape::SetTolerance (const double theToler)
4155 if (myToler != theToler) {
4157 SetShape(myShape, myType);
4161 double ElementsOnShape::GetTolerance() const
4166 void ElementsOnShape::SetAllNodes (bool theAllNodes)
4168 myAllNodesFlag = theAllNodes;
4171 void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
4173 myMeshModifTracer.SetMesh( theMesh );
4174 if ( myMeshModifTracer.IsMeshModified())
4176 size_t nbNodes = theMesh ? theMesh->NbNodes() : 0;
4177 if ( myNodeIsChecked.size() == nbNodes )
4179 std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4183 SMESHUtils::FreeVector( myNodeIsChecked );
4184 SMESHUtils::FreeVector( myNodeIsOut );
4185 myNodeIsChecked.resize( nbNodes, false );
4186 myNodeIsOut.resize( nbNodes );
4191 bool ElementsOnShape::getNodeIsOut( const SMDS_MeshNode* n, bool& isOut )
4193 if ( n->GetID() >= (int) myNodeIsChecked.size() ||
4194 !myNodeIsChecked[ n->GetID() ])
4197 isOut = myNodeIsOut[ n->GetID() ];
4201 void ElementsOnShape::setNodeIsOut( const SMDS_MeshNode* n, bool isOut )
4203 if ( n->GetID() < (int) myNodeIsChecked.size() )
4205 myNodeIsChecked[ n->GetID() ] = true;
4206 myNodeIsOut [ n->GetID() ] = isOut;
4210 void ElementsOnShape::SetShape (const TopoDS_Shape& theShape,
4211 const SMDSAbs_ElementType theType)
4213 bool shapeChanges = ( myShape != theShape );
4216 if ( myShape.IsNull() ) return;
4220 TopTools_IndexedMapOfShape shapesMap;
4221 TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX };
4222 TopExp_Explorer sub;
4223 for ( int i = 0; i < 4; ++i )
4225 if ( shapesMap.IsEmpty() )
4226 for ( sub.Init( myShape, shapeTypes[i] ); sub.More(); sub.Next() )
4227 shapesMap.Add( sub.Current() );
4229 for ( sub.Init( myShape, shapeTypes[i], shapeTypes[i-1] ); sub.More(); sub.Next() )
4230 shapesMap.Add( sub.Current() );
4234 myClassifiers.resize( shapesMap.Extent() );
4235 for ( int i = 0; i < shapesMap.Extent(); ++i )
4236 myClassifiers[ i ] = new Classifier( shapesMap( i+1 ), myToler );
4239 if ( theType == SMDSAbs_Node )
4241 SMESHUtils::FreeVector( myNodeIsChecked );
4242 SMESHUtils::FreeVector( myNodeIsOut );
4246 std::fill( myNodeIsChecked.begin(), myNodeIsChecked.end(), false );
4250 void ElementsOnShape::clearClassifiers()
4252 for ( size_t i = 0; i < myClassifiers.size(); ++i )
4253 delete myClassifiers[ i ];
4254 myClassifiers.clear();
4260 bool ElementsOnShape::IsSatisfy (long elemId)
4262 const SMDS_Mesh* mesh = myMeshModifTracer.GetMesh();
4263 const SMDS_MeshElement* elem =
4264 ( myType == SMDSAbs_Node ? mesh->FindNode( elemId ) : mesh->FindElement( elemId ));
4265 if ( !elem || myClassifiers.empty() )
4268 bool isSatisfy = myAllNodesFlag, isNodeOut;
4270 gp_XYZ centerXYZ (0, 0, 0);
4272 if ( !myOctree && myClassifiers.size() > 5 )
4273 myOctree = new OctreeClassifier( myClassifiers );
4275 std::vector< Classifier* >& classifiers = myOctree ? myWorkClassifiers : myClassifiers;
4277 SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator();
4278 while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
4280 SMESH_TNodeXYZ aPnt( aNodeItr->next() );
4284 if ( !getNodeIsOut( aPnt._node, isNodeOut ))
4288 myWorkClassifiers.clear();
4289 myOctree->GetClassifiersAtPoint( aPnt, myWorkClassifiers );
4291 for ( size_t i = 0; i < classifiers.size(); ++i )
4292 classifiers[i]->IsChecked() = false;
4294 for ( size_t i = 0; i < classifiers.size() && isNodeOut; ++i )
4295 if ( !classifiers[i]->IsChecked() )
4296 isNodeOut = classifiers[i]->IsOut( aPnt );
4298 setNodeIsOut( aPnt._node, isNodeOut );
4300 isSatisfy = !isNodeOut;
4303 // Check the center point for volumes MantisBug 0020168
4306 myClassifiers[0]->ShapeType() == TopAbs_SOLID)
4308 centerXYZ /= elem->NbNodes();
4310 for ( size_t i = 0; i < classifiers.size() && !isSatisfy; ++i )
4311 isSatisfy = ! classifiers[i]->IsOut( centerXYZ );
4317 void ElementsOnShape::Classifier::Init (const TopoDS_Shape& theShape, double theTol)
4322 bool isShapeBox = false;
4323 switch ( myShape.ShapeType() )
4327 if (( isShapeBox = isBox( theShape )))
4329 myIsOutFun = & ElementsOnShape::Classifier::isOutOfBox;
4333 mySolidClfr.Load(theShape);
4334 myIsOutFun = & ElementsOnShape::Classifier::isOutOfSolid;
4340 Standard_Real u1,u2,v1,v2;
4341 Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape ));
4342 surf->Bounds( u1,u2,v1,v2 );
4343 myProjFace.Init(surf, u1,u2, v1,v2, myTol );
4344 myIsOutFun = & ElementsOnShape::Classifier::isOutOfFace;
4349 Standard_Real u1, u2;
4350 Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( theShape ), u1, u2);
4351 myProjEdge.Init(curve, u1, u2);
4352 myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge;
4357 myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) );
4358 myIsOutFun = & ElementsOnShape::Classifier::isOutOfVertex;
4362 throw SALOME_Exception("Programmer error in usage of ElementsOnShape::Classifier");
4368 BRepBndLib::Add( myShape, box );
4370 myBox.Add( box.CornerMin() );
4371 myBox.Add( box.CornerMax() );
4372 gp_XYZ halfSize = 0.5 * ( box.CornerMax().XYZ() - box.CornerMin().XYZ() );
4373 for ( int iDim = 1; iDim <= 3; ++iDim )
4375 double x = halfSize.Coord( iDim );
4376 halfSize.SetCoord( iDim, x + Max( myTol, 1e-2 * x ));
4378 myBox.SetHSize( halfSize );
4382 bool ElementsOnShape::Classifier::isOutOfSolid (const gp_Pnt& p)
4384 mySolidClfr.Perform( p, myTol );
4385 return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON );
4388 bool ElementsOnShape::Classifier::isOutOfBox (const gp_Pnt& p)
4390 return myBox.IsOut( p.XYZ() );
4393 bool ElementsOnShape::Classifier::isOutOfFace (const gp_Pnt& p)
4395 myProjFace.Perform( p );
4396 if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol )
4398 // check relatively to the face
4399 Quantity_Parameter u, v;
4400 myProjFace.LowerDistanceParameters(u, v);
4401 gp_Pnt2d aProjPnt (u, v);
4402 BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol );
4403 if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON )
4409 bool ElementsOnShape::Classifier::isOutOfEdge (const gp_Pnt& p)
4411 myProjEdge.Perform( p );
4412 return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
4415 bool ElementsOnShape::Classifier::isOutOfVertex(const gp_Pnt& p)
4417 return ( myVertexXYZ.Distance( p ) > myTol );
4420 bool ElementsOnShape::Classifier::isBox (const TopoDS_Shape& theShape)
4422 TopTools_IndexedMapOfShape vMap;
4423 TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
4424 if ( vMap.Extent() != 8 )
4428 for ( int i = 1; i <= 8; ++i )
4429 myBox.Add( BRep_Tool::Pnt( TopoDS::Vertex( vMap( i ))).XYZ() );
4431 gp_XYZ pMin = myBox.CornerMin(), pMax = myBox.CornerMax();
4432 for ( int i = 1; i <= 8; ++i )
4434 gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vMap( i )));
4435 for ( int iC = 1; iC <= 3; ++ iC )
4437 double d1 = Abs( pMin.Coord( iC ) - p.Coord( iC ));
4438 double d2 = Abs( pMax.Coord( iC ) - p.Coord( iC ));
4439 if ( Min( d1, d2 ) > myTol )
4443 myBox.Enlarge( myTol );
4448 OctreeClassifier::OctreeClassifier( const std::vector< ElementsOnShape::Classifier* >& classifiers )
4449 :SMESH_Octree( new SMESH_TreeLimit )
4451 myClassifiers = classifiers;
4455 void ElementsOnShape::
4456 OctreeClassifier::GetClassifiersAtPoint( const gp_XYZ& point,
4457 std::vector< ElementsOnShape::Classifier* >& result )
4459 if ( getBox()->IsOut( point ))
4464 for ( size_t i = 0; i < myClassifiers.size(); ++i )
4465 if ( !myClassifiers[i]->GetBndBox()->IsOut( point ))
4466 result.push_back( myClassifiers[i] );
4470 for (int i = 0; i < nbChildren(); i++)
4471 ((OctreeClassifier*) myChildren[i])->GetClassifiersAtPoint( point, result );
4475 void ElementsOnShape::OctreeClassifier::buildChildrenData()
4477 // distribute myClassifiers among myChildren
4478 for ( size_t i = 0; i < myClassifiers.size(); ++i )
4480 for (int j = 0; j < nbChildren(); j++)
4482 if ( !myClassifiers[i]->GetBndBox()->IsOut( *myChildren[j]->getBox() ))
4484 ((OctreeClassifier*)myChildren[j])->myClassifiers.push_back( myClassifiers[i]);
4488 SMESHUtils::FreeVector( myClassifiers );
4490 // define if a child isLeaf()
4491 for (int i = 0; i < nbChildren(); i++)
4493 OctreeClassifier* child = static_cast<OctreeClassifier*>( myChildren[i]);
4494 child->myIsLeaf = ( child->myClassifiers.size() <= 5 );
4496 if ( child->myClassifiers.capacity() - child->myClassifiers.size() > 100 )
4497 SMESHUtils::CompactVector( child->myClassifiers );
4501 Bnd_B3d* ElementsOnShape::OctreeClassifier::buildRootBox()
4503 Bnd_B3d* box = new Bnd_B3d;
4504 for ( size_t i = 0; i < myClassifiers.size(); ++i )
4505 box->Add( *myClassifiers[i]->GetBndBox() );
4510 Class : BelongToGeom
4511 Description : Predicate for verifying whether entity belongs to
4512 specified geometrical support
4515 BelongToGeom::BelongToGeom()
4517 myType(SMDSAbs_All),
4518 myIsSubshape(false),
4519 myTolerance(Precision::Confusion())
4522 void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh )
4524 myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
4528 void BelongToGeom::SetGeom( const TopoDS_Shape& theShape )
4534 static bool IsSubShape (const TopTools_IndexedMapOfShape& theMap,
4535 const TopoDS_Shape& theShape)
4537 if (theMap.Contains(theShape)) return true;
4539 if (theShape.ShapeType() == TopAbs_COMPOUND ||
4540 theShape.ShapeType() == TopAbs_COMPSOLID)
4542 TopoDS_Iterator anIt (theShape, Standard_True, Standard_True);
4543 for (; anIt.More(); anIt.Next())
4545 if (!IsSubShape(theMap, anIt.Value())) {
4555 void BelongToGeom::init()
4557 if (!myMeshDS || myShape.IsNull()) return;
4559 // is sub-shape of main shape?
4560 TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
4561 if (aMainShape.IsNull()) {
4562 myIsSubshape = false;
4565 TopTools_IndexedMapOfShape aMap;
4566 TopExp::MapShapes(aMainShape, aMap);
4567 myIsSubshape = IsSubShape(aMap, myShape);
4570 //if (!myIsSubshape) // to be always ready to check an element not bound to geometry
4572 if ( !myElementsOnShapePtr )
4573 myElementsOnShapePtr.reset( new ElementsOnShape() );
4574 myElementsOnShapePtr->SetTolerance( myTolerance );
4575 myElementsOnShapePtr->SetAllNodes( true ); // "belong", while false means "lays on"
4576 myElementsOnShapePtr->SetMesh( myMeshDS );
4577 myElementsOnShapePtr->SetShape( myShape, myType );
4581 static bool IsContains( const SMESHDS_Mesh* theMeshDS,
4582 const TopoDS_Shape& theShape,
4583 const SMDS_MeshElement* theElem,
4584 TopAbs_ShapeEnum theFindShapeEnum,
4585 TopAbs_ShapeEnum theAvoidShapeEnum = TopAbs_SHAPE )
4587 TopExp_Explorer anExp( theShape,theFindShapeEnum,theAvoidShapeEnum );
4589 while( anExp.More() )
4591 const TopoDS_Shape& aShape = anExp.Current();
4592 if( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape ) ){
4593 if( aSubMesh->Contains( theElem ) )
4601 bool BelongToGeom::IsSatisfy (long theId)
4603 if (myMeshDS == 0 || myShape.IsNull())
4608 return myElementsOnShapePtr->IsSatisfy(theId);
4613 if (myType == SMDSAbs_Node)
4615 if( const SMDS_MeshNode* aNode = myMeshDS->FindNode( theId ) )
4617 if ( aNode->getshapeId() < 1 )
4618 return myElementsOnShapePtr->IsSatisfy(theId);
4620 const SMDS_PositionPtr& aPosition = aNode->GetPosition();
4621 SMDS_TypeOfPosition aTypeOfPosition = aPosition->GetTypeOfPosition();
4622 switch( aTypeOfPosition )
4624 case SMDS_TOP_VERTEX : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_VERTEX ));
4625 case SMDS_TOP_EDGE : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_EDGE ));
4626 case SMDS_TOP_FACE : return ( IsContains( myMeshDS,myShape,aNode,TopAbs_FACE ));
4627 case SMDS_TOP_3DSPACE: return ( IsContains( myMeshDS,myShape,aNode,TopAbs_SOLID ) ||
4628 IsContains( myMeshDS,myShape,aNode,TopAbs_SHELL ));
4635 if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ))
4637 if ( anElem->getshapeId() < 1 )
4638 return myElementsOnShapePtr->IsSatisfy(theId);
4640 if( myType == SMDSAbs_All )
4642 return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE ) ||
4643 IsContains( myMeshDS,myShape,anElem,TopAbs_FACE ) ||
4644 IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )||
4645 IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL ));
4647 else if( myType == anElem->GetType() )
4651 case SMDSAbs_Edge : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_EDGE ));
4652 case SMDSAbs_Face : return ( IsContains( myMeshDS,myShape,anElem,TopAbs_FACE ));
4653 case SMDSAbs_Volume: return ( IsContains( myMeshDS,myShape,anElem,TopAbs_SOLID )||
4654 IsContains( myMeshDS,myShape,anElem,TopAbs_SHELL ));
4664 void BelongToGeom::SetType (SMDSAbs_ElementType theType)
4670 SMDSAbs_ElementType BelongToGeom::GetType() const
4675 TopoDS_Shape BelongToGeom::GetShape()
4680 const SMESHDS_Mesh* BelongToGeom::GetMeshDS() const
4685 void BelongToGeom::SetTolerance (double theTolerance)
4687 myTolerance = theTolerance;
4692 double BelongToGeom::GetTolerance()
4699 Description : Predicate for verifying whether entiy lying or partially lying on
4700 specified geometrical support
4703 LyingOnGeom::LyingOnGeom()
4705 myType(SMDSAbs_All),
4706 myIsSubshape(false),
4707 myTolerance(Precision::Confusion())
4710 void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh )
4712 myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
4716 void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape )
4722 void LyingOnGeom::init()
4724 if (!myMeshDS || myShape.IsNull()) return;
4726 // is sub-shape of main shape?
4727 TopoDS_Shape aMainShape = myMeshDS->ShapeToMesh();
4728 if (aMainShape.IsNull()) {
4729 myIsSubshape = false;
4732 myIsSubshape = myMeshDS->IsGroupOfSubShapes( myShape );
4737 TopTools_IndexedMapOfShape shapes;
4738 TopExp::MapShapes( myShape, shapes );
4739 mySubShapesIDs.Clear();
4740 for ( int i = 1; i <= shapes.Extent(); ++i )
4742 int subID = myMeshDS->ShapeToIndex( shapes( i ));
4744 mySubShapesIDs.Add( subID );
4747 // else // to be always ready to check an element not bound to geometry
4749 if ( !myElementsOnShapePtr )
4750 myElementsOnShapePtr.reset( new ElementsOnShape() );
4751 myElementsOnShapePtr->SetTolerance( myTolerance );
4752 myElementsOnShapePtr->SetAllNodes( false ); // lays on, while true means "belong"
4753 myElementsOnShapePtr->SetMesh( myMeshDS );
4754 myElementsOnShapePtr->SetShape( myShape, myType );
4758 bool LyingOnGeom::IsSatisfy( long theId )
4760 if ( myMeshDS == 0 || myShape.IsNull() )
4765 return myElementsOnShapePtr->IsSatisfy(theId);
4770 const SMDS_MeshElement* elem =
4771 ( myType == SMDSAbs_Node ) ? myMeshDS->FindNode( theId ) : myMeshDS->FindElement( theId );
4773 if ( mySubShapesIDs.Contains( elem->getshapeId() ))
4776 if ( elem->GetType() != SMDSAbs_Node )
4778 SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator();
4779 while ( nodeItr->more() )
4781 const SMDS_MeshElement* aNode = nodeItr->next();
4782 if ( mySubShapesIDs.Contains( aNode->getshapeId() ))
4790 void LyingOnGeom::SetType( SMDSAbs_ElementType theType )
4796 SMDSAbs_ElementType LyingOnGeom::GetType() const
4801 TopoDS_Shape LyingOnGeom::GetShape()
4806 const SMESHDS_Mesh* LyingOnGeom::GetMeshDS() const
4811 void LyingOnGeom::SetTolerance (double theTolerance)
4813 myTolerance = theTolerance;
4818 double LyingOnGeom::GetTolerance()
4823 bool LyingOnGeom::Contains( const SMESHDS_Mesh* theMeshDS,
4824 const TopoDS_Shape& theShape,
4825 const SMDS_MeshElement* theElem,
4826 TopAbs_ShapeEnum theFindShapeEnum,
4827 TopAbs_ShapeEnum theAvoidShapeEnum )
4829 // if (IsContains(theMeshDS, theShape, theElem, theFindShapeEnum, theAvoidShapeEnum))
4832 // TopTools_MapOfShape aSubShapes;
4833 // TopExp_Explorer exp( theShape, theFindShapeEnum, theAvoidShapeEnum );
4834 // for ( ; exp.More(); exp.Next() )
4836 // const TopoDS_Shape& aShape = exp.Current();
4837 // if ( !aSubShapes.Add( aShape )) continue;
4839 // if ( SMESHDS_SubMesh* aSubMesh = theMeshDS->MeshElements( aShape ))
4841 // if ( aSubMesh->Contains( theElem ))
4844 // SMDS_ElemIteratorPtr nodeItr = theElem->nodesIterator();
4845 // while ( nodeItr->more() )
4847 // const SMDS_MeshElement* aNode = nodeItr->next();
4848 // if ( aSubMesh->Contains( aNode ))
4856 TSequenceOfXYZ::TSequenceOfXYZ(): myElem(0)
4859 TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n), myElem(0)
4862 TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t), myElem(0)
4865 TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray), myElem(theSequenceOfXYZ.myElem)
4868 template <class InputIterator>
4869 TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd), myElem(0)
4872 TSequenceOfXYZ::~TSequenceOfXYZ()
4875 TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
4877 myArray = theSequenceOfXYZ.myArray;
4878 myElem = theSequenceOfXYZ.myElem;
4882 gp_XYZ& TSequenceOfXYZ::operator()(size_type n)
4884 return myArray[n-1];
4887 const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const
4889 return myArray[n-1];
4892 void TSequenceOfXYZ::clear()
4897 void TSequenceOfXYZ::reserve(size_type n)
4902 void TSequenceOfXYZ::push_back(const gp_XYZ& v)
4904 myArray.push_back(v);
4907 TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
4909 return myArray.size();
4912 SMDSAbs_EntityType TSequenceOfXYZ::getElementEntity() const
4914 return myElem ? myElem->GetEntityType() : SMDSEntity_Last;
4917 TMeshModifTracer::TMeshModifTracer():
4918 myMeshModifTime(0), myMesh(0)
4921 void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh )
4923 if ( theMesh != myMesh )
4924 myMeshModifTime = 0;
4927 bool TMeshModifTracer::IsMeshModified()
4929 bool modified = false;
4932 modified = ( myMeshModifTime != myMesh->GetMTime() );
4933 myMeshModifTime = myMesh->GetMTime();