src/GMSHPlugin/GMSHPlugin_Mesher.cxx

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