src/GMSHPlugin/GMSHPlugin_Mesher.cxx

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