src/GHS3DPRLPlugin/GHS3DPRLPlugin_GHS3DPRL.cxx

   1 // Copyright (C) 2007-2015  CEA/DEN, EDF 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, or (at your option) any later version.
   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 // ---
  21 // File   : GHS3DPRLPlugin_GHS3DPRL.cxx
  22 // Author : Christian VAN WAMBEKE (CEA) (from Hexotic plugin Lioka RAZAFINDRAZAKA)
  23 // ---
  24 //
  25 #include "GHS3DPRLPlugin_GHS3DPRL.hxx"
  26 #include "GHS3DPRLPlugin_Hypothesis.hxx"
  27
  28 #include <SMDS_MeshElement.hxx>
  29 #include <SMDS_MeshNode.hxx>
  30 #include <SMESH_subMesh.hxx>
  31
  32 #include <TopExp_Explorer.hxx>
  33 #include <OSD_File.hxx>
  34
  35 #include "utilities.h"
  36
  37 #ifndef WIN32
  38 #include <sys/sysinfo.h>
  39 #endif
  40
  41 #ifdef _DEBUG_
  42 #define DUMP(txt) \
  43 //  cout << txt
  44 #else
  45 #define DUMP(txt)
  46 #endif
  47
  48 #include <SMESH_Gen.hxx>
  49 #include <SMESHDS_Mesh.hxx>
  50 #include <SMESH_ControlsDef.hxx>
  51
  52 #include <list>
  53 #include <Standard_ProgramError.hxx>
  54 #include <TCollection_AsciiString.hxx>
  55 #include <TopTools_MapOfShape.hxx>
  56 #include <TopoDS.hxx>
  57 #include <BRepGProp.hxx>
  58 #include <GProp_GProps.hxx>
  59
  60 /*
  61 #include <med.h>
  62 //#include <med_config.h>
  63 #include <med_utils.h>
  64 //#include <med_misc.h>
  65 #include <stdlib.h>
  66 using namespace med_2_2;*/
  67
  68 static void removeFile( const TCollection_AsciiString& fileName )
  69 {
  70   try {
  71     OSD_File( fileName ).Remove();
  72   }
  73   catch ( Standard_ProgramError ) {
  74     MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
  75   }
  76 }
  77
  78 //=============================================================================
  79 GHS3DPRLPlugin_GHS3DPRL::GHS3DPRLPlugin_GHS3DPRL(int hypId, int studyId, SMESH_Gen* gen)
  80   : SMESH_3D_Algo(hypId, studyId, gen)
  81 {
  82   MESSAGE("GHS3DPRLPlugin_GHS3DPRL::GHS3DPRLPlugin_GHS3DPRL");
  83   _name = "MG-Tetra Parallel";
  84   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  85   _countSubMesh=0;
  86   _nodeRefNumber=0;
  87   _compatibleHypothesis.push_back(GHS3DPRLPlugin_Hypothesis::GetHypType());
  88 }
  89
  90 //=============================================================================
  91 GHS3DPRLPlugin_GHS3DPRL::~GHS3DPRLPlugin_GHS3DPRL()
  92 {
  93   MESSAGE("GHS3DPRLPlugin_GHS3DPRL::~GHS3DPRLPlugin_GHS3DPRL");
  94 }
  95
  96 //=============================================================================
  97 bool GHS3DPRLPlugin_GHS3DPRL::CheckHypothesis
  98                          (SMESH_Mesh& aMesh,
  99                           const TopoDS_Shape& aShape,
 100                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 101 {
 102   //MESSAGE("GHS3DPRLPlugin_GHS3DPRL::CheckHypothesis");
 103   _hypothesis = NULL;
 104
 105   list<const SMESHDS_Hypothesis*>::const_iterator itl;
 106   const SMESHDS_Hypothesis* theHyp;
 107
 108   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
 109   int nbHyp = hyps.size();
 110   if (!nbHyp)
 111   {
 112     aStatus = SMESH_Hypothesis::HYP_OK;
 113     return true;  // can work with no hypothesis
 114   }
 115
 116   itl = hyps.begin();
 117   theHyp = (*itl); // use only the first hypothesis
 118
 119   string hypName = theHyp->GetName();
 120   if (hypName == GHS3DPRLPlugin_Hypothesis::GetHypType())
 121   {
 122     _hypothesis = static_cast<const GHS3DPRLPlugin_Hypothesis*> (theHyp);
 123     ASSERT(_hypothesis);
 124     aStatus = SMESH_Hypothesis::HYP_OK;
 125   }
 126   else
 127     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 128
 129   return aStatus == SMESH_Hypothesis::HYP_OK;
 130 }
 131
 132 //=======================================================================
 133 static bool writeGHS3DPRLFiles (const TCollection_AsciiString &  GHS3DPRL_In,
 134                                 SMESHDS_Mesh *                   theMesh,
 135                                 map <int,int> &                  theSmdsToGHS3DPRLIdMap,
 136                                 map <int,const SMDS_MeshNode*> & theGHS3DPRLIdToNodeMap)
 137 {
 138    bool Ok;
 139    int ifam=0;
 140    TCollection_AsciiString namefile(GHS3DPRL_In);
 141    namefile+=".points";
 142    removeFile(namefile);
 143    ofstream theFile;
 144    theFile.open(namefile.ToCString(),ios::out);
 145 #ifdef WIN32
 146    Ok = theFile.is_open();
 147 #else
 148    Ok=theFile.rdbuf()->is_open();
 149 #endif
 150    if (!Ok)
 151    {
 152       INFOS("Can't write into "<<namefile.ToCString());
 153       return false;
 154    }
 155    cout<<endl<<"writeGHS3DPRLFiles version 2.1 "<<endl;
 156    cout<<endl<<"Creating GHS3DPRL processed mesh file : "<<namefile<<endl;
 157
 158    int nbVertices=theMesh->NbNodes();
 159    int nbFaces=theMesh->NbFaces();        //triangles or quadrangles
 160    const char* space="  ";
 161    const int dummyint=1;                  //nrs,nsd,refnum=1 (for wrap)
 162
 163    // Writing SMESH points into GHS3DPRL File.points
 164    theFile<<nbVertices<<endl;
 165
 166    int aSmdsNodeID=1;
 167    const SMDS_MeshNode* node_2;
 168    SMDS_NodeIteratorPtr itOnNode=theMesh->nodesIterator();
 169    //int ifam=100;//test famille
 170    theFile.precision(15); theFile.setf(ios::scientific,ios::floatfield);
 171    //cout<<"set precision 15 on float\n";
 172    while (itOnNode->more())
 173    {
 174       node_2 = itOnNode->next();
 175       theSmdsToGHS3DPRLIdMap.insert(map <int,int>::value_type(node_2->GetID(),aSmdsNodeID));
 176       theGHS3DPRLIdToNodeMap.insert(map <int,const SMDS_MeshNode*>::value_type(aSmdsNodeID,node_2));
 177       theFile<<node_2->X()<<space<<node_2->Y()<<space<<node_2->Z()<<space<<ifam<<endl;
 178       aSmdsNodeID++;
 179       //if (aSmdsNodeID==11) ifam++;
 180    }
 181    //no specified points;
 182    theFile.close();
 183
 184    namefile=GHS3DPRL_In+".faces";
 185    removeFile(namefile);
 186    theFile.open(namefile.ToCString(),ios::out);
 187 #ifdef WIN32
 188    Ok=theFile.is_open();
 189 #else
 190    Ok=theFile.rdbuf()->is_open();
 191 #endif
 192    if (!Ok)
 193    {
 194       INFOS("Can't write into "<<namefile.ToCString());
 195       return false;
 196    }
 197    cout<<endl<<"Creating GHS3DPRL processed mesh file : "<<namefile<<endl;
 198
 199    // Writing SMESH faces into GHS3DPRL File.faces
 200    theFile<<nbFaces<<" 0"<<endl;   //NB_ELEMS DUMMY_INT
 201                                    //" 0" is a reserved parameter
 202
 203    const SMDS_MeshElement* aFace;
 204    map<int,int>::const_iterator itOnSmdsNode;
 205    SMDS_ElemIteratorPtr itOnFaceNode;
 206    SMDS_FaceIteratorPtr itOnSmdsFace = theMesh->facesIterator();
 207    long nbNoTriangles=0;
 208    int ifaces=0;
 209    //ifam=300;
 210    while (itOnSmdsFace->more())
 211    {
 212       aFace=itOnSmdsFace->next();
 213       itOnFaceNode=aFace->nodesIterator();
 214       const int nbNodes=aFace->NbNodes();
 215       if (nbNodes!=3) nbNoTriangles++;
 216       ifaces++;
 217       theFile<<nbNodes<<space;        // NB_NODES
 218       while (itOnFaceNode->more())
 219       {
 220           aSmdsNodeID=itOnFaceNode->next()->GetID();
 221           itOnSmdsNode=theSmdsToGHS3DPRLIdMap.find(aSmdsNodeID);
 222           ASSERT(itOnSmdsNode!=theSmdsToGHS3DPRLIdMap.end());
 223           theFile<<space<<(*itOnSmdsNode).second; //NODE_1 NODE_2 ...
 224       }
 225       //(NB_NODES+1) times: DUMMY_INT
 226       //if (ifaces==11) ifam++;
 227       theFile<<space<<ifam;
 228       for ( int i=1; i<=nbNodes; i++) theFile<<space<<200+i;
 229       theFile<<endl;
 230    }
 231    theFile.close();
 232
 233    cout<<"Processed mesh files created, they contains :\n";
 234    cout<<"    "<<nbVertices<<" vertices\n";
 235    if (nbNoTriangles==0)
 236       cout<<"    "<<nbFaces<<" faces\n\n";
 237    else
 238       cout<<"    "<<nbFaces<<" faces with "<<nbNoTriangles<<"faces no triangles\n\n";
 239    return true;
 240 }
 241
 242 //=======================================================================
 243 static bool getInt( int & theValue, char * & theLine )
 244 {
 245   char *ptr;
 246   theValue = strtol( theLine, &ptr, 10 );
 247   if ( ptr == theLine ||
 248       // there must not be neither '.' nor ',' nor 'E' ...
 249       (*ptr != ' ' && *ptr != '\n' && *ptr != '\0'))
 250     return false;
 251
 252   DUMP( "  " << theValue );
 253   theLine = ptr;
 254   return true;
 255 }
 256
 257 //=======================================================================
 258 static bool getDouble( double & theValue, char * & theLine )
 259 {
 260   char *ptr;
 261   theValue = strtod( theLine, &ptr );
 262   if ( ptr == theLine )
 263     return false;
 264
 265   DUMP( "   " << theValue );
 266   theLine = ptr;
 267   return true;
 268 }
 269
 270 //=======================================================================
 271
 272 #define GHS3DPRLPlugin_BUFLENGTH 256
 273 #define GHS3DPRLPlugin_ReadLine(aPtr,aBuf,aFile,aLineNb) \
 274 {  aPtr = fgets( aBuf, GHS3DPRLPlugin_BUFLENGTH - 2, aFile ); aLineNb++; DUMP(endl); }
 275
 276 //=======================================================================
 277 static bool readResult(FILE *                           theFile,
 278                        SMESHDS_Mesh *                   theMesh,
 279                        const TopoDS_Shape &             theShape,
 280                        map <int,const SMDS_MeshNode*> & theGHS3DPRLIdToNodeMap,
 281                        const TCollection_AsciiString &  GHS3DPRL_Out,
 282                        int &                            nodeRefNumber)
 283 {
 284   // ---------------------------------
 285   // Read generated elements and nodes
 286   // ---------------------------------
 287
 288   cout << "Reading GHS3DPRL output file : " << GHS3DPRL_Out << endl;
 289   cout << endl;
 290
 291   char aBuffer[ GHS3DPRLPlugin_BUFLENGTH ];
 292   char * aPtr;
 293   int aLineNb = 0;
 294   int shapeID = theMesh->ShapeToIndex( theShape );
 295
 296   int line = 1, EndOfFile = 0, nbElem = 0, nField = 10, nbRef = 0, aGHS3DPRLNodeID = 0;
 297   const char * theField;
 298
 299   vector<const char*> tabField = vector<const char*>(nField);
 300   vector<int> tabRef = vector<int>(nField);
 301
 302   tabField[0] = "MeshVersionFormatted";    tabRef[0] = 0;
 303   tabField[1] = "Dimension";               tabRef[1] = 0;
 304   tabField[2] = "Vertices";                tabRef[2] = 3;
 305   tabField[3] = "Edges";                   tabRef[3] = 2;
 306   tabField[4] = "Triangles";               tabRef[4] = 3;
 307   tabField[5] = "Quadrilaterals";          tabRef[5] = 4;
 308   tabField[6] = "Hexahedra";               tabRef[6] = 8;
 309   tabField[7] = "Corners";                 tabRef[7] = 1;
 310   tabField[8] = "Ridges";                  tabRef[0] = 1;
 311   tabField[9] = "End";                     tabRef[0] = 0;
 312
 313   nodeRefNumber += theMesh->NbNodes();
 314
 315   SMDS_NodeIteratorPtr itOnGHS3DPRLInputNode = theMesh->nodesIterator();
 316   while ( itOnGHS3DPRLInputNode->more() )
 317       theMesh->RemoveNode( itOnGHS3DPRLInputNode->next() );
 318
 319   while ( EndOfFile == 0  ) {
 320       GHS3DPRLPlugin_ReadLine( aPtr, aBuffer, theFile, aLineNb );
 321       for ( int iField = 0; iField < nField; iField++ ) {
 322           stringstream theMessage;
 323           theField = tabField[iField];
 324           if ( strncmp(aPtr, theField, strlen(theField)) == 0 ) {
 325               if ( strcmp(theField, "End") == 0 ) {
 326                   EndOfFile = 1;
 327                   theMessage << "End of GHS3DPRL output file has been reached";
 328               }
 329               else {
 330                   GHS3DPRLPlugin_ReadLine( aPtr, aBuffer, theFile, aLineNb );
 331                   line++;
 332                   getInt( nbElem, aPtr );
 333
 334                   if ( strcmp(theField, "MeshVersionFormatted") == 0 )
 335                       theMessage << "GHS3DPRL mesh descriptor : " << theField << " " << nbElem;
 336                   else if ( strcmp(theField, "Dimension") == 0 )
 337                       theMessage << "GHS3DPRL mesh of " << nbElem << "D dimension";
 338                   else if ( strcmp(theField, "Vertices")       == 0 ||
 339                             strcmp(theField, "Edges")          == 0 ||
 340                             strcmp(theField, "Quadrilaterals") == 0 ||
 341                             strcmp(theField, "Hexahedra")      == 0 ) {
 342                       nbRef = tabRef[iField];
 343                       GHS3DPRLPlugin_ReadLine( aPtr, aBuffer, theFile, aLineNb );           // read blank line
 344
 345                       if ( strcmp(theField, "Vertices") == 0 ) {
 346                           int aGHS3DPRLID;
 347                           vector<double> coord = vector<double>(nbRef);
 348                           SMDS_MeshNode * aGHS3DPRLNode;
 349
 350                           for ( int iElem = 0; iElem < nbElem; iElem++ ) {
 351                               aGHS3DPRLID = iElem + 1 + nodeRefNumber;
 352                               GHS3DPRLPlugin_ReadLine( aPtr, aBuffer, theFile, aLineNb );   // read file lines
 353                               for ( int iCoord = 0; iCoord < 3; iCoord++ )
 354                                   getDouble ( coord[ iCoord ], aPtr );
 355                               aGHS3DPRLNode = theMesh->AddNode(coord[0], coord[1], coord[2]);
 356                               theMesh->SetNodeInVolume( aGHS3DPRLNode, shapeID );
 357                               theGHS3DPRLIdToNodeMap[ aGHS3DPRLID ] = aGHS3DPRLNode;
 358                           }
 359                       }
 360                       else {
 361                           vector<const SMDS_MeshNode*> node = vector<const SMDS_MeshNode*>(nbRef);
 362                           SMDS_MeshElement* aGHS3DPRLElement;
 363                           map <int,const SMDS_MeshNode*>::iterator itOnGHS3DPRLNode;
 364
 365                           for ( int iElem = 0; iElem < nbElem; iElem++ ) {
 366                               GHS3DPRLPlugin_ReadLine( aPtr, aBuffer, theFile, aLineNb );   // read file lines
 367                               for ( int iRef = 0; iRef < nbRef; iRef++ ) {
 368                                   getInt ( aGHS3DPRLNodeID, aPtr );                         // read nbRef aGHS3DPRLNodeID
 369                                   aGHS3DPRLNodeID += nodeRefNumber;
 370                                   itOnGHS3DPRLNode = theGHS3DPRLIdToNodeMap.find( aGHS3DPRLNodeID );
 371                                   node[ iRef ] = itOnGHS3DPRLNode->second;
 372                               }
 373
 374                               if ( strcmp(theField, "Edges") == 0 )                        // create an element
 375                                   aGHS3DPRLElement = theMesh->AddEdge( node[0], node[1] );
 376                               else if ( strcmp(theField, "Quadrilaterals") == 0 )
 377                                   aGHS3DPRLElement = theMesh->AddFace( node[0], node[1], node[2], node[3] );
 378                               else if ( strcmp(theField, "Hexahedra") == 0 )
 379                                   aGHS3DPRLElement = theMesh->AddVolume( node[0], node[1], node[2], node[3], node[4], node[5], node[6], node[7] );
 380
 381                               theMesh->SetMeshElementOnShape( aGHS3DPRLElement, shapeID );
 382                           }
 383                       }
 384                       theMessage << nbElem << " " << theField << " created";
 385                   }
 386               }
 387               if ( theMessage.str().size() != 0 ) {
 388                   cout << theMessage.str() << endl;
 389                   break;
 390               }
 391           }
 392       }
 393   }
 394   cout << endl;
 395   return true;
 396 }
 397
 398 //=============================================================================
 399 // Pass parameters to GHS3DPRL
 400 void GHS3DPRLPlugin_GHS3DPRL::SetParameters(const GHS3DPRLPlugin_Hypothesis* hyp)
 401 {
 402   if (hyp) {
 403     MESSAGE("GHS3DPRLPlugin_GHS3DPRL::SetParameters");
 404     _MEDName = hyp->GetMEDName();  //"DOMAIN\0"
 405     _NbPart = hyp->GetNbPart();
 406     _KeepFiles = hyp->GetKeepFiles();
 407     _Background = hyp->GetBackground();
 408     _ToMergeSubdomains = hyp->GetToMergeSubdomains();
 409     _ToTagSubdomains = hyp->GetToTagSubdomains();
 410     _ToOutputInterfaces = hyp->GetToOutputInterfaces();
 411     _ToDiscardSubdomains = hyp->GetToDiscardSubdomains();
 412   }
 413 }
 414
 415 //=======================================================================
 416 //before launching salome
 417 //SALOME_TMP_DIR (for keep tepal intermediates files) could be set in user's directories
 418 static TCollection_AsciiString getTmpDir()
 419 {
 420   TCollection_AsciiString aTmpDir;
 421   char *Tmp_dir = getenv("SALOME_TMP_DIR");
 422   if (Tmp_dir == NULL) Tmp_dir = getenv("TMP");
 423   if(Tmp_dir != NULL)
 424   {
 425     aTmpDir = Tmp_dir;
 426 #ifdef WIN32
 427       if(aTmpDir.Value(aTmpDir.Length()) != '\\') aTmpDir+='\\';
 428 #else
 429       if(aTmpDir.Value(aTmpDir.Length()) != '/') aTmpDir+='/';
 430 #endif
 431   }
 432   else
 433   {
 434 #ifdef WIN32
 435       aTmpDir = TCollection_AsciiString("C:\\");
 436 #else
 437       aTmpDir = TCollection_AsciiString("/tmp/");
 438 #endif
 439   }
 440   return aTmpDir;
 441 }
 442
 443 //=============================================================================
 444 // Here we are going to use the GHS3DPRL mesher for tetra-hpc (formerly tepal in v3 (2014))
 445 bool GHS3DPRLPlugin_GHS3DPRL::Compute(SMESH_Mesh& theMesh,
 446                                      const TopoDS_Shape& theShape)
 447 {
 448    bool Ok=false;
 449    TCollection_AsciiString pluginerror("ghs3dprl: ");
 450    SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
 451    //cout<<"GetMeshDS done\n";
 452    if (_countSubMesh==0){
 453       MESSAGE("GHS3DPRLPlugin_GHS3DPRL::Compute for tetra-hpc");
 454       _countTotal=0;
 455       TopExp_Explorer expf(meshDS->ShapeToMesh(), TopAbs_SOLID);
 456       for ( ; expf.More(); expf.Next() ) _countTotal++;
 457    }
 458    _countSubMesh++;
 459    //cout<<"Compute _countSubMesh "<<_countSubMesh<<endl;
 460    //no stuff if multiples submesh, multiple call compute
 461    //mesh all in one pass tepal (the last)
 462    if (_countSubMesh != _countTotal ) return true;
 463
 464    //stuff on last call
 465    if (_hypothesis==NULL){
 466       pluginerror += "No existing parameters/hypothesis for GHS3DPRL";
 467       cout <<"\n"<<pluginerror<<"\n\n";
 468       error(COMPERR_ALGO_FAILED, pluginerror.ToCString() );
 469       return false;
 470    }
 471    SetParameters(_hypothesis);
 472    cout << "\n" << _name << " parameters :\n" << endl;
 473    cout << "     generic path/name of MED Files = " << _MEDName << endl;
 474    cout << "     number of partitions = " << _NbPart << endl;
 475    cout << "     merge subdomains = " << _ToMergeSubdomains << endl;
 476    cout << "     tag subdomains = " << _ToTagSubdomains << endl;
 477    cout << "     output interfaces = " << _ToOutputInterfaces << endl;
 478    cout << "     discard dubdomains = " << _ToDiscardSubdomains << endl;
 479    cout << "     keep intermediates files (from tetra-hpc) = " << _KeepFiles << endl;
 480    cout << "     background (from tetra-hpc) = " << _Background << "\n\n";
 481
 482       //string tmpDir=getTmpDir_new();
 483       TCollection_AsciiString
 484          tmpDir=getTmpDir(),
 485          GHS3DPRL_In,GHS3DPRL_Out,GHS3DPRL_Outxml,
 486          run_GHS3DPRL("tetrahpc2med "),rm("rm -f "),run_nokeep_files,
 487          NbPart,fileskinmed(""),fileskinmesh(""),path,casenamemed;  //_MEDName.c_str());
 488
 489       casenamemed += (char *)_MEDName.c_str();
 490       int n=casenamemed.SearchFromEnd('/');
 491       if (n>0) {
 492          path=casenamemed.SubString(1,n);
 493          casenamemed=casenamemed.SubString(n+1,casenamemed.Length());
 494       }
 495       else
 496          path=tmpDir;
 497
 498       if (casenamemed.Length()>20){
 499          casenamemed=casenamemed.SubString(1,20);
 500          cerr<<"MEDName truncated (no more 20 characters) = "<<casenamemed<<endl;
 501       }
 502       cout<<"path="<<path<<endl;
 503       cout<<"casenamemed="<<casenamemed<<endl;
 504
 505       map <int,int> aSmdsToGHS3DPRLIdMap;
 506       map <int,const SMDS_MeshNode*> aGHS3DPRLIdToNodeMap;
 507       GHS3DPRL_In = path + "GHS3DPRL";
 508       GHS3DPRL_Out = path + casenamemed;
 509       GHS3DPRL_Outxml = path + casenamemed + ".xml"; //master file
 510       NbPart=_NbPart;
 511       //tepal2med --casename=/home/whoami/tmp/GHS3DPRL --number=5 --medname=DOMAIN
 512       //          --limitswap=1000 --verbose=0 --test=no --menu=no --launchtepal=yes;
 513
 514       run_GHS3DPRL = run_GHS3DPRL +
 515                      " --casename=" + GHS3DPRL_In +
 516                      " --number=" + NbPart +
 517                      " --medname=" + GHS3DPRL_Out +
 518                      " --launchtetra=yes";
 519       //no more meshhole option
 520       //if (_ToMeshHoles) run_GHS3DPRL += " --meshholes=yes"; else run_GHS3DPRL += " --meshholes=no";
 521       if (_ToMergeSubdomains) run_GHS3DPRL += " --merge_subdomains=yes"; else run_GHS3DPRL += " --merge_subdomains=no";
 522       if (_ToTagSubdomains) run_GHS3DPRL += " --tag_subdomains=yes"; else run_GHS3DPRL += " --tag_subdomains=no";
 523       if (_ToOutputInterfaces) run_GHS3DPRL += " --output_interfaces=yes"; else run_GHS3DPRL += " --output_interfaces=no";
 524       if (_ToDiscardSubdomains) run_GHS3DPRL += " --discard_subdomains=yes"; else run_GHS3DPRL += " --discard_subdomains=no";
 525       if (_Background) run_GHS3DPRL += " --background=yes"; else run_GHS3DPRL += " --background=no";
 526       run_nokeep_files = rm +GHS3DPRL_In + "* " + path + "tetrahpc.log";
 527       system( run_nokeep_files.ToCString() ); //clean files
 528       run_nokeep_files = rm + GHS3DPRL_In + "* ";
 529
 530       cout<<"GHS3DPRL command : "<<run_GHS3DPRL.ToCString()<<endl;
 531       fileskinmesh=path + "GHS3DPRL.mesh";
 532       cout<<"Write input file for tetra_hpc.exe "<<fileskinmesh<<"...";
 533       GHS3DPRL_Out = path + casenamemed;
 534       removeFile( GHS3DPRL_Outxml ); //only the master xml file
 535       //Ok=writeGHS3DPRLFiles(GHS3DPRL_In, meshDS, aSmdsToGHS3DPRLIdMap, aGHS3DPRLIdToNodeMap);
 536       bool toCreateGroups = false;
 537       theMesh.ExportGMF(fileskinmesh.ToCString(), meshDS, toCreateGroups );
 538       cout<<" ...done\n";
 539       //else {
 540       //   cout<<" ...NOT done\n";
 541       //   pluginerror = pluginerror + "problem writing input tepal files " + GHS3DPRL_In + ".mesh";
 542
 543       //Ecriture dans un fichier MED ?v2.? meme si not Ok
 544       //create empty file -> avoid warning message
 545       //med_idt fid=MEDouvrir((const char *)fileskinmed.ToCString(),MED_CREATION);
 546       //med_err ret=MEDfermer(fid);
 547       //fileskinmed=fileskinmed + "cp /home/wambeke/empty.med "+ path + "GHS3DPRL_skin.med";
 548       //system( fileskinmed.ToCString() );
 549       fileskinmed=path + "GHS3DPRL_skin.med";
 550       cout<<"Write file "<<fileskinmed<<"...";
 551       theMesh.ExportMED(fileskinmed.ToCString(),"SKIN_INITIAL",true,1);
 552       cout<<" ...done\n";
 553
 554       /*
 555       if (!Ok) {
 556          error(COMPERR_ALGO_FAILED, pluginerror.ToCString());
 557          return false; //pb abandonne
 558       }*/
 559
 560       //sometimes it is better to wait flushing files on slow filesystem...
 561       system( "sleep 3" );
 562       system( run_GHS3DPRL.ToCString() );
 563       system( "sleep 3" );
 564
 565       //cout<<"!!!reprise plugin!!!"<<endl;
 566
 567       if (_Background) {
 568          pluginerror = pluginerror + "backgrounding... plugin is not waiting for output files "+ casenamemed + "_*.med";
 569          cout<<pluginerror<<endl;
 570          error(COMPERR_ALGO_FAILED, pluginerror.ToCString());
 571          return false; //but it is not a problem but if true my message is overwritten
 572          //return true; //but it is not a problem,
 573       }
 574
 575       // read a result, GHS3DPRL_Out is the name of master file (previous xml format)
 576       FILE * aResultFile = fopen( GHS3DPRL_Outxml.ToCString(), "r" );
 577       if (aResultFile){
 578           //Ok = readResult( aResultFile, meshDS, theShape, aGHS3DPRLIdToNodeMap, GHS3DPRL_Out, _nodeRefNumber );
 579           Ok = true;
 580           Ok = false; //but it is not a problem but if true my message is overwritten
 581           fclose(aResultFile);
 582           cout<<"GHS3DPRL OK output master file "<<casenamemed<<".xml exist !\n\n";
 583           pluginerror = pluginerror + "new tetraedra not in memory, but stored in files "+ casenamemed + "_*.med";
 584           cout<<pluginerror<<endl;
 585           error(COMPERR_ALGO_FAILED, pluginerror.ToCString());
 586           if (!_KeepFiles) system( run_nokeep_files.ToCString() );
 587       }
 588       else{
 589           Ok = false; //it is a problem AND my message is NOT overwritten
 590           pluginerror = pluginerror + "output master file " + casenamemed + ".xml do not exist";
 591           cout<<pluginerror<<endl;
 592           error(COMPERR_ALGO_FAILED, pluginerror.ToCString() );
 593           cout<<"GHS3DPRL KO output files "<<GHS3DPRL_Out<<" do not exist ! see intermediates files keeped:\n";
 594           TCollection_AsciiString run_list_files("ls -alt ");
 595           run_list_files +=  GHS3DPRL_Out + "* " + GHS3DPRL_In + "* " + path + "tetrahpc.log";
 596           system( run_list_files.ToCString() );
 597           cout<<endl;
 598       }
 599       _countSubMesh=0;
 600
 601
 602     return Ok;
 603 }
 604
 605 //=============================================================================
 606 // Here we are going to use the GHS3DPRL mesher (old obsolescent version tepal in v1 & v2 (before 2014))
 607 bool GHS3DPRLPlugin_GHS3DPRL::ComputeForTepal(SMESH_Mesh& theMesh,
 608                                      const TopoDS_Shape& theShape)
 609 {
 610    bool Ok;
 611    TCollection_AsciiString pluginerror("ghs3dprl: ");
 612    SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
 613    //cout<<"GetMeshDS done\n";
 614    if (_countSubMesh==0){
 615       MESSAGE("GHS3DPRLPlugin_GHS3DPRL::Compute");
 616       _countTotal=0;
 617       TopExp_Explorer expf(meshDS->ShapeToMesh(), TopAbs_SOLID);
 618       for ( ; expf.More(); expf.Next() ) _countTotal++;
 619    }
 620    _countSubMesh++;
 621    //cout<<"Compute _countSubMesh "<<_countSubMesh<<endl;
 622    //no stuff if multiples submesh, multiple call compute
 623    //mesh all in one pass tepal (the last)
 624    if (_countSubMesh != _countTotal ) return true;
 625
 626    //stuff on last call
 627    if (_hypothesis==NULL){
 628       pluginerror += "No existing parameters/hypothesis for GHS3DPRL";
 629       cout <<"\n"<<pluginerror<<"\n\n";
 630       error(COMPERR_ALGO_FAILED, pluginerror.ToCString() );
 631       return false;
 632    }
 633    SetParameters(_hypothesis);
 634    cout << "\n" << _name << " parameters :\n" << endl;
 635    cout << "     generic path/name of MED Files = " << _MEDName << endl;
 636    cout << "     number of partitions = " << _NbPart << endl;
 637    //cout << "     mesh holes = " << _ToMeshHoles << endl;
 638    cout << "     mesh holes = " << 1 << endl;
 639    cout << "     keep intermediates files (from tepal) = " << _KeepFiles << endl;
 640    cout << "     background (from tepal) = " << _Background << "\n\n";
 641
 642       //string tmpDir=getTmpDir_new();
 643       TCollection_AsciiString
 644          tmpDir=getTmpDir(),
 645          GHS3DPRL_In,GHS3DPRL_Out,GHS3DPRL_Outxml,
 646          run_GHS3DPRL("tepal2med "),rm("rm "),run_nokeep_files,
 647          NbPart,fileskinmed(""),path,casenamemed;  //_MEDName.c_str());
 648
 649       casenamemed += (char *)_MEDName.c_str();
 650       int n=casenamemed.SearchFromEnd('/');
 651       if (n>0) {
 652          path=casenamemed.SubString(1,n);
 653          casenamemed=casenamemed.SubString(n+1,casenamemed.Length());
 654       }
 655       else
 656          path=tmpDir;
 657
 658       if (casenamemed.Length()>20){
 659          casenamemed=casenamemed.SubString(1,20);
 660          cerr<<"MEDName truncated (no more 20 characters) = "<<casenamemed<<endl;
 661       }
 662       cout<<"path="<<path<<endl;
 663       cout<<"casenamemed="<<casenamemed<<endl;
 664
 665       map <int,int> aSmdsToGHS3DPRLIdMap;
 666       map <int,const SMDS_MeshNode*> aGHS3DPRLIdToNodeMap;
 667       GHS3DPRL_In = path + "GHS3DPRL";
 668       GHS3DPRL_Out = path + casenamemed;
 669       GHS3DPRL_Outxml = path + casenamemed + ".xml"; //master file
 670       NbPart=_NbPart;
 671       //tepal2med --casename=/home/whoami/tmp/GHS3DPRL --number=5 --medname=DOMAIN
 672       //          --limitswap=1000 --verbose=0 --test=no --menu=no --launchtepal=yes;
 673
 674       run_GHS3DPRL = run_GHS3DPRL +
 675                      " --casename=" + GHS3DPRL_In +
 676                      " --number=" + NbPart +
 677                      " --medname=" + GHS3DPRL_Out +
 678                      " --launchtepal=yes";
 679       //if (_ToMeshHoles) run_GHS3DPRL += " --meshholes=yes"; else run_GHS3DPRL += " --meshholes=no";
 680       if (1) run_GHS3DPRL += " --meshholes=yes"; else run_GHS3DPRL += " --meshholes=no";
 681       if (_Background) run_GHS3DPRL += " --background=yes"; else run_GHS3DPRL += " --background=no";
 682       run_nokeep_files = rm +GHS3DPRL_In + "* " + path + "tepal.log";
 683       system( run_nokeep_files.ToCString() ); //clean files
 684       run_nokeep_files = rm + GHS3DPRL_In + "* ";
 685
 686       cout<<"GHS3DPRL command : "<<run_GHS3DPRL.ToCString()<<endl;
 687       cout<<"Write files .faces .point ...";
 688       GHS3DPRL_Out = path + casenamemed;
 689       removeFile( GHS3DPRL_Outxml ); //only the master xml file
 690       Ok=writeGHS3DPRLFiles(GHS3DPRL_In, meshDS, aSmdsToGHS3DPRLIdMap, aGHS3DPRLIdToNodeMap);
 691       if (Ok) {cout<<" ...done\n";}
 692       else {
 693          cout<<" ...NOT done\n";
 694          pluginerror = pluginerror + "problem writing input tepal files " + GHS3DPRL_In + "[.faces|.points]";
 695       }
 696
 697       //Ecriture dans un fichier MED ?v2.? meme si not Ok
 698       //create empty file -> avoid warning message
 699       //med_idt fid=MEDouvrir((const char *)fileskinmed.ToCString(),MED_CREATION);
 700       //med_err ret=MEDfermer(fid);
 701       //fileskinmed=fileskinmed + "cp /home/wambeke/empty.med "+ path + "GHS3DPRL_skin.med";
 702       //system( fileskinmed.ToCString() );
 703       fileskinmed=path + "GHS3DPRL_skin.med";
 704       cout<<"Write file "<<fileskinmed<<"...";
 705       theMesh.ExportMED(fileskinmed.ToCString(),"SKIN_INITIAL",true,1);
 706       cout<<" ...done\n";
 707
 708       if (!Ok) {
 709          error(COMPERR_ALGO_FAILED, pluginerror.ToCString());
 710          return false; //pb abandonne
 711       }
 712
 713       //sometimes it is better to wait flushing files on slow filesystem...
 714       system( "sleep 3" );
 715       system( run_GHS3DPRL.ToCString() );
 716       system( "sleep 3" );
 717
 718       //cout<<"!!!reprise plugin!!!"<<endl;
 719
 720       if (_Background) {
 721          pluginerror = pluginerror + "backgrounding... plugin is not waiting for output files "+ casenamemed + "_*.med";
 722          cout<<pluginerror<<endl;
 723          error(COMPERR_ALGO_FAILED, pluginerror.ToCString());
 724          return false; //but it is not a problem but if true my message is overwritten
 725          //return true; //but it is not a problem,
 726       }
 727
 728       // read a result, GHS3DPRL_Out is the name of master file (previous xml format)
 729       FILE * aResultFile = fopen( GHS3DPRL_Outxml.ToCString(), "r" );
 730       if (aResultFile){
 731           //Ok = readResult( aResultFile, meshDS, theShape, aGHS3DPRLIdToNodeMap, GHS3DPRL_Out, _nodeRefNumber );
 732           Ok = true;
 733           Ok = false; //but it is not a problem but if true my message is overwritten
 734           fclose(aResultFile);
 735           cout<<"GHS3DPRL OK output master file "<<casenamemed<<".xml exist !\n\n";
 736           pluginerror = pluginerror + "new tetraedra not in memory, but stored in files "+ casenamemed + "_*.med";
 737           cout<<pluginerror<<endl;
 738           error(COMPERR_ALGO_FAILED, pluginerror.ToCString());
 739           if (!_KeepFiles) system( run_nokeep_files.ToCString() );
 740       }
 741       else{
 742           Ok = false; //it is a problem AND my message is NOT overwritten
 743           pluginerror = pluginerror + "output master file " + casenamemed + ".xml do not exist";
 744           cout<<pluginerror<<endl;
 745           error(COMPERR_ALGO_FAILED, pluginerror.ToCString() );
 746           cout<<"GHS3DPRL KO output files "<<GHS3DPRL_Out<<" do not exist ! see intermediates files keeped:\n";
 747           TCollection_AsciiString run_list_files("ls -alt ");
 748           run_list_files +=  GHS3DPRL_Out + "* " + GHS3DPRL_In + "* " + path + "tepal.log";
 749           system( run_list_files.ToCString() );
 750           cout<<endl;
 751       }
 752       _countSubMesh=0;
 753
 754
 755     return Ok;
 756    /*pid_t pid = fork();
 757    if (pid > 0) {
 758       //Processus pere
 759       cout<<"le pere est la\n";
 760       system("echo le pere > lepere.tmp");
 761       system("sleep 10");
 762       system("ps -edf |grep wambeke > pslepere.tmp");
 763       cout<<"le pere return 0\n";
 764       return 0; //ok
 765       //exit(0);
 766    } else if (pid == 0) {
 767       //Processus fils
 768       cout<<"le fils est la\n";
 769       //On rend le fils independant de tout terminal
 770       setsid();
 771       system("sleep 20");
 772       system("echo le fils > lefils.tmp");
 773       system("sleep 20");
 774       system("ps -edf |grep wambeke > pslefils.tmp");
 775       cout<<"le fils return 0\n";
 776       return 0; //ok
 777    } else {
 778       //Traitement d'erreur
 779       cout<<"ya probleme sur fork()\n";
 780       return 1; //ko
 781    }*/
 782 }
 783
 784 //=============================================================================
 785 ostream & GHS3DPRLPlugin_GHS3DPRL::SaveTo(ostream & save)
 786 {
 787   return save;
 788 }
 789
 790 //=============================================================================
 791 istream & GHS3DPRLPlugin_GHS3DPRL::LoadFrom(istream & load)
 792 {
 793   return load;
 794 }
 795
 796 //=============================================================================
 797 ostream & operator << (ostream & save, GHS3DPRLPlugin_GHS3DPRL & hyp)
 798 {
 799   return hyp.SaveTo( save );
 800 }
 801
 802 //=============================================================================
 803 istream & operator >> (istream & load, GHS3DPRLPlugin_GHS3DPRL & hyp)
 804 {
 805   return hyp.LoadFrom( load );
 806 }
 807
 808
 809 //=============================================================================
 810 /*!
 811  *
 812  */
 813 //=============================================================================
 814 bool GHS3DPRLPlugin_GHS3DPRL::Evaluate(SMESH_Mesh& aMesh,
 815                                        const TopoDS_Shape& aShape,
 816                                        MapShapeNbElems& aResMap)
 817 {
 818   int nbtri = 0, nbqua = 0;
 819   double fullArea = 0.0;
 820   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
 821     TopoDS_Face F = TopoDS::Face( exp.Current() );
 822     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
 823     MapShapeNbElemsItr anIt = aResMap.find(sm);
 824     if( anIt==aResMap.end() ) {
 825       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 826       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
 827                                             "Submesh can not be evaluated",this));
 828       return false;
 829     }
 830     std::vector<int> aVec = (*anIt).second;
 831     nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
 832     nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 833     GProp_GProps G;
 834     BRepGProp::SurfaceProperties(F,G);
 835     double anArea = G.Mass();
 836     fullArea += anArea;
 837   }
 838
 839   // collect info from edges
 840   int nb0d_e = 0, nb1d_e = 0;
 841   bool IsQuadratic = false;
 842   bool IsFirst = true;
 843   TopTools_MapOfShape tmpMap;
 844   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
 845     TopoDS_Edge E = TopoDS::Edge(exp.Current());
 846     if( tmpMap.Contains(E) )
 847       continue;
 848     tmpMap.Add(E);
 849     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
 850     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
 851     std::vector<int> aVec = (*anIt).second;
 852     nb0d_e += aVec[SMDSEntity_Node];
 853     nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
 854     if(IsFirst) {
 855       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 856       IsFirst = false;
 857     }
 858   }
 859   tmpMap.Clear();
 860
 861   double ELen = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
 862
 863   GProp_GProps G;
 864   BRepGProp::VolumeProperties(aShape,G);
 865   double aVolume = G.Mass();
 866   double tetrVol = 0.1179*ELen*ELen*ELen;
 867   double CoeffQuality = 0.9;
 868   int nbVols = (int)aVolume/tetrVol/CoeffQuality;
 869   int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
 870   int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5;
 871   std::vector<int> aVec(SMDSEntity_Last);
 872   for(int i=0; i<SMDSEntity_Last; i++) aVec[i]=0;
 873   if( IsQuadratic ) {
 874     aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
 875     aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
 876     aVec[SMDSEntity_Quad_Pyramid] = nbqua;
 877   }
 878   else {
 879     aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
 880     aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
 881     aVec[SMDSEntity_Pyramid] = nbqua;
 882   }
 883   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
 884   aResMap.insert(std::make_pair(sm,aVec));
 885
 886   return true;
 887 }