src/GMSHPlugin/GMSHPlugin_Mesher.cxx

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