src/GHS3DPlugin_GHS3D.cxx

   1 // Copyright (C) 2005  CEA/DEN, EDF R&D, OPEN CASCADE, PRINCIPIA R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 //=============================================================================
  20 // File      : GHS3DPlugin_GHS3D.cxx
  21 // Created   :
  22 // Author    : Edward AGAPOV, modified by Lioka RAZAFINDRAZAKA (CEA) 09/02/2007
  23 // Project   : SALOME
  24 // Copyright : CEA 2003
  25 // $Header$
  26 //=============================================================================
  27 using namespace std;
  28
  29 #include "GHS3DPlugin_GHS3D.hxx"
  30 #include "SMESH_Gen.hxx"
  31 #include "SMESH_Mesh.hxx"
  32
  33 #include "SMDS_MeshElement.hxx"
  34 #include "SMDS_MeshNode.hxx"
  35
  36 #include <TopExp_Explorer.hxx>
  37 #include <OSD_File.hxx>
  38
  39 #include "utilities.h"
  40
  41 #ifndef WIN32
  42 #include <sys/sysinfo.h>
  43 #endif
  44
  45 //#include <Standard_Stream.hxx>
  46
  47 #include <BRepGProp.hxx>
  48 #include <BRepBndLib.hxx>
  49 #include <BRepClass_FaceClassifier.hxx>
  50 #include <BRepClass3d_SolidClassifier.hxx>
  51 #include <TopAbs.hxx>
  52 #include <Bnd_Box.hxx>
  53 #include <GProp_GProps.hxx>
  54 #include <Precision.hxx>
  55
  56 #ifdef _DEBUG_
  57 #define DUMP(txt) \
  58 //  cout << txt
  59 #else
  60 #define DUMP(txt)
  61 #endif
  62
  63 // include for mmap
  64 #include "HDFconvert.hxx"
  65
  66 //=============================================================================
  67 /*!
  68  *
  69  */
  70 //=============================================================================
  71
  72 GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D(int hypId, int studyId, SMESH_Gen* gen)
  73   : SMESH_3D_Algo(hypId, studyId, gen)
  74 {
  75   MESSAGE("GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D");
  76   _name = "GHS3D_3D";
  77   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  78   _iShape=0;
  79   _nbShape=0;
  80 }
  81
  82 //=============================================================================
  83 /*!
  84  *
  85  */
  86 //=============================================================================
  87
  88 GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D()
  89 {
  90   MESSAGE("GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D");
  91 }
  92
  93 //=============================================================================
  94 /*!
  95  *
  96  */
  97 //=============================================================================
  98
  99 bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh&                          aMesh,
 100                                           const TopoDS_Shape&                  aShape,
 101                                           SMESH_Hypothesis::Hypothesis_Status& aStatus )
 102 {
 103 //  MESSAGE("GHS3DPlugin_GHS3D::CheckHypothesis");
 104   aStatus = SMESH_Hypothesis::HYP_OK;
 105   return true;
 106 }
 107
 108 //=======================================================================
 109 //function : writeFaces
 110 //purpose  :
 111 //=======================================================================
 112
 113 static bool writeFaces (ofstream &            theFile,
 114                         SMESHDS_Mesh *        theMesh,
 115                         const map <int,int> & theSmdsToGhs3dIdMap)
 116 {
 117   // record structure:
 118   //
 119   // NB_ELEMS DUMMY_INT
 120   // Loop from 1 to NB_ELEMS
 121   // NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
 122
 123   // get all faces bound to theShape
 124
 125   int nbFaces = 0;
 126   TopoDS_Shape theShape = theMesh->ShapeToMesh();
 127   list< const SMDS_MeshElement* > faces;
 128   TopExp_Explorer fExp( theShape, TopAbs_FACE );
 129   SMESHDS_SubMesh* sm;
 130   SMDS_ElemIteratorPtr eIt;
 131
 132   const char* space    = "  ";
 133   const int   dummyint = 0;
 134
 135   list< const SMDS_MeshElement* >::iterator f;
 136   map<int,int>::const_iterator it;
 137   SMDS_ElemIteratorPtr nodeIt;
 138   const SMDS_MeshElement* elem;
 139   int nbNodes;
 140   int aSmdsID;
 141
 142   for ( ; fExp.More(); fExp.Next() ) {
 143     sm = theMesh->MeshElements( fExp.Current() );
 144     if ( sm ) {
 145       eIt = sm->GetElements();
 146       while ( eIt->more() ) {
 147         faces.push_back( eIt->next() );
 148         nbFaces++;
 149       }
 150     }
 151   }
 152
 153   if ( nbFaces == 0 )
 154     return false;
 155
 156   cout << "  " << nbFaces << " triangles" << endl;
 157   cout << endl;
 158
 159   // NB_ELEMS DUMMY_INT
 160   theFile << space << nbFaces << space << dummyint << endl;
 161
 162   // Loop from 1 to NB_ELEMS
 163
 164   f = faces.begin();
 165   for ( ; f != faces.end(); ++f )
 166   {
 167     // NB_NODES PER FACE
 168     elem = *f;
 169     nbNodes = elem->NbNodes();
 170     theFile << space << nbNodes;
 171
 172     // NODE_NB_1 NODE_NB_2 ...
 173     nodeIt = elem->nodesIterator();
 174     while ( nodeIt->more() )
 175     {
 176       // find GHS3D ID
 177       aSmdsID = nodeIt->next()->GetID();
 178       it = theSmdsToGhs3dIdMap.find( aSmdsID );
 179       ASSERT( it != theSmdsToGhs3dIdMap.end() );
 180       theFile << space << (*it).second;
 181     }
 182
 183     // (NB_NODES + 1) times: DUMMY_INT
 184     for ( int i=0; i<=nbNodes; i++)
 185       theFile << space << dummyint;
 186
 187     theFile << endl;
 188   }
 189
 190   return true;
 191 }
 192
 193 //=======================================================================
 194 //function : writePoints
 195 //purpose  :
 196 //=======================================================================
 197
 198 static bool writePoints (ofstream &                       theFile,
 199                          SMESHDS_Mesh *                   theMesh,
 200                          map <int,int> &                  theSmdsToGhs3dIdMap,
 201                          map <int,const SMDS_MeshNode*> & theGhs3dIdToNodeMap)
 202 {
 203   // record structure:
 204   //
 205   // NB_NODES
 206   // Loop from 1 to NB_NODES
 207   //   X Y Z DUMMY_INT
 208
 209   int nbNodes = theMesh->NbNodes();
 210   if ( nbNodes == 0 )
 211     return false;
 212
 213   const char* space    = "  ";
 214   const int   dummyint = 0;
 215
 216   int aGhs3dID = 1;
 217   SMDS_NodeIteratorPtr it = theMesh->nodesIterator();
 218   const SMDS_MeshNode* node;
 219
 220   // NB_NODES
 221   theFile << space << nbNodes << endl;
 222   cout << "The 2D mesh contains :" << endl;
 223   cout << "  " << nbNodes << " nodes" << endl;
 224
 225   // Loop from 1 to NB_NODES
 226
 227   while ( it->more() )
 228   {
 229     node = it->next();
 230     theSmdsToGhs3dIdMap.insert( map <int,int>::value_type( node->GetID(), aGhs3dID ));
 231     theGhs3dIdToNodeMap.insert (map <int,const SMDS_MeshNode*>::value_type( aGhs3dID, node ));
 232     aGhs3dID++;
 233
 234     // X Y Z DUMMY_INT
 235     theFile
 236       << space << node->X()
 237       << space << node->Y()
 238       << space << node->Z()
 239       << space << dummyint;
 240
 241     theFile << endl;
 242   }
 243
 244   return true;
 245 }
 246
 247 //=======================================================================
 248 //function : findSolid
 249 //purpose  :
 250 //=======================================================================
 251
 252 static TopoDS_Shape findSolid(const SMDS_MeshNode *aNode[],
 253                               TopoDS_Shape        aSolid,
 254                               const TopoDS_Shape  shape[],
 255                               const double        box[][6],
 256                               const int           nShape) {
 257
 258   Standard_Real PX, PY, PZ;
 259   int iShape;
 260
 261   PX = ( aNode[0]->X() + aNode[1]->X() + aNode[2]->X() + aNode[3]->X() ) / 4.0;
 262   PY = ( aNode[0]->Y() + aNode[1]->Y() + aNode[2]->Y() + aNode[3]->Y() ) / 4.0;
 263   PZ = ( aNode[0]->Z() + aNode[1]->Z() + aNode[2]->Z() + aNode[3]->Z() ) / 4.0;
 264   gp_Pnt aPnt(PX, PY, PZ);
 265
 266   BRepClass3d_SolidClassifier SC (aSolid, aPnt, Precision::Confusion());
 267   if ( not(SC.State() == TopAbs_IN) ) {
 268     for (iShape = 0; iShape < nShape; iShape++) {
 269       aSolid = shape[iShape];
 270       if ( not( PX < box[iShape][0] || box[iShape][1] < PX ||
 271                 PY < box[iShape][2] || box[iShape][3] < PY ||
 272                 PZ < box[iShape][4] || box[iShape][5] < PZ) ) {
 273         BRepClass3d_SolidClassifier SC (aSolid, aPnt, Precision::Confusion());
 274         if (SC.State() == TopAbs_IN)
 275           break;
 276       }
 277     }
 278   }
 279   return aSolid;
 280 }
 281
 282 //=======================================================================
 283 //function : readMapIntLine
 284 //purpose  :
 285 //=======================================================================
 286
 287 static char* readMapIntLine(char* ptr, int tab[]) {
 288   char* ptrRet = "\n";
 289   long int intVal;
 290   int i = 1;
 291   cout << endl;
 292
 293   while ( *ptr != *ptrRet ) {
 294     intVal = strtol(ptr, &ptr, 10);
 295     if ( i < 4 )
 296       tab[i-1] = intVal;
 297     i++;
 298   }
 299   return ptr;
 300 }
 301
 302 //=======================================================================
 303 //function : readResultFile
 304 //purpose  :
 305 //=======================================================================
 306
 307 static bool readResultFile(const int                       fileOpen,
 308                            SMESHDS_Mesh*                   theMeshDS,
 309                            TopoDS_Shape                    tabShape[],
 310                            double                          tabBox[][6],
 311                            const int                       nShape,
 312                            map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap) {
 313
 314   struct stat  status;
 315   size_t       length;
 316
 317   char *ptr;
 318   char *tetraPtr;
 319   char *shapePtr;
 320
 321   int fileStat, fileClose;
 322   int nbElems, nbNodes, nbInputNodes;
 323   int nodeId, triangleId;
 324   int tab[3], tabID[nShape];
 325   int nbTriangle;
 326   int ID, shapeID, ghs3dShapeID;
 327
 328   double coord [3];
 329
 330   TopoDS_Shape aSolid;
 331   SMDS_MeshNode * aNewNode;
 332   const SMDS_MeshNode * node[4];
 333   map <int,const SMDS_MeshNode*>::iterator IdNode;
 334   SMDS_MeshElement* aTet;
 335
 336   for (int i=0; i<nShape; i++)
 337     tabID[i] = 0;
 338
 339   // Read the file state
 340   fileStat = fstat(fileOpen, &status);
 341   length   = status.st_size;
 342
 343   // Mapping the result file into memory
 344   ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
 345   fileClose = close(fileOpen);
 346
 347   ptr      = readMapIntLine(ptr, tab);
 348   tetraPtr = ptr;
 349
 350   nbElems      = tab[0];
 351   nbNodes      = tab[1];
 352   nbInputNodes = tab[2];
 353
 354   // Reading the nodeId
 355   for (int i=0; i < 4*nbElems; i++)
 356     nodeId = strtol(ptr, &ptr, 10);
 357
 358   // Reading the nodeCoor and update the nodeMap
 359   for (int iNode=0; iNode < nbNodes; iNode++) {
 360     for (int iCoor=0; iCoor < 3; iCoor++)
 361       coord[ iCoor ] = strtod(ptr, &ptr);
 362     if ((iNode+1) > nbInputNodes) {
 363       aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
 364       theGhs3dIdToNodeMap.insert(make_pair( (iNode+1), aNewNode ));
 365     }
 366   }
 367
 368   // Reading the triangles
 369   nbTriangle = strtol(ptr, &ptr, 10);
 370
 371   for (int i=0; i < 3*nbTriangle; i++)
 372     triangleId = strtol(ptr, &ptr, 10);
 373
 374   shapePtr = ptr;
 375
 376   // Associating the tetrahedrons with the shapes
 377   for (int iElem = 0; iElem < nbElems; iElem++) {
 378     for (int iNode = 0; iNode < 4; iNode++) {
 379       ID = strtol(tetraPtr, &tetraPtr, 10);
 380       IdNode = theGhs3dIdToNodeMap.find(ID);
 381       node[ iNode ] = IdNode->second;
 382     }
 383     aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
 384     ghs3dShapeID = strtol(shapePtr, &shapePtr, 10);
 385     if ( tabID[ ghs3dShapeID - 1 ] == 0 ) {
 386       if (iElem == 0)
 387         aSolid = tabShape[0];
 388       aSolid = findSolid(node, aSolid, tabShape, tabBox, nbTriangle);
 389       shapeID = theMeshDS->ShapeToIndex( aSolid );
 390       tabID[ ghs3dShapeID - 1] = shapeID;
 391     }
 392     else
 393       shapeID = tabID[ ghs3dShapeID - 1];
 394     theMeshDS->SetMeshElementOnShape( aTet, shapeID );
 395   }
 396   return true;
 397 }
 398
 399 //=======================================================================
 400 //function : getTmpDir
 401 //purpose  :
 402 //=======================================================================
 403
 404 static TCollection_AsciiString getTmpDir()
 405 {
 406   TCollection_AsciiString aTmpDir;
 407
 408   char *Tmp_dir = getenv("SALOME_TMP_DIR");
 409   if(Tmp_dir != NULL) {
 410     aTmpDir = Tmp_dir;
 411 #ifdef WIN32
 412     if(aTmpDir.Value(aTmpDir.Length()) != '\\') aTmpDir+='\\';
 413 #else
 414     if(aTmpDir.Value(aTmpDir.Length()) != '/') aTmpDir+='/';
 415 #endif
 416   }
 417   else {
 418 #ifdef WIN32
 419     aTmpDir = TCollection_AsciiString("C:\\");
 420 #else
 421     aTmpDir = TCollection_AsciiString("/tmp/");
 422 #endif
 423   }
 424   return aTmpDir;
 425 }
 426
 427 //=============================================================================
 428 /*!
 429  *Here we are going to use the GHS3D mesher
 430  */
 431 //=============================================================================
 432
 433 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
 434                                 const TopoDS_Shape& theShape)
 435 {
 436   bool Ok;
 437   SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
 438
 439   if (_iShape == 0 && _nbShape == 0) {
 440     cout << endl;
 441     cout << "Ghs3d execution..." << endl;
 442     cout << endl;
 443
 444     TopExp_Explorer exp (meshDS->ShapeToMesh(), TopAbs_SOLID);
 445     for (; exp.More(); exp.Next())
 446       _nbShape++;
 447   }
 448
 449   _iShape++;
 450
 451   if ( _iShape == _nbShape ) {
 452
 453     // create bounding box for every shape
 454
 455     int iShape = 0;
 456     TopoDS_Shape tabShape[_nbShape];
 457     double tabBox[_nbShape][6];
 458     Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
 459
 460     TopExp_Explorer expBox (meshDS->ShapeToMesh(), TopAbs_SOLID);
 461     for (; expBox.More(); expBox.Next()) {
 462       tabShape[iShape] = expBox.Current();
 463       Bnd_Box BoundingBox;
 464       BRepBndLib::Add(expBox.Current(), BoundingBox);
 465       BoundingBox.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
 466       tabBox[iShape][0] = Xmin; tabBox[iShape][1] = Xmax;
 467       tabBox[iShape][2] = Ymin; tabBox[iShape][3] = Ymax;
 468       tabBox[iShape][4] = Zmin; tabBox[iShape][5] = Zmax;
 469       iShape++;
 470     }
 471
 472     // make a unique working file name
 473     // to avoid access to the same files by eg different users
 474
 475     TCollection_AsciiString aGenericName, aTmpDir = getTmpDir();
 476     aGenericName = aTmpDir + "GHS3D_";
 477 #ifdef WIN32
 478     aGenericName += GetCurrentProcessId();
 479 #else
 480     aGenericName += getpid();
 481 #endif
 482     aGenericName += "_";
 483     aGenericName += meshDS->ShapeToIndex( theShape );
 484
 485     TCollection_AsciiString aFacesFileName, aPointsFileName, aResultFileName;
 486     TCollection_AsciiString aBadResFileName, aBbResFileName, aLogFileName;
 487     aFacesFileName  = aGenericName + ".faces";  // in faces
 488     aPointsFileName = aGenericName + ".points"; // in points
 489     aResultFileName = aGenericName + ".noboite";// out points and volumes
 490     aBadResFileName = aGenericName + ".boite";  // out bad result
 491     aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
 492     aLogFileName    = aGenericName + ".log";    // log
 493
 494     // -----------------
 495     // make input files
 496     // -----------------
 497
 498     ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
 499     ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out);
 500     // bool Ok =
 501     Ok =
 502 #ifdef WIN32
 503       aFacesFile->is_open() && aPointsFile->is_open();
 504 #else
 505       aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
 506 #endif
 507     if (!Ok) {
 508       INFOS( "Can't write into " << aTmpDir.ToCString());
 509       return false;
 510     }
 511     map <int,int> aSmdsToGhs3dIdMap;
 512     map <int,const SMDS_MeshNode*> aGhs3dIdToNodeMap;
 513
 514     Ok = writePoints( aPointsFile, meshDS, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap ) &&
 515          writeFaces ( aFacesFile,  meshDS, aSmdsToGhs3dIdMap );
 516
 517     aFacesFile.close();
 518     aPointsFile.close();
 519
 520     if ( ! Ok ) {
 521       if ( !getenv("GHS3D_KEEP_FILES") ) {
 522         OSD_File( aFacesFileName ).Remove();
 523         OSD_File( aPointsFileName ).Remove();
 524       }
 525       return false;
 526     }
 527
 528     // -----------------
 529     // run ghs3d mesher              WIN32???
 530     // -----------------
 531
 532     // ghs3d need to know amount of memory it may use (MB).
 533     // Default memory is defined at ghs3d installation but it may be not enough,
 534     // so allow to use about all available memory
 535
 536     TCollection_AsciiString memory;
 537 #ifndef WIN32
 538     struct sysinfo si;
 539     int err = sysinfo( &si );
 540     if ( err == 0 ) {
 541       int freeMem = si.totalram * si.mem_unit / 1024 / 1024;
 542       memory = "-m ";
 543       memory += int( 0.7 * freeMem );
 544     }
 545 #endif
 546
 547     MESSAGE("GHS3DPlugin_GHS3D::Compute");
 548     TCollection_AsciiString cmd( "ghs3d " ); // command to run
 549     cmd +=
 550       memory +                     // memory
 551       " -c0 -f " + aGenericName +  // file to read
 552            " 1>" + aLogFileName;   // dump into file
 553
 554     system( cmd.ToCString() ); // run
 555
 556     cout << endl;
 557     cout << "End of Ghs3d execution !" << endl;
 558
 559     // --------------
 560     // read a result
 561     // --------------
 562
 563     // Mapping the result file
 564
 565     int fileOpen;
 566     fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
 567     if ( fileOpen < 0 ) {
 568       cout << endl;
 569       cout << "Error when opening the " << aResultFileName.ToCString() << " file" << endl;
 570       cout << endl;
 571       Ok = false;
 572     }
 573     else
 574       Ok = readResultFile( fileOpen, meshDS, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap );
 575
 576     // ---------------------
 577     // remove working files
 578     // ---------------------
 579
 580     if ( Ok ) {
 581       OSD_File( aLogFileName ).Remove();
 582     }
 583     else if ( OSD_File( aLogFileName ).Size() > 0 ) {
 584       INFOS( "GHS3D Error: see " << aLogFileName.ToCString() );
 585     }
 586     else {
 587       OSD_File( aLogFileName ).Remove();
 588       INFOS( "GHS3D Error: command '" << cmd.ToCString() << "' failed" );
 589     }
 590
 591     if ( !getenv("GHS3D_KEEP_FILES") ) {
 592       OSD_File( aFacesFileName ).Remove();
 593       OSD_File( aPointsFileName ).Remove();
 594       OSD_File( aResultFileName ).Remove();
 595       OSD_File( aBadResFileName ).Remove();
 596       OSD_File( aBbResFileName ).Remove();
 597     }
 598     if ( _iShape == _nbShape ) {
 599       cout << "Output file " << aResultFileName.ToCString();
 600       if ( !Ok )
 601         cout << " not treated !" << endl;
 602       else
 603         cout << " treated !" << endl;
 604       cout << endl;
 605     }
 606   }
 607   return Ok;
 608 }
 609
 610
 611 //=============================================================================
 612 /*!
 613  *
 614  */
 615 //=============================================================================
 616
 617 ostream & GHS3DPlugin_GHS3D::SaveTo(ostream & save)
 618 {
 619   return save;
 620 }
 621
 622 //=============================================================================
 623 /*!
 624  *
 625  */
 626 //=============================================================================
 627
 628 istream & GHS3DPlugin_GHS3D::LoadFrom(istream & load)
 629 {
 630   return load;
 631 }
 632
 633 //=============================================================================
 634 /*!
 635  *
 636  */
 637 //=============================================================================
 638
 639 ostream & operator << (ostream & save, GHS3DPlugin_GHS3D & hyp)
 640 {
 641   return hyp.SaveTo( save );
 642 }
 643
 644 //=============================================================================
 645 /*!
 646  *
 647  */
 648 //=============================================================================
 649
 650 istream & operator >> (istream & load, GHS3DPlugin_GHS3D & hyp)
 651 {
 652   return hyp.LoadFrom( load );
 653 }