src/GMSHPlugin/GMSHPlugin_GMSH_3D_Remote.cxx

   1 // Copyright (C) 2007-2023  CEA, EDF, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 //=============================================================================
  24 // File      : GMSHPlugin_GMSH_3D_Remote.hxx
  25 // Created   : 09 Septembre 2023
  26 // Author    : Cesar Conopoima (OCC)
  27 // Project   : SALOME
  28 //=============================================================================
  29 //
  30 //
  31
  32 #include "GMSHPlugin_GMSH_3D_Remote.hxx"
  33 #include "GMSHPlugin_Hypothesis.hxx"
  34 #include "Utils_SALOME_Exception.hxx"
  35
  36 #include <SMESH_Gen.hxx>
  37 #include <SMESH_Mesh.hxx>
  38 #include <SMESH_ParallelMesh.hxx>
  39 #include <SMESH_MesherHelper.hxx>
  40 #include <SMESH_DriverShape.hxx>
  41 #include <SMESH_DriverMesh.hxx>
  42 #include <SMESHDS_Mesh.hxx>
  43 #include <SMESH_MeshLocker.hxx>
  44 #include <SMESH_ProxyMesh.hxx>
  45
  46 #include <TopoDS.hxx>
  47 #include <TopExp.hxx>
  48 #include <TopExp_Explorer.hxx>
  49
  50 #include <QString>
  51 #include <QProcess>
  52
  53 #ifdef WIN32
  54 #include <filesystem>
  55 namespace fs = std::filesystem;
  56 #else
  57 #include <boost/filesystem.hpp>
  58 namespace fs = boost::filesystem;
  59 #endif
  60
  61 //=============================================================================
  62 /*!
  63  * Constructor
  64  */
  65 //=============================================================================
  66
  67 GMSHPlugin_GMSH_3D_Remote::GMSHPlugin_GMSH_3D_Remote(int hypId, SMESH_Gen * gen)
  68   : GMSHPlugin_GMSH_3D(hypId, gen)
  69 {
  70   _name = "GMSH_3D_Remote";
  71 }
  72
  73 //=============================================================================
  74 /*!
  75  * Destructor
  76  */
  77 //=============================================================================
  78
  79 GMSHPlugin_GMSH_3D_Remote::~GMSHPlugin_GMSH_3D_Remote()
  80 {
  81 }
  82
  83 /**
  84  * @brief Fill the structure netgen_param with the information from the hypothesis
  85  *
  86  * @param param_file name of the file to saven the gmsh parameter
  87  * @param hyp the hypothesis
  88  */
  89 void GMSHPlugin_GMSH_3D_Remote::exportGmshParams( const std::string param_file, const SMESHDS_Hypothesis* hyp )
  90 {
  91   std::ofstream myfile(param_file);
  92   if ( const GMSHPlugin_Hypothesis* gmshHypo = dynamic_cast<const GMSHPlugin_Hypothesis*>(hyp) )
  93   {
  94     myfile << 1 << std::endl; // Mark existence of correct hypothesis
  95     myfile << (int) gmshHypo->Get2DAlgo()       << std::endl;
  96     myfile << (int) gmshHypo->Get3DAlgo()       << std::endl;
  97     myfile << (int) gmshHypo->GetRecomb2DAlgo() << std::endl;
  98     myfile << (int) gmshHypo->GetRecombineAll() << std::endl;
  99     myfile << (int) gmshHypo->GetSubdivAlgo()   << std::endl;
 100     myfile << (int) gmshHypo->GetRemeshAlgo()   << std::endl;
 101     myfile << (double) gmshHypo->GetSmouthSteps()  << std::endl;
 102     myfile << (double) gmshHypo->GetSizeFactor()   << std::endl;
 103 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=10
 104     myfile << (double) gmshHypo->GetMeshCurvatureSize()   << std::endl;
 105 #elif
 106     myfile << -1.0   << std::endl; // dummy value writte for conformity
 107 #endif
 108     myfile << (double) gmshHypo->GetMaxSize()     << std::endl;
 109     myfile << (double) gmshHypo->GetMinSize()     << std::endl;
 110     myfile << (int) gmshHypo->GetSecondOrder()    << std::endl;
 111     myfile << (int) gmshHypo->GetUseIncomplElem() << std::endl;
 112     if ( gmshHypo->GetCompoundOnEntries().size() > 0 )
 113     {
 114       TCompound defCompounds =  gmshHypo->GetCompoundOnEntries();
 115       for (TCompound::const_iterator it = defCompounds.begin();  it != defCompounds.end(); ++it )
 116       {
 117         std::string token = (it !=--defCompounds.end()) ? "," : "";
 118         myfile << *it << token;
 119       }
 120       myfile << std::endl;
 121     }
 122     else
 123     {
 124       myfile << 0 << std::endl; // mark the absence of compounds
 125     }
 126   }
 127   else
 128   {
 129     //Mark tha absence of parameters in the file
 130     myfile << 0 << std::endl;
 131   }
 132   myfile.close();
 133 }
 134
 135 //
 136
 137 /**
 138  * @brief write in a binary file the orientation for each surface element of the mesh
 139  *
 140  * @param aMesh The mesh
 141  * @param aShape the shape associated to the mesh
 142  * @param output_file name of the binary file
 143  */
 144 void GMSHPlugin_GMSH_3D_Remote::exportElementOrientation(SMESH_Mesh& aMesh,
 145                                                           const TopoDS_Shape& aShape,
 146                                                           const std::string output_file)
 147 {
 148   SMESH_MesherHelper helper(aMesh);
 149   SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
 150   std::map<vtkIdType, bool> elemOrientation;
 151
 152   for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next())
 153   {
 154     const TopoDS_Shape& aShapeFace = exFa.Current();
 155     int faceID = aMesh.GetMeshDS()->ShapeToIndex( aShapeFace );
 156     bool isRev = false;
 157     if ( helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
 158       // IsReversedSubMesh() can work wrong on strongly curved faces,
 159       // so we use it as less as possible
 160       isRev = helper.IsReversedSubMesh( TopoDS::Face( aShapeFace ));
 161
 162     const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace );
 163     if ( !aSubMeshDSFace ) continue;
 164
 165     SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 166
 167     while ( iteratorElem->more() ) // loop on elements on a geom face
 168     {
 169       // check mesh face
 170       const SMDS_MeshElement* elem = iteratorElem->next();
 171       if ( !elem )
 172         error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 173       if ( elem->NbCornerNodes() != 3 )
 174         error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
 175       elemOrientation[elem->GetID()] = isRev;
 176     } // loop on elements on a face
 177   } // loop on faces of a SOLID or SHELL
 178
 179   {
 180     std::ofstream df(output_file, ios::out|ios::binary);
 181     int size = elemOrientation.size();
 182     df.write((char*)&size, sizeof(int));
 183     for(auto const& [id, orient]:elemOrientation)
 184     {
 185       df.write((char*)&id, sizeof(vtkIdType));
 186       df.write((char*)&orient, sizeof(bool));
 187     }
 188     df.close();
 189   }
 190 }
 191
 192 /**
 193  * @brief Compute mesh associate to shape
 194  *
 195  * @param aMesh The mesh
 196  * @param aShape The shape
 197  * @return true fi there are some error
 198  */
 199 bool GMSHPlugin_GMSH_3D_Remote::Compute(SMESH_Mesh&         aMesh,
 200                                            const TopoDS_Shape& aShape)
 201 {
 202   {
 203     SMESH_MeshLocker myLocker(&aMesh);
 204     SMESH_Hypothesis::Hypothesis_Status hypStatus;
 205     GMSHPlugin_GMSH_3D::CheckHypothesis(aMesh, aShape, hypStatus); //in this call the _hypothesis is defined!
 206   }
 207
 208   SMESH_ParallelMesh& aParMesh = dynamic_cast<SMESH_ParallelMesh&>(aMesh);
 209   // Temporary folder for run
 210 #ifdef WIN32
 211   // On windows mesh does not have GetTmpFolder
 212   fs::path tmp_folder = aParMesh.GetTmpFolder() / fs::path("Volume-%%%%-%%%%");
 213 #else
 214   fs::path tmp_folder = aParMesh.GetTmpFolder() / fs::unique_path(fs::path("Volume-%%%%-%%%%"));
 215 #endif
 216   fs::create_directories(tmp_folder);
 217   // Using MESH2D generated after all triangles where created.
 218   fs::path mesh_file= aParMesh.GetTmpFolder() / fs::path("Mesh2D.med");
 219   fs::path element_orientation_file=tmp_folder / fs::path("element_orientation.dat");
 220   fs::path new_element_file=tmp_folder / fs::path("new_elements.dat");
 221   //fs::path tmp_mesh_file=tmp_folder / fs::path("tmp_mesh.med");
 222   // Not used kept for debug
 223   // fs::path output_mesh_file=tmp_folder / fs::path("output_mesh.med");
 224   fs::path shape_file=tmp_folder / fs::path("shape.brep");
 225   fs::path param_file=tmp_folder / fs::path("gmsh_param.txt");
 226   fs::path log_file=tmp_folder / fs::path("run.log");
 227   fs::path cmd_file=tmp_folder / fs::path("cmd.txt");
 228
 229   std::string mesh_name = "MESH";
 230
 231   {
 232     SMESH_MeshLocker myLocker(&aMesh);
 233     //Writing Shape
 234     SMESH_DriverShape::exportShape(shape_file.string(), aShape);
 235
 236     //Writing hypothesis to file
 237     exportGmshParams(param_file.string(), _hypothesis);
 238
 239     // Exporting element orientation
 240     exportElementOrientation(aMesh, aShape, element_orientation_file.string());
 241   }
 242
 243   // Calling run_mesher
 244   // Path to mesher script
 245   fs::path mesher_launcher = fs::path(std::getenv("SMESH_ROOT_DIR"))/
 246        fs::path("bin")/
 247        fs::path("salome")/
 248        fs::path("mesher_launcher.py");
 249
 250   std::string s_program="python3";
 251   std::list<std::string> params;
 252   params.push_back(mesher_launcher.string());
 253   params.push_back("GMSH3D");
 254   params.push_back(mesh_file.string());
 255   params.push_back(shape_file.string());
 256   params.push_back(param_file.string());
 257   params.push_back("--elem-orient-file=" + element_orientation_file.string());
 258   params.push_back("--new-element-file=" + new_element_file.string());
 259   // params.push_back("--output-mesh-file=" + output_mesh_file.string());
 260
 261    // Parallelism method parameters
 262   int method = aParMesh.GetParallelismMethod();
 263   if(method == ParallelismMethod::MultiThread){
 264     params.push_back("--method=local");
 265   } else if (method == ParallelismMethod::MultiNode){
 266     params.push_back("--method=cluster");
 267     params.push_back("--resource="+aParMesh.GetResource());
 268     params.push_back("--wc-key="+aParMesh.GetWcKey());
 269     params.push_back("--nb-proc=1");
 270     params.push_back("--nb-proc-per-node="+std::to_string(aParMesh.GetNbProcPerNode()));
 271     params.push_back("--nb-node="+std::to_string(aParMesh.GetNbNode()));
 272     params.push_back("--walltime="+aParMesh.GetWalltime());
 273   } else {
 274     throw SALOME_Exception("Unknown parallelism method "+method);
 275   }
 276
 277   std::string cmd = "";
 278   cmd += s_program;
 279   for(auto arg: params){
 280     cmd += " " + arg;
 281   }
 282   MESSAGE("Running command: ");
 283   MESSAGE(cmd);
 284   // Writing command in cmd.log
 285   {
 286     std::ofstream flog(cmd_file.string());
 287     flog << cmd << endl;
 288   }
 289
 290    // Building arguments for QProcess
 291   QString program = QString::fromStdString(s_program);
 292   QStringList arguments;
 293   for(auto arg : params){
 294     arguments << arg.c_str();
 295   }
 296
 297   QString out_file = log_file.string().c_str();
 298   QProcess myProcess;
 299   // myProcess.setProcessChannelMode(QProcess::MergedChannels);
 300   myProcess.setProcessChannelMode(QProcess::ForwardedChannels);
 301   myProcess.setStandardOutputFile(out_file);
 302
 303   myProcess.start(program, arguments);
 304   // Waiting for process to finish (argument -1 make it wait until the end of
 305   // the process otherwise it just waits 30 seconds)
 306   bool finished = myProcess.waitForFinished(-1);
 307   int ret = myProcess.exitCode();
 308   if(ret != 0 || !finished){
 309     // Run crahed
 310     std::string msg = "Issue with mesh_launcher: \n";
 311     msg += "See log for more details: " + log_file.string() + "\n";
 312     msg += cmd + "\n";
 313     throw SALOME_Exception(msg);
 314   }
 315   {
 316     SMESH_MeshLocker myLocker(&aMesh);
 317     // Binary file written from SA version of the mesher
 318     std::ifstream df(new_element_file.string(), ios::binary);
 319
 320     int numberOfNodes;
 321     int totalNumberOfNodes;
 322     int numberOfVolumes;
 323     double meshNodes[3];
 324     int    volNodes[4];
 325     int nodeID;
 326
 327     SMESH_MesherHelper helper(aMesh);
 328     // This function is mandatory for setElementsOnShape to work
 329     helper.IsQuadraticSubMesh(aShape);
 330     helper.SetElementsOnShape( true );
 331
 332     // Number of nodes in intial mesh
 333     df.read((char*) &numberOfNodes, sizeof(int));
 334     // Number of nodes added by netgen
 335     df.read((char*) &totalNumberOfNodes, sizeof(int));
 336     // Filling nodevec (correspondence netgen numbering mesh numbering)
 337     std::vector< const SMDS_MeshNode* > nodeVec ( totalNumberOfNodes + 1 );
 338     SMESHDS_Mesh * meshDS = helper.GetMeshDS();
 339     for (int nodeIndex = 1 ; nodeIndex <= numberOfNodes; ++nodeIndex )
 340     {
 341       //Id of the point
 342       df.read((char*) &nodeID, sizeof(int));
 343       nodeVec.at(nodeIndex) = meshDS->FindNode(nodeID);
 344     }
 345
 346     // Add new points and update nodeVec
 347     for (int nodeIndex = numberOfNodes +1; nodeIndex <= totalNumberOfNodes; ++nodeIndex )
 348     {
 349       df.read((char *) &meshNodes, sizeof(double)*3);
 350       nodeVec.at(nodeIndex) = helper.AddNode(meshNodes[0], meshNodes[1], meshNodes[2]);
 351     }
 352
 353     // Add tetrahedrons
 354     df.read((char*) &numberOfVolumes, sizeof(int));
 355     for ( int elemIndex = 1; elemIndex <= numberOfVolumes; ++elemIndex )
 356     {
 357       df.read((char*) &volNodes, sizeof(int)*4);
 358       auto n0 = meshDS->FindNode(nodeVec[volNodes[0]]->GetID());
 359       auto n1 = meshDS->FindNode(nodeVec[volNodes[1]]->GetID());
 360       auto n2 = meshDS->FindNode(nodeVec[volNodes[2]]->GetID());
 361       auto n3 = meshDS->FindNode(nodeVec[volNodes[3]]->GetID());
 362       if ( n0 && n1 && n2 && n3 )
 363         helper.AddVolume( n0, n2, n1, n3 );
 364     }
 365   }
 366
 367   return true;
 368 }
 369
 370 /**
 371  * @brief Assign submeshes to compute
 372  *
 373  * @param aSubMesh submesh to add
 374  */
 375 void GMSHPlugin_GMSH_3D_Remote::setSubMeshesToCompute(SMESH_subMesh * aSubMesh)
 376 {
 377   SMESH_MeshLocker myLocker(aSubMesh->GetFather());
 378   SMESH_Algo::setSubMeshesToCompute(aSubMesh);
 379 }