+ return chn.size() > 1;
+ }
+ //================================================================================
+ /*!
+ * \brief Try to heal a polygon whose ends are not connected
+ */
+ bool Hexahedron::closePolygon( _Face* polygon, vector<_Node*>& chainNodes ) const
+ {
+ int i = -1, nbLinks = polygon->_links.size();
+ if ( nbLinks < 3 )
+ return false;
+ vector< _OrientedLink > newLinks;
+ // find a node lying on the same FACE as the last one
+ _Node* node = polygon->_links.back().LastNode();
+ int avoidFace = node->IsLinked( polygon->_links.back().FirstNode()->_intPoint );
+ for ( i = nbLinks - 2; i >= 0; --i )
+ if ( node->IsLinked( polygon->_links[i].FirstNode()->_intPoint, avoidFace ))
+ break;
+ if ( i >= 0 )
+ {
+ for ( ; i < nbLinks; ++i )
+ newLinks.push_back( polygon->_links[i] );
+ }
+ else
+ {
+ // find a node lying on the same FACE as the first one
+ node = polygon->_links[0].FirstNode();
+ avoidFace = node->IsLinked( polygon->_links[0].LastNode()->_intPoint );
+ for ( i = 1; i < nbLinks; ++i )
+ if ( node->IsLinked( polygon->_links[i].LastNode()->_intPoint, avoidFace ))
+ break;
+ if ( i < nbLinks )
+ for ( nbLinks = i + 1, i = 0; i < nbLinks; ++i )
+ newLinks.push_back( polygon->_links[i] );
+ }
+ if ( newLinks.size() > 1 )
+ {
+ polygon->_links.swap( newLinks );
+ chainNodes.clear();
+ chainNodes.push_back( polygon->_links.back().LastNode() );
+ chainNodes.push_back( polygon->_links[0].FirstNode() );
+ return true;
+ }
+ return false;
+ }
+ //================================================================================
+ /*!
+ * \brief Finds nodes on the same EDGE as the first node of avoidSplit.
+ *
+ * This function is for a case where an EDGE lies on a quad which lies on a FACE
+ * so that a part of quad in ON and another part in IN
+ */
+ bool Hexahedron::findChainOnEdge( const vector< _OrientedLink >& splits,
+ const _OrientedLink& prevSplit,
+ const _OrientedLink& avoidSplit,
+ size_t & iS,
+ _Face& quad,
+ vector<_Node*>& chn )
+ {
+ if ( !isImplementEdges() )
+ return false;
+
+ _Node* pn1 = prevSplit.FirstNode();
+ _Node* pn2 = prevSplit.LastNode();
+ int avoidFace = pn1->IsLinked( pn2->_intPoint ); // FACE under the quad
+ if ( avoidFace < 1 && pn1->_intPoint )
+ return false;
+
+ _Node* n, *stopNode = avoidSplit.LastNode();
+
+ chn.clear();
+ if ( !quad._eIntNodes.empty() )
+ {
+ chn.push_back( pn2 );
+ bool found;
+ do
+ {
+ found = false;
+ for ( size_t iP = 0; iP < quad._eIntNodes.size(); ++iP )
+ if (( !quad._eIntNodes[ iP ]->IsUsedInFace( &quad )) &&
+ ( chn.back()->IsLinked( quad._eIntNodes[ iP ]->_intPoint, avoidFace )) &&
+ ( !avoidFace || quad._eIntNodes[ iP ]->IsOnFace( avoidFace )))
+ {
+ chn.push_back( quad._eIntNodes[ iP ]);
+ found = ( quad._eIntNodes[ iP ]->_usedInFace = &quad );
+ break;
+ }
+ } while ( found );
+ pn2 = chn.back();
+ }
+
+ int i;
+ for ( i = splits.size()-1; i >= 0; --i )
+ {
+ if ( !splits[i] )
+ continue;
+
+ n = splits[i].LastNode();
+ if ( n == stopNode )
+ break;
+ if (( n != pn1 ) &&
+ ( n->IsLinked( pn2->_intPoint, avoidFace )) &&
+ ( !avoidFace || n->IsOnFace( avoidFace )))
+ break;
+
+ n = splits[i].FirstNode();
+ if ( n == stopNode )
+ break;
+ if (( n->IsLinked( pn2->_intPoint, avoidFace )) &&
+ ( !avoidFace || n->IsOnFace( avoidFace )))
+ break;
+ n = 0;
+ }
+ if ( n && n != stopNode)
+ {
+ if ( chn.empty() )
+ chn.push_back( pn2 );
+ chn.push_back( n );
+ iS = i-1;
+ return true;
+ }
+ return false;
+ }
+ //================================================================================
+ /*!
+ * \brief Checks transition at the ginen intersection node of a link
+ */
+ bool Hexahedron::isOutPoint( _Link& link, int iP, SMESH_MesherHelper& helper ) const
+ {
+ bool isOut = false;
+
+ const bool moreIntPoints = ( iP+1 < (int) link._fIntPoints.size() );
+
+ // get 2 _Node's
+ _Node* n1 = link._fIntNodes[ iP ];
+ if ( !n1->Node() )
+ n1 = link._nodes[0];
+ _Node* n2 = moreIntPoints ? link._fIntNodes[ iP+1 ] : 0;
+ if ( !n2 || !n2->Node() )
+ n2 = link._nodes[1];
+ if ( !n2->Node() )
+ return true;
+
+ // get all FACEs under n1 and n2
+ set< TGeomID > faceIDs;
+ if ( moreIntPoints ) faceIDs.insert( link._fIntPoints[iP+1]->_faceIDs.begin(),
+ link._fIntPoints[iP+1]->_faceIDs.end() );
+ if ( n2->_intPoint ) faceIDs.insert( n2->_intPoint->_faceIDs.begin(),
+ n2->_intPoint->_faceIDs.end() );
+ if ( faceIDs.empty() )
+ return false; // n2 is inside
+ if ( n1->_intPoint ) faceIDs.insert( n1->_intPoint->_faceIDs.begin(),
+ n1->_intPoint->_faceIDs.end() );
+ faceIDs.insert( link._fIntPoints[iP]->_faceIDs.begin(),
+ link._fIntPoints[iP]->_faceIDs.end() );
+
+ // get a point between 2 nodes
+ gp_Pnt p1 = n1->Point();
+ gp_Pnt p2 = n2->Point();
+ gp_Pnt pOnLink = 0.8 * p1.XYZ() + 0.2 * p2.XYZ();
+
+ TopLoc_Location loc;
+
+ set< TGeomID >::iterator faceID = faceIDs.begin();
+ for ( ; faceID != faceIDs.end(); ++faceID )
+ {
+ // project pOnLink on a FACE
+ if ( *faceID < 1 ) continue;
+ const TopoDS_Face& face = TopoDS::Face( _grid->_shapes( *faceID ));
+ GeomAPI_ProjectPointOnSurf& proj =
+ helper.GetProjector( face, loc, 0.1*_grid->_tol );
+ gp_Pnt testPnt = pOnLink.Transformed( loc.Transformation().Inverted() );
+ proj.Perform( testPnt );
+ if ( proj.IsDone() && proj.NbPoints() > 0 )
+ {
+ Standard_Real u,v;
+ proj.LowerDistanceParameters( u,v );
+
+ if ( proj.LowerDistance() <= 0.1 * _grid->_tol )
+ {
+ isOut = false;
+ }
+ else
+ {
+ // find isOut by normals
+ gp_Dir normal;
+ if ( GeomLib::NormEstim( BRep_Tool::Surface( face, loc ),
+ gp_Pnt2d( u,v ),
+ 0.1*_grid->_tol,
+ normal ) < 3 )
+ {
+ if ( face.Orientation() == TopAbs_REVERSED )
+ normal.Reverse();
+ gp_Vec v( proj.NearestPoint(), testPnt );
+ isOut = ( v * normal > 0 );
+ }
+ }
+ if ( !isOut )
+ {
+ // classify a projection
+ if ( !n1->IsOnFace( *faceID ) || !n2->IsOnFace( *faceID ))
+ {
+ BRepTopAdaptor_FClass2d cls( face, Precision::Confusion() );
+ TopAbs_State state = cls.Perform( gp_Pnt2d( u,v ));
+ if ( state == TopAbs_OUT )
+ {
+ isOut = true;
+ continue;
+ }
+ }
+ return false;
+ }
+ }
+ }
+ return isOut;
+ }
+ //================================================================================
+ /*!
+ * \brief Sort nodes on a FACE
+ */
+ void Hexahedron::sortVertexNodes(vector<_Node*>& nodes, _Node* curNode, TGeomID faceID)
+ {
+ if ( nodes.size() > 20 ) return;
+
+ // get shapes under nodes
+ TGeomID nShapeIds[20], *nShapeIdsEnd = &nShapeIds[0] + nodes.size();
+ for ( size_t i = 0; i < nodes.size(); ++i )
+ if ( !( nShapeIds[i] = nodes[i]->ShapeID() ))
+ return;
+
+ // get shapes of the FACE
+ const TopoDS_Face& face = TopoDS::Face( _grid->_shapes( faceID ));
+ list< TopoDS_Edge > edges;
+ list< int > nbEdges;
+ int nbW = SMESH_Block::GetOrderedEdges (face, edges, nbEdges);
+ if ( nbW > 1 ) {
+ // select a WIRE - remove EDGEs of irrelevant WIREs from edges
+ list< TopoDS_Edge >::iterator e = edges.begin(), eEnd = e;
+ list< int >::iterator nE = nbEdges.begin();
+ for ( ; nbW > 0; ++nE, --nbW )
+ {
+ std::advance( eEnd, *nE );
+ for ( ; e != eEnd; ++e )
+ for ( int i = 0; i < 2; ++i )
+ {
+ TGeomID id = i==0 ?
+ _grid->_shapes.FindIndex( *e ) :
+ _grid->_shapes.FindIndex( SMESH_MesherHelper::IthVertex( 0, *e ));
+ if (( id > 0 ) &&
+ ( std::find( &nShapeIds[0], nShapeIdsEnd, id ) != nShapeIdsEnd ))
+ {
+ edges.erase( eEnd, edges.end() ); // remove rest wires
+ e = eEnd = edges.end();
+ --e;
+ nbW = 0;
+ break;
+ }
+ }
+ if ( nbW > 0 )
+ edges.erase( edges.begin(), eEnd ); // remove a current irrelevant wire
+ }
+ }
+ // rotate edges to have the first one at least partially out of the hexa
+ list< TopoDS_Edge >::iterator e = edges.begin(), eMidOut = edges.end();
+ for ( ; e != edges.end(); ++e )
+ {
+ if ( !_grid->_shapes.FindIndex( *e ))
+ continue;
+ bool isOut = false;
+ gp_Pnt p;
+ double uvw[3], f,l;
+ for ( int i = 0; i < 2 && !isOut; ++i )
+ {
+ if ( i == 0 )
+ {
+ TopoDS_Vertex v = SMESH_MesherHelper::IthVertex( 0, *e );
+ p = BRep_Tool::Pnt( v );
+ }
+ else if ( eMidOut == edges.end() )
+ {
+ TopLoc_Location loc;
+ Handle(Geom_Curve) c = BRep_Tool::Curve( *e, loc, f, l);
+ if ( c.IsNull() ) break;
+ p = c->Value( 0.5 * ( f + l )).Transformed( loc );
+ }
+ else
+ {
+ continue;
+ }
+
+ _grid->ComputeUVW( p.XYZ(), uvw );
+ if ( isOutParam( uvw ))
+ {
+ if ( i == 0 )
+ isOut = true;
+ else
+ eMidOut = e;
+ }
+ }
+ if ( isOut )
+ break;
+ }
+ if ( e != edges.end() )
+ edges.splice( edges.end(), edges, edges.begin(), e );
+ else if ( eMidOut != edges.end() )
+ edges.splice( edges.end(), edges, edges.begin(), eMidOut );
+
+ // sort nodes according to the order of edges
+ _Node* orderNodes [20];
+ //TGeomID orderShapeIDs[20];
+ size_t nbN = 0;
+ TGeomID id, *pID = 0;
+ for ( e = edges.begin(); e != edges.end(); ++e )
+ {
+ if (( id = _grid->_shapes.FindIndex( SMESH_MesherHelper::IthVertex( 0, *e ))) &&
+ (( pID = std::find( &nShapeIds[0], nShapeIdsEnd, id )) != nShapeIdsEnd ))
+ {
+ //orderShapeIDs[ nbN ] = id;
+ orderNodes [ nbN++ ] = nodes[ pID - &nShapeIds[0] ];
+ *pID = -1;
+ }
+ if (( id = _grid->_shapes.FindIndex( *e )) &&
+ (( pID = std::find( &nShapeIds[0], nShapeIdsEnd, id )) != nShapeIdsEnd ))
+ {
+ //orderShapeIDs[ nbN ] = id;
+ orderNodes [ nbN++ ] = nodes[ pID - &nShapeIds[0] ];
+ *pID = -1;
+ }
+ }
+ if ( nbN != nodes.size() )
+ return;
+
+ bool reverse = ( orderNodes[0 ]->Point().SquareDistance( curNode->Point() ) >
+ orderNodes[nbN-1]->Point().SquareDistance( curNode->Point() ));
+
+ for ( size_t i = 0; i < nodes.size(); ++i )
+ nodes[ i ] = orderNodes[ reverse ? nbN-1-i : i ];