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