src/GMSHPlugin/GMSHPlugin_Mesher.cxx

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