src/SMESHUtils/SMESH_Triangulate.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_Triangulate.cxx
  23 // Created   : Thu Jan 18 18:00:13 2018
  24 // Author    : Edward AGAPOV (eap)
  25
  26 // Extracted from ../DriverSTL/DriverSTL_W_SMDS_Mesh.cxx
  27
  28 #include "SMESH_MeshAlgos.hxx"
  29
  30 #include <Standard_ErrorHandler.hxx>
  31 #include <Standard_Failure.hxx>
  32 #include <gp_Ax2.hxx>
  33
  34 using namespace SMESH_MeshAlgos;
  35
  36 //================================================================================
  37 /*!
  38  * \brief Initialization
  39  */
  40 //================================================================================
  41
  42 void Triangulate::PolyVertex::SetNodeAndNext( const SMDS_MeshNode* n,
  43                                               PolyVertex&          v )
  44 {
  45   _nxyz.Set( n );
  46   _next = &v;
  47   v._prev = this;
  48 }
  49 //================================================================================
  50 /*!
  51  * \brief Remove self from a polygon
  52  */
  53 //================================================================================
  54
  55 Triangulate::PolyVertex* Triangulate::PolyVertex::Delete()
  56 {
  57   _prev->_next = _next;
  58   _next->_prev = _prev;
  59   return _next;
  60 }
  61
  62 //================================================================================
  63 /*!
  64  * \brief Return nodes of a triangle
  65  */
  66 //================================================================================
  67
  68 void Triangulate::PolyVertex::GetTriaNodes( const SMDS_MeshNode** nodes) const
  69 {
  70   nodes[0] = _prev->_nxyz._node;
  71   nodes[1] =  this->_nxyz._node;
  72   nodes[2] = _next->_nxyz._node;
  73 }
  74
  75 //================================================================================
  76 /*!
  77  * \brief Compute triangle area
  78  */
  79 //================================================================================
  80
  81 inline static double Area( const gp_XY& xy0, const gp_XY& xy1, const gp_XY& xy2 )
  82 {
  83   gp_XY vPrev = xy0 - xy1;
  84   gp_XY vNext = xy2 - xy1;
  85   return vNext ^ vPrev;
  86 }
  87
  88 //================================================================================
  89 /*!
  90  * \brief Compute triangle area
  91  */
  92 //================================================================================
  93
  94 double Triangulate::PolyVertex::TriaArea() const
  95 {
  96   return Area( _prev->_xy, this->_xy, _next->_xy );
  97 }
  98
  99 //================================================================================
 100 /*!
 101  * \brief Check if a vertex is inside a triangle
 102  */
 103 //================================================================================
 104
 105 bool Triangulate::PolyVertex::IsInsideTria( const PolyVertex* v )
 106 {
 107   if ( this ->_nxyz == v->_nxyz ||
 108        _prev->_nxyz == v->_nxyz ||
 109        _next->_nxyz == v->_nxyz )
 110     return false;
 111
 112   gp_XY p = _prev->_xy - v->_xy;
 113   gp_XY t =  this->_xy - v->_xy;
 114   gp_XY n = _next->_xy - v->_xy;
 115   const double tol = -1e-12;
 116   return (( p ^ t ) >= tol &&
 117           ( t ^ n ) >= tol &&
 118           ( n ^ p ) >= tol );
 119   // return ( Area( _prev, this, v ) > 0 &&
 120   //          Area( this, _next, v ) > 0 &&
 121   //          Area( _next, _prev, v ) > 0 );
 122 }
 123
 124 //================================================================================
 125 /*!
 126  * \brief Triangulate a polygon. Assure correct orientation for concave polygons
 127  */
 128 //================================================================================
 129
 130 bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
 131                                const size_t                        nbNodes )
 132 {
 133   // connect nodes into a ring
 134   _pv.resize( nbNodes );
 135   for ( size_t i = 1; i < nbNodes; ++i )
 136     _pv[i-1].SetNodeAndNext( nodes[i-1], _pv[i] );
 137   _pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], _pv[0] );
 138
 139   // get a polygon normal
 140   gp_XYZ normal(0,0,0), p0,v01,v02;
 141   p0  = _pv[0]._nxyz;
 142   v01 = _pv[1]._nxyz - p0;
 143   for ( size_t i = 2; i < nbNodes; ++i )
 144   {
 145     v02 = _pv[i]._nxyz - p0;
 146     normal += v01 ^ v02;
 147     v01 = v02;
 148   }
 149   // project nodes to the found plane
 150   gp_Ax2 axes;
 151   try {
 152     axes = gp_Ax2( p0, normal, v01 );
 153   }
 154   catch ( Standard_Failure ) {
 155     return false;
 156   }
 157   for ( size_t i = 0; i < nbNodes; ++i )
 158   {
 159     gp_XYZ p = _pv[i]._nxyz - p0;
 160     _pv[i]._xy.SetX( axes.XDirection().XYZ() * p );
 161     _pv[i]._xy.SetY( axes.YDirection().XYZ() * p );
 162   }
 163
 164   // in a loop, find triangles with positive area and having no vertices inside
 165   int iN = 0, nbTria = nbNodes - 2;
 166   nodes.reserve( nbTria * 3 );
 167   const double minArea = 1e-6;
 168   PolyVertex* v = &_pv[0], *vi;
 169   int nbVertices = nbNodes, nbBadTria = 0, isGoodTria;
 170   while ( nbBadTria < nbVertices )
 171   {
 172     if (( isGoodTria = v->TriaArea() > minArea ))
 173     {
 174       for ( vi = v->_next->_next;
 175             vi != v->_prev;
 176             vi = vi->_next )
 177       {
 178         if ( v->IsInsideTria( vi ))
 179           break;
 180       }
 181       isGoodTria = ( vi == v->_prev );
 182     }
 183     if ( isGoodTria )
 184     {
 185       v->GetTriaNodes( &nodes[ iN ] );
 186       iN += 3;
 187       v = v->Delete();
 188       if ( --nbVertices == 3 )
 189       {
 190         // last triangle remains
 191         v->GetTriaNodes( &nodes[ iN ] );
 192         return true;
 193       }
 194       nbBadTria = 0;
 195     }
 196     else
 197     {
 198       v = v->_next;
 199       ++nbBadTria;
 200     }
 201   }
 202
 203   // the polygon is invalid; add triangles with positive area
 204   nbBadTria = 0;
 205   while ( nbBadTria < nbVertices )
 206   {
 207     isGoodTria = v->TriaArea() > minArea;
 208     if ( isGoodTria )
 209     {
 210       v->GetTriaNodes( &nodes[ iN ] );
 211       iN += 3;
 212       v = v->Delete();
 213       if ( --nbVertices == 3 )
 214       {
 215         // last triangle remains
 216         v->GetTriaNodes( &nodes[ iN ] );
 217         return true;
 218       }
 219       nbBadTria = 0;
 220     }
 221     else
 222     {
 223       v = v->_next;
 224       ++nbBadTria;
 225     }
 226   }
 227
 228   // add all the rest triangles
 229   while ( nbVertices >= 3 )
 230   {
 231     v->GetTriaNodes( &nodes[ iN ] );
 232     iN += 3;
 233     v = v->Delete();
 234     --nbVertices;
 235   }
 236
 237   return true;
 238
 239 } // triangulate()
 240
 241 //================================================================================
 242 /*!
 243  * \brief Return nb triangles in a decomposed mesh face
 244  *  \retval int - number of triangles
 245  */
 246 //================================================================================
 247
 248 int Triangulate::GetNbTriangles( const SMDS_MeshElement* face )
 249 {
 250   // WARNING: counting triangles must be coherent with GetTriangles()
 251   switch ( face->GetEntityType() )
 252   {
 253   case SMDSEntity_BiQuad_Triangle:
 254   case SMDSEntity_BiQuad_Quadrangle:
 255     return face->NbNodes() - 1;
 256     // case SMDSEntity_Triangle:
 257     // case SMDSEntity_Quad_Triangle:
 258     // case SMDSEntity_Quadrangle:
 259     // case SMDSEntity_Quad_Quadrangle:
 260     // case SMDSEntity_Polygon:
 261     // case SMDSEntity_Quad_Polygon:
 262   default:
 263     return face->NbNodes() - 2;
 264   }
 265   return 0;
 266 }
 267
 268 //================================================================================
 269 /*!
 270  * \brief Decompose a mesh face into triangles
 271  *  \retval int - number of triangles
 272  */
 273 //================================================================================
 274
 275 int Triangulate::GetTriangles( const SMDS_MeshElement*             face,
 276                                std::vector< const SMDS_MeshNode*>& nodes)
 277 {
 278   if ( face->GetType() != SMDSAbs_Face )
 279     return 0;
 280
 281   // WARNING: decomposing into triangles must be coherent with getNbTriangles()
 282   int nbTria, i = 0, nbNodes = face->NbNodes();
 283   SMDS_NodeIteratorPtr nIt = face->interlacedNodesIterator();
 284   nodes.resize( nbNodes * 3 );
 285   nodes[ i++ ] = nIt->next();
 286   nodes[ i++ ] = nIt->next();
 287
 288   const SMDSAbs_EntityType type = face->GetEntityType();
 289   switch ( type )
 290   {
 291   case SMDSEntity_BiQuad_Triangle:
 292   case SMDSEntity_BiQuad_Quadrangle:
 293
 294     nbTria = ( type == SMDSEntity_BiQuad_Triangle ) ? 6 : 8;
 295     nodes[ i++ ] = face->GetNode( nbTria );
 296     for ( i = 3; i < 3*(nbTria-1); i += 3 )
 297     {
 298       nodes[ i+0 ] = nodes[ i-2 ];
 299       nodes[ i+1 ] = nIt->next();
 300       nodes[ i+2 ] = nodes[ 2 ];
 301     }
 302     nodes[ i+0 ] = nodes[ i-2 ];
 303     nodes[ i+1 ] = nodes[ 0 ];
 304     nodes[ i+2 ] = nodes[ 2 ];
 305     break;
 306
 307   case SMDSEntity_Triangle:
 308
 309     nbTria = 1;
 310     nodes[ i++ ] = nIt->next();
 311     break;
 312
 313   default:
 314
 315     // case SMDSEntity_Quad_Triangle:
 316     // case SMDSEntity_Quadrangle:
 317     // case SMDSEntity_Quad_Quadrangle:
 318     // case SMDSEntity_Polygon:
 319     // case SMDSEntity_Quad_Polygon:
 320
 321     nbTria = nbNodes - 2;
 322     while ( nIt->more() )
 323       nodes[ i++ ] = nIt->next();
 324
 325     if ( nbTria > 1 && !triangulate( nodes, nbNodes ))
 326     {
 327       nIt = face->interlacedNodesIterator();
 328       nodes[ 0 ] = nIt->next();
 329       nodes[ 1 ] = nIt->next();
 330       nodes[ 2 ] = nIt->next();
 331       for ( i = 3; i < 3*nbTria; i += 3 )
 332       {
 333         nodes[ i+0 ] = nodes[ 0 ];
 334         nodes[ i+1 ] = nodes[ i-1 ];
 335         nodes[ i+2 ] = nIt->next();
 336       }
 337     }
 338   }
 339
 340   return nbTria;
 341 }