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