src/GMSHPlugin/GMSHPlugin_Mesher.cxx

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