src/DriverSTL/DriverSTL_W_SMDS_Mesh.cxx

   1 // Copyright (C) 2007-2015  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
  23 #include "DriverSTL_W_SMDS_Mesh.h"
  24
  25 #ifdef WIN32
  26 #define NOMINMAX
  27 #endif
  28
  29 #include <Basics_Utils.hxx>
  30
  31 #include "SMDS_FaceOfNodes.hxx"
  32 #include "SMDS_IteratorOnIterators.hxx"
  33 #include "SMDS_Mesh.hxx"
  34 #include "SMDS_MeshElement.hxx"
  35 #include "SMDS_MeshNode.hxx"
  36 #include "SMDS_PolygonalFaceOfNodes.hxx"
  37 #include "SMDS_SetIterator.hxx"
  38 #include "SMDS_VolumeTool.hxx"
  39 #include "SMESH_File.hxx"
  40 #include "SMESH_TypeDefs.hxx"
  41
  42 #include <Standard_ErrorHandler.hxx>
  43 #include <Standard_Failure.hxx>
  44 #include <gp_Ax2.hxx>
  45
  46 #include <limits>
  47
  48
  49 // definition des constantes
  50 static const int LABEL_SIZE = 80;
  51
  52 DriverSTL_W_SMDS_Mesh::DriverSTL_W_SMDS_Mesh()
  53 {
  54   myIsAscii = false;
  55 }
  56
  57 void DriverSTL_W_SMDS_Mesh::SetIsAscii( const bool theIsAscii )
  58 {
  59   myIsAscii = theIsAscii;
  60 }
  61
  62 Driver_Mesh::Status DriverSTL_W_SMDS_Mesh::Perform()
  63 {
  64   Kernel_Utils::Localizer loc;
  65
  66   Status aResult = DRS_OK;
  67
  68   if ( !myMesh ) {
  69     fprintf(stderr, ">> ERROR : Mesh is null \n");
  70     return DRS_FAIL;
  71   }
  72   findVolumeTriangles();
  73   if ( myIsAscii )
  74     aResult = writeAscii();
  75   else
  76     aResult = writeBinary();
  77
  78   return aResult;
  79 }
  80
  81 //================================================================================
  82 /*!
  83  * \brief Destructor deletes temporary faces
  84  */
  85 //================================================================================
  86
  87 DriverSTL_W_SMDS_Mesh::~DriverSTL_W_SMDS_Mesh()
  88 {
  89   for ( unsigned i = 0; i < myVolumeFacets.size(); ++i )
  90     delete myVolumeFacets[i];
  91 }
  92
  93 //================================================================================
  94 /*!
  95  * \brief Finds free facets of volumes for which faces are missing in the mesh
  96  */
  97 //================================================================================
  98
  99 void DriverSTL_W_SMDS_Mesh::findVolumeTriangles()
 100 {
 101   myNbVolumeTrias = 0;
 102
 103   SMDS_VolumeTool theVolume;
 104   SMDS_VolumeIteratorPtr vIt = myMesh->volumesIterator();
 105   std::vector< const SMDS_MeshNode*> nodes;
 106   while ( vIt->more() )
 107   {
 108     theVolume.Set( vIt->next(), /*ignoreCentralNodes=*/false );
 109     for ( int iF = 0; iF < theVolume.NbFaces(); ++iF )
 110       if ( theVolume.IsFreeFace( iF ))
 111       {
 112         const SMDS_MeshNode** n = theVolume.GetFaceNodes(iF);
 113         int                 nbN = theVolume.NbFaceNodes(iF);
 114         nodes.assign( n, n+nbN );
 115         if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*noMedium=*/false))
 116         {
 117           if (( nbN == 9 || nbN == 7 ) &&
 118               ( !theVolume.IsPoly() )) // facet is bi-quaratic
 119           {
 120             int nbTria = nbN - 1;
 121             for ( int iT = 0; iT < nbTria; ++iT )
 122               myVolumeFacets.push_back( new SMDS_FaceOfNodes( n[8], n[0+iT], n[1+iT] ));
 123             myNbVolumeTrias += nbTria;
 124           }
 125           else
 126           {
 127             myVolumeFacets.push_back( new SMDS_PolygonalFaceOfNodes( nodes ));
 128             myNbVolumeTrias += nbN - 2;
 129           }
 130         }
 131       }
 132   }
 133 }
 134
 135 //================================================================================
 136 /*!
 137  * \brief Return iterator on both faces in the mesh and on temporary faces
 138  */
 139 //================================================================================
 140
 141 SMDS_ElemIteratorPtr DriverSTL_W_SMDS_Mesh::getFaces() const
 142 {
 143   SMDS_ElemIteratorPtr facesIter = myMesh->elementsIterator(SMDSAbs_Face);
 144   SMDS_ElemIteratorPtr tmpTriaIter( new SMDS_ElementVectorIterator( myVolumeFacets.begin(),
 145                                                                     myVolumeFacets.end()));
 146   typedef std::vector< SMDS_ElemIteratorPtr > TElemIterVector;
 147   TElemIterVector iters(2);
 148   iters[0] = facesIter;
 149   iters[1] = tmpTriaIter;
 150
 151   typedef SMDS_IteratorOnIterators<const SMDS_MeshElement *, TElemIterVector> TItersIter;
 152   return SMDS_ElemIteratorPtr( new TItersIter( iters ));
 153 }
 154
 155 // static methods
 156
 157 static void writeInteger( const Standard_Integer& theVal, SMESH_File& ofile )
 158 {
 159   union {
 160     Standard_Integer i;
 161     char c[4];
 162   } u;
 163
 164   u.i = theVal;
 165
 166   Standard_Integer entier;
 167   entier  =  u.c[0] & 0xFF;
 168   entier |= (u.c[1] & 0xFF) << 0x08;
 169   entier |= (u.c[2] & 0xFF) << 0x10;
 170   entier |= (u.c[3] & 0xFF) << 0x18;
 171
 172   ofile.write( entier );
 173 }
 174
 175 static void writeFloat( const Standard_ShortReal& theVal, SMESH_File& ofile)
 176 {
 177   union {
 178     Standard_ShortReal f;
 179     char c[4];
 180   } u;
 181
 182   u.f = theVal;
 183
 184   Standard_Integer entier;
 185
 186   entier  =  u.c[0] & 0xFF;
 187   entier |= (u.c[1] & 0xFF) << 0x08;
 188   entier |= (u.c[2] & 0xFF) << 0x10;
 189   entier |= (u.c[3] & 0xFF) << 0x18;
 190
 191   ofile.write( entier );
 192 }
 193
 194 static gp_XYZ getNormale( const SMDS_MeshNode* n1,
 195                           const SMDS_MeshNode* n2,
 196                           const SMDS_MeshNode* n3)
 197 {
 198   SMESH_TNodeXYZ xyz1( n1 );
 199   SMESH_TNodeXYZ xyz2( n2 );
 200   SMESH_TNodeXYZ xyz3( n3 );
 201   gp_XYZ q1 = xyz2 - xyz1;
 202   gp_XYZ q2 = xyz3 - xyz1;
 203   gp_XYZ n  = q1 ^ q2;
 204   double len = n.Modulus();
 205   if ( len > std::numeric_limits<double>::min() )
 206     n /= len;
 207
 208   return n;
 209 }
 210
 211 namespace
 212 {
 213   /*!
 214    * \brief Vertex of a polygon. Together with 2 neighbor Vertices represents a triangle
 215    */
 216   struct PolyVertex
 217   {
 218     SMESH_TNodeXYZ _nxyz;
 219     gp_XY          _xy;
 220     PolyVertex*    _prev;
 221     PolyVertex*    _next;
 222
 223     void SetNodeAndNext( const SMDS_MeshNode* n, PolyVertex& v )
 224     {
 225       _nxyz.Set( n );
 226       _next = &v;
 227       v._prev = this;
 228     }
 229     PolyVertex* Delete()
 230     {
 231       _prev->_next = _next;
 232       _next->_prev = _prev;
 233       return _next;
 234     }
 235     void GetTriaNodes( const SMDS_MeshNode** nodes) const
 236     {
 237       nodes[0] = _prev->_nxyz._node;
 238       nodes[1] =  this->_nxyz._node;
 239       nodes[2] = _next->_nxyz._node;
 240     }
 241
 242     inline static double Area( const PolyVertex* v0, const PolyVertex* v1, const PolyVertex* v2 )
 243     {
 244       gp_XY vPrev = v0->_xy - v1->_xy;
 245       gp_XY vNext = v2->_xy - v1->_xy;
 246       return vNext ^ vPrev;
 247     }
 248     double TriaArea() const { return Area( _prev, this, _next ); }
 249
 250     bool IsInsideTria( const PolyVertex* v )
 251     {
 252       gp_XY p = _prev->_xy - v->_xy;
 253       gp_XY t =  this->_xy - v->_xy;
 254       gp_XY n = _next->_xy - v->_xy;
 255       return (( p ^ t ) > 0 &&
 256               ( t ^ n ) > 0 &&
 257               ( n ^ p ) > 0 );
 258       // return ( Area( _prev, this, v ) > 0 &&
 259       //          Area( this, _next, v ) > 0 &&
 260       //          Area( _next, _prev, v ) > 0 );
 261     }
 262   };
 263
 264   //================================================================================
 265   /*!
 266    * \brief Triangulate a polygon. Assure correct orientation for concave polygons
 267    */
 268   //================================================================================
 269
 270   bool triangulate( std::vector< const SMDS_MeshNode*>& nodes, const size_t nbNodes )
 271   {
 272     // connect nodes into a ring
 273     std::vector< PolyVertex > pv( nbNodes );
 274     for ( size_t i = 1; i < nbNodes; ++i )
 275       pv[i-1].SetNodeAndNext( nodes[i-1], pv[i] );
 276     pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], pv[0] );
 277
 278     // get a polygon normal
 279     gp_XYZ normal(0,0,0), p0,v01,v02;
 280     p0  = pv[0]._nxyz;
 281     v01 = pv[1]._nxyz - p0;
 282     for ( size_t i = 2; i < nbNodes; ++i )
 283     {
 284       v02 = pv[i]._nxyz - p0;
 285       normal += v01 ^ v02;
 286       v01 = v02;
 287     }
 288     // project nodes to the found plane
 289     gp_Ax2 axes;
 290     try {
 291       axes = gp_Ax2( p0, normal, v01 );
 292     }
 293     catch ( Standard_Failure ) {
 294       return false;
 295     }
 296     for ( size_t i = 0; i < nbNodes; ++i )
 297     {
 298       gp_XYZ p = pv[i]._nxyz - p0;
 299       pv[i]._xy.SetX( axes.XDirection().XYZ() * p );
 300       pv[i]._xy.SetY( axes.YDirection().XYZ() * p );
 301     }
 302
 303     // in a loop, find triangles with positive area and having no vertices inside
 304     int iN = 0, nbTria = nbNodes - 2;
 305     nodes.reserve( nbTria * 3 );
 306     const double minArea = 1e-6;
 307     PolyVertex* v = &pv[0], *vi;
 308     int nbVertices = nbNodes, nbBadTria = 0, isGoodTria;
 309     while ( nbBadTria < nbVertices )
 310     {
 311       if (( isGoodTria = v->TriaArea() > minArea ))
 312       {
 313         for ( vi = v->_next->_next;
 314               vi != v->_prev;
 315               vi = vi->_next )
 316         {
 317           if ( v->IsInsideTria( vi ))
 318             break;
 319         }
 320         isGoodTria = ( vi == v->_prev );
 321       }
 322       if ( isGoodTria )
 323       {
 324         v->GetTriaNodes( &nodes[ iN ] );
 325         iN += 3;
 326         v = v->Delete();
 327         if ( --nbVertices == 3 )
 328         {
 329           // last triangle remains
 330           v->GetTriaNodes( &nodes[ iN ] );
 331           return true;
 332         }
 333         nbBadTria = 0;
 334       }
 335       else
 336       {
 337         v = v->_next;
 338         ++nbBadTria;
 339       }
 340     }
 341
 342     // the polygon is invalid; add triangles with positive area
 343     nbBadTria = 0;
 344     while ( nbBadTria < nbVertices )
 345     {
 346       isGoodTria = v->TriaArea() > minArea;
 347       if ( isGoodTria )
 348       {
 349         v->GetTriaNodes( &nodes[ iN ] );
 350         iN += 3;
 351         v = v->Delete();
 352         if ( --nbVertices == 3 )
 353         {
 354           // last triangle remains
 355           v->GetTriaNodes( &nodes[ iN ] );
 356           return true;
 357         }
 358         nbBadTria = 0;
 359       }
 360       else
 361       {
 362         v = v->_next;
 363         ++nbBadTria;
 364       }
 365     }
 366
 367     // add all the rest triangles
 368     while ( nbVertices >= 3 )
 369     {
 370       v->GetTriaNodes( &nodes[ iN ] );
 371       iN += 3;
 372       v = v->Delete();
 373     }
 374
 375     return true;
 376
 377   } // triangulate()
 378 } // namespace
 379
 380 //================================================================================
 381 /*!
 382  * \brief Return nb triangles in a decomposed mesh face
 383  *  \retval int - number of triangles
 384  */
 385 //================================================================================
 386
 387 static int getNbTriangles( const SMDS_MeshElement* face)
 388 {
 389   // WARNING: counting triangles must be coherent with getTriangles()
 390   switch ( face->GetEntityType() )
 391   {
 392   case SMDSEntity_BiQuad_Triangle:
 393   case SMDSEntity_BiQuad_Quadrangle:
 394     return face->NbNodes() - 1;
 395   // case SMDSEntity_Triangle:
 396   // case SMDSEntity_Quad_Triangle:
 397   // case SMDSEntity_Quadrangle:
 398   // case SMDSEntity_Quad_Quadrangle:
 399   // case SMDSEntity_Polygon:
 400   // case SMDSEntity_Quad_Polygon:
 401   default:
 402     return face->NbNodes() - 2;
 403   }
 404   return 0;
 405 }
 406
 407 //================================================================================
 408 /*!
 409  * \brief Decompose a mesh face into triangles
 410  *  \retval int - number of triangles
 411  */
 412 //================================================================================
 413
 414 static int getTriangles( const SMDS_MeshElement*             face,
 415                          std::vector< const SMDS_MeshNode*>& nodes)
 416 {
 417   // WARNING: decomposing into triangles must be coherent with getNbTriangles()
 418   int nbTria, i = 0, nbNodes = face->NbNodes();
 419   SMDS_NodeIteratorPtr nIt = face->interlacedNodesIterator();
 420   nodes.resize( nbNodes * 3 );
 421   nodes[ i++ ] = nIt->next();
 422   nodes[ i++ ] = nIt->next();
 423
 424   const SMDSAbs_EntityType type = face->GetEntityType();
 425   switch ( type )
 426   {
 427   case SMDSEntity_BiQuad_Triangle:
 428   case SMDSEntity_BiQuad_Quadrangle:
 429     nbTria = ( type == SMDSEntity_BiQuad_Triangle ) ? 6 : 8;
 430     nodes[ i++ ] = face->GetNode( nbTria );
 431     for ( i = 3; i < 3*(nbTria-1); i += 3 )
 432     {
 433       nodes[ i+0 ] = nodes[ i-2 ];
 434       nodes[ i+1 ] = nIt->next();
 435       nodes[ i+2 ] = nodes[ 2 ];
 436     }
 437     nodes[ i+0 ] = nodes[ i-2 ];
 438     nodes[ i+1 ] = nodes[ 0 ];
 439     nodes[ i+2 ] = nodes[ 2 ];
 440     break;
 441   case SMDSEntity_Triangle:
 442     nbTria = 1;
 443     nodes[ i++ ] = nIt->next();
 444     break;
 445   default:
 446     // case SMDSEntity_Quad_Triangle:
 447     // case SMDSEntity_Quadrangle:
 448     // case SMDSEntity_Quad_Quadrangle:
 449     // case SMDSEntity_Polygon:
 450     // case SMDSEntity_Quad_Polygon:
 451     nbTria = nbNodes - 2;
 452     while ( nIt->more() )
 453       nodes[ i++ ] = nIt->next();
 454
 455     if ( !triangulate( nodes, nbNodes ))
 456     {
 457       nIt = face->interlacedNodesIterator();
 458       nodes[ 0 ] = nIt->next();
 459       nodes[ 1 ] = nIt->next();
 460       nodes[ 2 ] = nIt->next();
 461       for ( i = 3; i < 3*nbTria; i += 3 )
 462       {
 463         nodes[ i+0 ] = nodes[ 0 ];
 464         nodes[ i+1 ] = nodes[ i-1 ];
 465         nodes[ i+2 ] = nIt->next();
 466       }
 467     }
 468     break;
 469   }
 470   return nbTria;
 471 }
 472
 473 // private methods
 474
 475 Driver_Mesh::Status DriverSTL_W_SMDS_Mesh::writeAscii() const
 476 {
 477   Status aResult = DRS_OK;
 478   if ( myFile.empty() ) {
 479     fprintf(stderr, ">> ERREOR : invalid file name \n");
 480     return DRS_FAIL;
 481   }
 482
 483   SMESH_File aFile( myFile, /*openForReading=*/false );
 484   aFile.openForWriting();
 485
 486   std::string buf("solid\n");
 487   aFile.writeRaw( buf.c_str(), buf.size() );
 488
 489   char sval[128];
 490   std::vector< const SMDS_MeshNode* > triaNodes;
 491
 492   SMDS_ElemIteratorPtr itFaces = getFaces();
 493   while ( itFaces->more() )
 494   {
 495     const SMDS_MeshElement* aFace = itFaces->next();
 496     int nbTria = getTriangles( aFace, triaNodes );
 497
 498     for ( int iT = 0, iN = 0; iT < nbTria; ++iT )
 499     {
 500       gp_XYZ normale = getNormale( triaNodes[iN],
 501                                    triaNodes[iN+1],
 502                                    triaNodes[iN+2] );
 503       sprintf (sval,
 504                " facet normal % 12e % 12e % 12e\n"
 505                "   outer loop\n" ,
 506                normale.X(), normale.Y(), normale.Z());
 507       aFile.writeRaw ( sval, 70 + strlen( sval + 70 )); // at least 70 but can be more (WIN)
 508
 509       for ( int jN = 0; jN < 3; ++jN, ++iN )
 510       {
 511         SMESH_TNodeXYZ node = triaNodes[iN];
 512         sprintf (sval,
 513                  "     vertex % 12e % 12e % 12e\n",
 514                  node.X(), node.Y(), node.Z() );
 515         aFile.writeRaw ( sval, 54 + strlen( sval + 54 ));
 516       }
 517       aFile.writeRaw ("   endloop\n"
 518                       " endfacet\n", 21 );
 519     }
 520   }
 521   aFile.writeRaw ("endsolid\n" , 9 );
 522
 523   return aResult;
 524 }
 525
 526 //================================================================================
 527 /*!
 528  * \brief Writes all triangles in binary format
 529  *  \return Driver_Mesh::Status - DRS_FAIL if no file name is provided
 530  */
 531 //================================================================================
 532
 533 Driver_Mesh::Status DriverSTL_W_SMDS_Mesh::writeBinary() const
 534 {
 535   Status aResult = DRS_OK;
 536
 537   if ( myFile.empty() ) {
 538     fprintf(stderr, ">> ERREOR : invalid filename \n");
 539     return DRS_FAIL;
 540   }
 541
 542   SMESH_File aFile( myFile );
 543   aFile.openForWriting();
 544
 545   // we first count the number of triangles
 546   int nbTri = myNbVolumeTrias;
 547   {
 548     SMDS_FaceIteratorPtr itFaces = myMesh->facesIterator();
 549     while ( itFaces->more() ) {
 550       const SMDS_MeshElement* aFace = itFaces->next();
 551       nbTri += getNbTriangles( aFace );
 552     }
 553   }
 554   std::string sval( LABEL_SIZE, ' ' );
 555   aFile.write( sval.c_str(), LABEL_SIZE );
 556
 557   // write number of triangles
 558   writeInteger( nbTri, aFile );
 559
 560   // Loop writing nodes
 561
 562   int dum=0;
 563
 564   std::vector< const SMDS_MeshNode* > triaNodes;
 565
 566   SMDS_ElemIteratorPtr itFaces = getFaces();
 567   while ( itFaces->more() )
 568   {
 569     const SMDS_MeshElement* aFace = itFaces->next();
 570     int nbTria = getTriangles( aFace, triaNodes );
 571
 572     for ( int iT = 0, iN = 0; iT < nbTria; ++iT )
 573     {
 574       gp_XYZ normale = getNormale( triaNodes[iN],
 575                                    triaNodes[iN+1],
 576                                    triaNodes[iN+2] );
 577       writeFloat(normale.X(),aFile);
 578       writeFloat(normale.Y(),aFile);
 579       writeFloat(normale.Z(),aFile);
 580
 581       for ( int jN = 0; jN < 3; ++jN, ++iN )
 582       {
 583         const SMDS_MeshNode* node = triaNodes[iN];
 584         writeFloat(node->X(),aFile);
 585         writeFloat(node->Y(),aFile);
 586         writeFloat(node->Z(),aFile);
 587       }
 588       aFile.writeRaw ( &dum, 2 );
 589     }
 590   }
 591
 592   return aResult;
 593 }