src/GMSHPlugin/GMSHPlugin_Mesher.cxx

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