src/GMSHPlugin/GMSHPlugin_GMSH_3D_Remote.cxx

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