src/GMSHPlugin/GMSHPlugin_Mesher.cxx

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