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 #include "SMESH_Comment.hxx"
  33 #include "SMESH_MesherHelper.hxx"
  34
  35 #include "SMDS_MeshElement.hxx"
  36 #include "SMDS_MeshNode.hxx"
  37
  38 #include <TopExp.hxx>
  39 #include <TopExp_Explorer.hxx>
  40 #include <OSD_File.hxx>
  41
  42 #include "utilities.h"
  43
  44 #ifndef WIN32
  45 #include <sys/sysinfo.h>
  46 #endif
  47
  48 //#include <Standard_Stream.hxx>
  49
  50 #include <BRepGProp.hxx>
  51 #include <BRepBndLib.hxx>
  52 #include <BRepClass_FaceClassifier.hxx>
  53 #include <BRepClass3d_SolidClassifier.hxx>
  54 #include <TopAbs.hxx>
  55 #include <Bnd_Box.hxx>
  56 #include <GProp_GProps.hxx>
  57 #include <Precision.hxx>
  58
  59 #define castToNode(n) static_cast<const SMDS_MeshNode *>( n );
  60
  61 #ifdef _DEBUG_
  62 #define DUMP(txt) \
  63 //  cout << txt
  64 #else
  65 #define DUMP(txt)
  66 #endif
  67
  68 extern "C"
  69 {
  70 #ifndef WNT
  71 #include <unistd.h>
  72 #include <sys/mman.h>
  73 #endif
  74 #include <sys/stat.h>
  75 #include <fcntl.h>
  76 }
  77
  78 //=============================================================================
  79 /*!
  80  *
  81  */
  82 //=============================================================================
  83
  84 GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D(int hypId, int studyId, SMESH_Gen* gen)
  85   : SMESH_3D_Algo(hypId, studyId, gen)
  86 {
  87   MESSAGE("GHS3DPlugin_GHS3D::GHS3DPlugin_GHS3D");
  88   _name = "GHS3D_3D";
  89   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  90   _onlyUnaryInput = false; // Compute() will be called on a compound of solids
  91   _iShape=0;
  92   _nbShape=0;
  93 }
  94
  95 //=============================================================================
  96 /*!
  97  *
  98  */
  99 //=============================================================================
 100
 101 GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D()
 102 {
 103   MESSAGE("GHS3DPlugin_GHS3D::~GHS3DPlugin_GHS3D");
 104 }
 105
 106 //=============================================================================
 107 /*!
 108  *
 109  */
 110 //=============================================================================
 111
 112 bool GHS3DPlugin_GHS3D::CheckHypothesis ( SMESH_Mesh&                          aMesh,
 113                                           const TopoDS_Shape&                  aShape,
 114                                           SMESH_Hypothesis::Hypothesis_Status& aStatus )
 115 {
 116 //  MESSAGE("GHS3DPlugin_GHS3D::CheckHypothesis");
 117   aStatus = SMESH_Hypothesis::HYP_OK;
 118   return true;
 119 }
 120
 121 //=======================================================================
 122 //function : findShape
 123 //purpose  :
 124 //=======================================================================
 125
 126 static TopoDS_Shape findShape(const SMDS_MeshNode *aNode[],
 127                               TopoDS_Shape        aShape,
 128                               const TopoDS_Shape  shape[],
 129                               double**            box,
 130                               const int           nShape) {
 131   double *pntCoor;
 132   int iShape, nbNode = 4;
 133
 134   pntCoor = new double[3];
 135   for ( int i=0; i<3; i++ ) {
 136     pntCoor[i] = 0;
 137     for ( int j=0; j<nbNode; j++ ) {
 138       if ( i == 0) pntCoor[i] += aNode[j]->X();
 139       if ( i == 1) pntCoor[i] += aNode[j]->Y();
 140       if ( i == 2) pntCoor[i] += aNode[j]->Z();
 141     }
 142     pntCoor[i] /= nbNode;
 143   }
 144
 145   gp_Pnt aPnt(pntCoor[0], pntCoor[1], pntCoor[2]);
 146   BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
 147   if ( not(SC.State() == TopAbs_IN) ) {
 148     for (iShape = 0; iShape < nShape; iShape++) {
 149       aShape = shape[iShape];
 150       if ( not( pntCoor[0] < box[iShape][0] || box[iShape][1] < pntCoor[0] ||
 151                 pntCoor[1] < box[iShape][2] || box[iShape][3] < pntCoor[1] ||
 152                 pntCoor[2] < box[iShape][4] || box[iShape][5] < pntCoor[2]) ) {
 153         BRepClass3d_SolidClassifier SC (aShape, aPnt, Precision::Confusion());
 154         if (SC.State() == TopAbs_IN)
 155           break;
 156       }
 157     }
 158   }
 159   return aShape;
 160 }
 161
 162 //=======================================================================
 163 //function : readMapIntLine
 164 //purpose  :
 165 //=======================================================================
 166
 167 static char* readMapIntLine(char* ptr, int tab[]) {
 168   long int intVal;
 169   cout << endl;
 170
 171   for ( int i=0; i<17; i++ ) {
 172     intVal = strtol(ptr, &ptr, 10);
 173     if ( i < 3 )
 174       tab[i] = intVal;
 175   }
 176   return ptr;
 177 }
 178
 179 //=======================================================================
 180 //function : readLine
 181 //purpose  :
 182 //=======================================================================
 183
 184 #define GHS3DPlugin_BUFLENGTH 256
 185 #define GHS3DPlugin_ReadLine(aPtr,aBuf,aFile,aLineNb) \
 186 {  aPtr = fgets( aBuf, GHS3DPlugin_BUFLENGTH - 2, aFile ); aLineNb++; DUMP(endl); }
 187
 188 //=======================================================================
 189 //function : countShape
 190 //purpose  :
 191 //=======================================================================
 192
 193 template < class Mesh, class Shape >
 194 static int countShape( Mesh* mesh, Shape shape ) {
 195   TopExp_Explorer expShape ( mesh->ShapeToMesh(), shape );
 196   int nbShape = 0;
 197   for ( ; expShape.More(); expShape.Next() ) {
 198       nbShape++;
 199   }
 200   return nbShape;
 201 }
 202
 203 //=======================================================================
 204 //function : writeFaces
 205 //purpose  :
 206 //=======================================================================
 207
 208 static bool writeFaces (ofstream &            theFile,
 209                         SMESHDS_Mesh *        theMesh,
 210                         const map <int,int> & theSmdsToGhs3dIdMap)
 211 {
 212   // record structure:
 213   //
 214   // NB_ELEMS DUMMY_INT
 215   // Loop from 1 to NB_ELEMS
 216   // NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
 217
 218   int nbShape = countShape( theMesh, TopAbs_FACE );
 219
 220   int *tabID;             tabID    = new int[nbShape];
 221   TopoDS_Shape *tabShape; tabShape = new TopoDS_Shape[nbShape];
 222   TopoDS_Shape aShape;
 223   SMESHDS_SubMesh* theSubMesh;
 224   const SMDS_MeshElement* aFace;
 225   const char* space    = "  ";
 226   const int   dummyint = 0;
 227   map<int,int>::const_iterator itOnMap;
 228   SMDS_ElemIteratorPtr itOnSubMesh, itOnSubFace;
 229   int shapeID, nbNodes, aSmdsID;
 230   bool idFound;
 231
 232   cout << "    " << theMesh->NbFaces() << " shapes of 2D dimension" << endl;
 233   cout << endl;
 234
 235   theFile << space << theMesh->NbFaces() << space << dummyint << endl;
 236
 237   TopExp_Explorer expface( theMesh->ShapeToMesh(), TopAbs_FACE );
 238   for ( int i = 0; expface.More(); expface.Next(), i++ ) {
 239     tabID[i] = 0;
 240     aShape   = expface.Current();
 241     shapeID  = theMesh->ShapeToIndex( aShape );
 242     idFound  = false;
 243     for ( int j=0; j<=i; j++) {
 244       if ( shapeID == tabID[j] ) {
 245         idFound = true;
 246         break;
 247       }
 248     }
 249     if ( not idFound ) {
 250       tabID[i]    = shapeID;
 251       tabShape[i] = aShape;
 252     }
 253   }
 254   for ( int i =0; i < nbShape; i++ ) {
 255     if ( not (tabID[i] == 0) ) {
 256       aShape      = tabShape[i];
 257       shapeID     = tabID[i];
 258       theSubMesh  = theMesh->MeshElements( aShape );
 259       itOnSubMesh = theSubMesh->GetElements();
 260       while ( itOnSubMesh->more() ) {
 261         aFace   = itOnSubMesh->next();
 262         nbNodes = aFace->NbNodes();
 263
 264         theFile << space << nbNodes;
 265
 266         itOnSubFace = aFace->nodesIterator();
 267         while ( itOnSubFace->more() ) {
 268           // find GHS3D ID
 269           aSmdsID = itOnSubFace->next()->GetID();
 270           itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
 271           ASSERT( itOnMap != theSmdsToGhs3dIdMap.end() );
 272
 273           theFile << space << (*itOnMap).second;
 274         }
 275
 276         // (NB_NODES + 1) times: DUMMY_INT
 277         for ( int j=0; j<=nbNodes; j++)
 278           theFile << space << dummyint;
 279
 280         theFile << endl;
 281       }
 282     }
 283   }
 284
 285   delete [] tabID;
 286   delete [] tabShape;
 287
 288   return true;
 289 }
 290
 291 //=======================================================================
 292 //function : writeFaces
 293 //purpose  : Write Faces in case if generate 3D mesh w/o geometry
 294 //=======================================================================
 295
 296 static bool writeFaces (ofstream &            theFile,
 297                         SMESHDS_Mesh *        theMesh,
 298                         vector <const SMDS_MeshNode*> & theNodeByGhs3dId)
 299 {
 300   // record structure:
 301   //
 302   // NB_ELEMS DUMMY_INT
 303   // Loop from 1 to NB_ELEMS
 304   //   NB_NODES NODE_NB_1 NODE_NB_2 ... (NB_NODES + 1) times: DUMMY_INT
 305
 306
 307   int nbFaces = 0;
 308   list< const SMDS_MeshElement* > faces;
 309   list< const SMDS_MeshElement* >::iterator f;
 310   map< const SMDS_MeshNode*,int >::iterator it;
 311   SMDS_ElemIteratorPtr nodeIt;
 312   const SMDS_MeshElement* elem;
 313   int nbNodes;
 314
 315   const char* space    = "  ";
 316   const int   dummyint = 0;
 317
 318   //get all faces from mesh
 319   SMDS_FaceIteratorPtr eIt = theMesh->facesIterator();
 320   while ( eIt->more() ) {
 321     const SMDS_MeshElement* elem = eIt->next();
 322     if ( !elem )
 323       return false;
 324     faces.push_back( elem );
 325     nbFaces++;
 326   }
 327
 328   if ( nbFaces == 0 )
 329     return false;
 330
 331   cout << "The initial 2D mesh contains " << nbFaces << " faces and ";
 332
 333   // NB_ELEMS DUMMY_INT
 334   theFile << space << nbFaces << space << dummyint << endl;
 335
 336   // Loop from 1 to NB_ELEMS
 337
 338   map<const SMDS_MeshNode*,int> aNodeToGhs3dIdMap;
 339   f = faces.begin();
 340   for ( ; f != faces.end(); ++f )
 341   {
 342     // NB_NODES PER FACE
 343     elem = *f;
 344     nbNodes = elem->NbNodes();
 345     theFile << space << nbNodes;
 346
 347     // NODE_NB_1 NODE_NB_2 ...
 348     nodeIt = elem->nodesIterator();
 349     while ( nodeIt->more() )
 350     {
 351       // find GHS3D ID
 352       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
 353       int newId = aNodeToGhs3dIdMap.size() + 1; // ghs3d ids count from 1
 354       it = aNodeToGhs3dIdMap.insert( make_pair( node, newId )).first;
 355       theFile << space << it->second;
 356     }
 357
 358     // (NB_NODES + 1) times: DUMMY_INT
 359     for ( int i=0; i<=nbNodes; i++)
 360       theFile << space << dummyint;
 361
 362     theFile << endl;
 363   }
 364
 365   // put nodes to theNodeByGhs3dId vector
 366   theNodeByGhs3dId.resize( aNodeToGhs3dIdMap.size() );
 367   map<const SMDS_MeshNode*,int>::const_iterator n2id = aNodeToGhs3dIdMap.begin();
 368   for ( ; n2id != aNodeToGhs3dIdMap.end(); ++ n2id)
 369   {
 370     theNodeByGhs3dId[ n2id->second - 1 ] = n2id->first; // ghs3d ids count from 1
 371   }
 372
 373   return true;
 374 }
 375
 376 //=======================================================================
 377 //function : writePoints
 378 //purpose  :
 379 //=======================================================================
 380
 381 static bool writePoints (ofstream &                       theFile,
 382                          SMESHDS_Mesh *                   theMesh,
 383                          map <int,int> &                  theSmdsToGhs3dIdMap,
 384                          map <int,const SMDS_MeshNode*> & theGhs3dIdToNodeMap)
 385 {
 386   // record structure:
 387   //
 388   // NB_NODES
 389   // Loop from 1 to NB_NODES
 390   //   X Y Z DUMMY_INT
 391
 392   int nbNodes = theMesh->NbNodes();
 393   if ( nbNodes == 0 )
 394     return false;
 395
 396   const char* space    = "  ";
 397   const int   dummyint = 0;
 398
 399   int aGhs3dID = 1;
 400   SMDS_NodeIteratorPtr it = theMesh->nodesIterator();
 401   const SMDS_MeshNode* node;
 402
 403   // NB_NODES
 404   theFile << space << nbNodes << endl;
 405   cout << endl;
 406   cout << "The initial 2D mesh contains :" << endl;
 407   cout << "    " << nbNodes << " vertices" << endl;
 408
 409   // Loop from 1 to NB_NODES
 410
 411   while ( it->more() )
 412   {
 413     node = it->next();
 414     theSmdsToGhs3dIdMap.insert( map <int,int>::value_type( node->GetID(), aGhs3dID ));
 415     theGhs3dIdToNodeMap.insert (map <int,const SMDS_MeshNode*>::value_type( aGhs3dID, node ));
 416     aGhs3dID++;
 417
 418     // X Y Z DUMMY_INT
 419     theFile
 420       << space << node->X()
 421       << space << node->Y()
 422       << space << node->Z()
 423       << space << dummyint;
 424
 425     theFile << endl;
 426   }
 427
 428   return true;
 429 }
 430
 431 //=======================================================================
 432 //function : writePoints
 433 //purpose  :
 434 //=======================================================================
 435
 436 static bool writePoints (ofstream &                            theFile,
 437                          SMESHDS_Mesh *                        theMesh,
 438                          const vector <const SMDS_MeshNode*> & theNodeByGhs3dId)
 439 {
 440   // record structure:
 441   //
 442   // NB_NODES
 443   // Loop from 1 to NB_NODES
 444   //   X Y Z DUMMY_INT
 445
 446   //int nbNodes = theMesh->NbNodes();
 447   int nbNodes = theNodeByGhs3dId.size();
 448   if ( nbNodes == 0 )
 449     return false;
 450
 451   const char* space    = "  ";
 452   const int   dummyint = 0;
 453
 454   const SMDS_MeshNode* node;
 455
 456   // NB_NODES
 457   theFile << space << nbNodes << endl;
 458   cout << nbNodes << " nodes" << endl;
 459
 460   // Loop from 1 to NB_NODES
 461
 462   vector<const SMDS_MeshNode*>::const_iterator nodeIt = theNodeByGhs3dId.begin();
 463   vector<const SMDS_MeshNode*>::const_iterator after  = theNodeByGhs3dId.end();
 464   for ( ; nodeIt != after; ++nodeIt )
 465   {
 466     node = *nodeIt;
 467
 468     // X Y Z DUMMY_INT
 469     theFile
 470       << space << node->X()
 471       << space << node->Y()
 472       << space << node->Z()
 473       << space << dummyint;
 474
 475     theFile << endl;
 476   }
 477
 478   return true;
 479 }
 480
 481 //=======================================================================
 482 //function : readResultFile
 483 //purpose  :
 484 //=======================================================================
 485
 486 static bool readResultFile(const int                       fileOpen,
 487                            SMESHDS_Mesh*                   theMeshDS,
 488                            TopoDS_Shape                    tabShape[],
 489                            double**                        tabBox,
 490                            const int                       nbShape,
 491                            map <int,const SMDS_MeshNode*>& theGhs3dIdToNodeMap) {
 492
 493   struct stat status;
 494   size_t      length;
 495
 496   char *ptr, *mapPtr;
 497   char *tetraPtr;
 498   char *shapePtr;
 499
 500   int fileStat;
 501   int nbElems, nbNodes, nbInputNodes;
 502   int nodeId, triangleId;
 503   int nbTriangle;
 504   int ID, shapeID, ghs3dShapeID;
 505   int IdShapeRef = 1;
 506   int compoundID = theMeshDS->ShapeToIndex( theMeshDS->ShapeToMesh() );
 507
 508   int *tab, *tabID, *nodeID, *nodeAssigne;
 509   double *coord;
 510   const SMDS_MeshNode **node;
 511
 512   tab    = new int[3];
 513   tabID  = new int[nbShape];
 514   nodeID = new int[4];
 515   coord  = new double[3];
 516   node   = new const SMDS_MeshNode*[4];
 517
 518   TopoDS_Shape aSolid;
 519   SMDS_MeshNode * aNewNode;
 520   map <int,const SMDS_MeshNode*>::iterator itOnNode;
 521   SMDS_MeshElement* aTet;
 522
 523   for (int i=0; i<nbShape; i++)
 524     tabID[i] = 0;
 525
 526   // Read the file state
 527   fileStat = fstat(fileOpen, &status);
 528   length   = status.st_size;
 529
 530   // Mapping the result file into memory
 531   ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
 532   mapPtr = ptr;
 533
 534   ptr      = readMapIntLine(ptr, tab);
 535   tetraPtr = ptr;
 536
 537   nbElems      = tab[0];
 538   nbNodes      = tab[1];
 539   nbInputNodes = tab[2];
 540
 541   nodeAssigne = new int[ nbNodes+1 ];
 542
 543   // Reading the nodeId
 544   for (int i=0; i < 4*nbElems; i++)
 545     nodeId = strtol(ptr, &ptr, 10);
 546
 547   // Reading the nodeCoor and update the nodeMap
 548   for (int iNode=0; iNode < nbNodes; iNode++) {
 549     for (int iCoor=0; iCoor < 3; iCoor++)
 550       coord[ iCoor ] = strtod(ptr, &ptr);
 551     nodeAssigne[ iNode+1 ] = 1;
 552     if ( (iNode+1) > nbInputNodes ) {
 553       nodeAssigne[ iNode+1 ] = 0;
 554       aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
 555       theGhs3dIdToNodeMap.insert(make_pair( (iNode+1), aNewNode ));
 556     }
 557   }
 558
 559   // Reading the number of triangles which corresponds to the number of shapes
 560   nbTriangle = strtol(ptr, &ptr, 10);
 561
 562   for (int i=0; i < 3*nbShape; i++)
 563     triangleId = strtol(ptr, &ptr, 10);
 564
 565   shapePtr = ptr;
 566
 567   // Associating the tetrahedrons to the shapes
 568   shapeID = compoundID;
 569   for (int iElem = 0; iElem < nbElems; iElem++) {
 570     for (int iNode = 0; iNode < 4; iNode++) {
 571       ID = strtol(tetraPtr, &tetraPtr, 10);
 572       itOnNode = theGhs3dIdToNodeMap.find(ID);
 573       node[ iNode ] = itOnNode->second;
 574       nodeID[ iNode ] = ID;
 575     }
 576     aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
 577     if ( nbShape > 1 ) {
 578       ghs3dShapeID = strtol(shapePtr, &shapePtr, 10) - IdShapeRef;
 579       if ( tabID[ ghs3dShapeID ] == 0 ) {
 580         if (iElem == 0)
 581           aSolid = tabShape[0];
 582         aSolid = findShape(node, aSolid, tabShape, tabBox, nbShape);
 583         shapeID = theMeshDS->ShapeToIndex( aSolid );
 584         tabID[ ghs3dShapeID ] = shapeID;
 585       }
 586       else
 587         shapeID = tabID[ ghs3dShapeID ];
 588     }
 589     // set new nodes and tetrahedron on to the shape
 590     for ( int i=0; i<4; i++ ) {
 591       if ( nodeAssigne[ nodeID[i] ] == 0 )
 592         theMeshDS->SetNodeInVolume( node[i], shapeID );
 593     }
 594     theMeshDS->SetMeshElementOnShape( aTet, shapeID );
 595     if ( (iElem + 1) == nbElems )
 596       cout << nbElems << " tetrahedrons have been associated to " << nbShape << " shapes" << endl;
 597   }
 598   munmap(mapPtr, length);
 599   close(fileOpen);
 600
 601   delete [] tab;
 602   delete [] tabID;
 603   delete [] nodeID;
 604   delete [] coord;
 605   delete [] node;
 606
 607   return true;
 608 }
 609
 610 //=======================================================================
 611 //function : readResultFile
 612 //purpose  :
 613 //=======================================================================
 614
 615 static bool readResultFile(const int                      fileOpen,
 616                            SMESHDS_Mesh*                  theMeshDS,
 617                            TopoDS_Shape                   aSolid,
 618                            vector <const SMDS_MeshNode*>& theNodeByGhs3dId) {
 619
 620   struct stat  status;
 621   size_t       length;
 622
 623   char *ptr, *mapPtr;
 624   char *tetraPtr;
 625   char *shapePtr;
 626
 627   int fileStat;
 628   int nbElems, nbNodes, nbInputNodes;
 629   int nodeId, triangleId;
 630   int nbTriangle;
 631   int ID, shapeID;
 632
 633   int *tab;
 634   double *coord;
 635   const SMDS_MeshNode **node;
 636
 637   tab   = new int[3];
 638   coord = new double[3];
 639   node  = new const SMDS_MeshNode*[4];
 640
 641   SMDS_MeshNode * aNewNode;
 642   map <int,const SMDS_MeshNode*>::iterator IdNode;
 643   SMDS_MeshElement* aTet;
 644
 645   // Read the file state
 646   fileStat = fstat(fileOpen, &status);
 647   length   = status.st_size;
 648
 649   // Mapping the result file into memory
 650   ptr = (char *) mmap(0,length,PROT_READ,MAP_PRIVATE,fileOpen,0);
 651   mapPtr = ptr;
 652
 653   ptr      = readMapIntLine(ptr, tab);
 654   tetraPtr = ptr;
 655
 656   nbElems      = tab[0];
 657   nbNodes      = tab[1];
 658   nbInputNodes = tab[2];
 659
 660   theNodeByGhs3dId.resize( nbNodes );
 661
 662   // Reading the nodeId
 663   for (int i=0; i < 4*nbElems; i++)
 664     nodeId = strtol(ptr, &ptr, 10);
 665
 666   // Reading the nodeCoor and update the nodeMap
 667   shapeID = theMeshDS->ShapeToIndex( aSolid );
 668   for (int iNode=0; iNode < nbNodes; iNode++) {
 669     for (int iCoor=0; iCoor < 3; iCoor++)
 670       coord[ iCoor ] = strtod(ptr, &ptr);
 671     if ((iNode+1) > nbInputNodes) {
 672       aNewNode = theMeshDS->AddNode( coord[0],coord[1],coord[2] );
 673       theMeshDS->SetNodeInVolume( aNewNode, shapeID );
 674       theNodeByGhs3dId[ iNode ] = aNewNode;
 675     }
 676   }
 677
 678   // Reading the triangles
 679   nbTriangle = strtol(ptr, &ptr, 10);
 680
 681   for (int i=0; i < 3*nbTriangle; i++)
 682     triangleId = strtol(ptr, &ptr, 10);
 683
 684   shapePtr = ptr;
 685
 686   // Associating the tetrahedrons to the shapes
 687   for (int iElem = 0; iElem < nbElems; iElem++) {
 688     for (int iNode = 0; iNode < 4; iNode++) {
 689       ID = strtol(tetraPtr, &tetraPtr, 10);
 690       node[ iNode ] = theNodeByGhs3dId[ ID-1 ];
 691     }
 692     aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
 693     shapeID = theMeshDS->ShapeToIndex( aSolid );
 694     theMeshDS->SetMeshElementOnShape( aTet, shapeID );
 695   }
 696   if ( nbElems )
 697     cout << nbElems << " tetrahedrons have been associated to " << nbTriangle << " shapes" << endl;
 698   munmap(mapPtr, length);
 699   close(fileOpen);
 700
 701   delete [] tab;
 702   delete [] coord;
 703   delete [] node;
 704
 705   return true;
 706 }
 707
 708 //=======================================================================
 709 //function : getTmpDir
 710 //purpose  :
 711 //=======================================================================
 712
 713 static TCollection_AsciiString getTmpDir()
 714 {
 715   TCollection_AsciiString aTmpDir;
 716
 717   char *Tmp_dir = getenv("SALOME_TMP_DIR");
 718   if(Tmp_dir != NULL) {
 719     aTmpDir = Tmp_dir;
 720 #ifdef WIN32
 721     if(aTmpDir.Value(aTmpDir.Length()) != '\\') aTmpDir+='\\';
 722 #else
 723     if(aTmpDir.Value(aTmpDir.Length()) != '/') aTmpDir+='/';
 724 #endif
 725   }
 726   else {
 727 #ifdef WIN32
 728     aTmpDir = TCollection_AsciiString("C:\\");
 729 #else
 730     aTmpDir = TCollection_AsciiString("/tmp/");
 731 #endif
 732   }
 733   return aTmpDir;
 734 }
 735
 736 //================================================================================
 737 /*!
 738  * \brief Look for a line containing a text in a file
 739   * \retval bool - true if the line is found
 740  */
 741 //================================================================================
 742
 743 static bool findLineContaing(const TCollection_AsciiString& theText,
 744                              const TCollection_AsciiString& theFile,
 745                              TCollection_AsciiString &      theFoundLine)
 746 {
 747   bool found = false;
 748   if ( FILE * aFile = fopen( theFile.ToCString(), "r" ))
 749   {
 750     char * aPtr;
 751     char aBuffer[ GHS3DPlugin_BUFLENGTH ];
 752     int aLineNb = 0;
 753     do {
 754       GHS3DPlugin_ReadLine( aPtr, aBuffer, aFile, aLineNb );
 755       if ( aPtr ) {
 756         theFoundLine = aPtr;
 757         found = theFoundLine.Search( theText ) >= 0;
 758       }
 759     } while ( aPtr && !found );
 760
 761     fclose( aFile );
 762   }
 763   return found;
 764 }
 765
 766 //=============================================================================
 767 /*!
 768  *Here we are going to use the GHS3D mesher
 769  */
 770 //=============================================================================
 771
 772 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
 773                                 const TopoDS_Shape& theShape)
 774 {
 775   bool Ok(false);
 776   SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
 777
 778   _nbShape = countShape( meshDS, TopAbs_SOLID ); // we count the number of shapes
 779
 780   // create bounding box for every shape inside the compound
 781
 782   int iShape = 0;
 783   TopoDS_Shape* tabShape;
 784   double**      tabBox;
 785   tabShape = new TopoDS_Shape[_nbShape];
 786   tabBox   = new double*[_nbShape];
 787   for (int i=0; i<_nbShape; i++)
 788     tabBox[i] = new double[6];
 789   Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
 790
 791   TopExp_Explorer expBox (meshDS->ShapeToMesh(), TopAbs_SOLID);
 792   for (; expBox.More(); expBox.Next()) {
 793     tabShape[iShape] = expBox.Current();
 794     Bnd_Box BoundingBox;
 795     BRepBndLib::Add(expBox.Current(), BoundingBox);
 796     BoundingBox.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
 797     tabBox[iShape][0] = Xmin; tabBox[iShape][1] = Xmax;
 798     tabBox[iShape][2] = Ymin; tabBox[iShape][3] = Ymax;
 799     tabBox[iShape][4] = Zmin; tabBox[iShape][5] = Zmax;
 800     iShape++;
 801   }
 802
 803   cout << endl;
 804   cout << "Ghs3d execution..." << endl;
 805
 806   // make a unique working file name
 807   // to avoid access to the same files by eg different users
 808
 809   TCollection_AsciiString aGenericName, aTmpDir = getTmpDir();
 810   aGenericName = aTmpDir + "GHS3D_";
 811 #ifdef WIN32
 812   aGenericName += GetCurrentProcessId();
 813 #else
 814   aGenericName += getpid();
 815 #endif
 816   aGenericName += "_";
 817   aGenericName += meshDS->ShapeToIndex( theShape );
 818
 819   TCollection_AsciiString aFacesFileName, aPointsFileName, aResultFileName;
 820   TCollection_AsciiString aBadResFileName, aBbResFileName, aLogFileName;
 821   aFacesFileName  = aGenericName + ".faces";  // in faces
 822   aPointsFileName = aGenericName + ".points"; // in points
 823   aResultFileName = aGenericName + ".noboite";// out points and volumes
 824   aBadResFileName = aGenericName + ".boite";  // out bad result
 825   aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
 826   aLogFileName    = aGenericName + ".log";    // log
 827
 828   // -----------------
 829   // make input files
 830   // -----------------
 831
 832   ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
 833   ofstream aPointsFile ( aPointsFileName.ToCString() , ios::out);
 834
 835   Ok =
 836 #ifdef WIN32
 837     aFacesFile->is_open() && aPointsFile->is_open();
 838 #else
 839     aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
 840 #endif
 841   if (!Ok) {
 842     INFOS( "Can't write into " << aTmpDir.ToCString());
 843     return error(SMESH_Comment("Can't write into ") << aTmpDir);
 844   }
 845   map <int,int> aSmdsToGhs3dIdMap;
 846   map <int,const SMDS_MeshNode*> aGhs3dIdToNodeMap;
 847
 848   Ok = writePoints( aPointsFile, meshDS, aSmdsToGhs3dIdMap, aGhs3dIdToNodeMap ) &&
 849        writeFaces ( aFacesFile,  meshDS, aSmdsToGhs3dIdMap );
 850
 851   aFacesFile.close();
 852   aPointsFile.close();
 853
 854   if ( ! Ok ) {
 855     if ( !getenv("GHS3D_KEEP_FILES") ) {
 856       OSD_File( aFacesFileName ).Remove();
 857       OSD_File( aPointsFileName ).Remove();
 858     }
 859     return error(COMPERR_BAD_INPUT_MESH);
 860   }
 861
 862   // -----------------
 863   // run ghs3d mesher              WIN32???
 864   // -----------------
 865
 866   // ghs3d need to know amount of memory it may use (MB).
 867   // Default memory is defined at ghs3d installation but it may be not enough,
 868   // so allow to use about all available memory
 869
 870   TCollection_AsciiString memory;
 871 #ifndef WIN32
 872   struct sysinfo si;
 873   int err = sysinfo( &si );
 874   if ( err == 0 ) {
 875     int freeMem = si.totalram * si.mem_unit / 1024 / 1024;
 876     memory = "-m ";
 877     memory += int( 0.7 * freeMem );
 878   }
 879 #endif
 880
 881   MESSAGE("GHS3DPlugin_GHS3D::Compute");
 882   TCollection_AsciiString cmd( "ghs3d " ); // command to run
 883   cmd +=
 884     memory +                     // memory
 885     " -c0 -f " + aGenericName +  // file to read
 886          " 1>" + aLogFileName;   // dump into file
 887
 888   system( cmd.ToCString() ); // run
 889
 890   cout << endl;
 891   cout << "End of Ghs3d execution !" << endl;
 892
 893   // --------------
 894   // read a result
 895   // --------------
 896
 897   // Mapping the result file
 898
 899   int fileOpen;
 900   fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
 901   if ( fileOpen < 0 ) {
 902     cout << endl;
 903     cout << "Can't open the " << aResultFileName.ToCString() << " GHS3D output file" << endl;
 904     Ok = false;
 905   }
 906   else
 907     Ok = readResultFile( fileOpen, meshDS, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap );
 908
 909   // ---------------------
 910   // remove working files
 911   // ---------------------
 912
 913   if ( Ok )
 914   {
 915     OSD_File( aLogFileName ).Remove();
 916   }
 917   else if ( OSD_File( aLogFileName ).Size() > 0 )
 918   {
 919     INFOS( "GHS3D Error, see the " << aLogFileName.ToCString() << " file" );
 920
 921     // get problem description from the log file
 922     SMESH_Comment comment;
 923     TCollection_AsciiString foundLine;
 924     if ( findLineContaing( "has expired",aLogFileName,foundLine) &&
 925          foundLine.Search("Licence") >= 0)
 926     {
 927       foundLine.LeftAdjust();
 928       comment << foundLine;
 929     }
 930     if ( findLineContaing( "%% ERROR",aLogFileName,foundLine))
 931     {
 932       foundLine.LeftAdjust();
 933       comment << foundLine;
 934     }
 935     if ( findLineContaing( "%% NO SAVING OPERATION",aLogFileName,foundLine))
 936     {
 937       comment << "Too many elements generated for a trial version.\n";
 938     }
 939     if ( comment.empty() )
 940       comment << "See " << aLogFileName << " for problem description";
 941     else
 942       comment << "See " << aLogFileName << " for more information";
 943     error(COMPERR_ALGO_FAILED, comment);
 944   }
 945   else
 946   {
 947     // the log file is empty
 948     OSD_File( aLogFileName ).Remove();
 949     INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" );
 950     error(COMPERR_ALGO_FAILED, "ghs3d: command not found" );
 951   }
 952
 953   if ( !getenv("GHS3D_KEEP_FILES") ) {
 954     OSD_File( aFacesFileName ).Remove();
 955     OSD_File( aPointsFileName ).Remove();
 956     OSD_File( aResultFileName ).Remove();
 957     OSD_File( aBadResFileName ).Remove();
 958     OSD_File( aBbResFileName ).Remove();
 959   }
 960   cout << "<" << aResultFileName.ToCString() << "> GHS3D output file ";
 961   if ( !Ok )
 962     cout << "not ";
 963   cout << "treated !" << endl;
 964   cout << endl;
 965
 966   _nbShape = 0;    // re-initializing _nbShape for the next Compute() method call
 967   delete [] tabShape;
 968   delete [] tabBox;
 969
 970   return Ok;
 971 }
 972
 973 //=============================================================================
 974 /*!
 975  *Here we are going to use the GHS3D mesher w/o geometry
 976  */
 977 //=============================================================================
 978 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
 979                                 SMESH_MesherHelper* aHelper)
 980 {
 981   MESSAGE("GHS3DPlugin_GHS3D::Compute()");
 982
 983   SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
 984   TopoDS_Shape theShape = aHelper->GetSubShape();
 985
 986   // make a unique working file name
 987   // to avoid access to the same files by eg different users
 988
 989   TCollection_AsciiString aGenericName, aTmpDir = getTmpDir();
 990   aGenericName = aTmpDir + "GHS3D_";
 991 #ifdef WIN32
 992   aGenericName += GetCurrentProcessId();
 993 #else
 994   aGenericName += getpid();
 995 #endif
 996   aGenericName += "_";
 997   aGenericName += meshDS->ShapeToIndex( theShape );
 998
 999   TCollection_AsciiString aFacesFileName, aPointsFileName, aResultFileName;
1000   TCollection_AsciiString aBadResFileName, aBbResFileName, aLogFileName;
1001   aFacesFileName  = aGenericName + ".faces";  // in faces
1002   aPointsFileName = aGenericName + ".points"; // in points
1003   aResultFileName = aGenericName + ".noboite";// out points and volumes
1004   aBadResFileName = aGenericName + ".boite";  // out bad result
1005   aBbResFileName  = aGenericName + ".bb";     // out vertex stepsize
1006   aLogFileName    = aGenericName + ".log";    // log
1007
1008   // -----------------
1009   // make input files
1010   // -----------------
1011
1012   ofstream aFacesFile  ( aFacesFileName.ToCString()  , ios::out);
1013   ofstream aPointsFile  ( aPointsFileName.ToCString()  , ios::out);
1014   bool Ok =
1015 #ifdef WIN32
1016     aFacesFile->is_open() && aPointsFile->is_open();
1017 #else
1018     aFacesFile.rdbuf()->is_open() && aPointsFile.rdbuf()->is_open();
1019 #endif
1020
1021   if (!Ok)
1022     return error( SMESH_Comment("Can't write into ") << aTmpDir.ToCString());
1023
1024   vector <const SMDS_MeshNode*> aNodeByGhs3dId;
1025
1026   Ok = (writeFaces ( aFacesFile, meshDS, aNodeByGhs3dId ) &&
1027         writePoints( aPointsFile, meshDS, aNodeByGhs3dId));
1028
1029   aFacesFile.close();
1030   aPointsFile.close();
1031
1032   if ( ! Ok ) {
1033     if ( !getenv("GHS3D_KEEP_FILES") ) {
1034       OSD_File( aFacesFileName ).Remove();
1035       OSD_File( aPointsFileName ).Remove();
1036     }
1037     return error(COMPERR_BAD_INPUT_MESH);
1038   }
1039
1040   // -----------------
1041   // run ghs3d mesher              WIN32???
1042   // -----------------
1043
1044   // ghs3d need to know amount of memory it may use (MB).
1045   // Default memory is defined at ghs3d installation but it may be not enough,
1046   // so allow to use about all available memory
1047   TCollection_AsciiString memory;
1048 #ifdef WIN32
1049   // ????
1050 #else
1051   struct sysinfo si;
1052   int err = sysinfo( &si );
1053   if ( !err ) {
1054     int freeMem = si.totalram * si.mem_unit / 1024 / 1024;
1055     memory = "-m ";
1056     memory += int( 0.7 * freeMem );
1057   }
1058 #endif
1059
1060   TCollection_AsciiString cmd( "ghs3d " ); // command to run
1061   cmd +=
1062     memory +                 // memory
1063     " -f " + aGenericName +  // file to read
1064     " 1>" + aLogFileName;    // dump into file
1065
1066
1067
1068   system( cmd.ToCString() ); // run
1069
1070   // --------------
1071   // read a result
1072   // --------------
1073   int fileOpen;
1074   fileOpen = open( aResultFileName.ToCString(), O_RDONLY);
1075   if ( fileOpen < 0 ) {
1076     cout << endl;
1077     cout << "Error when opening the " << aResultFileName.ToCString() << " file" << endl;
1078     cout << endl;
1079     Ok = false;
1080   }
1081   else {
1082     Ok = readResultFile( fileOpen, meshDS, theShape ,aNodeByGhs3dId );
1083   }
1084
1085   // ---------------------
1086   // remove working files
1087   // ---------------------
1088
1089   if ( Ok )
1090   {
1091     OSD_File( aLogFileName ).Remove();
1092   }
1093   else if ( OSD_File( aLogFileName ).Size() > 0 )
1094   {
1095     INFOS( "GHS3D Error, see the " << aLogFileName.ToCString() << " file" );
1096
1097     // get problem description from the log file
1098     SMESH_Comment comment;
1099     TCollection_AsciiString foundLine;
1100     if ( findLineContaing( "has expired",aLogFileName,foundLine) &&
1101          foundLine.Search("Licence") >= 0)
1102     {
1103       foundLine.LeftAdjust();
1104       comment << foundLine;
1105     }
1106     if ( findLineContaing( "%% ERROR",aLogFileName,foundLine))
1107     {
1108       foundLine.LeftAdjust();
1109       comment << foundLine;
1110     }
1111     if ( findLineContaing( "%% NO SAVING OPERATION",aLogFileName,foundLine))
1112     {
1113       comment << "Too many elements generated for a trial version.\n";
1114     }
1115     if ( comment.empty() )
1116       comment << "See " << aLogFileName << " for problem description";
1117     else
1118       comment << "See " << aLogFileName << " for more information";
1119     error(COMPERR_ALGO_FAILED, comment);
1120   }
1121   else {
1122     // the log file is empty
1123     OSD_File( aLogFileName ).Remove();
1124     INFOS( "GHS3D Error, command '" << cmd.ToCString() << "' failed" );
1125     error(COMPERR_ALGO_FAILED, "ghs3d: command not found" );
1126   }
1127
1128   if ( !getenv("GHS3D_KEEP_FILES") )
1129   {
1130     OSD_File( aFacesFileName ).Remove();
1131     OSD_File( aPointsFileName ).Remove();
1132     OSD_File( aResultFileName ).Remove();
1133     OSD_File( aBadResFileName ).Remove();
1134     OSD_File( aBbResFileName ).Remove();
1135   }
1136
1137   return Ok;
1138 }
1139
1140 //=============================================================================
1141 /*!
1142  *
1143  */
1144 //=============================================================================
1145
1146 ostream & GHS3DPlugin_GHS3D::SaveTo(ostream & save)
1147 {
1148   return save;
1149 }
1150
1151 //=============================================================================
1152 /*!
1153  *
1154  */
1155 //=============================================================================
1156
1157 istream & GHS3DPlugin_GHS3D::LoadFrom(istream & load)
1158 {
1159   return load;
1160 }
1161
1162 //=============================================================================
1163 /*!
1164  *
1165  */
1166 //=============================================================================
1167
1168 ostream & operator << (ostream & save, GHS3DPlugin_GHS3D & hyp)
1169 {
1170   return hyp.SaveTo( save );
1171 }
1172
1173 //=============================================================================
1174 /*!
1175  *
1176  */
1177 //=============================================================================
1178
1179 istream & operator >> (istream & load, GHS3DPlugin_GHS3D & hyp)
1180 {
1181   return hyp.LoadFrom( load );
1182 }