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