src/DriverSTL/DriverSTL_W_SMDS_Mesh.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
  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       const double tol = -1e-12;
 256       return (( p ^ t ) >= tol &&
 257               ( t ^ n ) >= tol &&
 258               ( n ^ p ) >= tol );
 259       // return ( Area( _prev, this, v ) > 0 &&
 260       //          Area( this, _next, v ) > 0 &&
 261       //          Area( _next, _prev, v ) > 0 );
 262     }
 263   };
 264
 265   //================================================================================
 266   /*!
 267    * \brief Triangulate a polygon. Assure correct orientation for concave polygons
 268    */
 269   //================================================================================
 270
 271   bool triangulate( std::vector< const SMDS_MeshNode*>& nodes, const size_t nbNodes )
 272   {
 273     // connect nodes into a ring
 274     std::vector< PolyVertex > pv( nbNodes );
 275     for ( size_t i = 1; i < nbNodes; ++i )
 276       pv[i-1].SetNodeAndNext( nodes[i-1], pv[i] );
 277     pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], pv[0] );
 278
 279     // get a polygon normal
 280     gp_XYZ normal(0,0,0), p0,v01,v02;
 281     p0  = pv[0]._nxyz;
 282     v01 = pv[1]._nxyz - p0;
 283     for ( size_t i = 2; i < nbNodes; ++i )
 284     {
 285       v02 = pv[i]._nxyz - p0;
 286       normal += v01 ^ v02;
 287       v01 = v02;
 288     }
 289     // project nodes to the found plane
 290     gp_Ax2 axes;
 291     try {
 292       axes = gp_Ax2( p0, normal, v01 );
 293     }
 294     catch ( Standard_Failure ) {
 295       return false;
 296     }
 297     for ( size_t i = 0; i < nbNodes; ++i )
 298     {
 299       gp_XYZ p = pv[i]._nxyz - p0;
 300       pv[i]._xy.SetX( axes.XDirection().XYZ() * p );
 301       pv[i]._xy.SetY( axes.YDirection().XYZ() * p );
 302     }
 303
 304     // in a loop, find triangles with positive area and having no vertices inside
 305     int iN = 0, nbTria = nbNodes - 2;
 306     nodes.reserve( nbTria * 3 );
 307     const double minArea = 1e-6;
 308     PolyVertex* v = &pv[0], *vi;
 309     int nbVertices = nbNodes, nbBadTria = 0, isGoodTria;
 310     while ( nbBadTria < nbVertices )
 311     {
 312       if (( isGoodTria = v->TriaArea() > minArea ))
 313       {
 314         for ( vi = v->_next->_next;
 315               vi != v->_prev;
 316               vi = vi->_next )
 317         {
 318           if ( v->IsInsideTria( vi ))
 319             break;
 320         }
 321         isGoodTria = ( vi == v->_prev );
 322       }
 323       if ( isGoodTria )
 324       {
 325         v->GetTriaNodes( &nodes[ iN ] );
 326         iN += 3;
 327         v = v->Delete();
 328         if ( --nbVertices == 3 )
 329         {
 330           // last triangle remains
 331           v->GetTriaNodes( &nodes[ iN ] );
 332           return true;
 333         }
 334         nbBadTria = 0;
 335       }
 336       else
 337       {
 338         v = v->_next;
 339         ++nbBadTria;
 340       }
 341     }
 342
 343     // the polygon is invalid; add triangles with positive area
 344     nbBadTria = 0;
 345     while ( nbBadTria < nbVertices )
 346     {
 347       isGoodTria = v->TriaArea() > minArea;
 348       if ( isGoodTria )
 349       {
 350         v->GetTriaNodes( &nodes[ iN ] );
 351         iN += 3;
 352         v = v->Delete();
 353         if ( --nbVertices == 3 )
 354         {
 355           // last triangle remains
 356           v->GetTriaNodes( &nodes[ iN ] );
 357           return true;
 358         }
 359         nbBadTria = 0;
 360       }
 361       else
 362       {
 363         v = v->_next;
 364         ++nbBadTria;
 365       }
 366     }
 367
 368     // add all the rest triangles
 369     while ( nbVertices >= 3 )
 370     {
 371       v->GetTriaNodes( &nodes[ iN ] );
 372       iN += 3;
 373       v = v->Delete();
 374       --nbVertices;
 375     }
 376
 377     return true;
 378
 379   } // triangulate()
 380 } // namespace
 381
 382 //================================================================================
 383 /*!
 384  * \brief Return nb triangles in a decomposed mesh face
 385  *  \retval int - number of triangles
 386  */
 387 //================================================================================
 388
 389 static int getNbTriangles( const SMDS_MeshElement* face)
 390 {
 391   // WARNING: counting triangles must be coherent with getTriangles()
 392   switch ( face->GetEntityType() )
 393   {
 394   case SMDSEntity_BiQuad_Triangle:
 395   case SMDSEntity_BiQuad_Quadrangle:
 396     return face->NbNodes() - 1;
 397   // case SMDSEntity_Triangle:
 398   // case SMDSEntity_Quad_Triangle:
 399   // case SMDSEntity_Quadrangle:
 400   // case SMDSEntity_Quad_Quadrangle:
 401   // case SMDSEntity_Polygon:
 402   // case SMDSEntity_Quad_Polygon:
 403   default:
 404     return face->NbNodes() - 2;
 405   }
 406   return 0;
 407 }
 408
 409 //================================================================================
 410 /*!
 411  * \brief Decompose a mesh face into triangles
 412  *  \retval int - number of triangles
 413  */
 414 //================================================================================
 415
 416 static int getTriangles( const SMDS_MeshElement*             face,
 417                          std::vector< const SMDS_MeshNode*>& nodes)
 418 {
 419   // WARNING: decomposing into triangles must be coherent with getNbTriangles()
 420   int nbTria, i = 0, nbNodes = face->NbNodes();
 421   SMDS_NodeIteratorPtr nIt = face->interlacedNodesIterator();
 422   nodes.resize( nbNodes * 3 );
 423   nodes[ i++ ] = nIt->next();
 424   nodes[ i++ ] = nIt->next();
 425
 426   const SMDSAbs_EntityType type = face->GetEntityType();
 427   switch ( type )
 428   {
 429   case SMDSEntity_BiQuad_Triangle:
 430   case SMDSEntity_BiQuad_Quadrangle:
 431     nbTria = ( type == SMDSEntity_BiQuad_Triangle ) ? 6 : 8;
 432     nodes[ i++ ] = face->GetNode( nbTria );
 433     for ( i = 3; i < 3*(nbTria-1); i += 3 )
 434     {
 435       nodes[ i+0 ] = nodes[ i-2 ];
 436       nodes[ i+1 ] = nIt->next();
 437       nodes[ i+2 ] = nodes[ 2 ];
 438     }
 439     nodes[ i+0 ] = nodes[ i-2 ];
 440     nodes[ i+1 ] = nodes[ 0 ];
 441     nodes[ i+2 ] = nodes[ 2 ];
 442     break;
 443   case SMDSEntity_Triangle:
 444     nbTria = 1;
 445     nodes[ i++ ] = nIt->next();
 446     break;
 447   default:
 448     // case SMDSEntity_Quad_Triangle:
 449     // case SMDSEntity_Quadrangle:
 450     // case SMDSEntity_Quad_Quadrangle:
 451     // case SMDSEntity_Polygon:
 452     // case SMDSEntity_Quad_Polygon:
 453     nbTria = nbNodes - 2;
 454     while ( nIt->more() )
 455       nodes[ i++ ] = nIt->next();
 456
 457     if ( !triangulate( nodes, nbNodes ))
 458     {
 459       nIt = face->interlacedNodesIterator();
 460       nodes[ 0 ] = nIt->next();
 461       nodes[ 1 ] = nIt->next();
 462       nodes[ 2 ] = nIt->next();
 463       for ( i = 3; i < 3*nbTria; i += 3 )
 464       {
 465         nodes[ i+0 ] = nodes[ 0 ];
 466         nodes[ i+1 ] = nodes[ i-1 ];
 467         nodes[ i+2 ] = nIt->next();
 468       }
 469     }
 470     break;
 471   }
 472   return nbTria;
 473 }
 474
 475 // private methods
 476
 477 Driver_Mesh::Status DriverSTL_W_SMDS_Mesh::writeAscii() const
 478 {
 479   Status aResult = DRS_OK;
 480   if ( myFile.empty() ) {
 481     fprintf(stderr, ">> ERREOR : invalid file name \n");
 482     return DRS_FAIL;
 483   }
 484
 485   SMESH_File aFile( myFile, /*openForReading=*/false );
 486   aFile.openForWriting();
 487
 488   std::string buf("solid\n");
 489   aFile.writeRaw( buf.c_str(), buf.size() );
 490
 491   char sval[128];
 492   std::vector< const SMDS_MeshNode* > triaNodes;
 493
 494   SMDS_ElemIteratorPtr itFaces = getFaces();
 495   while ( itFaces->more() )
 496   {
 497     const SMDS_MeshElement* aFace = itFaces->next();
 498     int nbTria = getTriangles( aFace, triaNodes );
 499
 500     for ( int iT = 0, iN = 0; iT < nbTria; ++iT )
 501     {
 502       gp_XYZ normale = getNormale( triaNodes[iN],
 503                                    triaNodes[iN+1],
 504                                    triaNodes[iN+2] );
 505       sprintf (sval,
 506                " facet normal % 12e % 12e % 12e\n"
 507                "   outer loop\n" ,
 508                normale.X(), normale.Y(), normale.Z());
 509       aFile.writeRaw ( sval, 70 + strlen( sval + 70 )); // at least 70 but can be more (WIN)
 510
 511       for ( int jN = 0; jN < 3; ++jN, ++iN )
 512       {
 513         SMESH_TNodeXYZ node = triaNodes[iN];
 514         sprintf (sval,
 515                  "     vertex % 12e % 12e % 12e\n",
 516                  node.X(), node.Y(), node.Z() );
 517         aFile.writeRaw ( sval, 54 + strlen( sval + 54 ));
 518       }
 519       aFile.writeRaw ("   endloop\n"
 520                       " endfacet\n", 21 );
 521     }
 522   }
 523   aFile.writeRaw ("endsolid\n" , 9 );
 524
 525   return aResult;
 526 }
 527
 528 //================================================================================
 529 /*!
 530  * \brief Writes all triangles in binary format
 531  *  \return Driver_Mesh::Status - DRS_FAIL if no file name is provided
 532  */
 533 //================================================================================
 534
 535 Driver_Mesh::Status DriverSTL_W_SMDS_Mesh::writeBinary() const
 536 {
 537   Status aResult = DRS_OK;
 538
 539   if ( myFile.empty() ) {
 540     fprintf(stderr, ">> ERREOR : invalid filename \n");
 541     return DRS_FAIL;
 542   }
 543
 544   SMESH_File aFile( myFile );
 545   aFile.openForWriting();
 546
 547   // we first count the number of triangles
 548   int nbTri = myNbVolumeTrias;
 549   {
 550     SMDS_FaceIteratorPtr itFaces = myMesh->facesIterator();
 551     while ( itFaces->more() ) {
 552       const SMDS_MeshElement* aFace = itFaces->next();
 553       nbTri += getNbTriangles( aFace );
 554     }
 555   }
 556   std::string sval( LABEL_SIZE, ' ' );
 557   aFile.write( sval.c_str(), LABEL_SIZE );
 558
 559   // write number of triangles
 560   writeInteger( nbTri, aFile );
 561
 562   // Loop writing nodes
 563
 564   int dum=0;
 565
 566   std::vector< const SMDS_MeshNode* > triaNodes;
 567
 568   SMDS_ElemIteratorPtr itFaces = getFaces();
 569   while ( itFaces->more() )
 570   {
 571     const SMDS_MeshElement* aFace = itFaces->next();
 572     int nbTria = getTriangles( aFace, triaNodes );
 573
 574     for ( int iT = 0, iN = 0; iT < nbTria; ++iT )
 575     {
 576       gp_XYZ normale = getNormale( triaNodes[iN],
 577                                    triaNodes[iN+1],
 578                                    triaNodes[iN+2] );
 579       writeFloat(normale.X(),aFile);
 580       writeFloat(normale.Y(),aFile);
 581       writeFloat(normale.Z(),aFile);
 582
 583       for ( int jN = 0; jN < 3; ++jN, ++iN )
 584       {
 585         const SMDS_MeshNode* node = triaNodes[iN];
 586         writeFloat(node->X(),aFile);
 587         writeFloat(node->Y(),aFile);
 588         writeFloat(node->Z(),aFile);
 589       }
 590       aFile.writeRaw ( &dum, 2 );
 591     }
 592   }
 593
 594   return aResult;
 595 }