src/GHS3DPRLPlugin/GHS3DPRLPlugin_GHS3DPRL.cxx

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