+ //================================================================================
+ /*!
+ * \brief Implements geom edges into the mesh
+ */
+ void Hexahedron::addEdges(SMESH_MesherHelper& helper,
+ vector< Hexahedron* >& hexes,
+ const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap)
+ {
+ if ( edge2faceIDsMap.empty() ) return;
+
+ // Prepare planes for intersecting with EDGEs
+ GridPlanes pln[3];
+ {
+ for ( int iDirZ = 0; iDirZ < 3; ++iDirZ ) // iDirZ gives normal direction to planes
+ {
+ GridPlanes& planes = pln[ iDirZ ];
+ int iDirX = ( iDirZ + 1 ) % 3;
+ int iDirY = ( iDirZ + 2 ) % 3;
+ planes._zNorm = ( _grid->_axes[ iDirX ] ^ _grid->_axes[ iDirY ] ).Normalized();
+ planes._zProjs.resize ( _grid->_coords[ iDirZ ].size() );
+ planes._zProjs [0] = 0;
+ const double zFactor = _grid->_axes[ iDirZ ] * planes._zNorm;
+ const vector< double > & u = _grid->_coords[ iDirZ ];
+ for ( size_t i = 1; i < planes._zProjs.size(); ++i )
+ {
+ planes._zProjs [i] = zFactor * ( u[i] - u[0] );
+ }
+ }
+ }
+ const double deflection = _grid->_minCellSize / 20.;
+ const double tol = _grid->_tol;
+ E_IntersectPoint ip;
+
+ // Intersect EDGEs with the planes
+ map< TGeomID, vector< TGeomID > >::const_iterator e2fIt = edge2faceIDsMap.begin();
+ for ( ; e2fIt != edge2faceIDsMap.end(); ++e2fIt )
+ {
+ const TGeomID edgeID = e2fIt->first;
+ const TopoDS_Edge & E = TopoDS::Edge( _grid->_shapes( edgeID ));
+ BRepAdaptor_Curve curve( E );
+ TopoDS_Vertex v1 = helper.IthVertex( 0, E, false );
+ TopoDS_Vertex v2 = helper.IthVertex( 1, E, false );
+
+ ip._faceIDs = e2fIt->second;
+ ip._shapeID = edgeID;
+
+ // discretize the EGDE
+ GCPnts_UniformDeflection discret( curve, deflection, true );
+ if ( !discret.IsDone() || discret.NbPoints() < 2 )
+ continue;
+
+ // perform intersection
+ for ( int iDirZ = 0; iDirZ < 3; ++iDirZ )
+ {
+ GridPlanes& planes = pln[ iDirZ ];
+ int iDirX = ( iDirZ + 1 ) % 3;
+ int iDirY = ( iDirZ + 2 ) % 3;
+ double xLen = _grid->_coords[ iDirX ].back() - _grid->_coords[ iDirX ][0];
+ double yLen = _grid->_coords[ iDirY ].back() - _grid->_coords[ iDirY ][0];
+ double zLen = _grid->_coords[ iDirZ ].back() - _grid->_coords[ iDirZ ][0];
+ int dIJK[3], d000[3] = { 0,0,0 };
+ double o[3] = { _grid->_coords[0][0],
+ _grid->_coords[1][0],
+ _grid->_coords[2][0] };
+
+ // locate the 1st point of a segment within the grid
+ gp_XYZ p1 = discret.Value( 1 ).XYZ();
+ double u1 = discret.Parameter( 1 );
+ double zProj1 = planes._zNorm * ( p1 - _grid->_origin );
+
+ _grid->ComputeUVW( p1, ip._uvw );
+ int iX1 = int(( ip._uvw[iDirX] - o[iDirX]) / xLen * (_grid->_coords[ iDirX ].size() - 1));
+ int iY1 = int(( ip._uvw[iDirY] - o[iDirY]) / yLen * (_grid->_coords[ iDirY ].size() - 1));
+ int iZ1 = int(( ip._uvw[iDirZ] - o[iDirZ]) / zLen * (_grid->_coords[ iDirZ ].size() - 1));
+ locateValue( iX1, ip._uvw[iDirX], _grid->_coords[ iDirX ], dIJK[ iDirX ], tol );
+ locateValue( iY1, ip._uvw[iDirY], _grid->_coords[ iDirY ], dIJK[ iDirY ], tol );
+ locateValue( iZ1, ip._uvw[iDirZ], _grid->_coords[ iDirZ ], dIJK[ iDirZ ], tol );
+
+ int ijk[3]; // grid index where a segment intersect a plane
+ ijk[ iDirX ] = iX1;
+ ijk[ iDirY ] = iY1;
+ ijk[ iDirZ ] = iZ1;
+
+ // add the 1st vertex point to a hexahedron
+ if ( iDirZ == 0 )
+ {
+ ip._point = p1;
+ ip._shapeID = _grid->_shapes.Add( v1 );
+ _grid->_edgeIntP.push_back( ip );
+ if ( !addIntersection( _grid->_edgeIntP.back(), hexes, ijk, d000 ))
+ _grid->_edgeIntP.pop_back();
+ ip._shapeID = edgeID;
+ }
+ for ( int iP = 2; iP <= discret.NbPoints(); ++iP )
+ {
+ // locate the 2nd point of a segment within the grid
+ gp_XYZ p2 = discret.Value( iP ).XYZ();
+ double u2 = discret.Parameter( iP );
+ double zProj2 = planes._zNorm * ( p2 - _grid->_origin );
+ int iZ2 = iZ1;
+ if ( Abs( zProj2 - zProj1 ) > std::numeric_limits<double>::min() )
+ {
+ locateValue( iZ2, zProj2, planes._zProjs, dIJK[ iDirZ ], tol );
+
+ // treat intersections with planes between 2 end points of a segment
+ int dZ = ( iZ1 <= iZ2 ) ? +1 : -1;
+ int iZ = iZ1 + ( iZ1 < iZ2 );
+ for ( int i = 0, nb = Abs( iZ1 - iZ2 ); i < nb; ++i, iZ += dZ )
+ {
+ ip._point = findIntPoint( u1, zProj1, u2, zProj2,
+ planes._zProjs[ iZ ],
+ curve, planes._zNorm, _grid->_origin );
+ _grid->ComputeUVW( ip._point.XYZ(), ip._uvw );
+ locateValue( ijk[iDirX], ip._uvw[iDirX], _grid->_coords[iDirX], dIJK[iDirX], tol );
+ locateValue( ijk[iDirY], ip._uvw[iDirY], _grid->_coords[iDirY], dIJK[iDirY], tol );
+ ijk[ iDirZ ] = iZ;
+
+ // add ip to hex "above" the plane
+ _grid->_edgeIntP.push_back( ip );
+ dIJK[ iDirZ ] = 0;
+ bool added = addIntersection(_grid->_edgeIntP.back(), hexes, ijk, dIJK);
+
+ // add ip to hex "below" the plane
+ ijk[ iDirZ ] = iZ-1;
+ if ( !addIntersection( _grid->_edgeIntP.back(), hexes, ijk, dIJK ) &&
+ !added)
+ _grid->_edgeIntP.pop_back();
+ }
+ }
+ iZ1 = iZ2;
+ p1 = p2;
+ u1 = u2;
+ zProj1 = zProj2;
+ }
+ // add the 2nd vertex point to a hexahedron
+ if ( iDirZ == 0 )
+ {
+ ip._shapeID = _grid->_shapes.Add( v2 );
+ ip._point = p1;
+ _grid->ComputeUVW( p1, ip._uvw );
+ locateValue( ijk[iDirX], ip._uvw[iDirX], _grid->_coords[iDirX], dIJK[iDirX], tol );
+ locateValue( ijk[iDirY], ip._uvw[iDirY], _grid->_coords[iDirY], dIJK[iDirY], tol );
+ ijk[ iDirZ ] = iZ1;
+ _grid->_edgeIntP.push_back( ip );
+ if ( !addIntersection( _grid->_edgeIntP.back(), hexes, ijk, d000 ))
+ _grid->_edgeIntP.pop_back();
+ ip._shapeID = edgeID;
+ }
+ } // loop on 3 grid directions
+ } // loop on EDGEs
+
+ }
+
+ //================================================================================
+ /*!
+ * \brief Finds intersection of a curve with a plane
+ * \param [in] u1 - parameter of one curve point
+ * \param [in] proj1 - projection of the curve point to the plane normal
+ * \param [in] u2 - parameter of another curve point
+ * \param [in] proj2 - projection of the other curve point to the plane normal
+ * \param [in] proj - projection of a point where the curve intersects the plane
+ * \param [in] curve - the curve
+ * \param [in] axis - the plane normal
+ * \param [in] origin - the plane origin
+ * \return gp_Pnt - the found intersection point
+ */
+ gp_Pnt Hexahedron::findIntPoint( double u1, double proj1,
+ double u2, double proj2,
+ double proj,
+ BRepAdaptor_Curve& curve,
+ const gp_XYZ& axis,
+ const gp_XYZ& origin)
+ {
+ double r = (( proj - proj1 ) / ( proj2 - proj1 ));
+ double u = u1 * ( 1 - r ) + u2 * r;
+ gp_Pnt p = curve.Value( u );
+ double newProj = axis * ( p.XYZ() - origin );
+ if ( Abs( proj - newProj ) > _grid->_tol / 10. )
+ {
+ if ( r > 0.5 )
+ return findIntPoint( u2, proj2, u, newProj, proj, curve, axis, origin );
+ else
+ return findIntPoint( u1, proj2, u, newProj, proj, curve, axis, origin );
+ }
+ return p;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Returns indices of a hexahedron sub-entities holding a point
+ * \param [in] ip - intersection point
+ * \param [out] facets - 0-3 facets holding a point
+ * \param [out] sub - index of a vertex or an edge holding a point
+ * \return int - number of facets holding a point
+ */
+ int Hexahedron::getEntity( const E_IntersectPoint* ip, int* facets, int& sub )
+ {
+ enum { X = 1, Y = 2, Z = 4 }; // == 001, 010, 100
+ int nbFacets = 0;
+ int vertex = 0, egdeMask = 0;
+
+ if ( Abs( _grid->_coords[0][ _i ] - ip->_uvw[0] ) < _grid->_tol ) {
+ facets[ nbFacets++ ] = SMESH_Block::ID_F0yz;
+ egdeMask |= X;
+ }
+ else if ( Abs( _grid->_coords[0][ _i+1 ] - ip->_uvw[0] ) < _grid->_tol ) {
+ facets[ nbFacets++ ] = SMESH_Block::ID_F1yz;
+ vertex |= X;
+ egdeMask |= X;
+ }
+ if ( Abs( _grid->_coords[1][ _j ] - ip->_uvw[1] ) < _grid->_tol ) {
+ facets[ nbFacets++ ] = SMESH_Block::ID_Fx0z;
+ egdeMask |= Y;
+ }
+ else if ( Abs( _grid->_coords[1][ _j+1 ] - ip->_uvw[1] ) < _grid->_tol ) {
+ facets[ nbFacets++ ] = SMESH_Block::ID_Fx1z;
+ vertex |= Y;
+ egdeMask |= Y;
+ }
+ if ( Abs( _grid->_coords[2][ _k ] - ip->_uvw[2] ) < _grid->_tol ) {
+ facets[ nbFacets++ ] = SMESH_Block::ID_Fxy0;
+ egdeMask |= Z;
+ }
+ else if ( Abs( _grid->_coords[2][ _k+1 ] - ip->_uvw[2] ) < _grid->_tol ) {
+ facets[ nbFacets++ ] = SMESH_Block::ID_Fxy1;
+ vertex |= Z;
+ egdeMask |= Z;
+ }
+
+ switch ( nbFacets )
+ {
+ case 0: sub = 0; break;
+ case 1: sub = facets[0]; break;
+ case 2: {
+ const int edge [3][8] = {
+ { SMESH_Block::ID_E00z, SMESH_Block::ID_E10z,
+ SMESH_Block::ID_E01z, SMESH_Block::ID_E11z },
+ { SMESH_Block::ID_E0y0, SMESH_Block::ID_E1y0, 0, 0,
+ SMESH_Block::ID_E0y1, SMESH_Block::ID_E1y1 },
+ { SMESH_Block::ID_Ex00, 0, SMESH_Block::ID_Ex10, 0,
+ SMESH_Block::ID_Ex01, 0, SMESH_Block::ID_Ex11 }
+ };
+ switch ( egdeMask ) {
+ case X | Y: sub = edge[ 0 ][ vertex ]; break;
+ case X | Z: sub = edge[ 1 ][ vertex ]; break;
+ default: sub = edge[ 2 ][ vertex ];
+ }
+ break;
+ }
+ //case 3:
+ default:
+ sub = vertex + SMESH_Block::ID_FirstV;
+ }
+
+ return nbFacets;
+ }
+ //================================================================================
+ /*!
+ * \brief Adds intersection with an EDGE
+ */
+ bool Hexahedron::addIntersection( const E_IntersectPoint& ip,
+ vector< Hexahedron* >& hexes,
+ int ijk[], int dIJK[] )
+ {
+ bool added = false;
+
+ size_t hexIndex[4] = {
+ _grid->CellIndex( ijk[0], ijk[1], ijk[2] ),
+ dIJK[0] ? _grid->CellIndex( ijk[0]+dIJK[0], ijk[1], ijk[2] ) : -1,
+ dIJK[1] ? _grid->CellIndex( ijk[0], ijk[1]+dIJK[1], ijk[2] ) : -1,
+ dIJK[2] ? _grid->CellIndex( ijk[0], ijk[1], ijk[2]+dIJK[2] ) : -1
+ };
+ for ( int i = 0; i < 4; ++i )
+ {
+ if ( /*0 <= hexIndex[i] &&*/ hexIndex[i] < hexes.size() && hexes[ hexIndex[i] ] )
+ {
+ Hexahedron* h = hexes[ hexIndex[i] ];
+ // check if ip is really inside the hex
+#ifdef _DEBUG_
+ if ( h->isOutParam( ip._uvw ))
+ throw SALOME_Exception("ip outside a hex");
+#endif
+ h->_eIntPoints.push_back( & ip );
+ added = true;
+ }
+ }
+ return added;
+ }
+ //================================================================================
+ /*!
+ * \brief Finds nodes at a path from one node to another via intersections with EDGEs
+ */
+ bool Hexahedron::findChain( _Node* n1,
+ _Node* n2,
+ _Face& quad,
+ vector<_Node*>& chn )
+ {
+ chn.clear();
+ chn.push_back( n1 );
+ for ( size_t iP = 0; iP < quad._eIntNodes.size(); ++iP )
+ if ( !quad._eIntNodes[ iP ]->IsUsedInFace( &quad ) &&
+ n1->IsLinked( quad._eIntNodes[ iP ]->_intPoint ) &&
+ n2->IsLinked( quad._eIntNodes[ iP ]->_intPoint ))
+ {
+ chn.push_back( quad._eIntNodes[ iP ]);
+ chn.push_back( n2 );
+ quad._eIntNodes[ iP ]->_usedInFace = &quad;
+ return true;
+ }
+ 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 ))
+ {
+ chn.push_back( quad._eIntNodes[ iP ]);
+ found = ( quad._eIntNodes[ iP ]->_usedInFace = &quad );
+ break;
+ }
+ } while ( found && ! chn.back()->IsLinked( n2->_intPoint ) );
+
+ if ( chn.back() != n2 && chn.back()->IsLinked( n2->_intPoint ))
+ chn.push_back( n2 );
+
+ 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 )
+ {
+ Quantity_Parameter 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 accoring 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 ];
+ }
+