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