From: eap Date: Wed, 25 Dec 2013 15:46:26 +0000 (+0000) Subject: 22360: EDF SMESH: Body Fitting algorithm: incorporate edges X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=60e380a37ac09f6f5c147c146fff1acf5bec3307;p=modules%2Fsmesh.git 22360: EDF SMESH: Body Fitting algorithm: incorporate edges --- diff --git a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx index b94506f68..50ace723e 100644 --- a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx +++ b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx @@ -459,7 +459,14 @@ namespace bool IsLinked( const B_IntersectPoint* other ) const { // node inside a SOLID is considered "linked" with any other node - return ( !other || !_intPoint ) ? true : _intPoint->HasCommonFace( other ); + if ( !other || !_intPoint || _intPoint->HasCommonFace( other )) + return true; + const F_IntersectPoint* ip1 = dynamic_cast< const F_IntersectPoint* >( _intPoint ); + const F_IntersectPoint* ip2 = dynamic_cast< const F_IntersectPoint* >( other ); + if ( ip1 && ip2 ) + return (( ip1->_transition != ip2->_transition ) && + ( ip1->_transition + ip2->_transition ) == Trans_IN + Trans_OUT ); + return false; } bool IsOnFace( int faceID ) const // returns true if faceID is found { @@ -591,8 +598,13 @@ namespace void addEdges(SMESH_MesherHelper& helper, vector< Hexahedron* >& intersectedHex, const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap); + gp_Pnt findIntPoint( double u1, double proj1, double u2, double proj2, + double proj, BRepAdaptor_Curve& curve, + const gp_XYZ& axis, const gp_XYZ& origin ); int getEntity( const E_IntersectPoint* ip, int* facets, int& sub ); - void addIntersection( const E_IntersectPoint& ip ); + bool addIntersection( const E_IntersectPoint& ip, + vector< Hexahedron* >& hexes, + int ijk[], int dIJK[] ); bool findChain( _Node* n1, _Node* n2, _Face& quad, vector<_Node*>& chainNodes ); int addElements(SMESH_MesherHelper& helper); bool isInHole() const; @@ -651,7 +663,8 @@ namespace /*! * \brief adjust \a i to have \a val between values[i] and values[i+1] */ - inline void locateValue( int & i, double val, const vector& values ) + inline void locateValue( int & i, double val, const vector& values, + int& di, double tol ) { val += values[0]; // input \a val is measured from 0. if ( i > values.size()-2 ) @@ -661,6 +674,13 @@ namespace ++i; while ( i > 0 && val < values[ i ]) --i; + + if ( i > 0 && val - values[ i ] < tol ) + di = -1; + else if ( i+2 < values.size() && values[ i+1 ] - val < tol ) + di = 1; + else + di = 0; } //============================================================================= /* @@ -1636,8 +1656,6 @@ namespace } // switch( nbFacets ) } // loop on _edgeIntPnts - - //SMESHUtils::FreeVector( _edgeIntPnts ); } } //================================================================================ @@ -1679,14 +1697,14 @@ namespace _Link polyLink; - bool hasEdgeIntersections = !_vertexNodes.empty(); + bool hasEdgeIntersections = !_edgeIntPnts.empty(); for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron { _Face& quad = _hexQuads[ iF ] ; - if ( !quad._edgeNodes.empty() ) - hasEdgeIntersections = true; + // if ( !quad._edgeNodes.empty() ) + // hasEdgeIntersections = true; _polygons.resize( _polygons.size() + 1 ); _Face& polygon = _polygons.back(); @@ -1734,7 +1752,7 @@ namespace polygon._links.push_back( split ); } } - if ( polygon._links.size() > 1 ) + if ( !polygon._links.empty() && polygon._links.size() + quad._edgeNodes.size() > 1 ) { _Node* n1 = polygon._links.back().LastNode(); _Node* n2 = polygon._links.front().FirstNode(); @@ -1749,6 +1767,9 @@ namespace polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() )); } } + } + if ( polygon._links.size() >= 3 ) + { // add polygon to its links for ( size_t iL = 0; iL < polygon._links.size(); ++iL ) { @@ -2139,6 +2160,7 @@ namespace } } const double deflection = _grid->_minCellSize / 20.; + const double tol = _grid->_tol; // int facets[6] = { SMESH_Block::ID_F0yz, SMESH_Block::ID_F1yz, // SMESH_Block::ID_Fx0z, SMESH_Block::ID_Fx1z, // SMESH_Block::ID_Fxy0, SMESH_Block::ID_Fxy1 }; @@ -2170,7 +2192,7 @@ namespace double xLen = _grid->_coords[ iDirX ].back() - _grid->_coords[ iDirX ][0]; double yLen = _grid->_coords[ iDirY ].back() - _grid->_coords[ iDirY ][0]; double zFactor = _grid->_axes[ iDirZ ] * planes._zNorm; - //int* zFacets = facets + iDirZ * 2; + int dIJK[3], d000[3] = { 0,0,0 }; // locate the 1st point of a segment within the grid gp_XYZ p1 = discret.Value( 1 ).XYZ(); @@ -2181,9 +2203,9 @@ namespace int iX1 = int( uv.X() / xLen * ( _grid->_coords[ iDirX ].size() - 1. )); int iY1 = int( uv.Y() / yLen * ( _grid->_coords[ iDirY ].size() - 1. )); int iZ1 = int( zProj1 / planes._zProjs.back() * ( planes._zProjs.size() - 1. )); - locateValue( iX1, uv.X(), _grid->_coords[ iDirX ] ); - locateValue( iY1, uv.Y(), _grid->_coords[ iDirY ] ); - locateValue( iZ1, zProj1, planes._zProjs ); + locateValue( iX1, uv.X(), _grid->_coords[ iDirX ], dIJK[ iDirX ], tol ); + locateValue( iY1, uv.Y(), _grid->_coords[ iDirY ], dIJK[ iDirY ], tol ); + locateValue( iZ1, zProj1, planes._zProjs , dIJK[ iDirZ ], tol ); int ijk[3]; // grid index where a segment intersect a plane ijk[ iDirX ] = iX1; @@ -2194,15 +2216,14 @@ namespace ip._uvw[ iDirZ ] = zProj1 / zFactor + _grid->_coords[ iDirZ ][0]; // add the 1st vertex point to a hexahedron - size_t hexIndex = _grid->CellIndex( ijk[0], ijk[1], ijk[2] ); - if ( iDirZ == 0 && hexIndex < hexes.size() && hexes[ hexIndex ] ) + if ( iDirZ == 0 ) { //ip._shapeID = _grid->_shapes.Add( helper.IthVertex( 0, curve.Edge(),/*CumOri=*/false)); ip._point = p1; _grid->_edgeIntP.push_back( ip ); - hexes[ hexIndex ]->addIntersection( _grid->_edgeIntP.back() ); + if ( !addIntersection( _grid->_edgeIntP.back(), hexes, ijk, d000 )) + _grid->_edgeIntP.pop_back(); } - for ( int iP = 2; iP <= discret.NbPoints(); ++iP ) { // locate the 2nd point of a segment within the grid @@ -2210,35 +2231,34 @@ namespace double u2 = discret.Parameter( iP ); double zProj2 = planes._zNorm * ( p2 - planes._origins[0] ); int iZ2 = iZ1; - locateValue( iZ2, zProj2, planes._zProjs ); + locateValue( iZ2, zProj2, planes._zProjs, dIJK[ iDirZ ], tol ); - // treat intersections with planes between 2 points of a segment + // 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 ) { - double r = (( planes._zProjs[ iZ ] - zProj1 ) / ( zProj2 - zProj1 )); - ip._point = curve.Value( u1 * ( 1 - r ) + u2 * r ); + ip._point = findIntPoint( u1, zProj1, u2, zProj2, + planes._zProjs[ iZ ], + curve, planes._zNorm, planes._origins[0] ); gp_XY uv = planes.GetUV( ip._point, planes._origins[ iZ ]); - locateValue( ijk[ iDirX ], uv.X(), _grid->_coords[ iDirX ] ); - locateValue( ijk[ iDirY ], uv.Y(), _grid->_coords[ iDirY ] ); + locateValue( ijk[ iDirX ], uv.X(), _grid->_coords[ iDirX ], dIJK[ iDirX ], tol ); + locateValue( ijk[ iDirY ], uv.Y(), _grid->_coords[ iDirY ], dIJK[ iDirY ], tol ); ijk[ iDirZ ] = iZ; ip._uvw[ iDirX ] = uv.X() + _grid->_coords[ iDirX ][0]; ip._uvw[ iDirY ] = uv.Y() + _grid->_coords[ iDirY ][0]; ip._uvw[ iDirZ ] = planes._zProjs[ iZ ] / zFactor + _grid->_coords[ iDirZ ][0]; - _grid->_edgeIntP.push_back( ip ); - // add ip to hex "above" the plane - hexIndex = _grid->CellIndex( ijk[0], ijk[1], ijk[2] ); - if ( hexIndex < hexes.size() && hexes[ hexIndex ] ) - hexes[ hexIndex ]->addIntersection( _grid->_edgeIntP.back() ); + _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; - hexIndex = _grid->CellIndex( ijk[0], ijk[1], ijk[2] ); - if ( hexIndex < hexes.size() && hexes[ hexIndex ] ) - hexes[ hexIndex ]->addIntersection( _grid->_edgeIntP.back() ); + if ( !addIntersection( _grid->_edgeIntP.back(), hexes, ijk, dIJK ) && + !added) + _grid->_edgeIntP.pop_back(); } iZ1 = iZ2; p1 = p2; @@ -2250,19 +2270,16 @@ namespace { orig = planes._origins[0] + planes._zNorm * zProj1; uv = planes.GetUV( p1, orig ); - locateValue( ijk[ iDirX ], uv.X(), _grid->_coords[ iDirX ] ); - locateValue( ijk[ iDirY ], uv.Y(), _grid->_coords[ iDirY ] ); + locateValue( ijk[ iDirX ], uv.X(), _grid->_coords[ iDirX ], dIJK[ iDirX ], tol ); + locateValue( ijk[ iDirY ], uv.Y(), _grid->_coords[ iDirY ], dIJK[ iDirY ], tol ); ijk[ iDirZ ] = iZ1; - hexIndex = _grid->CellIndex( ijk[0], ijk[1], ijk[2] ); - if ( hexIndex < hexes.size() && hexes[ hexIndex ] ) - { - ip._uvw[ iDirX ] = uv.X() + _grid->_coords[ iDirX ][0]; - ip._uvw[ iDirY ] = uv.Y() + _grid->_coords[ iDirY ][0]; - ip._uvw[ iDirZ ] = zProj1 / zFactor + _grid->_coords[ iDirZ ][0]; - ip._point = p1; - _grid->_edgeIntP.push_back( ip ); - hexes[ hexIndex ]->addIntersection( _grid->_edgeIntP.back() ); - } + ip._uvw[ iDirX ] = uv.X() + _grid->_coords[ iDirX ][0]; + ip._uvw[ iDirY ] = uv.Y() + _grid->_coords[ iDirY ][0]; + ip._uvw[ iDirZ ] = zProj1 / zFactor + _grid->_coords[ iDirZ ][0]; + ip._point = p1; + _grid->_edgeIntP.push_back( ip ); + if ( !addIntersection( _grid->_edgeIntP.back(), hexes, ijk, d000 )) + _grid->_edgeIntP.pop_back(); } } // loop on 3 grid directions } // loop on EDGEs @@ -2291,6 +2308,42 @@ namespace // } } + //================================================================================ + /*! + * \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 index of a hexahedron sub-entities holding a point @@ -2301,7 +2354,7 @@ namespace */ int Hexahedron::getEntity( const E_IntersectPoint* ip, int* facets, int& sub ) { - enum { X = 4, Y = 2, Z = 1 }; // == 100, 010, 001 + enum { X = 1, Y = 2, Z = 4 }; // == 001, 010, 100 int nbFacets = 0; int vertex = 0, egdeMask = 0; @@ -2309,7 +2362,7 @@ namespace facets[ nbFacets++ ] = SMESH_Block::ID_F0yz; egdeMask |= X; } - if ( Abs( _grid->_coords[0][ _i+1 ] - ip->_uvw[0] ) < _grid->_tol ) { + else if ( Abs( _grid->_coords[0][ _i+1 ] - ip->_uvw[0] ) < _grid->_tol ) { facets[ nbFacets++ ] = SMESH_Block::ID_F1yz; vertex |= X; egdeMask |= X; @@ -2318,7 +2371,7 @@ namespace facets[ nbFacets++ ] = SMESH_Block::ID_Fx0z; egdeMask |= Y; } - if ( Abs( _grid->_coords[1][ _j+1 ] - ip->_uvw[1] ) < _grid->_tol ) { + else if ( Abs( _grid->_coords[1][ _j+1 ] - ip->_uvw[1] ) < _grid->_tol ) { facets[ nbFacets++ ] = SMESH_Block::ID_Fx1z; vertex |= Y; egdeMask |= Y; @@ -2327,57 +2380,75 @@ namespace facets[ nbFacets++ ] = SMESH_Block::ID_Fxy0; egdeMask |= Z; } - if ( Abs( _grid->_coords[2][ _k+1 ] - ip->_uvw[2] ) < _grid->_tol ) { + else if ( Abs( _grid->_coords[2][ _k+1 ] - ip->_uvw[2] ) < _grid->_tol ) { facets[ nbFacets++ ] = SMESH_Block::ID_Fxy1; vertex |= Z; egdeMask |= Z; } - if ( nbFacets == 3 ) - { - sub = vertex + SMESH_Block::ID_FirstV; - } - else if ( nbFacets == 2 ) + + switch ( nbFacets ) { + case 0: sub = 0; break; + case 1: sub = facets[0]; break; + case 2: { const int edge [3][8] = { - { SMESH_Block::ID_E00z, 0, SMESH_Block::ID_E01z, 0, - SMESH_Block::ID_E10z, 0, SMESH_Block::ID_E11z }, - { SMESH_Block::ID_E0y0, SMESH_Block::ID_E0y1, 0, 0, - SMESH_Block::ID_E1y0, SMESH_Block::ID_E1y1, 0, 0 }, - { SMESH_Block::ID_Ex00, SMESH_Block::ID_Ex01, - SMESH_Block::ID_Ex10, SMESH_Block::ID_Ex11, 0, 0, 0, 0 } + { 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; } - else if ( nbFacets == 1 ) - { - sub = facets[0]; - } - else // nbFacets == 0 - { + //case 3: + default: + sub = vertex + SMESH_Block::ID_FirstV; } + return nbFacets; } //================================================================================ /*! - * \brief Adds intersection with an EDGE, either to a _Face or inside a hex + * \brief Adds intersection with an EDGE */ - void Hexahedron::addIntersection( const E_IntersectPoint& ip ) + bool Hexahedron::addIntersection( const E_IntersectPoint& ip, + vector< Hexahedron* >& hexes, + int ijk[], int dIJK[] ) { - // check if ip is really inside the hex + 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][ _i ] - _grid->_tol > ip._uvw[0] ) || - ( _grid->_coords[0][ _i+1 ] + _grid->_tol < ip._uvw[0] ) || - ( _grid->_coords[1][ _j ] - _grid->_tol > ip._uvw[1] ) || - ( _grid->_coords[1][ _j+1 ] + _grid->_tol < ip._uvw[1] ) || - ( _grid->_coords[2][ _k ] - _grid->_tol > ip._uvw[2] ) || - ( _grid->_coords[2][ _k+1 ] + _grid->_tol < ip._uvw[2] )) - throw SALOME_Exception("ip outside a hex"); + 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 - _edgeIntPnts.push_back( & ip ); + h->_edgeIntPnts.push_back( & ip ); + added = true; + } + } + return added; } //================================================================================ /*!