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