src/GHS3DPlugin_GHS3D.cxx

   1 //=============================================================================
   2 // File      : GHS3DPlugin_GHS3D.cxx
   3 // Created   :
   4 // Author    : Edward AGAPOV
   5 // Project   : SALOME
   6 // Copyright : CEA 2003
   7 // $Header$
   8 //=============================================================================
   9 using namespace std;
  10
  11 #include "GHS3DPlugin_GHS3D.hxx"
  12 #include "SMESH_Gen.hxx"
  13 #include "SMESH_Mesh.hxx"
  14
  15 #include "SMDS_MeshElement.hxx"
  16 #include "SMDS_MeshNode.hxx"
  17
  18 #include <TopExp_Explorer.hxx>
  19
  20 #include "utilities.h"
  21
  22 #include <qdir.h>
  23
  24 #ifdef _DEBUG_
  25 #define DUMP(txt) \
  26 //  cout << txt
  27 #else
  28 #define DUMP(txt)
  29 #endif
  30
  31 //=============================================================================
  32 /*!
  33  *
  34  */
  35 //=============================================================================
  36
  37 GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D(int hypId, int studyId, SMESH_Gen* gen)
  38   : SMESH_3D_Algo(hypId, studyId, gen)
  39 {
  40   MESSAGE("GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D");
  41   _name = "GHS3D";
  42   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  43 }
  44
  45 //=============================================================================
  46 /*!
  47  *
  48  */
  49 //=============================================================================
  50
  51 GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D()
  52 {
  53   MESSAGE("GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D");
  54 }
  55
  56 //=============================================================================
  57 /*!
  58  *
  59  */
  60 //=============================================================================
  61
  62 bool GHS3DPlugin_GHS3D::CheckHypothesis
  63                          (SMESH_Mesh& aMesh,
  64                           const TopoDS_Shape& aShape,
  65                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
  66 {
  67 //  MESSAGE("GHS3DPlugin_GHS3D::CheckHypothesis");
  68   aStatus = SMESH_Hypothesis::HYP_OK;
  69   return true;
  70 }
  71
  72 //=======================================================================
  73 //function : writeFaces
  74 //purpose  :
  75 //=======================================================================
  76
  77 static bool writeFaces (ofstream &            theFile,
  78                         SMESHDS_Mesh *        theMesh,
  79                         const map <int,int> & theSmdsToGhs3dIdMap)
  80 {
  81   // structure:
  82   //
  83   // NB_ELEMS DUMMY_INT
  84   // Loop from 1 to NB_ELEMS
  85   //   NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
  86
  87   int nbFaces = theMesh->NbFaces();
  88   if ( nbFaces == 0 )
  89     return false;
  90
  91   const char* space    = "  ";
  92   const int   dummyint = 0;
  93
  94   // NB_ELEMS DUMMY_INT
  95   theFile << space << nbFaces << space << dummyint << endl;
  96
  97   // Loop from 1 to NB_ELEMS
  98   SMDS_FaceIteratorPtr it = theMesh->facesIterator();
  99   while ( it->more() )
 100   {
 101     // NB_NODES
 102     const SMDS_MeshElement* elem = it->next();
 103     const int nbNodes = elem->NbNodes();
 104     theFile << space << nbNodes;
 105
 106     // NODE_NB_1 NODE_NB_2 ...
 107     SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator();
 108     while ( nodeIt->more() )
 109     {
 110       // find GHS3D ID
 111       int aSmdsID = nodeIt->next()->GetID();
 112       map<int,int>::const_iterator it = theSmdsToGhs3dIdMap.find( aSmdsID );
 113       ASSERT( it != theSmdsToGhs3dIdMap.end() );
 114       theFile << space << (*it).second;
 115     }
 116
 117     // (NB_NODES + 1) times: DUMMY_INT
 118     for ( int i=0; i<=nbNodes; i++)
 119       theFile << space << dummyint;
 120
 121     theFile << endl;
 122   }
 123
 124   return true;
 125 }
 126
 127 //=======================================================================
 128 //function : writePoints
 129 //purpose  :
 130 //=======================================================================
 131
 132 static bool writePoints (ofstream &                       theFile,
 133                          SMESHDS_Mesh *                   theMesh,
 134                          map <int,int> &                  theSmdsToGhs3dIdMap,
 135                          map <int,const SMDS_MeshNode*> & theGhs3dIdToNodeMap)
 136 {
 137   // structure:
 138   //
 139   // NB_NODES
 140   // Loop from 1 to NB_NODES
 141   //   X Y Z DUMMY_INT
 142
 143   int nbNodes = theMesh->NbNodes();
 144   if ( nbNodes == 0 )
 145     return false;
 146
 147   const char* space    = "  ";
 148   const int   dummyint = 0;
 149
 150   // NB_NODES
 151   theFile << space << nbNodes << endl;
 152
 153   // Loop from 1 to NB_NODES
 154   int aGhs3dID = 1;
 155   SMDS_NodeIteratorPtr it = theMesh->nodesIterator();
 156   while ( it->more() )
 157   {
 158     const SMDS_MeshNode* node = it->next();
 159     theSmdsToGhs3dIdMap.insert( map <int,int>::value_type( node->GetID(), aGhs3dID ));
 160     theGhs3dIdToNodeMap.insert (map <int,const SMDS_MeshNode*>::value_type( aGhs3dID, node ));
 161     aGhs3dID++;
 162
 163     // X Y Z DUMMY_INT
 164     theFile
 165       << space << node->X()
 166       << space << node->Y()
 167       << space << node->Z()
 168       << space << dummyint;
 169
 170     theFile << endl;
 171   }
 172
 173   return true;
 174 }
 175
 176 //=======================================================================
 177 //function : getInt
 178 //purpose  :
 179 //=======================================================================
 180
 181 static bool getInt( int & theValue, char * & theLine )
 182 {
 183   char *ptr;
 184   theValue = strtol( theLine, &ptr, 10 );
 185   if ( ptr == theLine ||
 186       // there must not be neither '.' nor ',' nor 'E' ...
 187       (*ptr != ' ' && *ptr != '\n' && *ptr != '\0'))
 188     return false;
 189
 190   DUMP( "  " << theValue );
 191   theLine = ptr;
 192   return true;
 193 }
 194
 195 //=======================================================================
 196 //function : getDouble
 197 //purpose  :
 198 //=======================================================================
 199
 200 static bool getDouble( double & theValue, char * & theLine )
 201 {
 202   char *ptr;
 203   theValue = strtod( theLine, &ptr );
 204   if ( ptr == theLine )
 205     return false;
 206
 207   DUMP( "   " << theValue );
 208   theLine = ptr;
 209   return true;
 210 }
 211
 212 //=======================================================================
 213 //function : readLine
 214 //purpose  :
 215 //=======================================================================
 216
 217 #define GHS3DPlugin_BUFLENGTH 256
 218 #define GHS3DPlugin_ReadLine(aPtr,aBuf,aFile,aLineNb) \
 219 {  aPtr = fgets( aBuf, GHS3DPlugin_BUFLENGTH - 2, aFile ); aLineNb++; DUMP(endl); }
 220
 221 //=======================================================================
 222 //function : readResult
 223 //purpose  :
 224 //=======================================================================
 225
 226 static bool readResult(FILE *                          theFile,
 227                        SMESHDS_Mesh *                  theMesh,
 228                        const TopoDS_Shape&             theShape,
 229                        map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap)
 230 {
 231   // structure:
 232
 233   // record 1:
 234   //  NB_ELEMENTS NB_NODES NB_INPUT_NODES (14 DUMMY_INT)
 235   // record 2:
 236   //  (NB_ELEMENTS * 4) node nbs
 237   // record 3:
 238   //  (NB_NODES) node XYZ
 239
 240   char aBuffer[ GHS3DPlugin_BUFLENGTH ];
 241   char * aPtr;
 242   int aLineNb = 0;
 243
 244   // get shell to set nodes in
 245   TopExp_Explorer exp( theShape, TopAbs_SHELL );
 246   TopoDS_Shell aShell = TopoDS::Shell( exp.Current() );
 247
 248   // -------------------------------------
 249   // record 1:
 250   // read nb generated elements and nodes
 251   // -------------------------------------
 252   int nbElems = 0 , nbNodes = 0, nbInputNodes = 0;
 253   GHS3DPlugin_ReadLine( aPtr, aBuffer, theFile, aLineNb );
 254   if (!aPtr ||
 255       !getInt( nbElems, aPtr ) ||
 256       !getInt( nbNodes, aPtr ) ||
 257       !getInt( nbInputNodes, aPtr))
 258     return false;
 259
 260   // -------------------------------------------
 261   // record 2:
 262   // read element nodes and create tetrahedrons
 263   // -------------------------------------------
 264   GHS3DPlugin_ReadLine( aPtr, aBuffer, theFile, aLineNb );
 265   for (int iElem = 0; iElem < nbElems; iElem++)
 266   {
 267     // read 4 nodes
 268     const SMDS_MeshNode * node[4];
 269     for (int iNode = 0; iNode < 4; iNode++)
 270     {
 271       // read Ghs3d node ID
 272       int ID = 0;
 273       if (!aPtr || ! getInt ( ID, aPtr ))
 274       {
 275         GHS3DPlugin_ReadLine( aPtr, aBuffer, theFile, aLineNb );
 276         if (!aPtr || ! getInt ( ID, aPtr ))
 277         {
 278           MESSAGE( "Cant read " << (iNode+1) << "-th node on line " << aLineNb );
 279           return false;
 280         }
 281       }
 282       // find/create a node with ID
 283       map <int,const SMDS_MeshNode*>::iterator IdNode = theGhs3dIdToNodeMap.find( ID );
 284       if ( IdNode == theGhs3dIdToNodeMap.end())
 285       {
 286         // ID is not yet in theGhs3dIdToNodeMap
 287         ASSERT ( ID > nbInputNodes ); // it should be a new one
 288         SMDS_MeshNode * aNewNode = theMesh->AddNode( 0.,0.,0. ); // read XYZ later
 289         theMesh->SetNodeInVolume( aNewNode, aShell );
 290         theGhs3dIdToNodeMap.insert
 291           ( map <int,const SMDS_MeshNode*>::value_type( ID, aNewNode ));
 292         node[ iNode ] = aNewNode;
 293       }
 294       else
 295       {
 296         node[ iNode ] = (*IdNode).second;
 297       }
 298     }
 299     // create a tetrahedron
 300     SMDS_MeshElement* aTet = theMesh->AddVolume( node[0], node[1], node[2], node[3] );
 301     theMesh->SetMeshElementOnShape( aTet, theShape );
 302   }
 303
 304   // ------------------------
 305   // record 3:
 306   // read and set nodes' XYZ
 307   // ------------------------
 308   GHS3DPlugin_ReadLine( aPtr, aBuffer, theFile, aLineNb );
 309   for (int iNode = 0; iNode < nbNodes; iNode++)
 310   {
 311     // read 3 coordinates
 312     double coord [3];
 313     for (int iCoord = 0; iCoord < 3; iCoord++)
 314     {
 315       if (!aPtr || ! getDouble ( coord[ iCoord ], aPtr ))
 316       {
 317         GHS3DPlugin_ReadLine( aPtr, aBuffer, theFile, aLineNb );
 318         if (!aPtr || ! getDouble ( coord[ iCoord ], aPtr ))
 319         {
 320           MESSAGE( "Cant read " << (iCoord+1) << "-th node coord on line " << aLineNb );
 321           return false;
 322         }
 323       }
 324     }
 325     // do not move old nodes
 326     int ID = iNode + 1;
 327     if (ID <= nbInputNodes)
 328       continue;
 329     // find a node
 330     map <int,const SMDS_MeshNode*>::iterator IdNode = theGhs3dIdToNodeMap.find( ID );
 331     ASSERT ( IdNode != theGhs3dIdToNodeMap.end());
 332     SMDS_MeshNode* node = const_cast<SMDS_MeshNode*> ( (*IdNode).second );
 333
 334     // set XYZ
 335     theMesh->MoveNode( node, coord[0], coord[1], coord[2] );
 336   }
 337
 338   return true;
 339 }
 340
 341 //=============================================================================
 342 /*!
 343  *Here we are going to use the GHS3D mesher
 344  */
 345 //=============================================================================
 346
 347 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
 348                                 const TopoDS_Shape& theShape)
 349 {
 350   MESSAGE("GHS3DPlugin_GHS3D::Compute");
 351
 352   // working dir
 353   const char* wdName = "tmp";
 354   QDir wd = QDir::root();                      // "/"
 355   if ( !wd.cd( wdName ) ) {                    // "/tmp"
 356     MESSAGE( "Cannot find the " << wdName << " directory" );
 357     return false;
 358   }
 359
 360   const QString aGenericName    = wd.filePath( "GHS3DTMP" );
 361   const QString aFacesFileName  = aGenericName + ".faces";
 362   const QString aPointsFileName = aGenericName + ".points";
 363   const QString aResultFileName = aGenericName + ".noboite";
 364   const QString aErrorFileName  = aGenericName + ".error";
 365
 366   // remove old files
 367   QFile( aFacesFileName ).remove();
 368   QFile( aPointsFileName ).remove();
 369   QFile( aResultFileName ).remove();
 370   QFile( aErrorFileName ).remove();
 371
 372
 373   // -----------------
 374   // make input files
 375   // -----------------
 376
 377   ofstream aFacesFile  ( aFacesFileName.latin1()  , ios::out);
 378   ofstream aPointsFile ( aPointsFileName.latin1() , ios::out);
 379   if (!aFacesFile  || !aFacesFile.rdbuf()->is_open() ||
 380       !aPointsFile || !aPointsFile.rdbuf()->is_open())
 381   {
 382     MESSAGE( "Can't write into " << wdName << " directory");
 383     return false;
 384   }
 385   SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
 386   map <int,int> aSmdsToGhs3dIdMap;
 387   map <int,const SMDS_MeshNode*> aGhs3dIdToNodeMap;
 388
 389   bool Ok =
 390     (writePoints( aPointsFile, meshDS, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap ) &&
 391      writeFaces ( aFacesFile, meshDS, aSmdsToGhs3dIdMap ));
 392
 393   aFacesFile.close();
 394   aPointsFile.close();
 395
 396   if ( ! Ok )
 397     return false;
 398
 399   // -----------------
 400   // run ghs3d mesher
 401   // -----------------
 402
 403   QString cmd = "ghs3d "
 404     "-m 1000 "                 // memory: 1000 M
 405       "-f " + aGenericName +   // file to read
 406         " 1>" + aErrorFileName; // dump into file
 407   if (system(cmd.latin1()))
 408   {
 409     MESSAGE ("Failed: <" << cmd.latin1() << ">");
 410     return false;
 411   }
 412
 413   // --------------
 414   // read a result
 415   // --------------
 416
 417   FILE * aResultFile = fopen( aResultFileName.latin1(), "r" );
 418   if (!aResultFile)
 419   {
 420     MESSAGE( "GHS3D ERROR: see " << aErrorFileName.latin1() );
 421     return false;
 422   }
 423
 424   Ok = readResult( aResultFile, meshDS, theShape, aGhs3dIdToNodeMap );
 425   fclose(aResultFile);
 426
 427   return Ok;
 428 }
 429
 430
 431 //=============================================================================
 432 /*!
 433  *
 434  */
 435 //=============================================================================
 436
 437 ostream & GHS3DPlugin_GHS3D::SaveTo(ostream & save)
 438 {
 439   return save;
 440 }
 441
 442 //=============================================================================
 443 /*!
 444  *
 445  */
 446 //=============================================================================
 447
 448 istream & GHS3DPlugin_GHS3D::LoadFrom(istream & load)
 449 {
 450   return load;
 451 }
 452
 453 //=============================================================================
 454 /*!
 455  *
 456  */
 457 //=============================================================================
 458
 459 ostream & operator << (ostream & save, GHS3DPlugin_GHS3D & hyp)
 460 {
 461   return hyp.SaveTo( save );
 462 }
 463
 464 //=============================================================================
 465 /*!
 466  *
 467  */
 468 //=============================================================================
 469
 470 istream & operator >> (istream & load, GHS3DPlugin_GHS3D & hyp)
 471 {
 472   return hyp.LoadFrom( load );
 473 }