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
  28 #include <SMESHDS_Mesh.hxx>
  29 #include <SMESH_Gen.hxx>
  30 #include <SMESH_subMesh.hxx>
  31
  32 #include <TopExp_Explorer.hxx>
  33 #include <OSD_File.hxx>
  34
  35 #include "utilities.h"
  36
  37 #ifndef WIN32
  38 #include <sys/sysinfo.h>
  39 #endif
  40
  41 #ifdef _DEBUG_
  42 #define DUMP(txt) \
  43 //  cout << txt
  44 #else
  45 #define DUMP(txt)
  46 #endif
  47
  48 #include <list>
  49
  50 #include <BRepGProp.hxx>
  51 #include <GProp_GProps.hxx>
  52 #include <Standard_ProgramError.hxx>
  53 #include <TCollection_AsciiString.hxx>
  54 #include <TopTools_MapOfShape.hxx>
  55 #include <TopoDS.hxx>
  56 #include <TopoDS_Edge.hxx>
  57 #include <TopoDS_Face.hxx>
  58
  59 using namespace std;
  60
  61 static void removeFile( const TCollection_AsciiString& fileName )
  62 {
  63   try {
  64     OSD_File( fileName ).Remove();
  65   }
  66   catch ( Standard_ProgramError ) {
  67     MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
  68   }
  69 }
  70
  71 //=============================================================================
  72 GHS3DPRLPlugin_GHS3DPRL::GHS3DPRLPlugin_GHS3DPRL(int hypId, int studyId, SMESH_Gen* gen)
  73   : SMESH_3D_Algo(hypId, studyId, gen)
  74 {
  75   MESSAGE("GHS3DPRLPlugin_GHS3DPRL::GHS3DPRLPlugin_GHS3DPRL");
  76   _name = "MG-Tetra Parallel";
  77   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  78   _countSubMesh=0;
  79   _nodeRefNumber=0;
  80   _compatibleHypothesis.push_back(GHS3DPRLPlugin_Hypothesis::GetHypType());
  81 }
  82
  83 //=============================================================================
  84 GHS3DPRLPlugin_GHS3DPRL::~GHS3DPRLPlugin_GHS3DPRL()
  85 {
  86   MESSAGE("GHS3DPRLPlugin_GHS3DPRL::~GHS3DPRLPlugin_GHS3DPRL");
  87 }
  88
  89 //=============================================================================
  90 bool GHS3DPRLPlugin_GHS3DPRL::CheckHypothesis
  91                          (SMESH_Mesh& aMesh,
  92                           const TopoDS_Shape& aShape,
  93                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
  94 {
  95   //MESSAGE("GHS3DPRLPlugin_GHS3DPRL::CheckHypothesis");
  96   _hypothesis = NULL;
  97
  98   list<const SMESHDS_Hypothesis*>::const_iterator itl;
  99   const SMESHDS_Hypothesis* theHyp;
 100
 101   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
 102   int nbHyp = hyps.size();
 103   if (!nbHyp)
 104   {
 105     aStatus = SMESH_Hypothesis::HYP_OK;
 106     return true;  // can work with no hypothesis
 107   }
 108
 109   itl = hyps.begin();
 110   theHyp = (*itl); // use only the first hypothesis
 111
 112   string hypName = theHyp->GetName();
 113   if (hypName == GHS3DPRLPlugin_Hypothesis::GetHypType())
 114   {
 115     _hypothesis = static_cast<const GHS3DPRLPlugin_Hypothesis*> (theHyp);
 116     ASSERT(_hypothesis);
 117     aStatus = SMESH_Hypothesis::HYP_OK;
 118   }
 119   else
 120     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 121
 122   return aStatus == SMESH_Hypothesis::HYP_OK;
 123 }
 124
 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)
 130 // {
 131 //    bool Ok;
 132 //    int ifam=0;
 133 //    TCollection_AsciiString namefile(GHS3DPRL_In);
 134 //    namefile+=".points";
 135 //    removeFile(namefile);
 136 //    ofstream theFile;
 137 //    theFile.open(namefile.ToCString(),ios::out);
 138 // #ifdef WIN32
 139 //    Ok = theFile.is_open();
 140 // #else
 141 //    Ok=theFile.rdbuf()->is_open();
 142 // #endif
 143 //    if (!Ok)
 144 //    {
 145 //       INFOS("Can't write into "<<namefile.ToCString());
 146 //       return false;
 147 //    }
 148 //    cout<<endl<<"writeGHS3DPRLFiles version tetra_hpc v2.1.11 "<<endl;
 149 //    cout<<endl<<"Creating GHS3DPRL processed mesh file : "<<namefile<<endl;
 150
 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)
 155
 156 //    // Writing SMESH points into GHS3DPRL File.points
 157 //    theFile<<nbVertices<<endl;
 158
 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())
 166 //    {
 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;
 171 //       aSmdsNodeID++;
 172 //       //if (aSmdsNodeID==11) ifam++;
 173 //    }
 174 //    //no specified points;
 175 //    theFile.close();
 176
 177 //    namefile=GHS3DPRL_In+".faces";
 178 //    removeFile(namefile);
 179 //    theFile.open(namefile.ToCString(),ios::out);
 180 // #ifdef WIN32
 181 //    Ok=theFile.is_open();
 182 // #else
 183 //    Ok=theFile.rdbuf()->is_open();
 184 // #endif
 185 //    if (!Ok)
 186 //    {
 187 //       INFOS("Can't write into "<<namefile.ToCString());
 188 //       return false;
 189 //    }
 190 //    cout<<endl<<"Creating GHS3DPRL processed mesh file : "<<namefile<<endl;
 191
 192 //    // Writing SMESH faces into GHS3DPRL File.faces
 193 //    theFile<<nbFaces<<" 0"<<endl;   //NB_ELEMS DUMMY_INT
 194 //                                    //" 0" is a reserved parameter
 195
 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;
 201 //    int ifaces=0;
 202 //    //ifam=300;
 203 //    while (itOnSmdsFace->more())
 204 //    {
 205 //       aFace=itOnSmdsFace->next();
 206 //       itOnFaceNode=aFace->nodesIterator();
 207 //       const int nbNodes=aFace->NbNodes();
 208 //       if (nbNodes!=3) nbNoTriangles++;
 209 //       ifaces++;
 210 //       theFile<<nbNodes<<space;        // NB_NODES
 211 //       while (itOnFaceNode->more())
 212 //       {
 213 //           aSmdsNodeID=itOnFaceNode->next()->GetID();
 214 //           itOnSmdsNode=theSmdsToGHS3DPRLIdMap.find(aSmdsNodeID);
 215 //           ASSERT(itOnSmdsNode!=theSmdsToGHS3DPRLIdMap.end());
 216 //           theFile<<space<<(*itOnSmdsNode).second; //NODE_1 NODE_2 ...
 217 //       }
 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;
 222 //       theFile<<endl;
 223 //    }
 224 //    theFile.close();
 225
 226 //    cout<<"Processed mesh files created, they contains :\n";
 227 //    cout<<"    "<<nbVertices<<" vertices\n";
 228 //    if (nbNoTriangles==0)
 229 //       cout<<"    "<<nbFaces<<" faces\n\n";
 230 //    else
 231 //       cout<<"    "<<nbFaces<<" faces with "<<nbNoTriangles<<"faces no triangles\n\n";
 232 //    return true;
 233 // }
 234
 235
 236 //=======================================================================
 237
 238 #define GHS3DPRLPlugin_BUFLENGTH 256
 239 #define GHS3DPRLPlugin_ReadLine(aPtr,aBuf,aFile,aLineNb) \
 240 {  aPtr = fgets( aBuf, GHS3DPRLPlugin_BUFLENGTH - 2, aFile ); aLineNb++; DUMP(endl); }
 241
 242
 243 //=============================================================================
 244 // Pass parameters to GHS3DPRL
 245 void GHS3DPRLPlugin_GHS3DPRL::SetParameters(const GHS3DPRLPlugin_Hypothesis* hyp)
 246 {
 247   if (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();
 258   }
 259 }
 260
 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()
 265 {
 266   TCollection_AsciiString aTmpDir;
 267   char *Tmp_dir = getenv("SALOME_TMP_DIR");
 268   if (Tmp_dir == NULL) Tmp_dir = getenv("TMP");
 269   if(Tmp_dir != NULL)
 270   {
 271     aTmpDir = Tmp_dir;
 272 #ifdef WIN32
 273       if(aTmpDir.Value(aTmpDir.Length()) != '\\') aTmpDir+='\\';
 274 #else
 275       if(aTmpDir.Value(aTmpDir.Length()) != '/') aTmpDir+='/';
 276 #endif
 277   }
 278   else
 279   {
 280 #ifdef WIN32
 281       aTmpDir = TCollection_AsciiString("C:\\");
 282 #else
 283       aTmpDir = TCollection_AsciiString("/tmp/");
 284 #endif
 285   }
 286   return aTmpDir;
 287 }
 288
 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)
 293 {
 294   bool Ok=false;
 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");
 300     _countTotal=0;
 301     TopExp_Explorer expf(meshDS->ShapeToMesh(), TopAbs_SOLID);
 302     for ( ; expf.More(); expf.Next() ) _countTotal++;
 303   }
 304   _countSubMesh++;
 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;
 309
 310   //stuff on last call
 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() );
 315     return false;
 316   }
 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";
 328
 329   //string tmpDir=getTmpDir_new();
 330   TCollection_AsciiString
 331     tmpDir=getTmpDir(),
 332     GHS3DPRL_In,
 333     GHS3DPRL_Out,
 334     GHS3DPRL_Outxml,
 335     run_GHS3DPRL("tetrahpc2med "),
 336     rm("rm -f "),
 337     run_nokeep_files,
 338     NbPart,
 339     Gradation,
 340     MinSize,
 341     MaxSize,
 342     fileskinmed(""),
 343     fileskinmesh(""),
 344     path,
 345     casenamemed;  //_MEDName.c_str());
 346
 347   casenamemed += (char *)_MEDName.c_str();
 348   int n=casenamemed.SearchFromEnd('/');
 349   if (n>0) {
 350     path=casenamemed.SubString(1,n);
 351     casenamemed=casenamemed.SubString(n+1,casenamemed.Length());
 352   }
 353   else
 354     path=tmpDir;
 355
 356   if (casenamemed.Length()>20){
 357     casenamemed=casenamemed.SubString(1,20);
 358     cerr<<"MEDName truncated (no more 20 characters) = "<<casenamemed<<endl;
 359   }
 360   //cout<<"path="<<path<<endl;
 361   //cout<<"casenamemed="<<casenamemed<<endl;
 362
 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
 368   NbPart=_NbPart;
 369   Gradation=_Gradation;
 370   MinSize=_MinSize;
 371   MaxSize=_MaxSize;
 372
 373   //an example:
 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;
 377
 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 +
 386     " --verbose=3" +
 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 + "* ";
 393
 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 );
 402   cout<<" ...done\n";
 403   //else {
 404   //   cout<<" ...NOT done\n";
 405   //   pluginerror = pluginerror + "problem writing input tepal files " + GHS3DPRL_In + ".mesh";
 406
 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";
 417
 418
 419   //sometimes it is better to wait flushing files on slow filesystem...
 420   system( "sleep 3" );
 421   //launch tetrahpc2med which launch mg-tetra_hpc.py which launch mg-tetra_hpc(_mpi?).exe
 422   system( run_GHS3DPRL.ToCString() );
 423   system( "sleep 3" );
 424
 425   if (_Background) {
 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,
 431   }
 432
 433   // read a result, GHS3DPRL_Out is the name of master file (previous xml format)
 434   FILE * aResultFile = fopen( GHS3DPRL_Outxml.ToCString(), "r" );
 435   if (aResultFile){
 436     //Ok = readResult( aResultFile, meshDS, theShape, aGHS3DPRLIdToNodeMap, GHS3DPRL_Out, _nodeRefNumber );
 437     Ok = true;
 438     Ok = false; //but it is not a problem but if true my message is overwritten
 439     fclose(aResultFile);
 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() );
 445   }
 446   else{
 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() );
 455     cout<<endl;
 456   }
 457   _countSubMesh=0;
 458
 459
 460   return Ok;
 461 }
 462
 463
 464
 465 //=============================================================================
 466 /*!
 467  *
 468  */
 469 //=============================================================================
 470 bool GHS3DPRLPlugin_GHS3DPRL::Evaluate(SMESH_Mesh& aMesh,
 471                                        const TopoDS_Shape& aShape,
 472                                        MapShapeNbElems& aResMap)
 473 {
 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));
 484       return false;
 485     }
 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]);
 489     GProp_GProps G;
 490     BRepGProp::SurfaceProperties(F,G);
 491     double anArea = G.Mass();
 492     fullArea += anArea;
 493   }
 494
 495   // collect info from edges
 496   int nb0d_e = 0, nb1d_e = 0;
 497   bool IsQuadratic = false;
 498   bool IsFirst = true;
 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) )
 503       continue;
 504     tmpMap.Add(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]);
 510     if(IsFirst) {
 511       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 512       IsFirst = false;
 513     }
 514   }
 515   tmpMap.Clear();
 516
 517   double ELen = sqrt(2.* ( fullArea/(nbtri+nbqua*2) ) / sqrt(3.0) );
 518
 519   GProp_GProps G;
 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;
 529   if( IsQuadratic ) {
 530     aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
 531     aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
 532     aVec[SMDSEntity_Quad_Pyramid] = nbqua;
 533   }
 534   else {
 535     aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
 536     aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
 537     aVec[SMDSEntity_Pyramid] = nbqua;
 538   }
 539   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
 540   aResMap.insert(std::make_pair(sm,aVec));
 541
 542   return true;
 543 }