src/SMESHUtils/SMESH_Delaunay.cxx

   1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   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.
  10 //
  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.
  15 //
  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
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 // File      : SMESH_Delaunay.cxx
  23 // Created   : Wed Apr 19 15:41:15 2017
  24 // Author    : Edward AGAPOV (eap)
  25
  26 #include "SMESH_Delaunay.hxx"
  27 #include "SMESH_MeshAlgos.hxx"
  28
  29 #include <BRepAdaptor_Surface.hxx>
  30 #include <BRepMesh_Delaun.hxx>
  31
  32 //================================================================================
  33 /*!
  34  * \brief Construct a Delaunay triangulation of given boundary nodes
  35  *  \param [in] boundaryNodes - vector of nodes of a wire
  36  *  \param [in] face - the face
  37  *  \param [in] faceID - the face ID
  38  *  \param [in] nbNodesToVisit - nb of non-marked nodes on the face
  39  */
  40 //================================================================================
  41
  42 SMESH_Delaunay::SMESH_Delaunay(const std::vector< const UVPtStructVec* > & boundaryNodes,
  43                                const TopoDS_Face&                          face,
  44                                const int                                   faceID)
  45   : _face( face ), _faceID( faceID ), _scale( 1., 1. )
  46 {
  47   // compute _scale
  48   BRepAdaptor_Surface surf( face );
  49   if ( surf.GetType() != GeomAbs_Plane )
  50   {
  51     const int nbDiv = 100;
  52     const double uRange = surf.LastUParameter() - surf.FirstUParameter();
  53     const double vRange = surf.LastVParameter() - surf.FirstVParameter();
  54     const double dU = uRange / nbDiv;
  55     const double dV = vRange / nbDiv;
  56     double u = surf.FirstUParameter(), v = surf.FirstVParameter();
  57     gp_Pnt p0U = surf.Value( u, v ), p0V = p0U;
  58     double lenU = 0, lenV = 0;
  59     for ( ; u < surf.LastUParameter(); u += dU, v += dV )
  60     {
  61       gp_Pnt p1U = surf.Value( u, surf.FirstVParameter() );
  62       lenU += p1U.Distance( p0U );
  63       p0U = p1U;
  64       gp_Pnt p1V = surf.Value( surf.FirstUParameter(), v );
  65       lenV += p1V.Distance( p0V );
  66       p0V = p1V;
  67     }
  68     _scale.SetCoord( lenU / uRange, lenV / vRange );
  69   }
  70
  71   // count boundary points
  72   int iP = 1, nbP = 0;
  73   for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW ) // loop on wires
  74   {
  75     nbP += boundaryNodes[iW]->size();
  76     if ( boundaryNodes[iW]->front().node == boundaryNodes[iW]->back().node )
  77       --nbP; // 1st and last points coincide
  78   }
  79   _bndNodes.resize( nbP );
  80
  81   // fill boundary points
  82   BRepMesh::Array1OfVertexOfDelaun bndVert( 1, 1 + nbP );
  83   BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier );
  84   for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW )
  85   {
  86     const UVPtStructVec& bndPnt = *boundaryNodes[iW];
  87     int i = 0, nb = bndPnt.size();
  88     if ( bndPnt[0].node == bndPnt.back().node )
  89       --nb;
  90     for ( ; i < nb;  ++i, ++iP )
  91     {
  92       _bndNodes[ iP-1 ] = bndPnt[i].node;
  93       bndPnt[i].node->setIsMarked( true );
  94
  95       v.ChangeCoord() = bndPnt[i].UV().Multiplied( _scale );
  96       bndVert( iP )   = v;
  97     }
  98   }
  99
 100   // triangulate the srcFace in 2D
 101   BRepMesh_Delaun Delaunay( bndVert );
 102   _triaDS = Delaunay.Result();
 103 }
 104
 105 //================================================================================
 106 /*!
 107  * \brief Prepare to the exploration of nodes
 108  */
 109 //================================================================================
 110
 111 void SMESH_Delaunay::InitTraversal(const int nbNodesToVisit)
 112 {
 113   _nbNodesToVisit = (size_t) nbNodesToVisit;
 114   _nbVisitedNodes = _iBndNode = 0;
 115   _noTriQueue.clear();
 116 }
 117
 118 //================================================================================
 119 /*!
 120  * \brief Return a node with its Barycentric Coordinates within the triangle
 121  *        defined by its node indices (zero based)
 122  *  \param [out] bc - Barycentric Coordinates of the returned node
 123  *  \param [out] triaNodes - indices of triangle nodes
 124  *  \return const SMDS_MeshNode* - the next node or NULL
 125  */
 126 //================================================================================
 127
 128 const SMDS_MeshNode* SMESH_Delaunay::NextNode( double bc[3], int triaNodes[3] )
 129 {
 130   while ( _nbVisitedNodes < _nbNodesToVisit )
 131   {
 132     while ( !_noTriQueue.empty() )
 133     {
 134       const SMDS_MeshNode*     node = _noTriQueue.front().first;
 135       const BRepMesh_Triangle* tria = _noTriQueue.front().second;
 136       _noTriQueue.pop_front();
 137       if ( node->isMarked() )
 138         continue;
 139       ++_nbVisitedNodes;
 140       node->setIsMarked( true );
 141
 142       // find a Delaunay triangle containing the src node
 143       gp_XY uv = getNodeUV( _face, node );
 144       tria = FindTriangle( uv, tria, bc, triaNodes );
 145       if ( tria )
 146       {
 147         addCloseNodes( node, tria, _faceID, _noTriQueue );
 148         return node;
 149       }
 150     }
 151     for ( ; _iBndNode < _bndNodes.size() &&  _noTriQueue.empty();  ++_iBndNode )
 152     {
 153       if ( const BRepMesh_Triangle* tria = GetTriangleNear( _iBndNode ))
 154         addCloseNodes( _bndNodes[ _iBndNode ], tria, _faceID, _noTriQueue );
 155     }
 156     if ( _noTriQueue.empty() )
 157       break;
 158   }
 159
 160   // if ( _nbVisitedNodes < _nbNodesToVisit )
 161   //   _nbVisitedNodes = std::numeric_limits<int>::max();
 162   return NULL;
 163 }
 164
 165 //================================================================================
 166 /*!
 167  * \brief Find a Delaunay triangle containing a given 2D point and return
 168  *        barycentric coordinates within the found triangle
 169  */
 170 //================================================================================
 171
 172 const BRepMesh_Triangle* SMESH_Delaunay::FindTriangle( const gp_XY&             UV,
 173                                                        const BRepMesh_Triangle* tria,
 174                                                        double                   bc[3],
 175                                                        int                      triaNodes[3] )
 176 {
 177   int   nodeIDs[3];
 178   gp_XY nodeUVs[3];
 179   int   linkIDs[3];
 180   Standard_Boolean ori[3];
 181
 182   gp_XY uv = UV.Multiplied( _scale );
 183
 184   while ( tria )
 185   {
 186     // check if the uv is in tria
 187
 188     _triaDS->ElementNodes( *tria, nodeIDs );
 189     nodeUVs[0] = _triaDS->GetNode( nodeIDs[0] ).Coord();
 190     nodeUVs[1] = _triaDS->GetNode( nodeIDs[1] ).Coord();
 191     nodeUVs[2] = _triaDS->GetNode( nodeIDs[2] ).Coord();
 192
 193     SMESH_MeshAlgos::GetBarycentricCoords( uv,
 194                                            nodeUVs[0], nodeUVs[1], nodeUVs[2],
 195                                            bc[0], bc[1] );
 196     if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 )
 197     {
 198       bc[2] = 1 - bc[0] - bc[1];
 199       triaNodes[0] = nodeIDs[0] - 1;
 200       triaNodes[1] = nodeIDs[1] - 1;
 201       triaNodes[2] = nodeIDs[2] - 1;
 202       return tria;
 203     }
 204
 205     // look for a neighbor triangle, which is adjacent to a link intersected
 206     // by a segment( triangle center -> uv )
 207
 208     gp_XY gc = ( nodeUVs[0] + nodeUVs[1] + nodeUVs[2] ) / 3.;
 209     gp_XY seg = uv - gc;
 210
 211     tria->Edges( linkIDs, ori );
 212     int triaID = _triaDS->IndexOf( *tria );
 213     tria = 0;
 214
 215     for ( int i = 0; i < 3; ++i )
 216     {
 217       const BRepMesh_PairOfIndex & triIDs = _triaDS->ElementsConnectedTo( linkIDs[i] );
 218       if ( triIDs.Extent() < 2 )
 219         continue; // no neighbor triangle
 220
 221       // check if a link intersects gc2uv
 222       const BRepMesh_Edge & link = _triaDS->GetLink( linkIDs[i] );
 223       const BRepMesh_Vertex & n1 = _triaDS->GetNode( link.FirstNode() );
 224       const BRepMesh_Vertex & n2 = _triaDS->GetNode( link.LastNode() );
 225       gp_XY uv1 = n1.Coord();
 226       gp_XY lin = n2.Coord() - uv1; // link direction
 227
 228       double crossSegLin = seg ^ lin;
 229       if ( Abs( crossSegLin ) < std::numeric_limits<double>::min() )
 230         continue; // parallel
 231
 232       double uSeg = ( uv1 - gc ) ^ lin / crossSegLin;
 233       if ( 0. <= uSeg && uSeg <= 1. )
 234       {
 235         tria = & _triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID )));
 236         break;
 237       }
 238     }
 239   }
 240   return tria;
 241 }
 242
 243 //================================================================================
 244 /*!
 245  * \brief Return a triangle sharing a given boundary node
 246  *  \param [in] iBndNode - index of the boundary node
 247  *  \return const BRepMesh_Triangle* - a found triangle
 248  */
 249 //================================================================================
 250
 251 const BRepMesh_Triangle* SMESH_Delaunay::GetTriangleNear( int iBndNode )
 252 {
 253   const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndNode + 1 );
 254   const BRepMesh_PairOfIndex &    triaIds = _triaDS->ElementsConnectedTo( linkIds.First() );
 255   const BRepMesh_Triangle&           tria = _triaDS->GetElement( triaIds.Index(1) );
 256   return &tria;
 257 }
 258
 259 //================================================================================
 260 /*!
 261  * \brief Return UV of the i-th source boundary node (zero based)
 262  */
 263 //================================================================================
 264
 265 gp_XY SMESH_Delaunay::GetBndUV(const int iNode) const
 266 {
 267   return _triaDS->GetNode( iNode+1 ).Coord();
 268 }
 269
 270 //================================================================================
 271 /*!
 272  * \brief Add non-marked nodes surrounding a given one to a queue
 273  */
 274 //================================================================================
 275
 276 void SMESH_Delaunay::addCloseNodes( const SMDS_MeshNode*     node,
 277                                     const BRepMesh_Triangle* tria,
 278                                     const int                faceID,
 279                                     TNodeTriaList &          _noTriQueue )
 280 {
 281   // find in-FACE nodes
 282   SMDS_ElemIteratorPtr elems = node->GetInverseElementIterator(SMDSAbs_Face);
 283   while ( elems->more() )
 284   {
 285     const SMDS_MeshElement* elem = elems->next();
 286     if ( elem->getshapeId() == faceID )
 287     {
 288       for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i )
 289       {
 290         const SMDS_MeshNode* n = elem->GetNode( i );
 291         if ( !n->isMarked() /*&& n->getshapeId() == faceID*/ )
 292           _noTriQueue.push_back( std::make_pair( n, tria ));
 293       }
 294     }
 295   }
 296 }