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
  28 #include "SMESH_Comment.hxx"
  29 #include "SMESH_File.hxx"
  30 #include "SMESH_MeshAlgos.hxx"
  31
  32 #include <BRepAdaptor_Surface.hxx>
  33 #include <BRepMesh_Delaun.hxx>
  34
  35 //================================================================================
  36 /*!
  37  * \brief Construct a Delaunay triangulation of given boundary nodes
  38  *  \param [in] boundaryNodes - vector of nodes of a wire
  39  *  \param [in] face - the face
  40  *  \param [in] faceID - the face ID
  41  *  \param [in] nbNodesToVisit - nb of non-marked nodes on the face
  42  */
  43 //================================================================================
  44
  45 SMESH_Delaunay::SMESH_Delaunay(const std::vector< const UVPtStructVec* > & boundaryNodes,
  46                                const TopoDS_Face&                          face,
  47                                const int                                   faceID)
  48   : _face( face ), _faceID( faceID ), _scale( 1., 1. )
  49 {
  50   // compute _scale
  51   BRepAdaptor_Surface surf( face );
  52   if ( surf.GetType() != GeomAbs_Plane )
  53   {
  54     const int nbDiv = 100;
  55     const double uRange = surf.LastUParameter() - surf.FirstUParameter();
  56     const double vRange = surf.LastVParameter() - surf.FirstVParameter();
  57     const double dU = uRange / nbDiv;
  58     const double dV = vRange / nbDiv;
  59     double u = surf.FirstUParameter(), v = surf.FirstVParameter();
  60     gp_Pnt p0U = surf.Value( u, v ), p0V = p0U;
  61     double lenU = 0, lenV = 0;
  62     for ( ; u < surf.LastUParameter(); u += dU, v += dV )
  63     {
  64       gp_Pnt p1U = surf.Value( u, surf.FirstVParameter() );
  65       lenU += p1U.Distance( p0U );
  66       p0U = p1U;
  67       gp_Pnt p1V = surf.Value( surf.FirstUParameter(), v );
  68       lenV += p1V.Distance( p0V );
  69       p0V = p1V;
  70     }
  71     _scale.SetCoord( lenU / uRange, lenV / vRange );
  72   }
  73
  74   // count boundary points
  75   int iP = 1, nbP = 0;
  76   for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW ) // loop on wires
  77   {
  78     nbP += boundaryNodes[iW]->size();
  79     if ( boundaryNodes[iW]->front().node == boundaryNodes[iW]->back().node )
  80       --nbP; // 1st and last points coincide
  81   }
  82   _bndNodes.resize( nbP );
  83
  84   // fill boundary points
  85   BRepMesh::Array1OfVertexOfDelaun bndVert( 1, 1 + nbP );
  86   BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier );
  87   for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW )
  88   {
  89     const UVPtStructVec& bndPnt = *boundaryNodes[iW];
  90     int i = 0, nb = bndPnt.size();
  91     if ( bndPnt[0].node == bndPnt.back().node )
  92       --nb;
  93     for ( ; i < nb;  ++i, ++iP )
  94     {
  95       _bndNodes[ iP-1 ] = bndPnt[i].node;
  96       bndPnt[i].node->setIsMarked( true );
  97
  98       v.ChangeCoord() = bndPnt[i].UV().Multiplied( _scale );
  99       bndVert( iP )   = v;
 100     }
 101   }
 102
 103   // triangulate the srcFace in 2D
 104   BRepMesh_Delaun Delaunay( bndVert );
 105   _triaDS = Delaunay.Result();
 106 }
 107
 108 //================================================================================
 109 /*!
 110  * \brief Prepare to the exploration of nodes
 111  */
 112 //================================================================================
 113
 114 void SMESH_Delaunay::InitTraversal(const int nbNodesToVisit)
 115 {
 116   _nbNodesToVisit = (size_t) nbNodesToVisit;
 117   _nbVisitedNodes = _iBndNode = 0;
 118   _noTriQueue.clear();
 119 }
 120
 121 //================================================================================
 122 /*!
 123  * \brief Return a node with its Barycentric Coordinates within the triangle
 124  *        defined by its node indices (zero based)
 125  *  \param [out] bc - Barycentric Coordinates of the returned node
 126  *  \param [out] triaNodes - indices of triangle nodes
 127  *  \return const SMDS_MeshNode* - the next node or NULL
 128  */
 129 //================================================================================
 130
 131 const SMDS_MeshNode* SMESH_Delaunay::NextNode( double bc[3], int triaNodes[3] )
 132 {
 133   while ( _nbVisitedNodes < _nbNodesToVisit )
 134   {
 135     while ( !_noTriQueue.empty() )
 136     {
 137       const SMDS_MeshNode*     node = _noTriQueue.front().first;
 138       const BRepMesh_Triangle* tria = _noTriQueue.front().second;
 139       _noTriQueue.pop_front();
 140       if ( node->isMarked() )
 141         continue;
 142       ++_nbVisitedNodes;
 143       node->setIsMarked( true );
 144
 145       // find a Delaunay triangle containing the src node
 146       gp_XY uv = getNodeUV( _face, node );
 147       tria = FindTriangle( uv, tria, bc, triaNodes );
 148       if ( tria )
 149       {
 150         addCloseNodes( node, tria, _faceID, _noTriQueue );
 151         return node;
 152       }
 153     }
 154     for ( ; _iBndNode < _bndNodes.size() &&  _noTriQueue.empty();  ++_iBndNode )
 155     {
 156       if ( const BRepMesh_Triangle* tria = GetTriangleNear( _iBndNode ))
 157         addCloseNodes( _bndNodes[ _iBndNode ], tria, _faceID, _noTriQueue );
 158     }
 159     if ( _noTriQueue.empty() )
 160       break;
 161   }
 162
 163   // if ( _nbVisitedNodes < _nbNodesToVisit )
 164   //   _nbVisitedNodes = std::numeric_limits<int>::max();
 165   return NULL;
 166 }
 167
 168 //================================================================================
 169 /*!
 170  * \brief Find a Delaunay triangle containing a given 2D point and return
 171  *        barycentric coordinates within the found triangle
 172  */
 173 //================================================================================
 174
 175 const BRepMesh_Triangle* SMESH_Delaunay::FindTriangle( const gp_XY&             UV,
 176                                                        const BRepMesh_Triangle* tria,
 177                                                        double                   bc[3],
 178                                                        int                      triaNodes[3] )
 179 {
 180   int   nodeIDs[3];
 181   gp_XY nodeUVs[3];
 182   int   linkIDs[3];
 183   Standard_Boolean ori[3];
 184
 185   gp_XY uv = UV.Multiplied( _scale );
 186
 187   while ( tria )
 188   {
 189     // check if the uv is in tria
 190
 191     _triaDS->ElementNodes( *tria, nodeIDs );
 192     nodeUVs[0] = _triaDS->GetNode( nodeIDs[0] ).Coord();
 193     nodeUVs[1] = _triaDS->GetNode( nodeIDs[1] ).Coord();
 194     nodeUVs[2] = _triaDS->GetNode( nodeIDs[2] ).Coord();
 195
 196     SMESH_MeshAlgos::GetBarycentricCoords( uv,
 197                                            nodeUVs[0], nodeUVs[1], nodeUVs[2],
 198                                            bc[0], bc[1] );
 199     if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 )
 200     {
 201       if ( _triaDS->GetNode( nodeIDs[0] ).Movability() != BRepMesh_Frontier ||
 202            _triaDS->GetNode( nodeIDs[1] ).Movability() != BRepMesh_Frontier ||
 203            _triaDS->GetNode( nodeIDs[2] ).Movability() != BRepMesh_Frontier )
 204       {
 205         return 0;
 206       }
 207       bc[2] = 1 - bc[0] - bc[1];
 208       triaNodes[0] = nodeIDs[0] - 1;
 209       triaNodes[1] = nodeIDs[1] - 1;
 210       triaNodes[2] = nodeIDs[2] - 1;
 211       return tria;
 212     }
 213
 214     // look for a neighbor triangle, which is adjacent to a link intersected
 215     // by a segment( triangle center -> uv )
 216
 217     gp_XY gc = ( nodeUVs[0] + nodeUVs[1] + nodeUVs[2] ) / 3.;
 218     gp_XY seg = uv - gc;
 219
 220     tria->Edges( linkIDs, ori );
 221     int triaID = _triaDS->IndexOf( *tria );
 222     tria = 0;
 223
 224     for ( int i = 0; i < 3; ++i )
 225     {
 226       const BRepMesh_PairOfIndex & triIDs = _triaDS->ElementsConnectedTo( linkIDs[i] );
 227       if ( triIDs.Extent() < 2 )
 228         continue; // no neighbor triangle
 229
 230       // check if a link intersects gc2uv
 231       const BRepMesh_Edge & link = _triaDS->GetLink( linkIDs[i] );
 232       const BRepMesh_Vertex & n1 = _triaDS->GetNode( link.FirstNode() );
 233       const BRepMesh_Vertex & n2 = _triaDS->GetNode( link.LastNode() );
 234       gp_XY uv1 = n1.Coord();
 235       gp_XY lin = n2.Coord() - uv1; // link direction
 236
 237       double crossSegLin = seg ^ lin;
 238       if ( Abs( crossSegLin ) < std::numeric_limits<double>::min() )
 239         continue; // parallel
 240
 241       double uSeg = ( uv1 - gc ) ^ lin / crossSegLin;
 242       if ( 0. <= uSeg && uSeg <= 1. )
 243       {
 244         tria = & _triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID )));
 245         if ( tria->Movability() != BRepMesh_Deleted )
 246           break;
 247       }
 248     }
 249   }
 250   return tria;
 251 }
 252
 253 //================================================================================
 254 /*!
 255  * \brief Return a triangle sharing a given boundary node
 256  *  \param [in] iBndNode - index of the boundary node
 257  *  \return const BRepMesh_Triangle* - a found triangle
 258  */
 259 //================================================================================
 260
 261 const BRepMesh_Triangle* SMESH_Delaunay::GetTriangleNear( int iBndNode )
 262 {
 263   int nodeIDs[3];
 264   int nbNbNodes = _bndNodes.size();
 265   const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndNode + 1 );
 266   BRepMesh::ListOfInteger::const_iterator iLink = linkIds.cbegin();
 267   for ( ; iLink != linkIds.cend(); ++iLink )
 268   {
 269     const BRepMesh_PairOfIndex & triaIds = _triaDS->ElementsConnectedTo( *iLink );
 270     {
 271       const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(1) );
 272       if ( tria.Movability() != BRepMesh_Deleted )
 273       {
 274         _triaDS->ElementNodes( tria, nodeIDs );
 275         if ( nodeIDs[0]-1 < nbNbNodes &&
 276              nodeIDs[1]-1 < nbNbNodes &&
 277              nodeIDs[2]-1 < nbNbNodes )
 278           return &tria;
 279       }
 280     }
 281     if ( triaIds.Extent() > 1 )
 282     {
 283       const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(2) );
 284       if ( tria.Movability() != BRepMesh_Deleted )
 285       {
 286         _triaDS->ElementNodes( tria, nodeIDs );
 287         if ( nodeIDs[0]-1 < nbNbNodes &&
 288              nodeIDs[1]-1 < nbNbNodes &&
 289              nodeIDs[2]-1 < nbNbNodes )
 290           return &tria;
 291       }
 292     }
 293   }
 294   return 0;
 295 }
 296
 297 //================================================================================
 298 /*!
 299  * \brief Return UV of the i-th source boundary node (zero based)
 300  */
 301 //================================================================================
 302
 303 gp_XY SMESH_Delaunay::GetBndUV(const int iNode) const
 304 {
 305   return _triaDS->GetNode( iNode+1 ).Coord();
 306 }
 307
 308 //================================================================================
 309 /*!
 310  * \brief Add non-marked nodes surrounding a given one to a queue
 311  */
 312 //================================================================================
 313
 314 void SMESH_Delaunay::addCloseNodes( const SMDS_MeshNode*     node,
 315                                     const BRepMesh_Triangle* tria,
 316                                     const int                faceID,
 317                                     TNodeTriaList &          _noTriQueue )
 318 {
 319   // find in-FACE nodes
 320   SMDS_ElemIteratorPtr elems = node->GetInverseElementIterator(SMDSAbs_Face);
 321   while ( elems->more() )
 322   {
 323     const SMDS_MeshElement* elem = elems->next();
 324     if ( elem->getshapeId() == faceID )
 325     {
 326       for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i )
 327       {
 328         const SMDS_MeshNode* n = elem->GetNode( i );
 329         if ( !n->isMarked() /*&& n->getshapeId() == faceID*/ )
 330           _noTriQueue.push_back( std::make_pair( n, tria ));
 331       }
 332     }
 333   }
 334 }
 335
 336 //================================================================================
 337 /*!
 338  * \brief Write a python script that creates an equal mesh in Mesh module
 339  */
 340 //================================================================================
 341
 342 void SMESH_Delaunay::ToPython() const
 343 {
 344   SMESH_Comment text;
 345   text << "import salome, SMESH\n";
 346   text << "salome.salome_init()\n";
 347   text << "from salome.smesh import smeshBuilder\n";
 348   text << "smesh = smeshBuilder.New(salome.myStudy)\n";
 349   text << "mesh=smesh.Mesh()\n";
 350   const char* endl = "\n";
 351
 352   for ( int i = 0; i < _triaDS->NbNodes(); ++i )
 353   {
 354     const BRepMesh_Vertex& v = _triaDS->GetNode( i+1 );
 355     text << "mesh.AddNode( " << v.Coord().X() << ", " << v.Coord().Y() << ", 0 )" << endl;
 356   }
 357
 358   int nodeIDs[3];
 359   for ( int i = 0; i < _triaDS->NbElements(); ++i )
 360   {
 361     const BRepMesh_Triangle& t = _triaDS->GetElement( i+1 );
 362     if ( t.Movability() == BRepMesh_Deleted )
 363       continue;
 364     _triaDS->ElementNodes( t, nodeIDs );
 365     text << "mesh.AddFace([ " << nodeIDs[0] << ", " << nodeIDs[1] << ", " << nodeIDs[2] << " ])" << endl;
 366   }
 367
 368   const char* fileName = "/tmp/Delaunay.py";
 369   SMESH_File file( fileName, false );
 370   file.remove();
 371   file.openForWriting();
 372   file.write( text.c_str(), text.size() );
 373   cout << "execfile( '" << fileName << "')" << endl;
 374 }