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