src/GMSHPlugin/GMSHPlugin_Mesher.cxx

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