1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 // File : GHS3DPRLPlugin_GHS3DPRL.cxx
22 // Author : Christian VAN WAMBEKE (CEA) (from Hexotic plugin Lioka RAZAFINDRAZAKA)
25 #include "GHS3DPRLPlugin_GHS3DPRL.hxx"
26 #include "GHS3DPRLPlugin_Hypothesis.hxx"
28 #include <SMESHDS_Mesh.hxx>
29 #include <SMESH_Gen.hxx>
30 #include <SMESH_subMesh.hxx>
32 #include <TopExp_Explorer.hxx>
33 #include <OSD_File.hxx>
35 #include "utilities.h"
38 #include <sys/sysinfo.h>
50 #include <BRepGProp.hxx>
51 #include <GProp_GProps.hxx>
52 #include <Standard_ProgramError.hxx>
53 #include <TCollection_AsciiString.hxx>
54 #include <TopTools_MapOfShape.hxx>
56 #include <TopoDS_Edge.hxx>
57 #include <TopoDS_Face.hxx>
61 static void removeFile( const TCollection_AsciiString& fileName )
64 OSD_File( fileName ).Remove();
66 catch ( Standard_ProgramError ) {
67 MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
71 //=============================================================================
72 GHS3DPRLPlugin_GHS3DPRL::GHS3DPRLPlugin_GHS3DPRL(int hypId, int studyId, SMESH_Gen* gen)
73 : SMESH_3D_Algo(hypId, studyId, gen)
75 MESSAGE("GHS3DPRLPlugin_GHS3DPRL::GHS3DPRLPlugin_GHS3DPRL");
76 _name = "MG-Tetra Parallel";
77 _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
80 _compatibleHypothesis.push_back(GHS3DPRLPlugin_Hypothesis::GetHypType());
83 //=============================================================================
84 GHS3DPRLPlugin_GHS3DPRL::~GHS3DPRLPlugin_GHS3DPRL()
86 MESSAGE("GHS3DPRLPlugin_GHS3DPRL::~GHS3DPRLPlugin_GHS3DPRL");
89 //=============================================================================
90 bool GHS3DPRLPlugin_GHS3DPRL::CheckHypothesis
92 const TopoDS_Shape& aShape,
93 SMESH_Hypothesis::Hypothesis_Status& aStatus)
95 //MESSAGE("GHS3DPRLPlugin_GHS3DPRL::CheckHypothesis");
98 list<const SMESHDS_Hypothesis*>::const_iterator itl;
99 const SMESHDS_Hypothesis* theHyp;
101 const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
102 int nbHyp = hyps.size();
105 aStatus = SMESH_Hypothesis::HYP_OK;
106 return true; // can work with no hypothesis
110 theHyp = (*itl); // use only the first hypothesis
112 string hypName = theHyp->GetName();
113 if (hypName == GHS3DPRLPlugin_Hypothesis::GetHypType())
115 _hypothesis = static_cast<const GHS3DPRLPlugin_Hypothesis*> (theHyp);
117 aStatus = SMESH_Hypothesis::HYP_OK;
120 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
122 return aStatus == SMESH_Hypothesis::HYP_OK;
125 //=======================================================================
126 // static bool writeGHS3DPRLFiles (const TCollection_AsciiString & GHS3DPRL_In,
127 // SMESHDS_Mesh * theMesh,
128 // map <int,int> & theSmdsToGHS3DPRLIdMap,
129 // map <int,const SMDS_MeshNode*> & theGHS3DPRLIdToNodeMap)
133 // TCollection_AsciiString namefile(GHS3DPRL_In);
134 // namefile+=".points";
135 // removeFile(namefile);
137 // theFile.open(namefile.ToCString(),ios::out);
139 // Ok = theFile.is_open();
141 // Ok=theFile.rdbuf()->is_open();
145 // INFOS("Can't write into "<<namefile.ToCString());
148 // cout<<endl<<"writeGHS3DPRLFiles version tetra_hpc v2.1.11 "<<endl;
149 // cout<<endl<<"Creating GHS3DPRL processed mesh file : "<<namefile<<endl;
151 // int nbVertices=theMesh->NbNodes();
152 // int nbFaces=theMesh->NbFaces(); //triangles or quadrangles
153 // const char* space=" ";
154 // //const int dummyint=1; //nrs,nsd,refnum=1 (for wrap)
156 // // Writing SMESH points into GHS3DPRL File.points
157 // theFile<<nbVertices<<endl;
159 // int aSmdsNodeID=1;
160 // const SMDS_MeshNode* node_2;
161 // SMDS_NodeIteratorPtr itOnNode=theMesh->nodesIterator();
162 // //int ifam=100;//test famille
163 // theFile.precision(15); theFile.setf(ios::scientific,ios::floatfield);
164 // //cout<<"set precision 15 on float\n";
165 // while (itOnNode->more())
167 // node_2 = itOnNode->next();
168 // theSmdsToGHS3DPRLIdMap.insert(map <int,int>::value_type(node_2->GetID(),aSmdsNodeID));
169 // theGHS3DPRLIdToNodeMap.insert(map <int,const SMDS_MeshNode*>::value_type(aSmdsNodeID,node_2));
170 // theFile<<node_2->X()<<space<<node_2->Y()<<space<<node_2->Z()<<space<<ifam<<endl;
172 // //if (aSmdsNodeID==11) ifam++;
174 // //no specified points;
177 // namefile=GHS3DPRL_In+".faces";
178 // removeFile(namefile);
179 // theFile.open(namefile.ToCString(),ios::out);
181 // Ok=theFile.is_open();
183 // Ok=theFile.rdbuf()->is_open();
187 // INFOS("Can't write into "<<namefile.ToCString());
190 // cout<<endl<<"Creating GHS3DPRL processed mesh file : "<<namefile<<endl;
192 // // Writing SMESH faces into GHS3DPRL File.faces
193 // theFile<<nbFaces<<" 0"<<endl; //NB_ELEMS DUMMY_INT
194 // //" 0" is a reserved parameter
196 // const SMDS_MeshElement* aFace;
197 // map<int,int>::const_iterator itOnSmdsNode;
198 // SMDS_ElemIteratorPtr itOnFaceNode;
199 // SMDS_FaceIteratorPtr itOnSmdsFace = theMesh->facesIterator();
200 // long nbNoTriangles=0;
203 // while (itOnSmdsFace->more())
205 // aFace=itOnSmdsFace->next();
206 // itOnFaceNode=aFace->nodesIterator();
207 // const int nbNodes=aFace->NbNodes();
208 // if (nbNodes!=3) nbNoTriangles++;
210 // theFile<<nbNodes<<space; // NB_NODES
211 // while (itOnFaceNode->more())
213 // aSmdsNodeID=itOnFaceNode->next()->GetID();
214 // itOnSmdsNode=theSmdsToGHS3DPRLIdMap.find(aSmdsNodeID);
215 // ASSERT(itOnSmdsNode!=theSmdsToGHS3DPRLIdMap.end());
216 // theFile<<space<<(*itOnSmdsNode).second; //NODE_1 NODE_2 ...
218 // //(NB_NODES+1) times: DUMMY_INT
219 // //if (ifaces==11) ifam++;
220 // theFile<<space<<ifam;
221 // for ( int i=1; i<=nbNodes; i++) theFile<<space<<200+i;
226 // cout<<"Processed mesh files created, they contains :\n";
227 // cout<<" "<<nbVertices<<" vertices\n";
228 // if (nbNoTriangles==0)
229 // cout<<" "<<nbFaces<<" faces\n\n";
231 // cout<<" "<<nbFaces<<" faces with "<<nbNoTriangles<<"faces no triangles\n\n";
236 //=======================================================================
238 #define GHS3DPRLPlugin_BUFLENGTH 256
239 #define GHS3DPRLPlugin_ReadLine(aPtr,aBuf,aFile,aLineNb) \
240 { aPtr = fgets( aBuf, GHS3DPRLPlugin_BUFLENGTH - 2, aFile ); aLineNb++; DUMP(endl); }
243 //=============================================================================
244 // Pass parameters to GHS3DPRL
245 void GHS3DPRLPlugin_GHS3DPRL::SetParameters(const GHS3DPRLPlugin_Hypothesis* hyp)
248 MESSAGE("GHS3DPRLPlugin_GHS3DPRL::SetParameters");
249 _MEDName = hyp->GetMEDName(); //"DOMAIN\0"
250 _NbPart = hyp->GetNbPart();
251 _KeepFiles = hyp->GetKeepFiles();
252 _Background = hyp->GetBackground();
253 _Multithread = hyp->GetMultithread();
254 _Gradation = hyp->GetGradation();
255 _MinSize = hyp->GetMinSize();
256 _MaxSize = hyp->GetMaxSize();
257 _AdvOptions = hyp->GetAdvancedOption();
261 //=======================================================================
262 //before launching salome
263 //SALOME_TMP_DIR (for keep tepal intermediates files) could be set in user's directories
264 static TCollection_AsciiString getTmpDir()
266 TCollection_AsciiString aTmpDir;
267 char *Tmp_dir = getenv("SALOME_TMP_DIR");
268 if (Tmp_dir == NULL) Tmp_dir = getenv("TMP");
273 if(aTmpDir.Value(aTmpDir.Length()) != '\\') aTmpDir+='\\';
275 if(aTmpDir.Value(aTmpDir.Length()) != '/') aTmpDir+='/';
281 aTmpDir = TCollection_AsciiString("C:\\");
283 aTmpDir = TCollection_AsciiString("/tmp/");
289 //=============================================================================
290 // Here we are going to use the GHS3DPRL mesher for tetra-hpc (formerly tepal in v3 (2014))
291 bool GHS3DPRLPlugin_GHS3DPRL::Compute(SMESH_Mesh& theMesh,
292 const TopoDS_Shape& theShape)
295 TCollection_AsciiString pluginerror("ghs3dprl: ");
296 SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
297 //cout<<"GetMeshDS done\n";
298 if (_countSubMesh==0){
299 MESSAGE("GHS3DPRLPlugin_GHS3DPRL::Compute for tetra-hpc");
301 TopExp_Explorer expf(meshDS->ShapeToMesh(), TopAbs_SOLID);
302 for ( ; expf.More(); expf.Next() ) _countTotal++;
305 //cout<<"Compute _countSubMesh "<<_countSubMesh<<endl;
306 //no stuff if multiples submesh, multiple call compute
307 //mesh all in one pass tepal (the last)
308 if (_countSubMesh != _countTotal ) return true;
311 if (_hypothesis==NULL){
312 pluginerror += "No existing parameters/hypothesis for GHS3DPRL";
313 cout <<"\n"<<pluginerror<<"\n\n";
314 error(COMPERR_ALGO_FAILED, pluginerror.ToCString() );
317 SetParameters(_hypothesis);
318 cout << "\n" << _name << " parameters:" << endl;
319 cout << " generic path/name of MED Files = " << _MEDName << endl;
320 cout << " number of partitions = " << _NbPart << endl;
321 cout << " gradation = " << _Gradation << endl;
322 cout << " min_size = " << _MinSize << endl;
323 cout << " max_size = " << _MaxSize << endl;
324 cout << " keep intermediates files (from tetra-hpc) = " << _KeepFiles << endl;
325 cout << " background (from tetra-hpc) = " << _Background << "\n";
326 cout << " multithread = " << _Multithread << "\n\n";
327 cout << " adv. options = '" << _AdvOptions << "'\n\n";
329 //string tmpDir=getTmpDir_new();
330 TCollection_AsciiString
335 run_GHS3DPRL("tetrahpc2med "),
345 casenamemed; //_MEDName.c_str());
347 casenamemed += (char *)_MEDName.c_str();
348 int n=casenamemed.SearchFromEnd('/');
350 path=casenamemed.SubString(1,n);
351 casenamemed=casenamemed.SubString(n+1,casenamemed.Length());
356 if (casenamemed.Length()>20){
357 casenamemed=casenamemed.SubString(1,20);
358 cerr<<"MEDName truncated (no more 20 characters) = "<<casenamemed<<endl;
360 //cout<<"path="<<path<<endl;
361 //cout<<"casenamemed="<<casenamemed<<endl;
363 map <int,int> aSmdsToGHS3DPRLIdMap;
364 map <int,const SMDS_MeshNode*> aGHS3DPRLIdToNodeMap;
365 GHS3DPRL_In = path + "GHS3DPRL";
366 GHS3DPRL_Out = path + casenamemed;
367 GHS3DPRL_Outxml = path + casenamemed + ".xml"; //master file
369 Gradation=_Gradation;
374 //tetrahpc2med --casename=/home/whoami/tmp/GHS3DPRL --number=5 --medname=DOMAIN
375 // --gradation=1.05 --min_size=1e-3 --max_size=1e-2
376 // --verbose=0 --menu=no --launchtetra=yes;
378 run_GHS3DPRL = run_GHS3DPRL +
379 " --casename=" + GHS3DPRL_In +
380 " --number=" + NbPart +
381 " --medname=" + GHS3DPRL_Out +
382 " --launchtetra=yes" +
383 " --gradation=" + Gradation +
384 " --min_size=" + MinSize +
385 " --max_size=" + MaxSize +
387 " " + _AdvOptions.c_str();
388 if (_Background) run_GHS3DPRL += " --background=yes"; else run_GHS3DPRL += " --background=no";
389 if (_Multithread) run_GHS3DPRL += " --multithread=yes"; else run_GHS3DPRL += " --multithread=no";
390 run_nokeep_files = rm +GHS3DPRL_In + "* " + path + "tetrahpc.log";
391 system( run_nokeep_files.ToCString() ); //clean files
392 run_nokeep_files = rm + GHS3DPRL_In + "* ";
394 cout<<"GHS3DPRL command :\n "<<run_GHS3DPRL.ToCString()<<endl;
395 fileskinmesh=path + "GHS3DPRL.mesh";
396 cout<<" Write input file for tetra_hpc.exe "<<fileskinmesh<<"...";
397 GHS3DPRL_Out = path + casenamemed;
398 removeFile( GHS3DPRL_Outxml ); //only the master xml file
399 //Ok=writeGHS3DPRLFiles(GHS3DPRL_In, meshDS, aSmdsToGHS3DPRLIdMap, aGHS3DPRLIdToNodeMap);
400 bool toCreateGroups = false;
401 theMesh.ExportGMF(fileskinmesh.ToCString(), meshDS, toCreateGroups );
404 // cout<<" ...NOT done\n";
405 // pluginerror = pluginerror + "problem writing input tepal files " + GHS3DPRL_In + ".mesh";
407 //Ecriture dans un fichier MED ?v2.? meme si not Ok
408 //create empty file -> avoid warning message
409 //med_idt fid=MEDouvrir((const char *)fileskinmed.ToCString(),MED_CREATION);
410 //med_err ret=MEDfermer(fid);
411 //fileskinmed=fileskinmed + "cp /home/wambeke/empty.med "+ path + "GHS3DPRL_skin.med";
412 //system( fileskinmed.ToCString() );
413 fileskinmed=path + "GHS3DPRL_skin.med";
414 cout<<" Write file "<<fileskinmed<<"...";
415 theMesh.ExportMED(fileskinmed.ToCString(),"SKIN_INITIAL",true,1);
416 cout<<" ...done\n\n";
419 //sometimes it is better to wait flushing files on slow filesystem...
421 //launch tetrahpc2med which launch mg-tetra_hpc.py which launch mg-tetra_hpc(_mpi?).exe
422 system( run_GHS3DPRL.ToCString() );
426 pluginerror = pluginerror + "backgrounding... plugin is not waiting for output files "+ casenamemed + "_*.med";
427 cout<<pluginerror<<endl;
428 error(COMPERR_ALGO_FAILED, pluginerror.ToCString());
429 return false; //but it is not a problem but if true my message is overwritten
430 //return true; //but it is not a problem,
433 // read a result, GHS3DPRL_Out is the name of master file (previous xml format)
434 FILE * aResultFile = fopen( GHS3DPRL_Outxml.ToCString(), "r" );
436 //Ok = readResult( aResultFile, meshDS, theShape, aGHS3DPRLIdToNodeMap, GHS3DPRL_Out, _nodeRefNumber );
438 Ok = false; //but it is not a problem but if true my message is overwritten
440 cout<<"GHS3DPRL OK output master file "<<casenamemed<<".xml exist !\n\n";
441 pluginerror = pluginerror + "new tetraedra not in memory, but stored in files "+ casenamemed + "_*.med";
442 cout<<pluginerror<<endl;
443 error(COMPERR_ALGO_FAILED, pluginerror.ToCString());
444 if (!_KeepFiles) system( run_nokeep_files.ToCString() );
447 Ok = false; //it is a problem AND my message is NOT overwritten
448 pluginerror = pluginerror + "output master file " + casenamemed + ".xml do not exist";
449 cout<<pluginerror<<endl;
450 error(COMPERR_ALGO_FAILED, pluginerror.ToCString() );
451 cout<<"GHS3DPRL KO output files "<<GHS3DPRL_Out<<" do not exist ! see intermediates files keeped:\n";
452 TCollection_AsciiString run_list_files("ls -alt ");
453 run_list_files += GHS3DPRL_Out + "* " + GHS3DPRL_In + "* " + path + "tetrahpc.log";
454 system( run_list_files.ToCString() );
465 //=============================================================================
469 //=============================================================================
470 bool GHS3DPRLPlugin_GHS3DPRL::Evaluate(SMESH_Mesh& aMesh,
471 const TopoDS_Shape& aShape,
472 MapShapeNbElems& aResMap)
474 int nbtri = 0, nbqua = 0;
475 double fullArea = 0.0;
476 for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
477 TopoDS_Face F = TopoDS::Face( exp.Current() );
478 SMESH_subMesh *sm = aMesh.GetSubMesh(F);
479 MapShapeNbElemsItr anIt = aResMap.find(sm);
480 if( anIt==aResMap.end() ) {
481 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
482 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
483 "Submesh can not be evaluated",this));
486 std::vector<int> aVec = (*anIt).second;
487 nbtri += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
488 nbqua += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
490 BRepGProp::SurfaceProperties(F,G);
491 double anArea = G.Mass();
495 // collect info from edges
496 int nb0d_e = 0, nb1d_e = 0;
497 bool IsQuadratic = false;
499 TopTools_MapOfShape tmpMap;
500 for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
501 const TopoDS_Edge& E = TopoDS::Edge(exp.Current());
502 if( tmpMap.Contains(E) )
505 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
506 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
507 std::vector<int> aVec = (*anIt).second;
508 nb0d_e += aVec[SMDSEntity_Node];
509 nb1d_e += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
511 IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
517 double ELen = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
520 BRepGProp::VolumeProperties(aShape,G);
521 double aVolume = G.Mass();
522 double tetrVol = 0.1179*ELen*ELen*ELen;
523 double CoeffQuality = 0.9;
524 int nbVols = (int)aVolume/tetrVol/CoeffQuality;
525 int nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
526 int nb1d_in = (int) ( nbVols*6 - nb1d_e - nb1d_f ) / 5;
527 std::vector<int> aVec(SMDSEntity_Last);
528 for(int i=0; i<SMDSEntity_Last; i++) aVec[i]=0;
530 aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
531 aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
532 aVec[SMDSEntity_Quad_Pyramid] = nbqua;
535 aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
536 aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
537 aVec[SMDSEntity_Pyramid] = nbqua;
539 SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
540 aResMap.insert(std::make_pair(sm,aVec));