src/GHS3DPRLPlugin/GHS3DPRLPlugin_GHS3DPRL.cxx

   1 // Copyright (C) 2007-2016  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 #include "MG_TetraHPC_API.hxx"
  28
  29 #include <SMESHDS_Mesh.hxx>
  30 #include <SMESH_Gen.hxx>
  31 #include <SMESH_TypeDefs.hxx>
  32 #include <SMESH_subMesh.hxx>
  33
  34 #include "utilities.h"
  35
  36 #include <list>
  37
  38 #include <BRepGProp.hxx>
  39 #include <GProp_GProps.hxx>
  40 #include <OSD_File.hxx>
  41 #include <Standard_ProgramError.hxx>
  42 #include <TCollection_AsciiString.hxx>
  43 #include <TopExp_Explorer.hxx>
  44 #include <TopTools_MapOfShape.hxx>
  45 #include <TopoDS.hxx>
  46 #include <TopoDS_Edge.hxx>
  47 #include <TopoDS_Face.hxx>
  48
  49 #define GMFVERSION GmfDouble
  50 #define GMFDIMENSION 3
  51
  52 using namespace std;
  53
  54 static void removeFile( const TCollection_AsciiString& fileName )
  55 {
  56   try {
  57     OSD_File( fileName ).Remove();
  58   }
  59   catch ( Standard_ProgramError ) {
  60     MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
  61   }
  62 }
  63
  64 //=============================================================================
  65 GHS3DPRLPlugin_GHS3DPRL::GHS3DPRLPlugin_GHS3DPRL(int hypId, SMESH_Gen* gen)
  66   : SMESH_3D_Algo(hypId, gen)
  67 {
  68   MESSAGE("GHS3DPRLPlugin_GHS3DPRL::GHS3DPRLPlugin_GHS3DPRL");
  69   _name = "MG-Tetra Parallel";
  70   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  71   _onlyUnaryInput = false; // Compute() will be called on a compound of solids
  72   _countSubMesh=0;
  73   _nodeRefNumber=0;
  74   _compatibleHypothesis.push_back(GHS3DPRLPlugin_Hypothesis::GetHypType());
  75 }
  76
  77 //=============================================================================
  78 GHS3DPRLPlugin_GHS3DPRL::~GHS3DPRLPlugin_GHS3DPRL()
  79 {
  80   MESSAGE("GHS3DPRLPlugin_GHS3DPRL::~GHS3DPRLPlugin_GHS3DPRL");
  81 }
  82
  83 //=============================================================================
  84 bool GHS3DPRLPlugin_GHS3DPRL::CheckHypothesis
  85                          (SMESH_Mesh& aMesh,
  86                           const TopoDS_Shape& aShape,
  87                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
  88 {
  89   //MESSAGE("GHS3DPRLPlugin_GHS3DPRL::CheckHypothesis");
  90   _hypothesis = NULL;
  91
  92   list<const SMESHDS_Hypothesis*>::const_iterator itl;
  93   const SMESHDS_Hypothesis* theHyp;
  94
  95   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
  96   int nbHyp = hyps.size();
  97   if (!nbHyp)
  98   {
  99     aStatus = SMESH_Hypothesis::HYP_MISSING;
 100     return false;  // can't work with no hypothesis
 101   }
 102
 103   itl = hyps.begin();
 104   theHyp = (*itl); // use only the first hypothesis
 105
 106   string hypName = theHyp->GetName();
 107   if (hypName == GHS3DPRLPlugin_Hypothesis::GetHypType())
 108   {
 109     _hypothesis = static_cast<const GHS3DPRLPlugin_Hypothesis*> (theHyp);
 110     ASSERT(_hypothesis);
 111     aStatus = SMESH_Hypothesis::HYP_OK;
 112   }
 113   else
 114     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 115
 116   return aStatus == SMESH_Hypothesis::HYP_OK;
 117 }
 118
 119 //=============================================================================
 120 // Pass parameters to GHS3DPRL
 121 void GHS3DPRLPlugin_GHS3DPRL::SetParameters(const GHS3DPRLPlugin_Hypothesis* hyp)
 122 {
 123   if (hyp) {
 124     MESSAGE("GHS3DPRLPlugin_GHS3DPRL::SetParameters");
 125     _MEDName     = hyp->GetMEDName();  //"DOMAIN\0"
 126     _NbPart      = hyp->GetNbPart();
 127     _KeepFiles   = hyp->GetKeepFiles();
 128     _Background  = hyp->GetBackground();
 129     _Multithread = hyp->GetMultithread();
 130     _Gradation   = hyp->GetGradation();
 131     _MinSize     = hyp->GetMinSize();
 132     _MaxSize     = hyp->GetMaxSize();
 133     _AdvOptions  = hyp->GetAdvancedOption();
 134   }
 135 }
 136
 137 //=======================================================================
 138 //before launching salome
 139 //SALOME_TMP_DIR (for keep tepal intermediates files) could be set in user's directories
 140 static TCollection_AsciiString getTmpDir()
 141 {
 142   TCollection_AsciiString aTmpDir;
 143   char *Tmp_dir = getenv("SALOME_TMP_DIR");
 144   if (Tmp_dir == NULL) Tmp_dir = getenv("TMP");
 145   if(Tmp_dir != NULL)
 146   {
 147     aTmpDir = Tmp_dir;
 148 #ifdef WIN32
 149       if(aTmpDir.Value(aTmpDir.Length()) != '\\') aTmpDir+='\\';
 150 #else
 151       if(aTmpDir.Value(aTmpDir.Length()) != '/') aTmpDir+='/';
 152 #endif
 153   }
 154   else
 155   {
 156 #ifdef WIN32
 157       aTmpDir = TCollection_AsciiString("C:\\");
 158 #else
 159       aTmpDir = TCollection_AsciiString("/tmp/");
 160 #endif
 161   }
 162   return aTmpDir;
 163 }
 164
 165 //=============================================================================
 166 // Write a skin mesh into a GMF file or pass it to MG-TetraHPC API
 167 static void exportGMF(MG_TetraHPC_API*    theTetraInput,
 168                       const char*         theFile,
 169                       const SMESHDS_Mesh* theMeshDS)
 170 {
 171   int meshID = theTetraInput->GmfOpenMesh( theFile, GmfWrite, GMFVERSION, GMFDIMENSION);
 172
 173   // nodes
 174   int iN = 0, nbNodes = theMeshDS->NbNodes();
 175   theTetraInput->GmfSetKwd( meshID, GmfVertices, nbNodes );
 176   std::map< const SMDS_MeshNode*, int, TIDCompare > node2IdMap;
 177   SMDS_NodeIteratorPtr nodeIt = theMeshDS->nodesIterator();
 178   SMESH_TNodeXYZ n;
 179   while ( nodeIt->more() )
 180   {
 181     n.Set( nodeIt->next() );
 182     theTetraInput->GmfSetLin( meshID, GmfVertices, n.X(), n.Y(), n.Z(), n._node->getshapeId() );
 183     node2IdMap.insert( node2IdMap.end(), std::make_pair( n._node, ++iN ));
 184   }
 185
 186   // triangles
 187   SMDS_ElemIteratorPtr elemIt = theMeshDS->elementGeomIterator( SMDSGeom_TRIANGLE );
 188   if ( elemIt->more() )
 189   {
 190     int nbTria = theMeshDS->GetMeshInfo().NbElements( SMDSGeom_TRIANGLE );
 191     theTetraInput->GmfSetKwd(meshID, GmfTriangles, nbTria );
 192     for ( int gmfID = 1; elemIt->more(); ++gmfID )
 193     {
 194       const SMDS_MeshElement* tria = elemIt->next();
 195       theTetraInput->GmfSetLin(meshID, GmfTriangles,
 196                               node2IdMap[ tria->GetNode( 0 )],
 197                               node2IdMap[ tria->GetNode( 1 )],
 198                               node2IdMap[ tria->GetNode( 2 )],
 199                               tria->getshapeId() );
 200     }
 201   }
 202   theTetraInput->GmfCloseMesh( meshID );
 203 }
 204
 205 //=============================================================================
 206 // Here we are going to use the GHS3DPRL mesher for tetra-hpc (formerly tepal in v3 (2014))
 207 bool GHS3DPRLPlugin_GHS3DPRL::Compute(SMESH_Mesh&         theMesh,
 208                                       const TopoDS_Shape& theShape)
 209 {
 210   bool Ok=false;
 211   TCollection_AsciiString pluginerror("ghs3dprl: ");
 212   SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
 213
 214   if (_hypothesis==NULL){
 215     pluginerror += "No existing parameters/hypothesis for GHS3DPRL";
 216     cout <<"\n"<<pluginerror<<"\n\n";
 217     error(COMPERR_ALGO_FAILED, pluginerror.ToCString() );
 218     return false;
 219   }
 220   SetParameters(_hypothesis);
 221   cout << "\n" << _name << " parameters:" << endl;
 222   cout << "   generic path/name of MED Files = " << _MEDName << endl;
 223   cout << "   number of partitions = " << _NbPart << endl;
 224   cout << "   gradation = " << _Gradation << endl;
 225   cout << "   min_size = " << _MinSize << endl;
 226   cout << "   max_size = " << _MaxSize << endl;
 227   cout << "   keep intermediates files (from tetra-hpc) = " << _KeepFiles << endl;
 228   cout << "   background (from tetra-hpc) = " << _Background << "\n";
 229   cout << "   multithread = " << _Multithread << "\n\n";
 230   cout << "   adv. options = '" << _AdvOptions << "'\n\n";
 231
 232   TCollection_AsciiString
 233     tmpDir=getTmpDir(),
 234     GHS3DPRL_In,
 235     GHS3DPRL_Out,
 236     GHS3DPRL_Out_Mesh,
 237     GHS3DPRL_Outxml,
 238     logFileName,
 239     run_GHS3DPRL("tetrahpc2med "),
 240     rm("rm -f "),
 241     run_nokeep_files,
 242     NbPart,
 243     Gradation,
 244     MinSize,
 245     MaxSize,
 246     fileskinmed(""),
 247     fileskinmesh(""),
 248     path,
 249     casenamemed;  //_MEDName.c_str());
 250
 251   casenamemed += (char *)_MEDName.c_str();
 252   int n=casenamemed.SearchFromEnd('/');
 253   if (n>0) {
 254     path=casenamemed.SubString(1,n);
 255     casenamemed=casenamemed.SubString(n+1,casenamemed.Length());
 256   }
 257   else
 258     path=tmpDir;
 259
 260   if (casenamemed.Length()>20){
 261     casenamemed=casenamemed.SubString(1,20);
 262     cerr<<"MEDName truncated (no more 20 characters) = "<<casenamemed<<endl;
 263   }
 264
 265   map <int,int> aSmdsToGHS3DPRLIdMap;
 266   map <int,const SMDS_MeshNode*> aGHS3DPRLIdToNodeMap;
 267   GHS3DPRL_In = path + "GHS3DPRL";
 268   GHS3DPRL_Out = path + casenamemed;
 269   GHS3DPRL_Out_Mesh = path + casenamemed + "_out.mesh";
 270   GHS3DPRL_Outxml = path + casenamemed + ".xml"; //master file
 271   logFileName = path + casenamemed + ".log"; // MG library output
 272   NbPart=_NbPart;
 273   Gradation=_Gradation;
 274   MinSize=_MinSize;
 275   MaxSize=_MaxSize;
 276
 277   //an example:
 278   //tetrahpc2med --casename=/home/whoami/tmp/GHS3DPRL --number=5 --medname=DOMAIN
 279   //             --gradation=1.05 --min_size=1e-3 --max_size=1e-2
 280   //             --verbose=0 --menu=no --launchtetra=yes;
 281
 282   run_GHS3DPRL = run_GHS3DPRL +
 283     " --casename=" + GHS3DPRL_In +
 284     " --number=" + NbPart +
 285     " --medname=" + GHS3DPRL_Out +
 286     " --launchtetra=yes" +
 287     " --gradation=" + Gradation +
 288     " --min_size=" + MinSize +
 289     " --max_size=" + MaxSize +
 290     " --verbose=3" +
 291     " " + _AdvOptions.c_str();
 292   if (_Background) run_GHS3DPRL += " --background=yes"; else run_GHS3DPRL += " --background=no";
 293   if (_Multithread) run_GHS3DPRL += " --multithread=yes"; else run_GHS3DPRL += " --multithread=no";
 294   run_nokeep_files = rm +GHS3DPRL_In + "* " + path + "tetrahpc.log";
 295   system( run_nokeep_files.ToCString() ); //clean files
 296   run_nokeep_files = rm + GHS3DPRL_In + "* ";
 297
 298   fileskinmesh=path + "GHS3DPRL.mesh";
 299   GHS3DPRL_Out = path + casenamemed;
 300   removeFile( GHS3DPRL_Outxml ); //only the master xml file
 301
 302   MG_TetraHPC_API mgTetraHPC( _computeCanceled, _progress );
 303   bool useLib = ( mgTetraHPC.IsLibrary() && !_Background && _Multithread );
 304   if ( !useLib )
 305     mgTetraHPC.SetUseExecutable();
 306
 307   exportGMF( &mgTetraHPC, fileskinmesh.ToCString(), meshDS );
 308
 309   if ( useLib )
 310   {
 311     TCollection_AsciiString cmd = TCollection_AsciiString("mg-tetra_hpc.exe") +
 312       " --number_of_subdomains=" + NbPart +
 313       " --gradation=" + Gradation +
 314       " --min_size=" + MinSize +
 315       " --max_size=" + MaxSize +
 316       " --verbose=3" +
 317       " --out=" + GHS3DPRL_Out_Mesh +
 318       " " + _AdvOptions.c_str();
 319
 320     cout << endl
 321          << "  Run mg-tetra_hpc as library. Creating a mesh file " << GHS3DPRL_Out_Mesh << endl
 322          << "  Creating a log file : " << logFileName << endl << endl;
 323     mgTetraHPC.SetLogFile( logFileName.ToCString() );
 324
 325     mgTetraHPC.Compute( cmd.ToCString() );
 326
 327     std::string log = mgTetraHPC.GetLog();
 328     if ( log.find(" Dlim "   ) != std::string::npos ||
 329          log.find(" license ") != std::string::npos )
 330       return error("License problem");
 331
 332     // set --launchtetra=no
 333     int yesPos = run_GHS3DPRL.Search("launchtetra") + sizeof("launchtetra");
 334     run_GHS3DPRL.SetValue( yesPos+0, 'n' );
 335     run_GHS3DPRL.SetValue( yesPos+1, 'o' );
 336     run_GHS3DPRL.SetValue( yesPos+2, ' ' );
 337   }
 338   else
 339   {
 340     cout<<"  Write input file for mg-tetra_hpc "<<fileskinmesh<<"...";
 341     cout<<" ...done\n";
 342   }
 343   fileskinmed=path + "GHS3DPRL_skin.med";
 344   cout<<"  Write file "<<fileskinmed<<"...";
 345   theMesh.ExportMED(fileskinmed.ToCString(),"SKIN_INITIAL",true,1);
 346   cout<<" ...done\n\n";
 347
 348   cout<<"GHS3DPRL command :\n  "<<run_GHS3DPRL.ToCString()<<endl;
 349   //sometimes it is better to wait flushing files on slow filesystem...
 350   system( "sleep 3" );
 351   //launch tetrahpc2med which launch mg-tetra_hpc.py which launch mg-tetra_hpc(_mpi?).exe
 352   system( run_GHS3DPRL.ToCString() );
 353   system( "sleep 3" );
 354
 355   if (_Background) {
 356     pluginerror = pluginerror + "backgrounding... plugin is not waiting for output files "+ casenamemed + "_*.med";
 357     cout<<pluginerror<<endl;
 358     error(COMPERR_ALGO_FAILED, pluginerror.ToCString());
 359     return false; //but it is not a problem but if true my message is overwritten
 360     //return true; //but it is not a problem,
 361   }
 362
 363   // read a result, GHS3DPRL_Out is the name of master file (previous xml format)
 364   FILE * aResultFile = fopen( GHS3DPRL_Outxml.ToCString(), "r" );
 365   if (aResultFile){
 366     Ok = false; //but it is not a problem but if true my message is overwritten
 367     fclose(aResultFile);
 368     cout<<"GHS3DPRL OK output master file "<<casenamemed<<".xml exist !\n\n";
 369     pluginerror = pluginerror + "new tetraedra not in memory, but stored in files "+ casenamemed + "_*.med";
 370     cout<<pluginerror<<endl;
 371     error(COMPERR_ALGO_FAILED, pluginerror.ToCString());
 372     if (!_KeepFiles) system( run_nokeep_files.ToCString() );
 373   }
 374   else{
 375     Ok = false; //it is a problem AND my message is NOT overwritten
 376     pluginerror = pluginerror + "output master file " + casenamemed + ".xml do not exist";
 377     cout<<pluginerror<<endl;
 378     error(COMPERR_ALGO_FAILED, pluginerror.ToCString() );
 379     cout<<"GHS3DPRL KO output files "<<GHS3DPRL_Out<<" do not exist ! see intermediates files keeped:\n";
 380     TCollection_AsciiString run_list_files("ls -alt ");
 381     run_list_files +=  GHS3DPRL_Out + "* " + GHS3DPRL_In + "* " + path + "tetrahpc.log";
 382     system( run_list_files.ToCString() );
 383     cout<<endl;
 384   }
 385
 386   return Ok;
 387 }
 388
 389
 390
 391 //=============================================================================
 392 /*!
 393  *
 394  */
 395 //=============================================================================
 396 bool GHS3DPRLPlugin_GHS3DPRL::Evaluate(SMESH_Mesh& aMesh,
 397                                        const TopoDS_Shape& aShape,
 398                                        MapShapeNbElems& aResMap)
 399 {
 400   int nbtri = 0, nbqua = 0;
 401   double fullArea = 0.0;
 402   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
 403     TopoDS_Face F = TopoDS::Face( exp.Current() );
 404     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
 405     MapShapeNbElemsItr anIt = aResMap.find(sm);
 406     if( anIt==aResMap.end() ) {
 407       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 408       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
 409                                             "Submesh can not be evaluated",this));
 410       return false;
 411     }
 412     std::vector<int> aVec = (*anIt).second;
 413     nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
 414     nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 415     GProp_GProps G;
 416     BRepGProp::SurfaceProperties(F,G);
 417     double anArea = G.Mass();
 418     fullArea += anArea;
 419   }
 420
 421   // collect info from edges
 422   int nb0d_e = 0, nb1d_e = 0;
 423   bool IsQuadratic = false;
 424   bool IsFirst = true;
 425   TopTools_MapOfShape tmpMap;
 426   for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
 427     const TopoDS_Edge& E = TopoDS::Edge(exp.Current());
 428     if( tmpMap.Contains(E) )
 429       continue;
 430     tmpMap.Add(E);
 431     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
 432     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
 433     std::vector<int> aVec = (*anIt).second;
 434     nb0d_e += aVec[SMDSEntity_Node];
 435     nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
 436     if(IsFirst) {
 437       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 438       IsFirst = false;
 439     }
 440   }
 441   tmpMap.Clear();
 442
 443   double ELen = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
 444
 445   GProp_GProps G;
 446   BRepGProp::VolumeProperties(aShape,G);
 447   double aVolume = G.Mass();
 448   double tetrVol = 0.1179*ELen*ELen*ELen;
 449   double CoeffQuality = 0.9;
 450   int nbVols = (int)aVolume/tetrVol/CoeffQuality;
 451   int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
 452   int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5;
 453   std::vector<int> aVec(SMDSEntity_Last);
 454   for(int i=0; i<SMDSEntity_Last; i++) aVec[i]=0;
 455   if( IsQuadratic ) {
 456     aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
 457     aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
 458     aVec[SMDSEntity_Quad_Pyramid] = nbqua;
 459   }
 460   else {
 461     aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
 462     aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
 463     aVec[SMDSEntity_Pyramid] = nbqua;
 464   }
 465   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
 466   aResMap.insert(std::make_pair(sm,aVec));
 467
 468   return true;
 469 }