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