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_FaceOfNodes.hxx>
  25 #include <SMDS_MeshElement.hxx>
  26 #include <SMDS_MeshNode.hxx>
  27 #include <SMESHDS_Mesh.hxx>
  28 #include <SMESH_Block.hxx>
  29 #include <SMESH_Comment.hxx>
  30 #include <SMESH_ComputeError.hxx>
  31 #include <SMESH_File.hxx>
  32 #include <SMESH_Gen_i.hxx>
  33 #include <SMESH_Mesh.hxx>
  34 #include <SMESH_MesherHelper.hxx>
  35 #include <SMESH_subMesh.hxx>
  36 #include <utilities.h>
  37
  38 #include <vector>
  39 #include <limits>
  40
  41 #include <BRep_Tool.hxx>
  42 #include <Bnd_B3d.hxx>
  43 #include <GCPnts_AbscissaPoint.hxx>
  44 #include <GeomAdaptor_Curve.hxx>
  45 #include <NCollection_Map.hxx>
  46 #include <OSD_File.hxx>
  47 #include <OSD_Path.hxx>
  48 #include <Standard_ErrorHandler.hxx>
  49 #include <Standard_ProgramError.hxx>
  50 #include <TCollection_AsciiString.hxx>
  51 #include <TopExp.hxx>
  52 #include <TopExp_Explorer.hxx>
  53 #include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
  54 #include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
  55 #include <TopTools_DataMapOfShapeInteger.hxx>
  56 #include <TopTools_DataMapOfShapeShape.hxx>
  57 #include <TopTools_ListIteratorOfListOfShape.hxx>
  58 #include <TopTools_MapOfShape.hxx>
  59 #include <TopoDS.hxx>
  60
  61 #if GMSH_MAJOR_VERSION >=4
  62 #include <GmshGlobal.h>
  63 #include <gmsh/Context.h>
  64 #endif
  65
  66 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
  67 #include <omp.h>
  68 #endif
  69
  70 using namespace std;
  71
  72 namespace
  73 {
  74   struct ShapeBounds
  75   {
  76     SBoundingBox3d _bounds;
  77     TopoDS_Shape   _shape;
  78   };
  79
  80   //================================================================================
  81   /*!
  82    * \brief Retrieve ShapeBounds from a compound GEdge
  83    */
  84   //================================================================================
  85
  86   bool getBoundsOfShapes( GEdge*                       gEdge,
  87                           std::vector< ShapeBounds > & topoEdges )
  88   {
  89     topoEdges.clear();
  90 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
  91     for ( size_t i = 0; i < gEdge->compound.size(); ++i )
  92     {
  93       GEdge* gE = static_cast< GEdge* >( gEdge->compound[ i ]);
  94       topoEdges.push_back( ShapeBounds{ gE->bounds(), *((TopoDS_Edge*)gE->getNativePtr()) });
  95     }
  96 #else
  97     for ( size_t i = 0; i < gEdge->_compound.size(); ++i )
  98     {
  99       GEdge* gE = static_cast< GEdge* >( gEdge->_compound[ i ]);
 100       topoEdges.push_back( ShapeBounds{ gE->bounds(), *((TopoDS_Edge*)gE->getNativePtr()) });
 101     }
 102 #endif
 103     return topoEdges.size();
 104   }
 105
 106   //================================================================================
 107   /*!
 108    * \brief Retrieve ShapeBounds from a compound GFace
 109    */
 110   //================================================================================
 111
 112   bool getBoundsOfShapes( GFace*                       gFace,
 113                           std::vector< ShapeBounds > & topoFaces )
 114   {
 115     topoFaces.clear();
 116 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 117     for ( size_t i = 0; i < gFace->compound.size(); ++i )
 118     {
 119       GFace* gF = static_cast< GFace* >( gFace->compound[ i ]);
 120       topoFaces.push_back( ShapeBounds{ gF->bounds(), *((TopoDS_Face*)gF->getNativePtr()) });
 121     }
 122 #else
 123     for ( size_t i = 0; i < gFace->_compound.size(); ++i )
 124     {
 125       GFace* gF = static_cast< GFace* >( gFace->_compound[ i ]);
 126       topoFaces.push_back( ShapeBounds{ gF->bounds(), *((TopoDS_Face*)gF->getNativePtr()) });
 127     }
 128 #endif
 129     return topoFaces.size();
 130   }
 131   //================================================================================
 132   /*!
 133    * \brief Find a shape whose bounding box includes a given point
 134    */
 135   //================================================================================
 136
 137   TopoDS_Shape getShapeAtPoint( const SPoint3& point, const std::vector< ShapeBounds > & shapes )
 138   {
 139     TopoDS_Shape shape;
 140     float distmin = std::numeric_limits<float>::max();
 141     for ( size_t i = 0; i < shapes.size(); ++i )
 142     {
 143       float dist = GMSHPlugin_Mesher::DistBoundingBox( shapes[i]._bounds, point );
 144       if (dist < distmin)
 145       {
 146         shape = shapes[i]._shape;
 147         distmin = dist;
 148         if ( distmin == 0. )
 149           break;
 150       }
 151     }
 152     return shape;
 153   }
 154 }
 155
 156 //=============================================================================
 157 /*!
 158  *
 159  */
 160 //=============================================================================
 161
 162 GMSHPlugin_Mesher::GMSHPlugin_Mesher (SMESH_Mesh* mesh,
 163                                       const TopoDS_Shape& aShape)
 164   : _mesh    (mesh),
 165     _shape   (aShape)
 166 {
 167   // il faudra peut être mettre un truc par defaut si l'utilisateur ne rentre rien en para
 168   //defaultParameters();
 169 }
 170
 171 //void GMSHPlugin_Mesher::defaultParameters(){}
 172
 173 void GMSHPlugin_Mesher::SetParameters(const GMSHPlugin_Hypothesis* hyp)
 174 {
 175   if (hyp != NULL)
 176   {
 177     _algo2d          = hyp->Get2DAlgo();
 178     _algo3d          = hyp->Get3DAlgo();
 179     _recomb2DAlgo    = hyp->GetRecomb2DAlgo();
 180     _recombineAll    = hyp->GetRecombineAll();
 181     _subdivAlgo      = hyp->GetSubdivAlgo();
 182     _remeshAlgo      = hyp->GetRemeshAlgo();
 183     _remeshPara      = hyp->GetRemeshPara();
 184     _smouthSteps     = hyp->GetSmouthSteps();
 185     _sizeFactor      = hyp->GetSizeFactor();
 186     _minSize         = hyp->GetMinSize();
 187     _maxSize         = hyp->GetMaxSize();
 188     _secondOrder     = hyp->GetSecondOrder();
 189     _useIncomplElem  = hyp->GetUseIncomplElem();
 190     _is2d            = hyp->GetIs2d();
 191     _compounds       = hyp->GetCompoundOnEntries();
 192   }
 193   else
 194   {
 195     _algo2d          = 0;
 196     _algo3d          = 0;
 197     _recomb2DAlgo    = 0;
 198     _recombineAll    = false;
 199     _subdivAlgo      = 0;
 200     _remeshAlgo      = 0;
 201     _remeshPara      = 0;
 202     _smouthSteps     = 1;
 203     _sizeFactor      = 1;
 204     _minSize         = 0;
 205     _maxSize         = 1e22;
 206     _secondOrder     = false;
 207     _useIncomplElem  = true;
 208     _is2d            = false;
 209     _compounds.clear();
 210   }
 211 }
 212
 213
 214 //================================================================================
 215 /*!
 216  * \brief Set maximum number of threads to be used by Gmsh
 217  */
 218 //================================================================================
 219
 220 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 221 void GMSHPlugin_Mesher::SetMaxThreadsGmsh()
 222 {
 223   MESSAGE("GMSHPlugin_Mesher::SetMaxThreadsGmsh");
 224   if (_compounds.size() > 0)
 225     _maxThreads = 1;
 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
 265   std::map <int,double> mapAlgo3d;
 266   mapAlgo3d[0]=1; // Delaunay
 267   mapAlgo3d[1]=4; // Frontal
 268   mapAlgo3d[2]=7; // MMG3D
 269   mapAlgo3d[3]=9; // R-tree
 270   mapAlgo3d[4]=10;// HXT
 271
 272   int ok;
 273   ok = GmshSetOption("Mesh", "Algorithm"                , mapAlgo2d[_algo2d])    ;
 274   ASSERT(ok);
 275   if ( !_is2d)
 276   {
 277     ok = GmshSetOption("Mesh", "Algorithm3D"            , mapAlgo3d[_algo3d])    ;
 278     ASSERT(ok);
 279   }
 280   ok = GmshSetOption("Mesh", "RecombinationAlgorithm"   , (double)_recomb2DAlgo) ;
 281   ASSERT(ok);
 282   ok = GmshSetOption("Mesh", "RecombineAll"             , (_recombineAll)?1.:0.) ;
 283   ASSERT(ok);
 284   ok = GmshSetOption("Mesh", "SubdivisionAlgorithm"     , (double)_subdivAlgo)   ;
 285   ASSERT(ok);
 286   ok = GmshSetOption("Mesh", "RemeshAlgorithm"          , (double)_remeshAlgo)   ;
 287   //ASSERT(ok);
 288   ok = GmshSetOption("Mesh", "RemeshParametrization"    , (double)_remeshPara)   ;
 289   //ASSERT(ok);
 290   ok = GmshSetOption("Mesh", "Smoothing"                , (double)_smouthSteps)  ;
 291   //ASSERT(ok);
 292 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 293   ok = GmshSetOption("Mesh", "MeshSizeFactor"              , _sizeFactor)     ;
 294   ASSERT(ok);
 295   ok = GmshSetOption("Mesh", "MeshSizeMin"                 , _minSize)        ;
 296   ASSERT(ok);
 297   ok = GmshSetOption("Mesh", "MeshSizeMax"                 , _maxSize)        ;
 298   ASSERT(ok);
 299 #else
 300   ok = GmshSetOption("Mesh", "CharacteristicLengthFactor"  , _sizeFactor)     ;
 301   ASSERT(ok);
 302   ok = GmshSetOption("Mesh", "CharacteristicLengthMin"     , _minSize)        ;
 303   ASSERT(ok);
 304   ok = GmshSetOption("Mesh", "CharacteristicLengthMax"     , _maxSize)        ;
 305   ASSERT(ok);
 306 #endif
 307   ok = GmshSetOption("Mesh", "ElementOrder"             , (_secondOrder)?2.:1.)  ;
 308   ASSERT(ok);
 309   if (_secondOrder)
 310   {
 311     ok = GmshSetOption("Mesh", "SecondOrderIncomplete"  ,(_useIncomplElem)?1.:0.);
 312     ASSERT(ok);
 313   }
 314
 315 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 316 /*ok = GmshSetOption("Mesh", "MaxNumThreads1D"          , 0. )  ; // Coarse-grain algo threads
 317   ASSERT(ok);
 318   ok = GmshSetOption("Mesh", "MaxNumThreads2D"          , 0. )  ; // Coarse-grain algo threads
 319   ASSERT(ok);
 320   ok = GmshSetOption("Mesh", "MaxNumThreads3D"          , 0. )  ; // Fine-grain algo threads HXT
 321   ASSERT(ok);
 322 **/
 323   ok = GmshSetOption("General", "NumThreads"            , _maxThreads )  ; // system default i.e. OMP_NUM_THREADS
 324   ASSERT(ok);
 325 #endif
 326 }
 327
 328 //================================================================================
 329 /*!
 330  * \brief Create and add Compounds into GModel _gModel.
 331  */
 332 //================================================================================
 333
 334 void GMSHPlugin_Mesher::CreateGmshCompounds()
 335 {
 336   MESSAGE("GMSHPlugin_Mesher::CreateGmshCompounds");
 337
 338   SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen();
 339
 340   OCC_Internals* occgeo = _gModel->getOCCInternals();
 341   bool toSynchronize = false;
 342
 343   for(std::set<std::string>::const_iterator its = _compounds.begin();its != _compounds.end(); ++its )
 344   {
 345     GEOM::GEOM_Object_var aGeomObj;
 346     TopoDS_Shape geomShape = TopoDS_Shape();
 347     SALOMEDS::SObject_var aSObj = SMESH_Gen_i::GetSMESHGen()->getStudyServant()->FindObjectID( (*its).c_str() );
 348     SALOMEDS::GenericAttribute_var anAttr;
 349     if (!aSObj->_is_nil() && aSObj->FindAttribute(anAttr, "AttributeIOR"))
 350     {
 351       SALOMEDS::AttributeIOR_var anIOR = SALOMEDS::AttributeIOR::_narrow(anAttr);
 352       CORBA::String_var aVal = anIOR->Value();
 353       CORBA::Object_var obj = SMESH_Gen_i::GetSMESHGen()->getStudyServant()->ConvertIORToObject(aVal);
 354       aGeomObj = GEOM::GEOM_Object::_narrow(obj);
 355     }
 356     geomShape = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
 357     if ( geomShape.IsNull() )
 358       continue;
 359
 360     TopAbs_ShapeEnum geomType = geomShape.ShapeType();
 361     if ( geomType == TopAbs_COMPOUND)// voir s'il ne faut pas mettre une erreur dans le cas contraire
 362     {
 363       MESSAGE("shapeType == TopAbs_COMPOUND");
 364       TopoDS_Iterator it(geomShape);
 365       if ( !it.More() )
 366         continue;
 367       TopAbs_ShapeEnum shapeType = it.Value().ShapeType();
 368       std::vector< std::pair< int, int > > dimTags;
 369       for ( ; it.More(); it.Next())
 370       {
 371         const TopoDS_Shape& topoShape = it.Value();
 372         ASSERT(topoShape.ShapeType() == shapeType);
 373         if ( _mesh->GetMeshDS()->ShapeToIndex( topoShape ) > 0 )
 374           occgeo->importShapes( &topoShape, false, dimTags );
 375         else
 376         {
 377           TopoDS_Shape face = TopExp_Explorer( _shape, shapeType ).Current();
 378           SMESH_subMesh* sm = _mesh->GetSubMesh( face );
 379           sm->GetComputeError() =
 380             SMESH_ComputeError::New
 381             ( COMPERR_WARNING, "Compound shape does not belong to the main geometry. Ignored");
 382         }
 383       }
 384       std::vector<int> tags;
 385       int dim = ( shapeType == TopAbs_EDGE ) ? 1 : 2;
 386       for ( size_t i = 0; i < dimTags.size(); ++i )
 387       {
 388         if ( dimTags[i].first == dim )
 389           tags.push_back( dimTags[i].second );
 390       }
 391       if ( !tags.empty() )
 392       {
 393         _gModel->getGEOInternals()->setCompoundMesh( dim, tags );
 394         toSynchronize = true;
 395       }
 396       if ( toSynchronize )
 397         _gModel->getGEOInternals()->synchronize(_gModel);
 398     }
 399   }
 400 }
 401
 402 //================================================================================
 403 /*!
 404  * \brief For a compound mesh set the mesh components to be transmitted to SMESH
 405  */
 406 //================================================================================
 407
 408 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 409 void GMSHPlugin_Mesher::SetCompoundMeshVisibility()
 410 {
 411
 412   // Loop over all faces, if the face belongs to a compound entry then
 413   // for all (boundary) edges whithin the face visibility is set to 0,
 414   // if the face doesn't belong to a compound entry then visibility is
 415   // set to 1 for all its (boundary) edges. Later, in FillSMesh() func
 416   // getVisibility() (returns either 1 or 0) is used to decide weather
 417   // the mesh of edge should be transmitted  to SMESH or not.
 418
 419   for ( GModel::fiter itF = _gModel->firstFace(); itF != _gModel->lastFace(); ++itF )
 420   {
 421     std::vector< GEdge *> faceEdges = (*itF)->edges();
 422
 423     for ( auto itE = faceEdges.begin(); itE != faceEdges.end(); ++itE )
 424     {
 425       if ( ((*itF)->compound.size()) )
 426         (*itE)->setVisibility(0);
 427       else
 428         (*itE)->setVisibility(1);
 429     }
 430   }
 431
 432
 433   // Loop over all edges, if the  edge belongs to a compound entry then
 434   // for all (boundary) vertices whithin the  edge visibility is set to
 435   // 0, if the edge doesn't belong to a  compound entry then visibility
 436   // is set to 1 for all its (boundary) vertices. Later, in FillSMesh()
 437   // func getVisibility() (returns either 1 or 0) is used to decide we-
 438   // ather the mesh of vertices should be transmitted  to SMESH or not.
 439
 440   for ( GModel::eiter itE = _gModel->firstEdge(); itE != _gModel->lastEdge(); ++itE )
 441   {
 442     std::vector<GVertex *> bndVerticies = (*itE)->vertices();
 443
 444     for( auto itV = bndVerticies.begin(); itV != bndVerticies.end(); ++itV )
 445     {
 446             if(((*itE)->compound.size()))
 447                (*itV)->setVisibility(0);
 448             else
 449                (*itV)->setVisibility(1);
 450     }
 451   }
 452
 453 }
 454 #endif
 455
 456 //================================================================================
 457 /*!
 458  * \brief Write mesh from GModel instance to SMESH instance
 459  */
 460 //================================================================================
 461
 462 void GMSHPlugin_Mesher::FillSMesh()
 463 {
 464   SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
 465
 466   // ADD 0D ELEMENTS
 467   for ( GModel::viter it = _gModel->firstVertex(); it != _gModel->lastVertex(); ++it)
 468   {
 469     GVertex *gVertex = *it;
 470
 471     // GET topoVertex CORRESPONDING TO gVertex
 472     TopoDS_Vertex topoVertex = *((TopoDS_Vertex*)gVertex->getNativePtr());
 473
 474     if (gVertex->getVisibility() == 0) // belongs to a compound
 475     {
 476       SMESH_subMesh* sm = _mesh->GetSubMesh(topoVertex);
 477       sm->SetIsAlwaysComputed(true); // prevent from displaying errors
 478       continue;
 479     }
 480
 481     // FILL SMESH FOR topoVertex
 482     //nodes
 483     for(unsigned int i = 0; i < gVertex->mesh_vertices.size(); i++)
 484     {
 485       MVertex *v = gVertex->mesh_vertices[i];
 486       if(v->getIndex() >= 0)
 487       {
 488         SMDS_MeshNode *node = meshDS->AddNodeWithID(v->x(),v->y(),v->z(),v->getNum());
 489         meshDS->SetNodeOnVertex( node, topoVertex );
 490       }
 491     }
 492     // WE DON'T ADD 0D ELEMENTS because it does not follow the salome meshers philosophy
 493     //elements
 494     // for(unsigned int i = 0; i < gVertex->getNumMeshElements(); i++)
 495     // {
 496     //   MElement *e = gVertex->getMeshElement(i);
 497     //   std::vector<MVertex*> verts;
 498     //   e->getVertices(verts);
 499     //   ASSERT(verts.size()==1);
 500     //   SMDS_Mesh0DElement* zeroDElement = 0;
 501     //   zeroDElement = meshDS->Add0DElementWithID(verts[0]->getNum(),e->getNum());
 502     //   meshDS->SetMeshElementOnShape(zeroDElement, topoVertex);
 503     // }
 504   }
 505
 506   // ADD 1D ELEMENTS
 507   for(GModel::eiter it = _gModel->firstEdge(); it != _gModel->lastEdge(); ++it)
 508   {
 509     GEdge *gEdge = *it;
 510
 511     // GET topoEdge CORRESPONDING TO gEdge
 512     TopoDS_Edge topoEdge;
 513     std::vector< ShapeBounds > topoEdges;
 514 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 515     if(gEdge->haveParametrization())
 516 #else
 517     if ( gEdge->geomType() != GEntity::CompoundCurve )
 518 #endif
 519     {
 520       topoEdge = *((TopoDS_Edge*)gEdge->getNativePtr());
 521       if (gEdge->getVisibility() == 0) // belongs to a compound
 522       {
 523         SMESH_subMesh* sm = _mesh->GetSubMesh(topoEdge);
 524         sm->SetIsAlwaysComputed(true); // prevent from displaying errors
 525         continue;
 526       }
 527     }
 528     bool isCompound = getBoundsOfShapes( gEdge, topoEdges );
 529
 530     // FILL SMESH FOR topoEdge
 531     //nodes
 532     for ( size_t i = 0; i < gEdge->mesh_vertices.size(); i++ )
 533     {
 534       MVertex *v = gEdge->mesh_vertices[i];
 535       if ( v->getIndex() >= 0 )
 536       {
 537         SMDS_MeshNode *node = meshDS->AddNodeWithID(v->x(),v->y(),v->z(),v->getNum());
 538
 539         if ( isCompound )
 540           topoEdge = TopoDS::Edge( getShapeAtPoint( v->point(), topoEdges ));
 541
 542         meshDS->SetNodeOnEdge( node, topoEdge );
 543       }
 544     }
 545   }
 546
 547   for ( GModel::eiter it = _gModel->firstEdge(); it != _gModel->lastEdge(); ++it )
 548   {
 549     GEdge *gEdge = *it;
 550     if ( gEdge->getVisibility() == 0) // belongs to a compound
 551       continue;
 552
 553     TopoDS_Edge topoEdge;
 554     std::vector< ShapeBounds > topoEdges;
 555     bool isCompound = getBoundsOfShapes( gEdge, topoEdges );
 556     if ( !isCompound )
 557       topoEdge = *((TopoDS_Edge*)gEdge->getNativePtr());
 558
 559     //elements
 560     std::vector<MVertex*> verts(3);
 561     for ( size_t i = 0; i < gEdge->getNumMeshElements(); i++ )
 562     {
 563       MElement *e = gEdge->getMeshElement(i);
 564       verts.clear();
 565       e->getVertices(verts);
 566
 567       // if a node wasn't set, it is assigned here
 568       for ( size_t j = 0; j < verts.size(); j++ )
 569       {
 570         if ( verts[j]->onWhat()->getVisibility() == 0 )
 571         {
 572           SMDS_MeshNode *node = meshDS->AddNodeWithID(verts[i]->x(),verts[j]->y(),verts[j]->z(),verts[j]->getNum());
 573           meshDS->SetNodeOnEdge( node, topoEdge );
 574           verts[j]->setEntity(gEdge);
 575         }
 576       }
 577
 578       SMDS_MeshEdge* edge = 0;
 579       switch (verts.size())
 580       {
 581         case 2:
 582           edge = meshDS->AddEdgeWithID(verts[0]->getNum(),
 583                                        verts[1]->getNum(),e->getNum());
 584           break;
 585         case 3:
 586           edge = meshDS->AddEdgeWithID(verts[0]->getNum(),
 587                                        verts[1]->getNum(),
 588                                        verts[2]->getNum(),e->getNum());
 589           break;
 590         default:
 591           ASSERT(false);
 592           continue;
 593       }
 594       if ( isCompound )
 595         topoEdge = TopoDS::Edge( getShapeAtPoint( e->barycenter(), topoEdges ));
 596
 597       meshDS->SetMeshElementOnShape( edge, topoEdge );
 598     }
 599   }
 600
 601   // ADD 2D ELEMENTS
 602   for ( GModel::fiter it = _gModel->firstFace(); it != _gModel->lastFace(); ++it)
 603   {
 604     GFace *gFace = *it;
 605
 606 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 607     // Gmsh since version 4.3 is now producing extra surface and mesh when
 608     // compounds are involved. Since in Gmsh meshing procedure needs acess
 609     // to each of the original topology and the meshed topology. Hence  we
 610     // bypass the additional mesh in case of compounds. Note, similar cri-
 611     // teria also occus in the following 'for' loop.
 612     if ( _compounds.size() && gFace->geomType() == GEntity::DiscreteSurface )
 613       continue;
 614 #endif
 615
 616     // GET topoFace CORRESPONDING TO gFace
 617     TopoDS_Face topoFace;
 618     std::vector< ShapeBounds > topoFaces;
 619
 620 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
 621     if(gFace->haveParametrization())
 622 #else
 623     if ( gFace->geomType() != GEntity::CompoundSurface )
 624 #endif
 625     {
 626       topoFace = *((TopoDS_Face*)gFace->getNativePtr());
 627       if (gFace->getVisibility() == 0) // belongs to a compound
 628       {
 629         SMESH_subMesh* sm = _mesh->GetSubMesh(topoFace);
 630         sm->SetIsAlwaysComputed(true); // prevent from displaying errors
 631         continue;
 632       }
 633     }
 634     bool isCompound = getBoundsOfShapes( gFace, topoFaces );
 635
 636     // FILL SMESH FOR topoFace
 637     //nodes
 638     for ( size_t i = 0; i < gFace->mesh_vertices.size(); i++ )
 639     {
 640       MVertex *v = gFace->mesh_vertices[i];
 641       if ( v->getIndex() >= 0 )
 642       {
 643         SMDS_MeshNode *node = meshDS->AddNodeWithID(v->x(),v->y(),v->z(),v->getNum());
 644
 645         if ( isCompound )
 646           topoFace = TopoDS::Face( getShapeAtPoint( v->point(), topoFaces ));
 647
 648         meshDS->SetNodeOnFace( node, topoFace );
 649       }
 650     }
 651   }
 652
 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     if ( _compounds.size() && gFace->geomType() == GEntity::DiscreteSurface )
 659       continue;
 660
 661     bool isCompound = (!gFace->haveParametrization());
 662 #else
 663     bool isCompound = ( gFace->geomType() == GEntity::CompoundSurface );
 664 #endif
 665
 666     if ( !isCompound && gFace->getVisibility() == 0 )
 667       continue;  // belongs to a compound
 668
 669     TopoDS_Face topoFace;
 670     std::vector< ShapeBounds > topoFaces;
 671     if ( isCompound )
 672       getBoundsOfShapes( gFace, topoFaces );
 673     else
 674       topoFace = *((TopoDS_Face*)gFace->getNativePtr());
 675
 676     //elements
 677     std::vector<MVertex*> verts;
 678     for ( size_t i = 0; i < gFace->getNumMeshElements(); i++ )
 679     {
 680       MElement *e = gFace->getMeshElement(i);
 681       verts.clear();
 682       e->getVertices(verts);
 683       SMDS_MeshFace* face = 0;
 684
 685       // if a node wasn't set, it is assigned here
 686       for ( size_t j = 0; j < verts.size(); j++)
 687       {
 688         if(verts[j]->onWhat()->getVisibility() == 0)
 689         {
 690           SMDS_MeshNode *node = meshDS->AddNodeWithID(verts[j]->x(),verts[j]->y(),verts[j]->z(),verts[j]->getNum());
 691           meshDS->SetNodeOnFace( node, topoFace );
 692           verts[i]->setEntity(gFace);
 693         }
 694       }
 695       switch (verts.size())
 696       {
 697         case 3:
 698           face = meshDS->AddFaceWithID(verts[0]->getNum(),
 699                                        verts[1]->getNum(),
 700                                        verts[2]->getNum(),e->getNum());
 701           break;
 702         case 4:
 703           face = meshDS->AddFaceWithID(verts[0]->getNum(),
 704                                        verts[1]->getNum(),
 705                                        verts[2]->getNum(),
 706                                        verts[3]->getNum(),e->getNum());
 707           break;
 708         case 6:
 709           face = meshDS->AddFaceWithID(verts[0]->getNum(),
 710                                        verts[1]->getNum(),
 711                                        verts[2]->getNum(),
 712                                        verts[3]->getNum(),
 713                                        verts[4]->getNum(),
 714                                        verts[5]->getNum(),e->getNum());
 715           break;
 716         case 8:
 717           face = meshDS->AddFaceWithID(verts[0]->getNum(),
 718                                        verts[1]->getNum(),
 719                                        verts[2]->getNum(),
 720                                        verts[3]->getNum(),
 721                                        verts[4]->getNum(),
 722                                        verts[5]->getNum(),
 723                                        verts[6]->getNum(),
 724                                        verts[7]->getNum(),e->getNum());
 725           break;
 726         case 9:
 727           face = meshDS->AddFaceWithID(verts[0]->getNum(),
 728                                        verts[1]->getNum(),
 729                                        verts[2]->getNum(),
 730                                        verts[3]->getNum(),
 731                                        verts[4]->getNum(),
 732                                        verts[5]->getNum(),
 733                                        verts[6]->getNum(),
 734                                        verts[7]->getNum(),
 735                                        verts[8]->getNum(),e->getNum());
 736           break;
 737         default:
 738           ASSERT(false);
 739           continue;
 740       }
 741
 742       if ( isCompound )
 743         topoFace = TopoDS::Face( getShapeAtPoint( e->barycenter(), topoFaces ));
 744
 745       meshDS->SetMeshElementOnShape(face, topoFace);
 746     }
 747   }
 748
 749   // ADD 3D ELEMENTS
 750   for ( GModel::riter it = _gModel->firstRegion(); it != _gModel->lastRegion(); ++it)
 751   {
 752     GRegion *gRegion = *it;
 753     if (gRegion->getVisibility() == 0)
 754       continue;
 755
 756     // GET topoSolid CORRESPONDING TO gRegion
 757     TopoDS_Solid topoSolid = *((TopoDS_Solid*)gRegion->getNativePtr());
 758
 759     // FILL SMESH FOR topoSolid
 760
 761     //nodes
 762     for(unsigned int i = 0; i < gRegion->mesh_vertices.size(); i++)
 763     {
 764       MVertex *v = gRegion->mesh_vertices[i];
 765       if(v->getIndex() >= 0)
 766       {
 767         SMDS_MeshNode *node = meshDS->AddNodeWithID(v->x(),v->y(),v->z(),v->getNum());
 768         meshDS->SetNodeInVolume( node, topoSolid );
 769       }
 770     }
 771
 772     //elements
 773     std::vector<MVertex*> verts;
 774     for(unsigned int i = 0; i < gRegion->getNumMeshElements(); i++)
 775     {
 776       MElement *e = gRegion->getMeshElement(i);
 777       verts.clear();
 778       e->getVertices(verts);
 779       SMDS_MeshVolume* volume = 0;
 780       switch (verts.size()){
 781         case 4:
 782           volume = meshDS->AddVolumeWithID(verts[0]->getNum(),
 783                                            verts[2]->getNum(),
 784                                            verts[1]->getNum(),
 785                                            verts[3]->getNum(),e->getNum());
 786           break;
 787         case 5:
 788           volume = meshDS->AddVolumeWithID(verts[0]->getNum(),
 789                                            verts[3]->getNum(),
 790                                            verts[2]->getNum(),
 791                                            verts[1]->getNum(),
 792                                            verts[4]->getNum(),e->getNum());
 793           break;
 794         case 6:
 795           volume = meshDS->AddVolumeWithID(verts[0]->getNum(),
 796                                            verts[2]->getNum(),
 797                                            verts[1]->getNum(),
 798                                            verts[3]->getNum(),
 799                                            verts[5]->getNum(),
 800                                            verts[4]->getNum(),e->getNum());
 801           break;
 802         case 8:
 803           volume = meshDS->AddVolumeWithID(verts[0]->getNum(),
 804                                            verts[3]->getNum(),
 805                                            verts[2]->getNum(),
 806                                            verts[1]->getNum(),
 807                                            verts[4]->getNum(),
 808                                            verts[7]->getNum(),
 809                                            verts[6]->getNum(),
 810                                            verts[5]->getNum(),e->getNum());
 811           break;
 812         case 10:
 813           volume = meshDS->AddVolumeWithID(verts[0]->getNum(),
 814                                            verts[2]->getNum(),
 815                                            verts[1]->getNum(),
 816                                            verts[3]->getNum(),
 817                                            verts[6]->getNum(),
 818                                            verts[5]->getNum(),
 819                                            verts[4]->getNum(),
 820                                            verts[7]->getNum(),
 821                                            verts[8]->getNum(),
 822                                            verts[9]->getNum(),e->getNum());
 823           break;
 824         case 13:
 825           volume = meshDS->AddVolumeWithID(verts[0]->getNum(),
 826                                            verts[3]->getNum(),
 827                                            verts[2]->getNum(),
 828                                            verts[1]->getNum(),
 829                                            verts[4]->getNum(),
 830                                            verts[6]->getNum(),
 831                                            verts[10]->getNum(),
 832                                            verts[8]->getNum(),
 833                                            verts[5]->getNum(),
 834                                            verts[7]->getNum(),
 835                                            verts[12]->getNum(),
 836                                            verts[11]->getNum(),
 837                                            verts[9]->getNum(),e->getNum());
 838           break;
 839         case 14: // same as case 13, because no pyra14 in smesh
 840           volume = meshDS->AddVolumeWithID(verts[0]->getNum(),
 841                                            verts[3]->getNum(),
 842                                            verts[2]->getNum(),
 843                                            verts[1]->getNum(),
 844                                            verts[4]->getNum(),
 845                                            verts[6]->getNum(),
 846                                            verts[10]->getNum(),
 847                                            verts[8]->getNum(),
 848                                            verts[5]->getNum(),
 849                                            verts[7]->getNum(),
 850                                            verts[12]->getNum(),
 851                                            verts[11]->getNum(),
 852                                            verts[9]->getNum(),e->getNum());
 853           break;
 854         case 15:
 855           volume = meshDS->AddVolumeWithID(verts[0]->getNum(),
 856                                            verts[2]->getNum(),
 857                                            verts[1]->getNum(),
 858                                            verts[3]->getNum(),
 859                                            verts[5]->getNum(),
 860                                            verts[4]->getNum(),
 861                                            verts[7]->getNum(),
 862                                            verts[9]->getNum(),
 863                                            verts[6]->getNum(),
 864                                            verts[13]->getNum(),
 865                                            verts[14]->getNum(),
 866                                            verts[12]->getNum(),
 867                                            verts[8]->getNum(),
 868                                            verts[11]->getNum(),
 869                                            verts[10]->getNum(),e->getNum());
 870           break;
 871         case 18: // same as case 15, because no penta18 in smesh
 872           volume = meshDS->AddVolumeWithID(verts[0]->getNum(),
 873                                            verts[2]->getNum(),
 874                                            verts[1]->getNum(),
 875                                            verts[3]->getNum(),
 876                                            verts[5]->getNum(),
 877                                            verts[4]->getNum(),
 878                                            verts[7]->getNum(),
 879                                            verts[9]->getNum(),
 880                                            verts[6]->getNum(),
 881                                            verts[13]->getNum(),
 882                                            verts[14]->getNum(),
 883                                            verts[12]->getNum(),
 884                                            verts[8]->getNum(),
 885                                            verts[11]->getNum(),
 886                                            verts[10]->getNum(),e->getNum());
 887           break;
 888         case 20:
 889           volume = meshDS->AddVolumeWithID(verts[0]->getNum(),
 890                                            verts[3]->getNum(),
 891                                            verts[2]->getNum(),
 892                                            verts[1]->getNum(),
 893                                            verts[4]->getNum(),
 894                                            verts[7]->getNum(),
 895                                            verts[6]->getNum(),
 896                                            verts[5]->getNum(),
 897                                            verts[9]->getNum(),
 898                                            verts[13]->getNum(),
 899                                            verts[11]->getNum(),
 900                                            verts[8]->getNum(),
 901                                            verts[17]->getNum(),
 902                                            verts[19]->getNum(),
 903                                            verts[18]->getNum(),
 904                                            verts[16]->getNum(),
 905                                            verts[10]->getNum(),
 906                                            verts[15]->getNum(),
 907                                            verts[14]->getNum(),
 908                                            verts[12]->getNum(),e->getNum());
 909           break;
 910         case 27:
 911           volume = meshDS->AddVolumeWithID(verts[0]->getNum(),
 912                                            verts[3]->getNum(),
 913                                            verts[2]->getNum(),
 914                                            verts[1]->getNum(),
 915                                            verts[4]->getNum(),
 916                                            verts[7]->getNum(),
 917                                            verts[6]->getNum(),
 918                                            verts[5]->getNum(),
 919                                            verts[9]->getNum(),
 920                                            verts[13]->getNum(),
 921                                            verts[11]->getNum(),
 922                                            verts[8]->getNum(),
 923                                            verts[17]->getNum(),
 924                                            verts[19]->getNum(),
 925                                            verts[18]->getNum(),
 926                                            verts[16]->getNum(),
 927                                            verts[10]->getNum(),
 928                                            verts[15]->getNum(),
 929                                            verts[14]->getNum(),
 930                                            verts[12]->getNum(),
 931                                            verts[20]->getNum(),
 932                                            verts[22]->getNum(),
 933                                            verts[24]->getNum(),
 934                                            verts[23]->getNum(),
 935                                            verts[21]->getNum(),
 936                                            verts[25]->getNum(),
 937                                            verts[26]->getNum(),
 938                                            e->getNum());
 939           break;
 940         default:
 941           ASSERT(false);
 942           continue;
 943       }
 944       meshDS->SetMeshElementOnShape(volume, topoSolid);
 945     }
 946   }
 947
 948   //return 0;
 949 }
 950
 951 //================================================================================
 952 /*!
 953  * \brief Find if SPoint point is in SBoundingBox3d bounds
 954  */
 955 //================================================================================
 956
 957 float GMSHPlugin_Mesher::DistBoundingBox(const SBoundingBox3d& bounds, const SPoint3& point)
 958 {
 959   SPoint3 min = bounds.min();
 960   SPoint3 max = bounds.max();
 961
 962   float x,y,z;
 963
 964   if (point.x() < min.x())
 965     x = min.x()-point.x();
 966   else if (point.x() > max.x())
 967     x = point.x()-max.x();
 968   else
 969     x = 0.;
 970
 971   if (point.y() < min.y())
 972     y = min.y()-point.y();
 973   else if (point.y() > max.y())
 974     y = point.y()-max.y();
 975   else
 976     y = 0.;
 977
 978   if (point.z() < min.z())
 979     z = min.z()-point.z();
 980   else if (point.z() > max.z())
 981     z = point.z()-max.z();
 982   else
 983     z = 0.;
 984
 985   return x*x+y*y+z*z;
 986 }
 987 //================================================================================
 988 /*!
 989  * \brief Reimplemented GmshMessage call. Actions done if errors occurs
 990  *        during gmsh meshing. We define here what to display and what to do.
 991  */
 992 //================================================================================
 993 void  GMSHPlugin_Mesher::mymsg::operator()(std::string level, std::string msg)
 994 {
 995   //MESSAGE("level="<< level.c_str() << ", msg=" << msg.c_str()<< "\n");
 996   printf("level=%s msg=%s\n", level.c_str(), msg.c_str());
 997
 998   if(level == "Fatal" || level == "Error")
 999   {
1000     std::ostringstream oss;
1001     if (level == "Fatal")
1002       oss << "Fatal error during Generation of Gmsh Mesh\n";
1003     else
1004       oss << "Error during Generation of Gmsh Mesh\n";
1005     oss << "  " << msg.c_str() << "\n";
1006     GEntity *e = _gModel->getCurrentMeshEntity();
1007     if(e)
1008     {
1009       oss << "  error occurred while meshing entity:\n" <<
1010              "    tag=" << e->tag() << "\n" <<
1011              "    dimension=" << e->dim() << "\n" <<
1012              "    native pointer=" << e->getNativePtr();
1013       //if(e->geomType() != GEntity::CompoundCurve and e->geomType() != GEntity::CompoundSurface)
1014       //{
1015         //SMESH_subMesh *sm = _mesh->GetSubMesh(*((TopoDS_Shape*)e->getNativePtr()));
1016         //SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1017         //SMESH_Comment comment;
1018         //comment << SMESH_Comment(oss.str);
1019         //std::string str = oss.str();
1020         //smError.reset( new SMESH_ComputeError( str ));
1021
1022         // plutot que de faire de la merde ici, on pourait simplement
1023         // remplir une liste pour dire sur quelles entités gmsh se plante
1024         // (puis faire le fillsmesh)
1025         // puis faire une nouvelle routine qui réécrit les messages d'erreur
1026         // probleme : gmsh peut planté en Fatal, dans ce cas pas de fillsmesh
1027       //}
1028     }
1029     if (level == "Fatal")
1030     {
1031         CTX::instance()->lock = 0;
1032         throw oss.str();
1033     }
1034     else
1035         printf("%s\n", oss.str().c_str());
1036   }
1037 }
1038
1039 //=============================================================================
1040 /*!
1041  * Here we are going to use the GMSH mesher
1042  */
1043 //=============================================================================
1044
1045 bool GMSHPlugin_Mesher::Compute()
1046 {
1047   MESSAGE("GMSHPlugin_Mesher::Compute");
1048
1049   int err = 0;
1050
1051 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1052   SetMaxThreadsGmsh();
1053 #endif
1054   GmshInitialize();
1055   SetGmshOptions();
1056   _gModel = new GModel();
1057   mymsg msg(_gModel);
1058   GmshSetMessageHandler(&msg);
1059   _gModel->importOCCShape((void*)&_shape);
1060   if (_compounds.size() > 0) CreateGmshCompounds();
1061   MESSAGE("GModel::Mesh");
1062   try
1063   {
1064     _gModel->mesh((_is2d)?2:3);
1065 #ifdef WITH_SMESH_CANCEL_COMPUTE
1066
1067 #endif
1068   }
1069   catch (std::string& str)
1070   {
1071     err = 1;
1072     MESSAGE(str);
1073   }
1074   catch (...)
1075   {
1076     err = 1;
1077     MESSAGE("Unrecoverable error during Generation of Gmsh Mesh");
1078   }
1079
1080   if (!err)
1081   {
1082 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1083     if (_compounds.size() > 0)
1084       SetCompoundMeshVisibility();
1085 #endif
1086     FillSMesh();
1087   }
1088   GmshSetMessageHandler(nullptr);
1089   delete _gModel;
1090   GmshFinalize();
1091   MESSAGE("GMSHPlugin_Mesher::Compute:End");
1092   return !err;
1093 }