src/SMESHUtils/SMESH_Delaunay.cxx

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