+ //================================================================================
+ /*!
+ * \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._uNorm = ( _grid->_axes[ iDirY ] ^ _grid->_axes[ iDirZ ] ).Normalized();
+ // planes._vNorm = ( _grid->_axes[ iDirZ ] ^ _grid->_axes[ iDirX ] ).Normalized();
+ 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 ( int 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];
+ //double zFactor = _grid->_axes[ iDirZ ] * planes._zNorm;
+ 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;
+ 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
+
+ // Create nodes at found intersections
+ // const E_IntersectPoint* eip;
+ // for ( size_t i = 0; i < hexes.size(); ++i )
+ // {
+ // Hexahedron* h = hexes[i];
+ // if ( !h ) continue;
+ // for ( int iF = 0; iF < 6; ++iF )
+ // {
+ // _Face& quad = h->_hexQuads[ iF ];
+ // for ( size_t iP = 0; iP < quad._edgeNodes.size(); ++iP )
+ // if ( !quad._edgeNodes[ iP ]._node )
+ // if (( eip = quad._edgeNodes[ iP ].EdgeIntPnt() ))
+ // quad._edgeNodes[ iP ]._intPoint->_node = helper.AddNode( eip->_point.X(),
+ // eip->_point.Y(),
+ // eip->_point.Z() );
+ // }
+ // for ( size_t iP = 0; iP < hexes[i]->_vertexNodes.size(); ++iP )
+ // if (( eip = h->_vertexNodes[ iP ].EdgeIntPnt() ))
+ // h->_vertexNodes[ iP ]._intPoint->_node = helper.AddNode( eip->_point.X(),
+ // eip->_point.Y(),
+ // eip->_point.Z() );
+ // }
+ }
+
+ //================================================================================
+ /*!
+ * \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 (( _grid->_coords[0][ h->_i ] - _grid->_tol > ip._uvw[0] ) ||
+ ( _grid->_coords[0][ h->_i+1 ] + _grid->_tol < ip._uvw[0] ) ||
+ ( _grid->_coords[1][ h->_j ] - _grid->_tol > ip._uvw[1] ) ||
+ ( _grid->_coords[1][ h->_j+1 ] + _grid->_tol < ip._uvw[1] ) ||
+ ( _grid->_coords[2][ h->_k ] - _grid->_tol > ip._uvw[2] ) ||
+ ( _grid->_coords[2][ h->_k+1 ] + _grid->_tol < ip._uvw[2] ))
+ throw SALOME_Exception("ip outside a hex");
+#endif
+ h->_edgeIntPnts.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._edgeNodes.size(); ++iP )
+ if ( !quad._edgeNodes[ iP ]._isUsedInFace &&
+ n1->IsLinked( quad._edgeNodes[ iP ]._intPoint ) &&
+ n2->IsLinked( quad._edgeNodes[ iP ]._intPoint ))
+ {
+ chn.push_back( & quad._edgeNodes[ iP ]);
+ chn.push_back( n2 );
+ quad._edgeNodes[ iP ]._isUsedInFace = true;
+ return true;
+ }
+ bool found;
+ do
+ {
+ found = false;
+ for ( size_t iP = 0; iP < quad._edgeNodes.size(); ++iP )
+ if ( !quad._edgeNodes[ iP ]._isUsedInFace &&
+ chn.back()->IsLinked( quad._edgeNodes[ iP ]._intPoint ))
+ {
+ chn.push_back( & quad._edgeNodes[ iP ]);
+ found = quad._edgeNodes[ iP ]._isUsedInFace = true;
+ 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 Checks transition at the 1st node of a link
+ */
+ bool Hexahedron::is1stNodeOut( _Link& link /*int iLink*/ ) const
+ {
+ // new version is for the case: tangent transition at the 1st node
+ bool isOut = false;
+ if ( link._intNodes.size() > 1 )
+ {
+ // check transition at the next intersection
+ switch ( link._intNodes[1].FaceIntPnt()->_transition ) {
+ case Trans_OUT: return false;
+ case Trans_IN : return true;
+ default: ; // tangent transition
+ }
+ }
+ gp_Pnt p1 = link._nodes[0]->Point();
+ gp_Pnt p2 = link._nodes[1]->Point();
+ gp_Pnt testPnt = 0.8 * p1.XYZ() + 0.2 * p2.XYZ();
+
+ TGeomID faceID = link._intNodes[0]._intPoint->_faceIDs[0];
+ const TopoDS_Face& face = TopoDS::Face( _grid->_shapes( faceID ));
+ TopLoc_Location loc;
+ GeomAPI_ProjectPointOnSurf& proj =
+ _grid->_helper->GetProjector( face, loc, 0.1*_grid->_tol );
+ testPnt.Transform( loc );
+ proj.Perform( testPnt );
+ if ( proj.IsDone() &&
+ proj.NbPoints() > 0 &&
+ proj.LowerDistance() > _grid->_tol )
+ {
+ Quantity_Parameter u,v;
+ proj.LowerDistanceParameters( u,v );
+ 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 );
+ return v * normal > 0;
+ }
+ }
+ // if ( !_hexLinks[ iLink ]._nodes[0]->Node() ) // no node
+ // return true;
+ // if ( !_hexLinks[ iLink ]._nodes[0]->_intPoint ) // no intersection with geometry
+ // return false;
+ // switch ( _hexLinks[ iLink ]._nodes[0]->FaceIntPnt()->_transition ) {
+ // case Trans_OUT: return true;
+ // case Trans_IN : return false;
+ // default: ; // tangent transition
+ // }
+
+// // ijk of a GridLine corresponding to the link
+// int iDir = iLink / 4;
+// int indSub = iLink % 4;
+// LineIndexer li = _grid->GetLineIndexer( iDir );
+// li.SetIJK( _i,_j,_k );
+// size_t lineIndex[4] = { li.LineIndex (),
+// li.LineIndex10(),
+// li.LineIndex01(),
+// li.LineIndex11() };
+// GridLine& line = _grid->_lines[ iDir ][ lineIndex[ indSub ]];
+
+// // analyze transition of previous ip
+// bool isOut = true;
+// multiset< F_IntersectPoint >::const_iterator ip = line._intPoints.begin();
+// for ( ; ip != line._intPoints.end(); ++ip )
+// {
+// if ( &(*ip) == _hexLinks[ iLink ]._nodes[0]->_intPoint )
+// break;
+// switch ( ip->_transition ) {
+// case Trans_OUT: isOut = true;
+// case Trans_IN : isOut = false;
+// default:;
+// }
+// }
+// #ifdef _DEBUG_
+// if ( ip == line._intPoints.end() )
+// cout << "BUG: Wrong GridLine. IKJ = ( "<< _i << " " << _j << " " << _k << " )" << endl;
+// #endif
+ return isOut;
+ }