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     return false;
 342
 343   } // triangulate()
 344 } // namespace
 345
 346 //================================================================================
 347 /*!
 348  * \brief Return nb triangles in a decomposed mesh face
 349  *  \retval int - number of triangles
 350  */
 351 //================================================================================
 352
 353 static int getNbTriangles( const SMDS_MeshElement* face)
 354 {
 355   // WARNING: counting triangles must be coherent with getTriangles()
 356   switch ( face->GetEntityType() )
 357   {
 358   case SMDSEntity_BiQuad_Triangle:
 359   case SMDSEntity_BiQuad_Quadrangle:
 360     return face->NbNodes() - 1;
 361   // case SMDSEntity_Triangle:
 362   // case SMDSEntity_Quad_Triangle:
 363   // case SMDSEntity_Quadrangle:
 364   // case SMDSEntity_Quad_Quadrangle:
 365   // case SMDSEntity_Polygon:
 366   // case SMDSEntity_Quad_Polygon:
 367   default:
 368     return face->NbNodes() - 2;
 369   }
 370   return 0;
 371 }
 372
 373 //================================================================================
 374 /*!
 375  * \brief Decompose a mesh face into triangles
 376  *  \retval int - number of triangles
 377  */
 378 //================================================================================
 379
 380 static int getTriangles( const SMDS_MeshElement*             face,
 381                          std::vector< const SMDS_MeshNode*>& nodes)
 382 {
 383   // WARNING: decomposing into triangles must be coherent with getNbTriangles()
 384   int nbTria, i = 0, nbNodes = face->NbNodes();
 385   SMDS_NodeIteratorPtr nIt = face->interlacedNodesIterator();
 386   nodes.resize( nbNodes * 3 );
 387   nodes[ i++ ] = nIt->next();
 388   nodes[ i++ ] = nIt->next();
 389
 390   const SMDSAbs_EntityType type = face->GetEntityType();
 391   switch ( type )
 392   {
 393   case SMDSEntity_BiQuad_Triangle:
 394   case SMDSEntity_BiQuad_Quadrangle:
 395     nbTria = ( type == SMDSEntity_BiQuad_Triangle ) ? 6 : 8;
 396     nodes[ i++ ] = face->GetNode( nbTria );
 397     while ( i < 3*(nbTria-1) )
 398     {
 399       nodes[ i++ ] = nodes[ i-2 ];
 400       nodes[ i++ ] = nIt->next();
 401       nodes[ i++ ] = nodes[ 2 ];
 402     }
 403     nodes[ i++ ] = nodes[ i-2 ];
 404     nodes[ i++ ] = nodes[ 0 ];
 405     nodes[ i++ ] = nodes[ 2 ];
 406     break;
 407   case SMDSEntity_Triangle:
 408     nbTria = 1;
 409     nodes[ i++ ] = nIt->next();
 410     break;
 411   default:
 412     // case SMDSEntity_Quad_Triangle:
 413     // case SMDSEntity_Quadrangle:
 414     // case SMDSEntity_Quad_Quadrangle:
 415     // case SMDSEntity_Polygon:
 416     // case SMDSEntity_Quad_Polygon:
 417     nbTria = nbNodes - 2;
 418     while ( nIt->more() )
 419       nodes[ i++ ] = nIt->next();
 420
 421     if ( !triangulate( nodes, nbNodes ))
 422     {
 423       nIt = face->interlacedNodesIterator();
 424       nodes[ 0 ] = nIt->next();
 425       nodes[ 1 ] = nIt->next();
 426       nodes[ 2 ] = nIt->next();
 427       for ( i = 3; i < 3*nbTria; i += 3 )
 428       {
 429         nodes[ i+0 ] = nodes[ 0 ];
 430         nodes[ i+1 ] = nodes[ i-1 ];
 431         nodes[ i+2 ] = nIt->next();
 432       }
 433     }
 434     break;
 435   }
 436   return nbTria;
 437 }
 438
 439 // private methods
 440
 441 Driver_Mesh::Status DriverSTL_W_SMDS_Mesh::writeAscii() const
 442 {
 443   Status aResult = DRS_OK;
 444   if ( myFile.empty() ) {
 445     fprintf(stderr, ">> ERREOR : invalid file name \n");
 446     return DRS_FAIL;
 447   }
 448
 449   SMESH_File aFile( myFile, /*openForReading=*/false );
 450   aFile.openForWriting();
 451
 452   std::string buf("solid\n");
 453   aFile.writeRaw( buf.c_str(), buf.size() );
 454
 455   char sval[128];
 456   std::vector< const SMDS_MeshNode* > triaNodes;
 457
 458   SMDS_ElemIteratorPtr itFaces = getFaces();
 459   while ( itFaces->more() )
 460   {
 461     const SMDS_MeshElement* aFace = itFaces->next();
 462     int nbTria = getTriangles( aFace, triaNodes );
 463
 464     for ( int iT = 0, iN = 0; iT < nbTria; ++iT )
 465     {
 466       gp_XYZ normale = getNormale( triaNodes[iN],
 467                                    triaNodes[iN+1],
 468                                    triaNodes[iN+2] );
 469       sprintf (sval,
 470                " facet normal % 12e % 12e % 12e\n"
 471                "   outer loop\n" ,
 472                normale.X(), normale.Y(), normale.Z());
 473       aFile.writeRaw ( sval, 70 + strlen( sval + 70 )); // at least 70 but can be more (WIN)
 474
 475       for ( int jN = 0; jN < 3; ++jN, ++iN )
 476       {
 477         SMESH_TNodeXYZ node = triaNodes[iN];
 478         sprintf (sval,
 479                  "     vertex % 12e % 12e % 12e\n",
 480                  node.X(), node.Y(), node.Z() );
 481         aFile.writeRaw ( sval, 54 + strlen( sval + 54 ));
 482       }
 483       aFile.writeRaw ("   endloop\n"
 484                       " endfacet\n", 21 );
 485     }
 486   }
 487   aFile.writeRaw ("endsolid\n" , 9 );
 488
 489   return aResult;
 490 }
 491
 492 //================================================================================
 493 /*!
 494  * \brief Writes all triangles in binary format
 495  *  \return Driver_Mesh::Status - DRS_FAIL if no file name is provided
 496  */
 497 //================================================================================
 498
 499 Driver_Mesh::Status DriverSTL_W_SMDS_Mesh::writeBinary() const
 500 {
 501   Status aResult = DRS_OK;
 502
 503   if ( myFile.empty() ) {
 504     fprintf(stderr, ">> ERREOR : invalid filename \n");
 505     return DRS_FAIL;
 506   }
 507
 508   SMESH_File aFile( myFile );
 509   aFile.openForWriting();
 510
 511   // we first count the number of triangles
 512   int nbTri = myNbVolumeTrias;
 513   {
 514     SMDS_FaceIteratorPtr itFaces = myMesh->facesIterator();
 515     while ( itFaces->more() ) {
 516       const SMDS_MeshElement* aFace = itFaces->next();
 517       nbTri += getNbTriangles( aFace );
 518     }
 519   }
 520   std::string sval( LABEL_SIZE, ' ' );
 521   aFile.write( sval.c_str(), LABEL_SIZE );
 522
 523   // write number of triangles
 524   writeInteger( nbTri, aFile );
 525
 526   // Loop writing nodes
 527
 528   int dum=0;
 529
 530   std::vector< const SMDS_MeshNode* > triaNodes;
 531
 532   SMDS_ElemIteratorPtr itFaces = getFaces();
 533   while ( itFaces->more() )
 534   {
 535     const SMDS_MeshElement* aFace = itFaces->next();
 536     int nbTria = getTriangles( aFace, triaNodes );
 537
 538     for ( int iT = 0, iN = 0; iT < nbTria; ++iT )
 539     {
 540       gp_XYZ normale = getNormale( triaNodes[iN],
 541                                    triaNodes[iN+1],
 542                                    triaNodes[iN+2] );
 543       writeFloat(normale.X(),aFile);
 544       writeFloat(normale.Y(),aFile);
 545       writeFloat(normale.Z(),aFile);
 546
 547       for ( int jN = 0; jN < 3; ++jN, ++iN )
 548       {
 549         const SMDS_MeshNode* node = triaNodes[iN];
 550         writeFloat(node->X(),aFile);
 551         writeFloat(node->Y(),aFile);
 552         writeFloat(node->Z(),aFile);
 553       }
 554       aFile.writeRaw ( &dum, 2 );
 555     }
 556   }
 557
 558   return aResult;
 559 }