src/GMSHPlugin/GMSHPlugin_Mesher.cxx

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