src/GMSHPlugin/GMSHPlugin_Mesher.cxx

   1 // Copyright (C) 2012-2015  ALNEOS
   2 // Copyright (C) 2016-2023  EDF
   3 //
   4 // This library is free software; you can redistribute it and/or
   5 // modify it under the terms of the GNU Lesser General Public
   6 // License as published by the Free Software Foundation; either
   7 // version 2.1 of the License, or (at your option) any later version.
   8 //
   9 // This library is distributed in the hope that it will be useful,
  10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12 // Lesser General Public License for more details.
  13 //
  14 // You should have received a copy of the GNU Lesser General Public
  15 // License along with this library; if not, write to the Free Software
  16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  17 //
  18 // See http://www.alneos.com/ or email : contact@alneos.fr
  19 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 //
  21 #include "GMSHPlugin_Mesher.hxx"
  22 #include "GMSHPlugin_Hypothesis_2D.hxx"
  23
  24 #include <SMDS_MeshElement.hxx>
  25 #include <SMDS_MeshNode.hxx>
  26 #include <SMESHDS_Mesh.hxx>
  27 #include <SMESH_Comment.hxx>
  28 #include <SMESH_ComputeError.hxx>
  29 #include <SMESH_Gen_i.hxx>
  30 #include <SMESH_Mesh.hxx>
  31 #include <SMESH_MesherHelper.hxx>
  32 #include <SMESH_subMesh.hxx>
  33 #include <StdMeshers_FaceSide.hxx>
  34 #include <utilities.h>
  35 #include <StdMeshers_QuadToTriaAdaptor.hxx>
  36
  37 //CAS
  38 #include <BRep_Tool.hxx>
  39 #include <GeomAPI_ProjectPointOnCurve.hxx>
  40 #include <gp_Pnt.hxx>
  41
  42 #include <gmsh.h>
  43 #include <vector>
  44 #include <limits>
  45
  46 #include <TopExp_Explorer.hxx>
  47 #include <TopoDS.hxx>
  48
  49 #include <MLine.h>
  50 #include <MTriangle.h>
  51 #include <MQuadrangle.h>
  52 #if GMSH_MAJOR_VERSION >=4
  53 #include <GmshGlobal.h>
  54 #include <gmsh/Context.h>
  55 #endif
  56
  57 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
  58 #include <omp.h>
  59 #endif
  60
  61 using namespace std;
  62
  63 namespace
  64 {
  65   struct ShapeBounds
  66   {
  67     SBoundingBox3d _bounds;
  68     TopoDS_Shape   _shape;
  69   };
  70
  71   //================================================================================
  72   /*!
  73    * \brief Retrieve ShapeBounds from a compound GEdge
  74    */
  75   //================================================================================
  76
  77   bool getBoundsOfShapes( GEdge*                       gEdge,
  78                           std::vector< ShapeBounds > & topoEdges )
  79   {
  80     topoEdges.clear();
  81 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
  82     for ( size_t i = 0; i < gEdge->compound.size(); ++i )
  83     {
  84       GEdge* gE = static_cast< GEdge* >( gEdge->compound[ i ]);
  85       topoEdges.push_back( ShapeBounds{ gE->bounds(), *((TopoDS_Edge*)gE->getNativePtr()) });
  86     }
  87 #else
  88     for ( size_t i = 0; i < gEdge->_compound.size(); ++i )
  89     {
  90       GEdge* gE = static_cast< GEdge* >( gEdge->_compound[ i ]);
  91       topoEdges.push_back( ShapeBounds{ gE->bounds(), *((TopoDS_Edge*)gE->getNativePtr()) });
  92     }
  93 #endif
  94     return topoEdges.size();
  95   }
  96
  97   //================================================================================
  98   /*!
  99    * \brief Retrieve ShapeBounds from a compound GFace
 100    */
 101   //================================================================================
 102
 103   bool getBoundsOfShapes( GFace*                       gFace,
 104                           std::vector< ShapeBounds > & topoFaces )
 105   {
 106     topoFaces.clear();
 107 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 108     for ( size_t i = 0; i < gFace->compound.size(); ++i )
 109     {
 110       GFace* gF = static_cast< GFace* >( gFace->compound[ i ]);
 111       topoFaces.push_back( ShapeBounds{ gF->bounds(), *((TopoDS_Face*)gF->getNativePtr()) });
 112     }
 113 #else
 114     for ( size_t i = 0; i < gFace->_compound.size(); ++i )
 115     {
 116       GFace* gF = static_cast< GFace* >( gFace->_compound[ i ]);
 117       topoFaces.push_back( ShapeBounds{ gF->bounds(), *((TopoDS_Face*)gF->getNativePtr()) });
 118     }
 119 #endif
 120     return topoFaces.size();
 121   }
 122   //================================================================================
 123   /*!
 124    * \brief Find a shape whose bounding box includes a given point
 125    */
 126   //================================================================================
 127
 128   TopoDS_Shape getShapeAtPoint( const SPoint3& point, const std::vector< ShapeBounds > & shapes )
 129   {
 130     TopoDS_Shape shape;
 131     float distmin = std::numeric_limits<float>::max();
 132     for ( size_t i = 0; i < shapes.size(); ++i )
 133     {
 134       float dist = GMSHPlugin_Mesher::DistBoundingBox( shapes[i]._bounds, point );
 135       if (dist < distmin)
 136       {
 137         shape = shapes[i]._shape;
 138         distmin = dist;
 139         if ( distmin == 0. )
 140           break;
 141       }
 142     }
 143     return shape;
 144   }
 145
 146   double segmentSize( const UVPtStructVec& nodeParam, size_t i )
 147   {
 148     double l1 = SMESH_NodeXYZ( nodeParam[i].node ).Distance( nodeParam[i-1].node );
 149     double l2 = SMESH_NodeXYZ( nodeParam[i].node ).Distance( nodeParam[i+1].node );
 150     return 0.5 * ( l1 + l2 );
 151   }
 152 }
 153
 154 //=============================================================================
 155 /*!
 156  *
 157  */
 158 //=============================================================================
 159
 160 GMSHPlugin_Mesher::GMSHPlugin_Mesher (SMESH_Mesh*         mesh,
 161                                       const TopoDS_Shape& aShape,
 162                                       bool                is2D,
 163                                       bool                is3D)
 164   : _mesh    (mesh),
 165     _shape   (aShape),
 166     _is2d    (is2D),
 167     _is3d    (is3D)
 168 {
 169   // il faudra peut être mettre un truc par defaut si l'utilisateur ne rentre rien en para
 170   //defaultParameters();
 171 }
 172
 173 //void GMSHPlugin_Mesher::defaultParameters(){}
 174
 175 void GMSHPlugin_Mesher::SetParameters(const GMSHPlugin_Hypothesis* hyp)
 176 {
 177   if (hyp != NULL)
 178   {
 179     _algo2d          = hyp->Get2DAlgo();
 180     _algo3d          = hyp->Get3DAlgo();
 181     _recomb2DAlgo    = hyp->GetRecomb2DAlgo();
 182     _recombineAll    = hyp->GetRecombineAll();
 183     _subdivAlgo      = hyp->GetSubdivAlgo();
 184     _remeshAlgo      = hyp->GetRemeshAlgo();
 185     _remeshPara      = hyp->GetRemeshPara();
 186     _smouthSteps     = hyp->GetSmouthSteps();
 187     _sizeFactor      = hyp->GetSizeFactor();
 188 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=10
 189     _meshCurvatureSize = hyp->GetMeshCurvatureSize();
 190 #endif
 191     _minSize         = hyp->GetMinSize();
 192     _maxSize         = hyp->GetMaxSize();
 193     _secondOrder     = hyp->GetSecondOrder();
 194     _useIncomplElem  = hyp->GetUseIncomplElem();
 195     _verbLvl         = hyp->GetVerbosityLevel();
 196     _compounds       = hyp->GetCompoundOnEntries();
 197     // 6 in the enum corresponds to 99 in gmsh
 198     if(_verbLvl == 6)
 199       _verbLvl = 99;
 200   }
 201   else
 202   {
 203     _algo2d          = 0;
 204     _algo3d          = 0;
 205     _recomb2DAlgo    = 0;
 206     _recombineAll    = false;
 207     _subdivAlgo      = 0;
 208     _remeshAlgo      = 0;
 209     _remeshPara      = 0;
 210     _smouthSteps     = 1;
 211     _sizeFactor      = 1;
 212 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=10
 213     _meshCurvatureSize = 0;
 214 #endif
 215     _minSize         = 0;
 216     _maxSize         = 1e22;
 217     _secondOrder     = false;
 218     _useIncomplElem  = true;
 219     _compounds.clear();
 220   }
 221 }
 222
 223
 224 //================================================================================
 225 /*!
 226  * \brief Set maximum number of threads to be used by Gmsh
 227  */
 228 //================================================================================
 229
 230 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 231 void GMSHPlugin_Mesher::SetMaxThreadsGmsh()
 232 {
 233   MESSAGE("GMSHPlugin_Mesher::SetMaxThreadsGmsh");
 234   // compound meshing (_compounds.size() > 0) and quad meshing (_algo2d >= 5) will
 235   // not be multi-threaded
 236   if (_compounds.size() > 0 || _algo2d >= 5){
 237     _maxThreads = 1;
 238   }
 239   else
 240     _maxThreads = omp_get_max_threads();
 241 }
 242 #endif
 243
 244 //================================================================================
 245 /*!
 246  * \brief Initialize GMSH model
 247  */
 248  //================================================================================
 249 void GMSHPlugin_Mesher::FillGMSHMesh()
 250 {
 251   gmsh::initialize();
 252   gmsh::model::add("mesh");
 253
 254   SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
 255
 256   int aNbOfNodes = 0;
 257   int aNbOfCurves = 0;
 258   int aNbOfLines = 0;
 259   std::vector<int> aTrinagle(3, 0);
 260
 261   const int invalid_ID = -1;
 262
 263   // typedef for maps
 264   typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 265   typedef TNodeToIDMap::value_type                     TN2ID;
 266   typedef map<std::pair<int, int>, int> TLineToIDMap;
 267   TNodeToIDMap aNodeToID;
 268   TLineToIDMap aLineToID;
 269
 270   TopAbs_ShapeEnum aMainType = _mesh->GetShapeToMesh().ShapeType();
 271   bool aCheckReverse = (aMainType == TopAbs_COMPOUND || aMainType == TopAbs_COMPSOLID);
 272
 273   SMESH_MesherHelper aHelper(*_mesh);
 274   SMESH_ProxyMesh::Ptr aProxyMesh(new SMESH_ProxyMesh(*_mesh));
 275   if (_mesh->NbQuadrangles() > 0)
 276   {
 277     StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
 278     Adaptor->Compute(*_mesh, _shape, aProxyMesh.get());
 279     aProxyMesh.reset(Adaptor);
 280   }
 281
 282   std::map<const SMDS_MeshElement*, bool, TIDCompare> listElements;
 283   for (TopExp_Explorer exFa(_shape, TopAbs_FACE); exFa.More(); exFa.Next())
 284   {
 285     const TopoDS_Shape& aShapeFace = exFa.Current();
 286     int faceID = meshDS->ShapeToIndex(aShapeFace);
 287
 288     bool isRev = false;
 289     if (aCheckReverse && aHelper.NbAncestors(aShapeFace, *_mesh, _shape.ShapeType()) > 1)
 290       // IsReversedSubMesh() can work wrong on strongly curved faces,
 291       // so we use it as less as possible
 292       isRev = aHelper.IsReversedSubMesh(TopoDS::Face(aShapeFace));
 293
 294     const SMESHDS_SubMesh* aSubMeshDSFace = aProxyMesh->GetSubMesh(aShapeFace);
 295     if (!aSubMeshDSFace) continue;
 296
 297     SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 298     if (aHelper.IsQuadraticSubMesh(_shape) &&
 299       dynamic_cast<const SMESH_ProxyMesh::SubMesh*>(aSubMeshDSFace))
 300     {
 301       // add medium nodes of proxy triangles to helper
 302       while (iteratorElem->more())
 303         aHelper.AddTLinks(static_cast<const SMDS_MeshFace*>(iteratorElem->next()));
 304
 305       iteratorElem = aSubMeshDSFace->GetElements();
 306     }
 307     while (iteratorElem->more()) // loop on elements on a geom face
 308     {
 309       // check mesh face
 310       const SMDS_MeshElement* elem = iteratorElem->next();
 311       if (!elem)
 312         return;
 313       if (elem->NbCornerNodes() != 3)
 314         return;
 315       listElements[elem] = isRev;
 316     }
 317   }
 318
 319   for (auto const& ls : listElements)
 320   {
 321     // Add nodes of triangles and triangles them-selves to netgen mesh
 322     // add three nodes of
 323     bool hasDegen = false;
 324     for (int iN = 0; iN < 3; ++iN)
 325     {
 326       const SMDS_MeshNode* aNode = ls.first->GetNode(iN);
 327       const int shapeID = aNode->getshapeId();
 328       if (aNode->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
 329         aHelper.IsDegenShape(shapeID))
 330       {
 331         // ignore all nodes on degeneraged edge and use node on its vertex instead
 332         TopoDS_Shape vertex = TopoDS_Iterator(meshDS->IndexToShape(shapeID)).Value();
 333         aNode = SMESH_Algo::VertexNode(TopoDS::Vertex(vertex), meshDS);
 334         hasDegen = true;
 335       }
 336       int& ngID = aNodeToID.insert(TN2ID(aNode, invalid_ID)).first->second;
 337       if (ngID == invalid_ID)
 338       {
 339         ngID = ++aNbOfNodes;
 340         gmsh::model::occ::addPoint(aNode->X(), aNode->Y(), aNode->Z(), 1.e-2, ngID);
 341       }
 342       aTrinagle[ls.second ? 2 - iN : iN] = ngID;
 343     }
 344     // add triangle
 345     if (hasDegen && (aTrinagle[0] == aTrinagle[1] ||
 346       aTrinagle[0] == aTrinagle[2] ||
 347       aTrinagle[2] == aTrinagle[1]))
 348         continue;
 349
 350
 351     std::vector<int> LinesID(3, 0);
 352     for (int anIndex = 0; anIndex < 3; ++anIndex)
 353     {
 354       int aNextIndex = (anIndex + 1) % 3;
 355       if (aLineToID.find({ aTrinagle[anIndex], aTrinagle[aNextIndex] }) == aLineToID.end()
 356         && aLineToID.find({ aTrinagle[aNextIndex], aTrinagle[anIndex] }) == aLineToID.end())
 357       {
 358         LinesID[anIndex] = aLineToID.insert({ { aTrinagle[aNextIndex], aTrinagle[anIndex] }, ++aNbOfLines }).first->second;
 359         gmsh::model::occ::addLine(aTrinagle[anIndex], aTrinagle[aNextIndex], LinesID[anIndex]);
 360       }
 361       else
 362       {
 363         LinesID[anIndex] = aLineToID.find({ aTrinagle[anIndex], aTrinagle[aNextIndex] })->second;
 364         if (LinesID[anIndex] == 0)
 365           LinesID[anIndex] = aLineToID.find({ aTrinagle[aNextIndex], aTrinagle[anIndex] })->second;
 366
 367       }
 368     }
 369     if (!aProxyMesh->IsTemporary(ls.first))
 370       swap(aTrinagle[1], aTrinagle[2]);
 371
 372     int aTag = gmsh::model::occ::addCurveLoop(LinesID);
 373   }
 374
 375   // Generate 1D and 2D mesh
 376   _gModel->mesh( /*dim=*/ 1);
 377   Set1DSubMeshes(_gModel);
 378   _gModel->mesh( /*dim=*/ 2);
 379 }
 380
 381 //================================================================================
 382 /*!
 383  * \brief Set Gmsh Options
 384  */
 385  //================================================================================
 386 void GMSHPlugin_Mesher::SetGmshOptions()
 387 {
 388   MESSAGE("GMSHPlugin_Mesher::SetGmshOptions");
 389   /*
 390   printf("We chose _algo2d         %d \n", _algo2d        );
 391   printf("We chose _algo3d         %d \n", _algo3d        );
 392   printf("We chose _recomb2DAlgo   %d \n", _recomb2DAlgo  );
 393   printf("We chose _recombineAll   %d \n", (_recombineAll)?1:0);
 394   printf("We chose _subdivAlgo     %d \n", _subdivAlgo    );
 395   printf("We chose _remeshAlgo     %d \n", _remeshAlgo    );
 396   printf("We chose _remeshPara     %d \n", _remeshPara    );
 397   printf("We chose _smouthSteps    %e \n", _smouthSteps   );
 398   printf("We chose _sizeFactor     %e \n", _sizeFactor    );
 399   printf("We chose _minSize        %e \n", _minSize       );
 400   printf("We chose _maxSize        %e \n", _maxSize       );
 401   printf("We chose _secondOrder    %d \n", (_secondOrder)?1:0);
 402   printf("We chose _useIncomplElem %d \n", (_useIncomplElem)?1:0);
 403   printf("We are in dimension      %d \n", (_is2d)?2:3);
 404   //*/
 405
 406   std::map <int,double> mapAlgo2d;
 407   mapAlgo2d[0]=2; // Automatic
 408   mapAlgo2d[1]=1; // MeshAdapt
 409   mapAlgo2d[2]=5; // Delaunay
 410   mapAlgo2d[3]=6; // Frontal-Delaunay
 411   mapAlgo2d[4]=8; // DelQuad (Frontal-Delaunay for Quads)
 412   mapAlgo2d[5]=9; // Packing of parallelograms
 413 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=10
 414   mapAlgo2d[6]=11;// Quasistructured quads with cross-fields
 415 #endif
 416
 417   std::map <int,double> mapAlgo3d;
 418   mapAlgo3d[0]=1; // Delaunay
 419   mapAlgo3d[1]=4; // Frontal
 420   mapAlgo3d[2]=7; // MMG3D
 421   mapAlgo3d[3]=9; // R-tree
 422   mapAlgo3d[4]=10;// HXT
 423
 424   int ok;
 425   ok = GmshSetOption("Mesh", "Algorithm"                , mapAlgo2d[_algo2d])    ;
 426   ASSERT(ok);
 427   if ( !_is2d)
 428   {
 429     ok = GmshSetOption("Mesh", "Algorithm3D"            , mapAlgo3d[_algo3d])    ;
 430     ASSERT(ok);
 431   }
 432   ok = GmshSetOption("Mesh", "RecombinationAlgorithm"   , (double)_recomb2DAlgo) ;
 433   ASSERT(ok);
 434   ok = GmshSetOption("Mesh", "RecombineAll"             , (_recombineAll)?1.:0.) ;
 435   ASSERT(ok);
 436   ok = GmshSetOption("Mesh", "SubdivisionAlgorithm"     , (double)_subdivAlgo)   ;
 437   ASSERT(ok);
 438   ok = GmshSetOption("Mesh", "RemeshAlgorithm"          , (double)_remeshAlgo)   ;
 439   //ASSERT(ok);
 440   ok = GmshSetOption("Mesh", "RemeshParametrization"    , (double)_remeshPara)   ;
 441   //ASSERT(ok);
 442   ok = GmshSetOption("Mesh", "Smoothing"                , (double)_smouthSteps)  ;
 443   //ASSERT(ok);
 444 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=10
 445   ok = GmshSetOption("Mesh", "MeshSizeFromCurvature"       , _meshCurvatureSize) ;
 446   ASSERT(ok);
 447 #endif
 448 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 449   ok = GmshSetOption("Mesh", "MeshSizeFactor"              , _sizeFactor)     ;
 450   ASSERT(ok);
 451   ok = GmshSetOption("Mesh", "MeshSizeMin"                 , _minSize)        ;
 452   ASSERT(ok);
 453   ok = GmshSetOption("Mesh", "MeshSizeMax"                 , _maxSize)        ;
 454   ASSERT(ok);
 455 #else
 456   ok = GmshSetOption("Mesh", "CharacteristicLengthFactor"  , _sizeFactor)     ;
 457   ASSERT(ok);
 458   ok = GmshSetOption("Mesh", "CharacteristicLengthMin"     , _minSize)        ;
 459   ASSERT(ok);
 460   ok = GmshSetOption("Mesh", "CharacteristicLengthMax"     , _maxSize)        ;
 461   ASSERT(ok);
 462 #endif
 463   ok = GmshSetOption("Mesh", "ElementOrder"             , (_secondOrder)?2.:1.)  ;
 464   ASSERT(ok);
 465   if (_secondOrder)
 466   {
 467     ok = GmshSetOption("Mesh", "SecondOrderIncomplete"  ,(_useIncomplElem)?1.:0.);
 468     ASSERT(ok);
 469   }
 470
 471   ok = GmshSetOption("General", "Verbosity"            , (double) _verbLvl )  ; // Verbosity level
 472   ASSERT(ok);
 473
 474 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 475 /*ok = GmshSetOption("Mesh", "MaxNumThreads1D"          , 0. )  ; // Coarse-grain algo threads
 476   ASSERT(ok);
 477   ok = GmshSetOption("Mesh", "MaxNumThreads2D"          , 0. )  ; // Coarse-grain algo threads
 478   ASSERT(ok);
 479   ok = GmshSetOption("Mesh", "MaxNumThreads3D"          , 0. )  ; // Fine-grain algo threads HXT
 480   ASSERT(ok);
 481 **/
 482   ok = GmshSetOption("General", "NumThreads"            , _maxThreads )  ; // system default i.e. OMP_NUM_THREADS
 483   ASSERT(ok);
 484 #ifdef WIN32
 485   if ( GMSHPlugin_Hypothesis::Algo3D::hxt == _algo3d ){
 486     MESSAGE("GMSHPlugin_Mesher::SetGmshOptions: HXT algorithm is being used. Setting number of threads to 1.");
 487     ok = GmshSetOption("Mesh", "MaxNumThreads3D"       , 1. );
 488     ASSERT(ok);
 489   } // hxt
 490 #endif
 491 #endif
 492 }
 493
 494 //================================================================================
 495 /*!
 496  * \brief Create and add Compounds into GModel _gModel.
 497  */
 498 //================================================================================
 499
 500 void GMSHPlugin_Mesher::CreateGmshCompounds()
 501 {
 502   MESSAGE("GMSHPlugin_Mesher::CreateGmshCompounds");
 503
 504   SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen();
 505
 506   OCC_Internals* occgeo = _gModel->getOCCInternals();
 507   bool toSynchronize = false;
 508
 509   for(std::set<std::string>::const_iterator its = _compounds.begin();its != _compounds.end(); ++its )
 510   {
 511     GEOM::GEOM_Object_var aGeomObj;
 512     TopoDS_Shape geomShape = TopoDS_Shape();
 513     SALOMEDS::SObject_var aSObj = SMESH_Gen_i::GetSMESHGen()->getStudyServant()->FindObjectID( (*its).c_str() );
 514     SALOMEDS::GenericAttribute_var anAttr;
 515     if (!aSObj->_is_nil() && aSObj->FindAttribute(anAttr, "AttributeIOR"))
 516     {
 517       SALOMEDS::AttributeIOR_var anIOR = SALOMEDS::AttributeIOR::_narrow(anAttr);
 518       CORBA::String_var aVal = anIOR->Value();
 519       CORBA::Object_var obj = SMESH_Gen_i::GetSMESHGen()->getStudyServant()->ConvertIORToObject(aVal);
 520       aGeomObj = GEOM::GEOM_Object::_narrow(obj);
 521     }
 522     geomShape = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
 523     if ( geomShape.IsNull() )
 524       continue;
 525
 526     TopAbs_ShapeEnum geomType = geomShape.ShapeType();
 527     if ( geomType == TopAbs_COMPOUND)// voir s'il ne faut pas mettre une erreur dans le cas contraire
 528     {
 529       MESSAGE("shapeType == TopAbs_COMPOUND");
 530       TopoDS_Iterator it(geomShape);
 531       if ( !it.More() )
 532         continue;
 533       TopAbs_ShapeEnum shapeType = it.Value().ShapeType();
 534       std::vector< std::pair< int, int > > dimTags;
 535       for ( ; it.More(); it.Next())
 536       {
 537         const TopoDS_Shape& topoShape = it.Value();
 538         ASSERT(topoShape.ShapeType() == shapeType);
 539         if ( _mesh->GetMeshDS()->ShapeToIndex( topoShape ) > 0 )
 540           occgeo->importShapes( &topoShape, false, dimTags );
 541         else
 542         {
 543           TopoDS_Shape face = TopExp_Explorer( _shape, shapeType ).Current();
 544           SMESH_subMesh* sm = _mesh->GetSubMesh( face );
 545           sm->GetComputeError() =
 546             SMESH_ComputeError::New
 547             ( COMPERR_WARNING, "Compound shape does not belong to the main geometry. Ignored");
 548         }
 549       }
 550       std::vector<int> tags;
 551       int dim = ( shapeType == TopAbs_EDGE ) ? 1 : 2;
 552       for ( size_t i = 0; i < dimTags.size(); ++i )
 553       {
 554         if ( dimTags[i].first == dim )
 555           tags.push_back( dimTags[i].second );
 556       }
 557       if ( !tags.empty() )
 558       {
 559         _gModel->getGEOInternals()->setCompoundMesh( dim, tags );
 560         toSynchronize = true;
 561       }
 562       if ( toSynchronize )
 563         _gModel->getGEOInternals()->synchronize(_gModel);
 564     }
 565   }
 566 }
 567
 568 //================================================================================
 569 /*!
 570  * \brief For a compound mesh set the mesh components to be transmitted to SMESH
 571  */
 572 //================================================================================
 573
 574 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 575 void GMSHPlugin_Mesher::SetCompoundMeshVisibility()
 576 {
 577
 578   // Loop over all faces, if the face belongs to a compound entry then
 579   // for all (boundary) edges whithin the face visibility is set to 0,
 580   // if the face doesn't belong to a compound entry then visibility is
 581   // set to 1 for all its (boundary) edges. Later, in FillSMesh() func
 582   // getVisibility() (returns either 1 or 0) is used to decide weather
 583   // the mesh of edge should be transmitted  to SMESH or not.
 584
 585   for ( GModel::fiter itF = _gModel->firstFace(); itF != _gModel->lastFace(); ++itF )
 586   {
 587     std::vector< GEdge *> faceEdges = (*itF)->edges();
 588
 589     for ( auto itE = faceEdges.begin(); itE != faceEdges.end(); ++itE )
 590     {
 591       if ( ((*itF)->compound.size()) )
 592         (*itE)->setVisibility(0);
 593       else
 594         (*itE)->setVisibility(1);
 595     }
 596   }
 597
 598
 599   // Loop over all edges, if the  edge belongs to a compound entry then
 600   // for all (boundary) vertices whithin the  edge visibility is set to
 601   // 0, if the edge doesn't belong to a  compound entry then visibility
 602   // is set to 1 for all its (boundary) vertices. Later, in FillSMesh()
 603   // func getVisibility() (returns either 1 or 0) is used to decide we-
 604   // ather the mesh of vertices should be transmitted  to SMESH or not.
 605
 606   for ( GModel::eiter itE = _gModel->firstEdge(); itE != _gModel->lastEdge(); ++itE )
 607   {
 608     std::vector<GVertex *> bndVerticies = (*itE)->vertices();
 609
 610     for( auto itV = bndVerticies.begin(); itV != bndVerticies.end(); ++itV )
 611     {
 612       if(((*itE)->compound.size()))
 613         (*itV)->setVisibility(0);
 614       else
 615         (*itV)->setVisibility(1);
 616     }
 617   }
 618
 619 }
 620 #endif
 621
 622 //================================================================================
 623 /*!
 624  * \brief Get a node by a GMSH mesh vertex
 625  */
 626 //================================================================================
 627
 628 const SMDS_MeshNode* GMSHPlugin_Mesher::Node( const MVertex* v )
 629 {
 630   std::map< const MVertex *, const SMDS_MeshNode* >::iterator v2n = _nodeMap.find( v );
 631   if ( v2n != _nodeMap.end() )
 632     return v2n->second;
 633
 634   return nullptr;
 635 }
 636
 637 //================================================================================
 638 /*!
 639  * \brief Return a corresponding sub-mesh if a shape is meshed
 640  */
 641 //================================================================================
 642
 643 SMESHDS_SubMesh* GMSHPlugin_Mesher::HasSubMesh( const TopoDS_Shape& s )
 644 {
 645   if ( SMESHDS_SubMesh*  sm = _mesh->GetMeshDS()->MeshElements( s ))
 646   {
 647     if ( s.ShapeType() == TopAbs_VERTEX )
 648       return ( sm->NbNodes() > 0 ) ? sm : nullptr;
 649     else
 650       return ( sm->NbElements() > 0 ) ? sm : nullptr;
 651   }
 652   return nullptr;
 653 }
 654
 655 //================================================================================
 656 /*!
 657  * \brief Write mesh from GModel instance to SMESH instance
 658  */
 659 //================================================================================
 660
 661 void GMSHPlugin_Mesher::FillSMesh()
 662 {
 663   SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
 664
 665   // ADD 0D ELEMENTS
 666   for ( GModel::viter it = _gModel->firstVertex(); it != _gModel->lastVertex(); ++it)
 667   {
 668     GVertex *gVertex = *it;
 669
 670     // GET topoVertex CORRESPONDING TO gVertex
 671     TopoDS_Vertex topoVertex = *((TopoDS_Vertex*)gVertex->getNativePtr());
 672
 673     if (gVertex->getVisibility() == 0) // belongs to a compound
 674     {
 675       SMESH_subMesh* sm = _mesh->GetSubMesh(topoVertex);
 676       sm->SetIsAlwaysComputed(true); // prevent from displaying errors
 677       continue;
 678     }
 679
 680     // FILL SMESH FOR topoVertex
 681     //nodes
 682     for( size_t i = 0; i < gVertex->mesh_vertices.size(); i++)
 683     {
 684       MVertex *v = gVertex->mesh_vertices[i];
 685       if(v->getIndex() >= 0)
 686       {
 687         if ( SMESHDS_SubMesh* sm = HasSubMesh( topoVertex ))
 688         {
 689           const SMDS_MeshNode *node = sm->GetNodes()->next();
 690           _nodeMap.insert({ v, node });
 691         }
 692         else
 693         {
 694           SMDS_MeshNode *node = meshDS->AddNode( v->x(),v->y(),v->z() );
 695           meshDS->SetNodeOnVertex( node, topoVertex );
 696           _nodeMap.insert({ v, node });
 697         }
 698       }
 699     }
 700     // WE DON'T ADD 0D ELEMENTS because it does not follow the salome meshers philosophy
 701     //elements
 702     // for(unsigned int i = 0; i < gVertex->getNumMeshElements(); i++)
 703     // {
 704     //   MElement *e = gVertex->getMeshElement(i);
 705     //   std::vector<MVertex*> verts;
 706     //   e->getVertices(verts);
 707     //   ASSERT(verts.size()==1);
 708     //   SMDS_Mesh0DElement* zeroDElement = 0;
 709     //   zeroDElement = meshDS->Add0DElementWithID(verts[0]->getNum(),e->getNum());
 710     //   meshDS->SetMeshElementOnShape(zeroDElement, topoVertex);
 711     // }
 712   }
 713
 714   // ADD 1D ELEMENTS
 715   for(GModel::eiter it = _gModel->firstEdge(); it != _gModel->lastEdge(); ++it)
 716   {
 717     GEdge *gEdge = *it;
 718
 719     // GET topoEdge CORRESPONDING TO gEdge
 720     TopoDS_Edge topoEdge;
 721     std::vector< ShapeBounds > topoEdges;
 722 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 723     if(gEdge->haveParametrization())
 724 #else
 725     if ( gEdge->geomType() != GEntity::CompoundCurve )
 726 #endif
 727     {
 728       topoEdge = *((TopoDS_Edge*)gEdge->getNativePtr());
 729       if (gEdge->getVisibility() == 0) // belongs to a compound
 730       {
 731         SMESH_subMesh* sm = _mesh->GetSubMesh(topoEdge);
 732         sm->SetIsAlwaysComputed(true); // prevent from displaying errors
 733         continue;
 734       }
 735       if ( HasSubMesh( topoEdge ))
 736         continue; // a meshed sub-mesh
 737     }
 738     bool isCompound = getBoundsOfShapes( gEdge, topoEdges );
 739
 740     // FILL SMESH FOR topoEdge
 741     //nodes
 742     for ( size_t i = 0; i < gEdge->mesh_vertices.size(); i++ )
 743     {
 744       MVertex *v = gEdge->mesh_vertices[i];
 745       if ( v->getIndex() >= 0 )
 746       {
 747         SMDS_MeshNode *node = meshDS->AddNode( v->x(),v->y(),v->z() );
 748
 749         if ( isCompound )
 750           topoEdge = TopoDS::Edge( getShapeAtPoint( v->point(), topoEdges ));
 751
 752         // Based on BLSURFPlugin_BLSURF
 753         gp_Pnt point3D( v->x(),v->y(),v->z() );
 754         Standard_Real p0 = 0.0;
 755         Standard_Real p1 = 1.0;
 756         TopLoc_Location loc;
 757         Handle(Geom_Curve) curve = BRep_Tool::Curve(topoEdge, loc, p0, p1);
 758
 759         if ( !curve.IsNull() )
 760         {
 761           if ( !loc.IsIdentity() )
 762             point3D.Transform( loc.Transformation().Inverted() );
 763
 764           GeomAPI_ProjectPointOnCurve proj(point3D, curve, p0, p1);
 765
 766           double pa = 0.;
 767           if ( proj.NbPoints() > 0 )
 768             pa = (double)proj.LowerDistanceParameter();
 769
 770           meshDS->SetNodeOnEdge( node, topoEdge, pa );
 771         }
 772         else
 773         {
 774           meshDS->SetNodeOnEdge( node, topoEdge );
 775         }
 776         //END on BLSURFPlugin_BLSURF
 777
 778
 779         _nodeMap.insert({ v, node });
 780       }
 781     }
 782   }
 783
 784   for ( GModel::eiter it = _gModel->firstEdge(); it != _gModel->lastEdge(); ++it )
 785   {
 786     GEdge *gEdge = *it;
 787     if ( gEdge->getVisibility() == 0) // belongs to a compound
 788       continue;
 789
 790     TopoDS_Edge topoEdge;
 791     std::vector< ShapeBounds > topoEdges;
 792     bool isCompound = getBoundsOfShapes( gEdge, topoEdges );
 793     if ( !isCompound )
 794       topoEdge = *((TopoDS_Edge*)gEdge->getNativePtr());
 795
 796     if ( HasSubMesh( topoEdge ))
 797       continue; // a meshed sub-mesh
 798
 799     //elements
 800     std::vector<MVertex*> verts(3);
 801     for ( size_t i = 0; i < gEdge->getNumMeshElements(); i++ )
 802     {
 803       MElement *e = gEdge->getMeshElement(i);
 804       verts.clear();
 805       e->getVertices(verts);
 806
 807       // if a node wasn't set, it is assigned here
 808       for ( size_t j = 0; j < verts.size(); j++ )
 809       {
 810         if ( verts[j]->onWhat()->getVisibility() == 0 )
 811         {
 812           SMDS_MeshNode *node = meshDS->AddNode(verts[j]->x(),verts[j]->y(),verts[j]->z() );
 813
 814           gp_Pnt point3D( verts[j]->x(),verts[j]->y(),verts[j]->z() );
 815           Standard_Real p0 = 0.0;
 816           Standard_Real p1 = 1.0;
 817           TopLoc_Location loc;
 818           Handle(Geom_Curve) curve = BRep_Tool::Curve(topoEdge, loc, p0, p1);
 819
 820           if ( !curve.IsNull() )
 821           {
 822             if ( !loc.IsIdentity() )
 823               point3D.Transform( loc.Transformation().Inverted() );
 824
 825             GeomAPI_ProjectPointOnCurve proj(point3D, curve, p0, p1);
 826
 827             double pa = 0.;
 828             if ( proj.NbPoints() > 0 )
 829               pa = (double)proj.LowerDistanceParameter();
 830
 831             meshDS->SetNodeOnEdge( node, topoEdge, pa );
 832           }
 833           else
 834           {
 835             meshDS->SetNodeOnEdge( node, topoEdge );
 836           }
 837
 838           verts[j]->setEntity(gEdge);
 839           _nodeMap.insert({ verts[j], node });
 840         }
 841       }
 842
 843       SMDS_MeshEdge* edge = 0;
 844       switch (verts.size())
 845       {
 846         case 2:
 847           edge = meshDS->AddEdge(Node( verts[0]),
 848                                  Node( verts[1]));
 849           break;
 850         case 3:
 851           edge = meshDS->AddEdge(Node( verts[0]),
 852                                  Node( verts[1]),
 853                                  Node( verts[2]));
 854           break;
 855         default:
 856           ASSERT(false);
 857           continue;
 858       }
 859       if ( isCompound )
 860         topoEdge = TopoDS::Edge( getShapeAtPoint( e->barycenter(), topoEdges ));
 861
 862       meshDS->SetMeshElementOnShape( edge, topoEdge );
 863     }
 864   }
 865
 866   // ADD 2D ELEMENTS
 867   for ( GModel::fiter it = _gModel->firstFace(); it != _gModel->lastFace(); ++it)
 868   {
 869     GFace *gFace = *it;
 870
 871 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 872     // Gmsh since version 4.3 is now producing extra surface and mesh when
 873     // compounds are involved. Since in Gmsh meshing procedure needs acess
 874     // to each of the original topology and the meshed topology. Hence  we
 875     // bypass the additional mesh in case of compounds. Note, similar cri-
 876     // teria also occurs in the following 'for' loop.
 877     if ( _compounds.size() && gFace->geomType() == GEntity::DiscreteSurface )
 878       continue;
 879 #endif
 880
 881     // GET topoFace CORRESPONDING TO gFace
 882     TopoDS_Face topoFace;
 883     std::vector< ShapeBounds > topoFaces;
 884
 885 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 886     if(gFace->haveParametrization())
 887 #else
 888     if ( gFace->geomType() != GEntity::CompoundSurface )
 889 #endif
 890     {
 891       topoFace = *((TopoDS_Face*)gFace->getNativePtr());
 892       if (gFace->getVisibility() == 0) // belongs to a compound
 893       {
 894         SMESH_subMesh* sm = _mesh->GetSubMesh(topoFace);
 895         sm->SetIsAlwaysComputed(true); // prevent from displaying errors
 896         continue;
 897       }
 898       if ( HasSubMesh( topoFace ))
 899         continue; // a meshed sub-mesh
 900     }
 901     bool isCompound = getBoundsOfShapes( gFace, topoFaces );
 902
 903     // FILL SMESH FOR topoFace
 904     //nodes
 905     for ( size_t i = 0; i < gFace->mesh_vertices.size(); i++ )
 906     {
 907       MVertex *v = gFace->mesh_vertices[i];
 908       if ( v->getIndex() >= 0 )
 909       {
 910         SMDS_MeshNode *node = meshDS->AddNode( v->x(),v->y(),v->z() );
 911
 912         if ( isCompound )
 913           topoFace = TopoDS::Face( getShapeAtPoint( v->point(), topoFaces ));
 914
 915         meshDS->SetNodeOnFace( node, topoFace );
 916         _nodeMap.insert({ v, node });
 917       }
 918     }
 919   }
 920
 921   for ( GModel::fiter it = _gModel->firstFace(); it != _gModel->lastFace(); ++it)
 922   {
 923     GFace *gFace = *it;
 924
 925 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 926     if ( _compounds.size() && gFace->geomType() == GEntity::DiscreteSurface )
 927       continue;
 928
 929     bool isCompound = (!gFace->haveParametrization());
 930 #else
 931     bool isCompound = ( gFace->geomType() == GEntity::CompoundSurface );
 932 #endif
 933
 934     if ( !isCompound && gFace->getVisibility() == 0 )
 935       continue;  // belongs to a compound
 936
 937     TopoDS_Face topoFace;
 938     std::vector< ShapeBounds > topoFaces;
 939     if ( isCompound )
 940       getBoundsOfShapes( gFace, topoFaces );
 941     else
 942       topoFace = *((TopoDS_Face*)gFace->getNativePtr());
 943
 944     if ( HasSubMesh( topoFace ))
 945       continue; // a meshed sub-mesh
 946
 947     //elements
 948     std::vector<MVertex*> verts;
 949     for ( size_t i = 0; i < gFace->getNumMeshElements(); i++ )
 950     {
 951       MElement *e = gFace->getMeshElement(i);
 952       verts.clear();
 953       e->getVertices(verts);
 954       SMDS_MeshFace* face = 0;
 955
 956       // if a node wasn't set, it is assigned here
 957       for ( size_t j = 0; j < verts.size(); j++)
 958       {
 959         if(verts[j]->onWhat()->getVisibility() == 0)
 960         {
 961           SMDS_MeshNode *node = meshDS->AddNode(verts[j]->x(),verts[j]->y(),verts[j]->z());
 962           meshDS->SetNodeOnFace( node, topoFace );
 963           _nodeMap.insert({ verts[j], node });
 964           verts[j]->setEntity(gFace);
 965         }
 966       }
 967       switch (verts.size())
 968       {
 969         case 3:
 970           face = meshDS->AddFace(Node( verts[0]),
 971                                  Node( verts[1]),
 972                                  Node( verts[2]));
 973           break;
 974         case 4:
 975           face = meshDS->AddFace(Node( verts[0]),
 976                                  Node( verts[1]),
 977                                  Node( verts[2]),
 978                                  Node( verts[3]));
 979           break;
 980         case 6:
 981           face = meshDS->AddFace(Node( verts[0]),
 982                                  Node( verts[1]),
 983                                  Node( verts[2]),
 984                                  Node( verts[3]),
 985                                  Node( verts[4]),
 986                                  Node( verts[5]));
 987           break;
 988         case 8:
 989           face = meshDS->AddFace(Node( verts[0]),
 990                                  Node( verts[1]),
 991                                  Node( verts[2]),
 992                                  Node( verts[3]),
 993                                  Node( verts[4]),
 994                                  Node( verts[5]),
 995                                  Node( verts[6]),
 996                                  Node( verts[7]));
 997           break;
 998         case 9:
 999           face = meshDS->AddFace(Node( verts[0]),
1000                                  Node( verts[1]),
1001                                  Node( verts[2]),
1002                                  Node( verts[3]),
1003                                  Node( verts[4]),
1004                                  Node( verts[5]),
1005                                  Node( verts[6]),
1006                                  Node( verts[7]),
1007                                  Node( verts[8]));
1008           break;
1009         default:
1010           ASSERT(false);
1011           continue;
1012       }
1013
1014       if ( isCompound )
1015         topoFace = TopoDS::Face( getShapeAtPoint( e->barycenter(), topoFaces ));
1016
1017       meshDS->SetMeshElementOnShape(face, topoFace);
1018     }
1019   }
1020
1021   // ADD 3D ELEMENTS
1022   for ( GModel::riter it = _gModel->firstRegion(); it != _gModel->lastRegion(); ++it)
1023   {
1024     GRegion *gRegion = *it;
1025     if (gRegion->getVisibility() == 0)
1026       continue;
1027
1028     // GET topoSolid CORRESPONDING TO gRegion
1029     TopoDS_Solid topoSolid = *((TopoDS_Solid*)gRegion->getNativePtr());
1030
1031     // FILL SMESH FOR topoSolid
1032
1033     //nodes
1034     for( size_t i = 0; i < gRegion->mesh_vertices.size(); i++)
1035     {
1036       MVertex *v = gRegion->mesh_vertices[i];
1037       if(v->getIndex() >= 0)
1038       {
1039         SMDS_MeshNode *node = meshDS->AddNode( v->x(),v->y(),v->z() );
1040         meshDS->SetNodeInVolume( node, topoSolid );
1041         _nodeMap.insert({ v, node });
1042       }
1043     }
1044
1045     //elements
1046     std::vector<MVertex*> verts;
1047     for( size_t i = 0; i < gRegion->getNumMeshElements(); i++)
1048     {
1049       MElement *e = gRegion->getMeshElement(i);
1050       verts.clear();
1051       e->getVertices(verts);
1052       SMDS_MeshVolume* volume = 0;
1053       switch (verts.size()){
1054       case 4:
1055         volume = meshDS->AddVolume(Node( verts[0]),
1056                                    Node( verts[2]),
1057                                    Node( verts[1]),
1058                                    Node( verts[3]));
1059         break;
1060       case 5:
1061         volume = meshDS->AddVolume(Node( verts[0]),
1062                                    Node( verts[3]),
1063                                    Node( verts[2]),
1064                                    Node( verts[1]),
1065                                    Node( verts[4]));
1066         break;
1067       case 6:
1068         volume = meshDS->AddVolume(Node( verts[0]),
1069                                    Node( verts[2]),
1070                                    Node( verts[1]),
1071                                    Node( verts[3]),
1072                                    Node( verts[5]),
1073                                    Node( verts[4]));
1074         break;
1075       case 8:
1076         volume = meshDS->AddVolume(Node( verts[0]),
1077                                    Node( verts[3]),
1078                                    Node( verts[2]),
1079                                    Node( verts[1]),
1080                                    Node( verts[4]),
1081                                    Node( verts[7]),
1082                                    Node( verts[6]),
1083                                    Node( verts[5]));
1084         break;
1085       case 10:
1086         volume = meshDS->AddVolume(Node( verts[0]),
1087                                    Node( verts[2]),
1088                                    Node( verts[1]),
1089                                    Node( verts[3]),
1090                                    Node( verts[6]),
1091                                    Node( verts[5]),
1092                                    Node( verts[4]),
1093                                    Node( verts[7]),
1094                                    Node( verts[8]),
1095                                    Node( verts[9]));
1096         break;
1097       case 13:
1098         volume = meshDS->AddVolume(Node( verts[0] ),
1099                                    Node( verts[3] ),
1100                                    Node( verts[2] ),
1101                                    Node( verts[1] ),
1102                                    Node( verts[4] ),
1103                                    Node( verts[6] ),
1104                                    Node( verts[10] ),
1105                                    Node( verts[8] ),
1106                                    Node( verts[5] ),
1107                                    Node( verts[7] ),
1108                                    Node( verts[12] ),
1109                                    Node( verts[11] ),
1110                                    Node( verts[9]));
1111         break;
1112       case 14: // same as case 13, because no pyra14 in smesh
1113         volume = meshDS->AddVolume(Node( verts[0] ),
1114                                    Node( verts[3] ),
1115                                    Node( verts[2] ),
1116                                    Node( verts[1] ),
1117                                    Node( verts[4] ),
1118                                    Node( verts[6] ),
1119                                    Node( verts[10] ),
1120                                    Node( verts[8] ),
1121                                    Node( verts[5] ),
1122                                    Node( verts[7] ),
1123                                    Node( verts[12] ),
1124                                    Node( verts[11] ),
1125                                    Node( verts[9]));
1126         break;
1127       case 15:
1128         volume = meshDS->AddVolume(Node( verts[0] ),
1129                                    Node( verts[2] ),
1130                                    Node( verts[1] ),
1131                                    Node( verts[3] ),
1132                                    Node( verts[5] ),
1133                                    Node( verts[4] ),
1134                                    Node( verts[7] ),
1135                                    Node( verts[9] ),
1136                                    Node( verts[6] ),
1137                                    Node( verts[13] ),
1138                                    Node( verts[14] ),
1139                                    Node( verts[12] ),
1140                                    Node( verts[8] ),
1141                                    Node( verts[11] ),
1142                                    Node( verts[10]));
1143         break;
1144       case 18: // same as case 15, because no penta18 in smesh
1145         volume = meshDS->AddVolume(Node( verts[0] ),
1146                                    Node( verts[2] ),
1147                                    Node( verts[1] ),
1148                                    Node( verts[3] ),
1149                                    Node( verts[5] ),
1150                                    Node( verts[4] ),
1151                                    Node( verts[7] ),
1152                                    Node( verts[9] ),
1153                                    Node( verts[6] ),
1154                                    Node( verts[13] ),
1155                                    Node( verts[14] ),
1156                                    Node( verts[12] ),
1157                                    Node( verts[8] ),
1158                                    Node( verts[11] ),
1159                                    Node( verts[10]));
1160         break;
1161       case 20:
1162         volume = meshDS->AddVolume(Node( verts[0] ),
1163                                    Node( verts[3] ),
1164                                    Node( verts[2] ),
1165                                    Node( verts[1] ),
1166                                    Node( verts[4] ),
1167                                    Node( verts[7] ),
1168                                    Node( verts[6] ),
1169                                    Node( verts[5] ),
1170                                    Node( verts[9] ),
1171                                    Node( verts[13] ),
1172                                    Node( verts[11] ),
1173                                    Node( verts[8] ),
1174                                    Node( verts[17] ),
1175                                    Node( verts[19] ),
1176                                    Node( verts[18] ),
1177                                    Node( verts[16] ),
1178                                    Node( verts[10] ),
1179                                    Node( verts[15] ),
1180                                    Node( verts[14] ),
1181                                    Node( verts[12]));
1182         break;
1183       case 27:
1184         volume = meshDS->AddVolume(Node( verts[0] ),
1185                                    Node( verts[3] ),
1186                                    Node( verts[2] ),
1187                                    Node( verts[1] ),
1188                                    Node( verts[4] ),
1189                                    Node( verts[7] ),
1190                                    Node( verts[6] ),
1191                                    Node( verts[5] ),
1192                                    Node( verts[9] ),
1193                                    Node( verts[13] ),
1194                                    Node( verts[11] ),
1195                                    Node( verts[8] ),
1196                                    Node( verts[17] ),
1197                                    Node( verts[19] ),
1198                                    Node( verts[18] ),
1199                                    Node( verts[16] ),
1200                                    Node( verts[10] ),
1201                                    Node( verts[15] ),
1202                                    Node( verts[14] ),
1203                                    Node( verts[12] ),
1204                                    Node( verts[20] ),
1205                                    Node( verts[22] ),
1206                                    Node( verts[24] ),
1207                                    Node( verts[23] ),
1208                                    Node( verts[21] ),
1209                                    Node( verts[25] ),
1210                                    Node( verts[26] ));
1211         break;
1212       default:
1213         ASSERT(false);
1214         continue;
1215       }
1216       meshDS->SetMeshElementOnShape(volume, topoSolid);
1217     }
1218   }
1219
1220   //return 0;
1221 }
1222
1223 //================================================================================
1224 /*!
1225  * \brief Find if SPoint point is in SBoundingBox3d bounds
1226  */
1227 //================================================================================
1228
1229 float GMSHPlugin_Mesher::DistBoundingBox(const SBoundingBox3d& bounds, const SPoint3& point)
1230 {
1231   SPoint3 min = bounds.min();
1232   SPoint3 max = bounds.max();
1233
1234   float x,y,z;
1235
1236   if (point.x() < min.x())
1237     x = min.x()-point.x();
1238   else if (point.x() > max.x())
1239     x = point.x()-max.x();
1240   else
1241     x = 0.;
1242
1243   if (point.y() < min.y())
1244     y = min.y()-point.y();
1245   else if (point.y() > max.y())
1246     y = point.y()-max.y();
1247   else
1248     y = 0.;
1249
1250   if (point.z() < min.z())
1251     z = min.z()-point.z();
1252   else if (point.z() > max.z())
1253     z = point.z()-max.z();
1254   else
1255     z = 0.;
1256
1257   return x*x+y*y+z*z;
1258 }
1259 //================================================================================
1260 /*!
1261  * \brief Reimplemented GmshMessage call. Actions done if errors occurs
1262  *        during gmsh meshing. We define here what to display and what to do.
1263  */
1264 //================================================================================
1265 void  GMSHPlugin_Mesher::mymsg::operator()(std::string level, std::string msg)
1266 {
1267   //MESSAGE("level="<< level.c_str() << ", msg=" << msg.c_str()<< "\n");
1268   printf("level=%s msg=%s\n", level.c_str(), msg.c_str());
1269
1270   if(level == "Fatal" || level == "Error")
1271   {
1272     std::ostringstream oss;
1273     if (level == "Fatal")
1274       oss << "Fatal error during Generation of Gmsh Mesh\n";
1275     else
1276       oss << "Error during Generation of Gmsh Mesh\n";
1277     oss << "  " << msg.c_str() << "\n";
1278     GEntity *e = _gModel->getCurrentMeshEntity();
1279     if(e)
1280     {
1281       oss << "  error occurred while meshing entity:\n" <<
1282              "    tag=" << e->tag() << "\n" <<
1283              "    dimension=" << e->dim() << "\n" <<
1284              "    native pointer=" << e->getNativePtr();
1285       //if(e->geomType() != GEntity::CompoundCurve and e->geomType() != GEntity::CompoundSurface)
1286       //{
1287         //SMESH_subMesh *sm = _mesh->GetSubMesh(*((TopoDS_Shape*)e->getNativePtr()));
1288         //SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1289         //SMESH_Comment comment;
1290         //comment << SMESH_Comment(oss.str);
1291         //std::string str = oss.str();
1292         //smError.reset( new SMESH_ComputeError( str ));
1293
1294         // plutot que de faire de la merde ici, on pourait simplement
1295         // remplir une liste pour dire sur quelles entités gmsh se plante
1296         // (puis faire le fillsmesh)
1297         // puis faire une nouvelle routine qui réécrit les messages d'erreur
1298         // probleme : gmsh peut planté en Fatal, dans ce cas pas de fillsmesh
1299       //}
1300     }
1301     if (level == "Fatal")
1302     {
1303         CTX::instance()->lock = 0;
1304         throw oss.str();
1305     }
1306     else
1307         printf("%s\n", oss.str().c_str());
1308   }
1309 }
1310
1311 //=============================================================================
1312 /*!
1313  * Here we are going to use the GMSH mesher
1314  */
1315 //=============================================================================
1316
1317 bool GMSHPlugin_Mesher::Compute()
1318 {
1319   MESSAGE("GMSHPlugin_Mesher::Compute");
1320
1321   int err = 0;
1322
1323 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1324   SetMaxThreadsGmsh();
1325 #endif
1326   //RNV: to avoid modification of PATH and PYTHONPATH
1327   char* argv[] = {"-noenv"};
1328   GmshInitialize(1,argv);
1329   SetGmshOptions();
1330   _gModel = new GModel();
1331   mymsg msg(_gModel);
1332   GmshSetMessageHandler(&msg);
1333   _gModel->importOCCShape((void*)&_shape);
1334   if (_compounds.size() > 0) CreateGmshCompounds();
1335   try
1336   {
1337
1338     HideComputedEntities(_gModel);
1339     if (_is3d)
1340     {
1341       FillGMSHMesh();
1342
1343       Set2DSubMeshes(_gModel);
1344
1345       _gModel->mesh( /*dim=*/ 3);
1346     }
1347     else
1348     {
1349       //CTX::instance()->mesh.maxNumThreads1D=1;
1350
1351       _gModel->mesh( /*dim=*/ 1);
1352
1353       Set1DSubMeshes(_gModel);
1354
1355       //_gModel->writeUNV("/tmp/1D.unv", 1,0,0);
1356       //CTX::instance()->mesh.maxNumThreads2D=1;
1357
1358       _gModel->mesh( /*dim=*/ 2);
1359
1360       if (!_is2d)
1361       {
1362         Set2DSubMeshes(_gModel);
1363
1364         //CTX::instance()->mesh.maxNumThreads3D=1;
1365
1366         _gModel->mesh( /*dim=*/ 3);
1367       }
1368       RestoreVisibility(_gModel);
1369     }
1370 #ifdef WITH_SMESH_CANCEL_COMPUTE
1371
1372 #endif
1373   }
1374   catch (std::string& str)
1375   {
1376     err = 1;
1377     std::cerr << "GMSH: exception caught: " << str << std::endl;
1378     MESSAGE(str);
1379   }
1380   catch (...)
1381   {
1382     err = 1;
1383     std::cerr << "GMSH: Unknown exception caught: " << std::endl;
1384     MESSAGE("Unrecoverable error during Generation of Gmsh Mesh");
1385   }
1386
1387   if (!err)
1388   {
1389 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1390     if (_compounds.size() > 0)
1391       SetCompoundMeshVisibility();
1392 #endif
1393     FillSMesh();
1394   }
1395   GmshSetMessageHandler(nullptr);
1396   delete _gModel;
1397   GmshFinalize();
1398   MESSAGE("GMSHPlugin_Mesher::Compute:End");
1399   return !err;
1400 }
1401
1402 //================================================================================
1403 /*!
1404  * \brief Set 1D sub-meshes to GModel. GMSH 1D mesh is made by now.
1405  *  \param [inout] _gModel - GMSH model
1406  */
1407 //================================================================================
1408
1409 void GMSHPlugin_Mesher::Set1DSubMeshes( GModel* gModel )
1410 {
1411   SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
1412
1413   for(GModel::eiter it = gModel->firstEdge(); it != gModel->lastEdge(); ++it)
1414   {
1415     GEdge *gEdge = *it;
1416
1417 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1418     if ( !gEdge->haveParametrization())
1419 #else
1420     if ( gEdge->geomType() == GEntity::CompoundCurve )
1421 #endif
1422       continue;
1423
1424     TopoDS_Edge topoEdge = *((TopoDS_Edge*)gEdge->getNativePtr());
1425     if ( !HasSubMesh( topoEdge ))
1426       continue; // empty sub-mesh
1427
1428     gEdge->deleteMesh();
1429
1430     // get node parameters on topoEdge
1431     StdMeshers_FaceSide side( TopoDS_Face(), topoEdge, _mesh, /*fwd=*/true, /*skpMedium=*/true);
1432     const UVPtStructVec& nodeParam = side.GetUVPtStruct();
1433     if ( nodeParam.empty() )
1434       throw std::string("Pb with StdMeshers_FaceSide::GetUVPtStruct()");
1435
1436     // get GMSH mesh vertices on VERTEX'es
1437     std::vector<MVertex *> mVertices( nodeParam.size(), nullptr );
1438     GVertex * gV0 = gEdge->getBeginVertex(), *gV1 = gEdge->getEndVertex();
1439     mVertices[0]     = gV0->mesh_vertices[ 0 ];
1440     mVertices.back() = gV1->mesh_vertices[ 0 ];
1441     TopoDS_Vertex v01 = *((TopoDS_Vertex*) gV0->getNativePtr());
1442     TopoDS_Shape  v02 = SMESH_MesherHelper::GetSubShapeByNode( nodeParam[0].node, meshDS );
1443     bool      reverse = !v01.IsSame( v02 );
1444     if ( mVertices[0] == mVertices.back() )
1445       reverse = ( nodeParam[0].param > nodeParam.back().param );
1446     const SMDS_MeshNode* n0 = reverse ? nodeParam.back().node : nodeParam[0].node;
1447     const SMDS_MeshNode* n1 = reverse ? nodeParam[0].node : nodeParam.back().node;
1448     _nodeMap.insert({ mVertices[ 0 ],   n0 });
1449     _nodeMap.insert({ mVertices.back(), n1 });
1450
1451     // create GMSH mesh vertices on gEdge
1452     for ( size_t i = 1; i < nodeParam.size() - 1; ++i )
1453     {
1454       size_t iN = reverse ? ( nodeParam.size() - 1 - i ) : i;
1455       SMESH_NodeXYZ xyz = nodeParam[ iN ].node;
1456       double lc = segmentSize( nodeParam, iN );
1457       // SVector3 der = gEdge->firstDer(nodeParam[ iN ].param);
1458       // double lc = norm(der) / segmentSize( nodeParam, i );
1459
1460       mVertices[ i ] = new MEdgeVertex( xyz.X(), xyz.Y(), xyz.Z(),
1461                                         gEdge, nodeParam[ iN ].param, 0, lc);
1462       gEdge->mesh_vertices.push_back( mVertices[ i ]);
1463       _nodeMap.insert({ mVertices[ i ], nodeParam[ iN ].node });
1464     }
1465     // create GMSH mesh edges
1466     for ( size_t i = 1; i < mVertices.size(); ++i )
1467     {
1468       gEdge->lines.push_back( new MLine( mVertices[ i - 1 ],
1469                                          mVertices[ i ]));
1470     }
1471     /*{
1472       cout << endl << "EDGE " << gEdge->tag() <<
1473         ( topoEdge.Orientation() == TopAbs_FORWARD ? " F" : " R") << endl;
1474       MVertex* mv = gV0->mesh_vertices[ 0 ];
1475       cout << "V0: " << mv->x() << ", " << mv->y() << ", " << mv->z() << ", "<<endl;
1476       for ( size_t i = 0; i < gEdge->mesh_vertices.size(); ++i )
1477       {
1478         MEdgeVertex* mv = (MEdgeVertex*) gEdge->mesh_vertices[i];
1479         cout << i << ": " << mv->x() << ", " << mv->y() << ", " << mv->z() << ", ";
1480         double t;
1481         mv->getParameter(0, t );
1482         cout << ":\t t = "<< t << " lc = " << mv->getLc() << endl;
1483       }
1484       mv = gV1->mesh_vertices[ 0 ];
1485       cout << "V1: " << mv->x() << ", " << mv->y() << ", " << mv->z() << ", "<<endl;
1486       }*/
1487   }
1488   return;
1489 }
1490
1491 //================================================================================
1492 /*!
1493  * \brief Set 2D sub-meshes to GModel. GMSH 2D mesh is made by now.
1494  *  \param [inout] _gModel - GMSH model
1495  */
1496 //================================================================================
1497
1498 void GMSHPlugin_Mesher::Set2DSubMeshes( GModel* gModel )
1499 {
1500   if ( _nodeMap.empty() )
1501     return; // no sub-meshes
1502
1503   SMESH_MesherHelper helper( *_mesh );
1504
1505   std::map< const SMDS_MeshNode* , const MVertex * > nodes2mvertMap;
1506   for ( auto & v2n : _nodeMap )
1507     nodes2mvertMap.insert({ v2n.second, v2n.first });
1508
1509   std::vector<MVertex *> mVertices;
1510
1511   for(GModel::fiter it = gModel->firstFace(); it != gModel->lastFace(); ++it)
1512   {
1513     GFace *gFace = *it;
1514
1515 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1516     if ( !gFace->haveParametrization())
1517 #else
1518       if ( gFace->geomType() == GEntity::CompoundSurface )
1519 #endif
1520         continue;
1521
1522     TopoDS_Face topoFace = *((TopoDS_Face*)gFace->getNativePtr());
1523     SMESHDS_SubMesh*  sm = HasSubMesh( topoFace );
1524     if ( !sm )
1525       continue;
1526     //_gModel->writeUNV("/tmp/befDEl.unv", 1,0,0);
1527
1528     gFace->deleteMesh();
1529
1530     bool reverse = false;
1531     if ( gFace->getRegion(0) )
1532     {
1533       //GRegion * gRegion = gFace->getRegion(0);
1534       TopoDS_Shape topoSolid = *((TopoDS_Shape*)gFace->getNativePtr());
1535       TopAbs_Orientation faceOriInSolid = helper.GetSubShapeOri( topoSolid, topoFace );
1536       if ( faceOriInSolid >= 0 )
1537         reverse =
1538           helper.IsReversedSubMesh( TopoDS::Face( topoFace.Oriented( faceOriInSolid )));
1539     }
1540
1541     for ( SMDS_ElemIteratorPtr fIt = sm->GetElements(); fIt->more(); )
1542     {
1543       const SMDS_MeshElement* f = fIt->next();
1544
1545       int nbN = f->NbCornerNodes();
1546       if ( nbN > 4 )
1547         throw std::string("Polygon sub-meshes not supported");
1548
1549       mVertices.resize( nbN );
1550       for ( int i = 0; i < nbN; ++i )
1551       {
1552         const SMDS_MeshNode* n = f->GetNode( i );
1553         MVertex *           mv = nullptr;
1554         auto n2v = nodes2mvertMap.find( n );
1555         if ( n2v != nodes2mvertMap.end() )
1556         {
1557           mv = const_cast< MVertex*>( n2v->second );
1558         }
1559         else
1560         {
1561           if ( n->GetPosition()->GetDim() < 2 )
1562             throw std::string("Wrong mapping of edge nodes to GMSH nodes");
1563           SMESH_NodeXYZ xyz = n;
1564           bool ok = true;
1565           gp_XY uv = helper.GetNodeUV( topoFace, n, nullptr, &ok );
1566           mv = new MFaceVertex( xyz.X(), xyz.Y(), xyz.Z(), gFace, uv.X(), uv.Y() );
1567           gFace->mesh_vertices.push_back( mv );
1568           nodes2mvertMap.insert({ n, mv });
1569           _nodeMap.insert      ({ mv, n });
1570         }
1571         mVertices[ i ] = mv;
1572       }
1573       // create GMSH mesh faces
1574       switch ( nbN ) {
1575       case 3:
1576         if ( reverse )
1577           gFace->triangles.push_back (new MTriangle(mVertices[0], mVertices[2], mVertices[1]));
1578         else
1579           gFace->triangles.push_back (new MTriangle(mVertices[0], mVertices[1], mVertices[2]));
1580         break;
1581       case 4:
1582         if ( reverse )
1583           gFace->quadrangles.push_back (new MQuadrangle(mVertices[0], mVertices[3],
1584                                                         mVertices[2], mVertices[1]));
1585         else
1586           gFace->quadrangles.push_back (new MQuadrangle(mVertices[0], mVertices[1],
1587                                                         mVertices[2], mVertices[3]));
1588         break;
1589       default:;
1590       }
1591     }
1592   } // loop on GMSH faces
1593
1594   return;
1595 }
1596
1597 //================================================================================
1598 /*!
1599  * \brief Set visibility 0 to already computed geom entities
1600  *        to prevent their meshing
1601  */
1602 //================================================================================
1603
1604 void GMSHPlugin_Mesher::HideComputedEntities( GModel* gModel )
1605 {
1606   CTX::instance()->mesh.meshOnlyVisible = true;
1607
1608   for(GModel::eiter it = gModel->firstEdge(); it != gModel->lastEdge(); ++it)
1609   {
1610     GEdge *gEdge = *it;
1611
1612 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1613     if ( !gEdge->haveParametrization())
1614 #else
1615       if ( gEdge->geomType() == GEntity::CompoundCurve )
1616 #endif
1617         continue;
1618
1619     TopoDS_Edge topoEdge = *((TopoDS_Edge*)gEdge->getNativePtr());
1620     if ( HasSubMesh( topoEdge ))
1621       gEdge->setVisibility(0);
1622   }
1623
1624
1625   for(GModel::fiter it = gModel->firstFace(); it != gModel->lastFace(); ++it)
1626   {
1627     GFace *gFace = *it;
1628
1629 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1630     if ( !gFace->haveParametrization())
1631 #else
1632       if ( gFace->geomType() == GEntity::CompoundSurface )
1633 #endif
1634         continue;
1635
1636     TopoDS_Face topoFace = *((TopoDS_Face*)gFace->getNativePtr());
1637     if ( HasSubMesh( topoFace ))
1638       gFace->setVisibility(0);
1639   }
1640 }
1641
1642 //================================================================================
1643 /*!
1644  * \brief Restore visibility of all geom entities
1645  */
1646 //================================================================================
1647
1648 void GMSHPlugin_Mesher::RestoreVisibility( GModel* gModel )
1649 {
1650   for(GModel::eiter it = gModel->firstEdge(); it != gModel->lastEdge(); ++it)
1651   {
1652     GEdge *gEdge = *it;
1653     gEdge->setVisibility(1);
1654   }
1655   for(GModel::fiter it = gModel->firstFace(); it != gModel->lastFace(); ++it)
1656   {
1657     GFace *gFace = *it;
1658     gFace->setVisibility(1);
1659   }
1660 }
1661
1662 /*
1663   void GMSHPlugin_Mesher::toPython( GModel* )
1664   {
1665   const char*  pyFile = "/tmp/gMesh.py";
1666   ofstream outfile( pyFile, ios::out );
1667   if ( !outfile ) return;
1668
1669   outfile << "import salome, SMESH" << std::endl
1670           << "from salome.smesh import smeshBuilder" << std::endl
1671           << "smesh = smeshBuilder.New()" << std::endl
1672           << "mesh = smesh.Mesh()" << std::endl << std::endl;
1673
1674   outfile << "## VERTICES" << endl;
1675   for ( GModel::viter it = _gModel->firstVertex(); it != _gModel->lastVertex(); ++it)
1676   {
1677     GVertex *gVertex = *it;
1678
1679     for(unsigned int i = 0; i < gVertex->mesh_vertices.size(); i++)
1680     {
1681       MVertex *v = gVertex->mesh_vertices[i];
1682       if ( v->getIndex() >= 0)
1683       {
1684         outfile << "n" << v->getNum() << " = mesh.AddNode("
1685                 << v->x() << ", " << v->y() << ", " << v->z()<< ")"
1686                 << " ## tag = " << gVertex->tag() << endl;
1687       }
1688     }
1689   }
1690
1691   for(GModel::eiter it = _gModel->firstEdge(); it != _gModel->lastEdge(); ++it)
1692   {
1693     GEdge *gEdge = *it;
1694     outfile << "## GEdge " << gEdge->tag() << endl;
1695
1696 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1697     if(gEdge->haveParametrization())
1698 #else
1699       if ( gEdge->geomType() != GEntity::CompoundCurve )
1700 #endif
1701         for ( size_t i = 0; i < gEdge->mesh_vertices.size(); i++ )
1702         {
1703           MVertex *v = gEdge->mesh_vertices[i];
1704           if ( v->getIndex() >= 0 )
1705           {
1706             outfile << "n" << v->getNum() << " = mesh.AddNode("
1707                     << v->x() << ", " << v->y() << ", " << v->z()<< ")"
1708                     << " ## tag = " << gEdge->tag() << endl;
1709           }
1710         }
1711   }
1712
1713   for ( GModel::fiter it = _gModel->firstFace(); it != _gModel->lastFace(); ++it)
1714   {
1715     GFace *gFace = *it;
1716     if ( _compounds.size() && gFace->geomType() == GEntity::DiscreteSurface )
1717       continue;
1718     outfile << "## GFace " << gFace->tag() << endl;
1719
1720     for ( size_t i = 0; i < gFace->mesh_vertices.size(); i++ )
1721     {
1722       MVertex *v = gFace->mesh_vertices[i];
1723       if ( v->getIndex() >= 0 )
1724       {
1725         outfile << "n" << v->getNum() << " = mesh.AddNode("
1726                 << v->x() << ", " << v->y() << ", " << v->z()<< ")"
1727                 << " ## tag = " << gFace->tag() << endl;
1728       }
1729     }
1730   }
1731
1732   std::vector<MVertex*> verts(3);
1733   for(GModel::eiter it = _gModel->firstEdge(); it != _gModel->lastEdge(); ++it)
1734   {
1735     GEdge *gEdge = *it;
1736     outfile << "## GEdge " << gEdge->tag() << endl;
1737
1738 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1739     if(gEdge->haveParametrization())
1740 #else
1741       if ( gEdge->geomType() != GEntity::CompoundCurve )
1742 #endif
1743     for ( size_t i = 0; i < gEdge->getNumMeshElements(); i++ )
1744     {
1745       MElement *e = gEdge->getMeshElement(i);
1746       verts.clear();
1747       e->getVertices(verts);
1748
1749       outfile << "e" << e->getNum() << " = mesh.AddEdge(["
1750               << "n" << verts[0]->getNum() << ","
1751               << "n" << verts[1]->getNum();
1752       if ( verts.size() == 3 )
1753         outfile << "n" << verts[2]->getNum();
1754       outfile << "])"<< endl;
1755     }
1756   }
1757
1758   for ( GModel::fiter it = _gModel->firstFace(); it != _gModel->lastFace(); ++it)
1759   {
1760     GFace *gFace = *it;
1761     if ( _compounds.size() && gFace->geomType() == GEntity::DiscreteSurface )
1762       continue;
1763     outfile << "## GFace " << gFace->tag() << endl;
1764
1765     for ( size_t i = 0; i < gFace->getNumMeshElements(); i++ )
1766     {
1767       MElement *e = gFace->getMeshElement(i);
1768       verts.clear();
1769       e->getVertices(verts);
1770
1771       outfile << "f" << e->getNum() << " = mesh.AddFace([";
1772       for ( size_t j = 0; j < verts.size(); j++)
1773       {
1774         outfile << "n" << verts[j]->getNum();
1775         if ( j < verts.size()-1 )
1776           outfile << ", ";
1777       }
1778       outfile << "])" << endl;
1779     }
1780   }
1781   std::cout << "Write " << pyFile << std::endl;
1782 }
1783 */