src/NETGENPlugin/NETGENPlugin_Mesher.cxx

   1 // Copyright (C) 2007-2014  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 //  NETGENPlugin : C++ implementation
  24 // File      : NETGENPlugin_Mesher.cxx
  25 // Author    : Michael Sazonov (OCN)
  26 // Date      : 31/03/2006
  27 // Project   : SALOME
  28 //=============================================================================
  29
  30 #include "NETGENPlugin_Mesher.hxx"
  31 #include "NETGENPlugin_Hypothesis_2D.hxx"
  32 #include "NETGENPlugin_SimpleHypothesis_3D.hxx"
  33
  34 #include <SMDS_FaceOfNodes.hxx>
  35 #include <SMDS_MeshElement.hxx>
  36 #include <SMDS_MeshNode.hxx>
  37 #include <SMESHDS_Mesh.hxx>
  38 #include <SMESH_Block.hxx>
  39 #include <SMESH_Comment.hxx>
  40 #include <SMESH_ComputeError.hxx>
  41 #include <SMESH_File.hxx>
  42 #include <SMESH_Gen_i.hxx>
  43 #include <SMESH_Mesh.hxx>
  44 #include <SMESH_MesherHelper.hxx>
  45 #include <SMESH_subMesh.hxx>
  46 #include <StdMeshers_QuadToTriaAdaptor.hxx>
  47 #include <StdMeshers_ViscousLayers2D.hxx>
  48
  49 #include <SALOMEDS_Tool.hxx>
  50
  51 #include <utilities.h>
  52
  53 #include <BRepBuilderAPI_Copy.hxx>
  54 #include <BRep_Tool.hxx>
  55 #include <Bnd_B3d.hxx>
  56 #include <NCollection_Map.hxx>
  57 #include <Standard_ErrorHandler.hxx>
  58 #include <Standard_ProgramError.hxx>
  59 #include <TopExp.hxx>
  60 #include <TopExp_Explorer.hxx>
  61 #include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
  62 #include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
  63 #include <TopTools_DataMapOfShapeInteger.hxx>
  64 #include <TopTools_DataMapOfShapeShape.hxx>
  65 #include <TopTools_MapOfShape.hxx>
  66 #include <TopoDS.hxx>
  67 #include <OSD_File.hxx>
  68 #include <OSD_Path.hxx>
  69
  70 // Netgen include files
  71 #ifndef OCCGEOMETRY
  72 #define OCCGEOMETRY
  73 #endif
  74 #include <occgeom.hpp>
  75 #include <meshing.hpp>
  76 //#include <ngexception.hpp>
  77 namespace netgen {
  78 #ifdef NETGEN_V5
  79   extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, MeshingParameters&, int, int);
  80 #else
  81   extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
  82 #endif
  83   //extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh);
  84   extern MeshingParameters mparam;
  85   extern volatile multithreadt multithread;
  86   extern bool merge_solids;
  87 }
  88
  89 #include <vector>
  90 #include <limits>
  91
  92 #ifdef WIN32
  93 #include <process.h>
  94 #endif
  95 using namespace nglib;
  96 using namespace std;
  97
  98 #ifdef _DEBUG_
  99 #define nodeVec_ACCESS(index) ((SMDS_MeshNode*) nodeVec.at((index)))
 100 #else
 101 #define nodeVec_ACCESS(index) ((SMDS_MeshNode*) nodeVec[index])
 102 #endif
 103
 104 #define NGPOINT_COORDS(p) p(0),p(1),p(2)
 105
 106 // dump elements added to ng mesh
 107 //#define DUMP_SEGMENTS
 108 //#define DUMP_TRIANGLES
 109 //#define DUMP_TRIANGLES_SCRIPT "/tmp/trias.py" //!< debug AddIntVerticesInSolids()
 110
 111 TopTools_IndexedMapOfShape ShapesWithLocalSize;
 112 std::map<int,double> VertexId2LocalSize;
 113 std::map<int,double> EdgeId2LocalSize;
 114 std::map<int,double> FaceId2LocalSize;
 115
 116 //=============================================================================
 117 /*!
 118  *
 119  */
 120 //=============================================================================
 121
 122 NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh*         mesh,
 123                                           const TopoDS_Shape& aShape,
 124                                           const bool          isVolume)
 125   : _mesh    (mesh),
 126     _shape   (aShape),
 127     _isVolume(isVolume),
 128     _optimize(true),
 129     _fineness(NETGENPlugin_Hypothesis::GetDefaultFineness()),
 130     _isViscousLayers2D(false),
 131     _ngMesh(NULL),
 132     _occgeom(NULL),
 133     _curShapeIndex(-1),
 134     _progressTic(1),
 135     _totalTime(1.0),
 136     _simpleHyp(NULL),
 137     _ptrToMe(NULL)
 138 {
 139   SetDefaultParameters();
 140   ShapesWithLocalSize.Clear();
 141   VertexId2LocalSize.clear();
 142   EdgeId2LocalSize.clear();
 143   FaceId2LocalSize.clear();
 144 }
 145
 146 //================================================================================
 147 /*!
 148  * Destuctor
 149  */
 150 //================================================================================
 151
 152 NETGENPlugin_Mesher::~NETGENPlugin_Mesher()
 153 {
 154   if ( _ptrToMe )
 155     *_ptrToMe = NULL;
 156   _ptrToMe = 0;
 157   _ngMesh = NULL;
 158 }
 159
 160 //================================================================================
 161 /*!
 162  * Set pointer to NETGENPlugin_Mesher* field of the holder, that will be
 163  * nullified at destruction of this
 164  */
 165 //================================================================================
 166
 167 void NETGENPlugin_Mesher::SetSelfPointer( NETGENPlugin_Mesher ** ptr )
 168 {
 169   if ( _ptrToMe )
 170     *_ptrToMe = NULL;
 171
 172   _ptrToMe = ptr;
 173
 174   if ( _ptrToMe )
 175     *_ptrToMe = this;
 176 }
 177
 178 //================================================================================
 179 /*!
 180  * \brief Initialize global NETGEN parameters with default values
 181  */
 182 //================================================================================
 183
 184 void NETGENPlugin_Mesher::SetDefaultParameters()
 185 {
 186   netgen::MeshingParameters& mparams = netgen::mparam;
 187   // maximal mesh edge size
 188   mparams.maxh            = 0;//NETGENPlugin_Hypothesis::GetDefaultMaxSize();
 189   mparams.minh            = 0;
 190   // minimal number of segments per edge
 191   mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
 192   // rate of growth of size between elements
 193   mparams.grading         = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
 194   // safety factor for curvatures (elements per radius)
 195   mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
 196   // create elements of second order
 197   mparams.secondorder     = NETGENPlugin_Hypothesis::GetDefaultSecondOrder();
 198   // quad-dominated surface meshing
 199   if (_isVolume)
 200     mparams.quad          = 0;
 201   else
 202     mparams.quad          = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed();
 203   _fineness               = NETGENPlugin_Hypothesis::GetDefaultFineness();
 204   mparams.uselocalh       = NETGENPlugin_Hypothesis::GetDefaultSurfaceCurvature();
 205   netgen::merge_solids    = NETGENPlugin_Hypothesis::GetDefaultFuseEdges();
 206 }
 207
 208 //=============================================================================
 209 /*!
 210  *
 211  */
 212 //=============================================================================
 213 void SetLocalSize(TopoDS_Shape GeomShape, double LocalSize)
 214 {
 215   TopAbs_ShapeEnum GeomType = GeomShape.ShapeType();
 216   if (GeomType == TopAbs_COMPOUND) {
 217     for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()) {
 218       SetLocalSize(it.Value(), LocalSize);
 219     }
 220     return;
 221   }
 222   int key;
 223   if (! ShapesWithLocalSize.Contains(GeomShape))
 224     key = ShapesWithLocalSize.Add(GeomShape);
 225   else
 226     key = ShapesWithLocalSize.FindIndex(GeomShape);
 227   if (GeomType == TopAbs_VERTEX) {
 228     VertexId2LocalSize[key] = LocalSize;
 229   } else if (GeomType == TopAbs_EDGE) {
 230     EdgeId2LocalSize[key] = LocalSize;
 231   } else if (GeomType == TopAbs_FACE) {
 232     FaceId2LocalSize[key] = LocalSize;
 233   }
 234 }
 235
 236 //=============================================================================
 237 /*!
 238  * Pass parameters to NETGEN
 239  */
 240 //=============================================================================
 241 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
 242 {
 243   if (hyp)
 244   {
 245     netgen::MeshingParameters& mparams = netgen::mparam;
 246     // Initialize global NETGEN parameters:
 247     // maximal mesh segment size
 248     mparams.maxh            = hyp->GetMaxSize();
 249     // maximal mesh element linear size
 250     mparams.minh            = hyp->GetMinSize();
 251     // minimal number of segments per edge
 252     mparams.segmentsperedge = hyp->GetNbSegPerEdge();
 253     // rate of growth of size between elements
 254     mparams.grading         = hyp->GetGrowthRate();
 255     // safety factor for curvatures (elements per radius)
 256     mparams.curvaturesafety = hyp->GetNbSegPerRadius();
 257     // create elements of second order
 258     mparams.secondorder     = hyp->GetSecondOrder() ? 1 : 0;
 259     // quad-dominated surface meshing
 260     // only triangles are allowed for volumic mesh (before realizing IMP 0021676)
 261     //if (!_isVolume)
 262       mparams.quad          = hyp->GetQuadAllowed() ? 1 : 0;
 263     _optimize               = hyp->GetOptimize();
 264     _fineness               = hyp->GetFineness();
 265     mparams.uselocalh       = hyp->GetSurfaceCurvature();
 266     netgen::merge_solids    = hyp->GetFuseEdges();
 267     _simpleHyp = NULL;
 268
 269     SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen();
 270     CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager");
 271     SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject);
 272     SALOMEDS::Study_var myStudy = aStudyMgr->GetStudyByID(hyp->GetStudyId());
 273
 274     const NETGENPlugin_Hypothesis::TLocalSize localSizes = hyp->GetLocalSizesAndEntries();
 275     NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin();
 276     for (it ; it != localSizes.end() ; it++)
 277       {
 278         std::string entry = (*it).first;
 279         double val = (*it).second;
 280         // --
 281         GEOM::GEOM_Object_var aGeomObj;
 282         TopoDS_Shape S = TopoDS_Shape();
 283         SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
 284         if (!aSObj->_is_nil()) {
 285           CORBA::Object_var obj = aSObj->GetObject();
 286           aGeomObj = GEOM::GEOM_Object::_narrow(obj);
 287           aSObj->UnRegister();
 288         }
 289         if ( !aGeomObj->_is_nil() )
 290           S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
 291         // --
 292         SetLocalSize(S, val);
 293       }
 294   }
 295 }
 296
 297 //=============================================================================
 298 /*!
 299  * Pass simple parameters to NETGEN
 300  */
 301 //=============================================================================
 302
 303 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp)
 304 {
 305   _simpleHyp = hyp;
 306   if ( _simpleHyp )
 307     SetDefaultParameters();
 308 }
 309
 310 //=============================================================================
 311 /*!
 312  *  Link - a pair of integer numbers
 313  */
 314 //=============================================================================
 315 struct Link
 316 {
 317   int n1, n2;
 318   Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
 319   Link() : n1(0), n2(0) {}
 320 };
 321
 322 int HashCode(const Link& aLink, int aLimit)
 323 {
 324   return HashCode(aLink.n1 + aLink.n2, aLimit);
 325 }
 326
 327 Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
 328 {
 329   return (aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ||
 330           aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1);
 331 }
 332
 333 namespace
 334 {
 335   //================================================================================
 336   /*!
 337    * \brief return id of netgen point corresponding to SMDS node
 338    */
 339   //================================================================================
 340   typedef map< const SMDS_MeshNode*, int > TNode2IdMap;
 341
 342   int ngNodeId( const SMDS_MeshNode* node,
 343                 netgen::Mesh&        ngMesh,
 344                 TNode2IdMap&         nodeNgIdMap)
 345   {
 346     int newNgId = ngMesh.GetNP() + 1;
 347
 348     TNode2IdMap::iterator node_id = nodeNgIdMap.insert( make_pair( node, newNgId )).first;
 349
 350     if ( node_id->second == newNgId)
 351     {
 352 #if defined(DUMP_SEGMENTS) || defined(DUMP_TRIANGLES)
 353       cout << "Ng " << newNgId << " - " << node;
 354 #endif
 355       netgen::MeshPoint p( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
 356       ngMesh.AddPoint( p );
 357     }
 358     return node_id->second;
 359   }
 360
 361   //================================================================================
 362   /*!
 363    * \brief Return computed EDGEs connected to the given one
 364    */
 365   //================================================================================
 366
 367   list< TopoDS_Edge > getConnectedEdges( const TopoDS_Edge&                 edge,
 368                                          const TopoDS_Face&                 face,
 369                                          const set< SMESH_subMesh* > &      computedSM,
 370                                          const SMESH_MesherHelper&          helper,
 371                                          map< SMESH_subMesh*, set< int > >& addedEdgeSM2Faces)
 372   {
 373     // get ordered EDGEs
 374     list< TopoDS_Edge > edges;
 375     list< int > nbEdgesInWire;
 376     int nbWires = SMESH_Block::GetOrderedEdges( face, edges, nbEdgesInWire);
 377
 378     // find <edge> within <edges>
 379     list< TopoDS_Edge >::iterator eItFwd = edges.begin();
 380     for ( ; eItFwd != edges.end(); ++eItFwd )
 381       if ( edge.IsSame( *eItFwd ))
 382         break;
 383     if ( eItFwd == edges.end()) return list< TopoDS_Edge>();
 384
 385     if ( eItFwd->Orientation() >= TopAbs_INTERNAL )
 386     {
 387       // connected INTERNAL edges returned from GetOrderedEdges() are wrongly oriented
 388       // so treat each INTERNAL edge separately
 389       TopoDS_Edge e = *eItFwd;
 390       edges.clear();
 391       edges.push_back( e );
 392       return edges;
 393     }
 394
 395     // get all computed EDGEs connected to <edge>
 396
 397     list< TopoDS_Edge >::iterator eItBack = eItFwd, ePrev;
 398     TopoDS_Vertex vCommon;
 399     TopTools_MapOfShape eAdded; // map used not to add a seam edge twice to <edges>
 400     eAdded.Add( edge );
 401
 402     // put edges before <edge> to <edges> back
 403     while ( edges.begin() != eItFwd )
 404       edges.splice( edges.end(), edges, edges.begin() );
 405
 406     // search forward
 407     ePrev = eItFwd;
 408     while ( ++eItFwd != edges.end() )
 409     {
 410       SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( *eItFwd );
 411
 412       bool connected = TopExp::CommonVertex( *ePrev, *eItFwd, vCommon );
 413       bool computed  = sm->IsMeshComputed();
 414       bool added     = addedEdgeSM2Faces[ sm ].count( helper.GetSubShapeID() );
 415       bool doubled   = !eAdded.Add( *eItFwd );
 416       bool orientOK  = (( ePrev ->Orientation() < TopAbs_INTERNAL ) ==
 417                         ( eItFwd->Orientation() < TopAbs_INTERNAL )    );
 418       if ( !connected || !computed || !orientOK || added || doubled )
 419       {
 420         // stop advancement; move edges from tail to head
 421         while ( edges.back() != *ePrev )
 422           edges.splice( edges.begin(), edges, --edges.end() );
 423         break;
 424       }
 425       ePrev = eItFwd;
 426     }
 427     // search backward
 428     while ( eItBack != edges.begin() )
 429     {
 430       ePrev = eItBack;
 431       --eItBack;
 432       SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( *eItBack );
 433
 434       bool connected = TopExp::CommonVertex( *ePrev, *eItBack, vCommon );
 435       bool computed  = sm->IsMeshComputed();
 436       bool added     = addedEdgeSM2Faces[ sm ].count( helper.GetSubShapeID() );
 437       bool doubled   = !eAdded.Add( *eItBack );
 438       bool orientOK  = (( ePrev  ->Orientation() < TopAbs_INTERNAL ) ==
 439                         ( eItBack->Orientation() < TopAbs_INTERNAL )    );
 440       if ( !connected || !computed || !orientOK || added || doubled)
 441       {
 442         // stop advancement
 443         edges.erase( edges.begin(), ePrev );
 444         break;
 445       }
 446     }
 447     if ( edges.front() != edges.back() )
 448     {
 449       // assure that the 1st vertex is meshed
 450       TopoDS_Edge eLast = edges.back();
 451       while ( !SMESH_Algo::VertexNode( SMESH_MesherHelper::IthVertex( 0, edges.front()), helper.GetMeshDS())
 452               &&
 453               edges.front() != eLast )
 454         edges.splice( edges.end(), edges, edges.begin() );
 455     }
 456     return edges;
 457   }
 458
 459   //================================================================================
 460   /*!
 461    * \brief Make triangulation of a shape precise enough
 462    */
 463   //================================================================================
 464
 465   void updateTriangulation( const TopoDS_Shape& shape )
 466   {
 467     // static set< Poly_Triangulation* > updated;
 468
 469     // TopLoc_Location loc;
 470     // TopExp_Explorer fExp( shape, TopAbs_FACE );
 471     // for ( ; fExp.More(); fExp.Next() )
 472     // {
 473     //   Handle(Poly_Triangulation) triangulation =
 474     //     BRep_Tool::Triangulation ( TopoDS::Face( fExp.Current() ), loc);
 475     //   if ( triangulation.IsNull() ||
 476     //        updated.insert( triangulation.operator->() ).second )
 477     //   {
 478     //     BRepTools::Clean (shape);
 479         try {
 480           OCC_CATCH_SIGNALS;
 481           BRepMesh_IncrementalMesh e(shape, 0.01, true);
 482         }
 483         catch (Standard_Failure)
 484         {
 485         }
 486   //       updated.erase( triangulation.operator->() );
 487   //       triangulation = BRep_Tool::Triangulation ( TopoDS::Face( fExp.Current() ), loc);
 488   //       updated.insert( triangulation.operator->() );
 489   //     }
 490   //   }
 491   }
 492   //================================================================================
 493   /*!
 494    * \brief Returns a medium node either existing in SMESH of created by NETGEN
 495    *  \param [in] corner1 - corner node 1
 496    *  \param [in] corner2 - corner node 2
 497    *  \param [in] defaultMedium - the node created by NETGEN
 498    *  \param [in] helper - holder of medium nodes existing in SMESH
 499    *  \return const SMDS_MeshNode* - the result node
 500    */
 501   //================================================================================
 502
 503   const SMDS_MeshNode* mediumNode( const SMDS_MeshNode*      corner1,
 504                                    const SMDS_MeshNode*      corner2,
 505                                    const SMDS_MeshNode*      defaultMedium,
 506                                    const SMESH_MesherHelper* helper)
 507   {
 508     if ( helper )
 509     {
 510       TLinkNodeMap::const_iterator l2n =
 511         helper->GetTLinkNodeMap().find( SMESH_TLink( corner1, corner2 ));
 512       if ( l2n != helper->GetTLinkNodeMap().end() )
 513         defaultMedium = l2n->second;
 514     }
 515     return defaultMedium;
 516   }
 517
 518   //================================================================================
 519   /*!
 520    * \brief Assure that mesh on given shapes is quadratic
 521    */
 522   //================================================================================
 523
 524   void makeQuadratic( const TopTools_IndexedMapOfShape& shapes,
 525                       SMESH_Mesh*                       mesh )
 526   {
 527     for ( int i = 1; i <= shapes.Extent(); ++i )
 528     {
 529       SMESHDS_SubMesh* smDS = mesh->GetMeshDS()->MeshElements( shapes(i) );
 530       if ( !smDS ) continue;
 531       SMDS_ElemIteratorPtr elemIt = smDS->GetElements();
 532       if ( !elemIt->more() ) continue;
 533       const SMDS_MeshElement* e = elemIt->next();
 534       if ( !e || e->IsQuadratic() )
 535         continue;
 536
 537       TIDSortedElemSet elems;
 538       elems.insert( e );
 539       while ( elemIt->more() )
 540         elems.insert( elems.end(), elemIt->next() );
 541
 542       SMESH_MeshEditor( mesh ).ConvertToQuadratic( /*3d=*/false, elems, /*biQuad=*/false );
 543     }
 544   }
 545
 546 }
 547
 548 //================================================================================
 549 /*!
 550  * \brief Initialize netgen::OCCGeometry with OCCT shape
 551  */
 552 //================================================================================
 553
 554 void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry&     occgeo,
 555                                              const TopoDS_Shape&      shape,
 556                                              SMESH_Mesh&              mesh,
 557                                              list< SMESH_subMesh* > * meshedSM,
 558                                              NETGENPlugin_Internals*  intern)
 559 {
 560   updateTriangulation( shape );
 561
 562   Bnd_Box bb;
 563   BRepBndLib::Add (shape, bb);
 564   double x1,y1,z1,x2,y2,z2;
 565   bb.Get (x1,y1,z1,x2,y2,z2);
 566   MESSAGE("shape bounding box:\n" <<
 567           "(" << x1 << " " << y1 << " " << z1 << ") " <<
 568           "(" << x2 << " " << y2 << " " << z2 << ")");
 569   netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
 570   netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
 571   occgeo.boundingbox = netgen::Box<3> (p1,p2);
 572
 573   occgeo.shape = shape;
 574   occgeo.changed = 1;
 575
 576   // fill maps of shapes of occgeo with not yet meshed subshapes
 577
 578   // get root submeshes
 579   list< SMESH_subMesh* > rootSM;
 580   const int shapeID = mesh.GetMeshDS()->ShapeToIndex( shape );
 581   if ( shapeID > 0 ) { // SMESH_subMesh with ID 0 may exist, don't use it!
 582     rootSM.push_back( mesh.GetSubMesh( shape ));
 583   }
 584   else {
 585     for ( TopoDS_Iterator it( shape ); it.More(); it.Next() )
 586       rootSM.push_back( mesh.GetSubMesh( it.Value() ));
 587   }
 588
 589   // add subshapes of empty submeshes
 590   list< SMESH_subMesh* >::iterator rootIt = rootSM.begin(), rootEnd = rootSM.end();
 591   for ( ; rootIt != rootEnd; ++rootIt ) {
 592     SMESH_subMesh * root = *rootIt;
 593     SMESH_subMeshIteratorPtr smIt = root->getDependsOnIterator(/*includeSelf=*/true,
 594                                                                /*complexShapeFirst=*/true);
 595     // to find a right orientation of subshapes (PAL20462)
 596     TopTools_IndexedMapOfShape subShapes;
 597     TopExp::MapShapes(root->GetSubShape(), subShapes);
 598     while ( smIt->more() )
 599     {
 600       SMESH_subMesh* sm = smIt->next();
 601       TopoDS_Shape shape = sm->GetSubShape();
 602       if ( intern && intern->isShapeToPrecompute( shape ))
 603         continue;
 604       if ( !meshedSM || sm->IsEmpty() )
 605       {
 606         if ( shape.ShapeType() != TopAbs_VERTEX )
 607           shape = subShapes( subShapes.FindIndex( shape ));// shape -> index -> oriented shape
 608         if ( shape.Orientation() >= TopAbs_INTERNAL )
 609           shape.Orientation( TopAbs_FORWARD ); // isuue 0020676
 610         switch ( shape.ShapeType() ) {
 611         case TopAbs_FACE  : occgeo.fmap.Add( shape ); break;
 612         case TopAbs_EDGE  : occgeo.emap.Add( shape ); break;
 613         case TopAbs_VERTEX: occgeo.vmap.Add( shape ); break;
 614         case TopAbs_SOLID :occgeo.somap.Add( shape ); break;
 615         default:;
 616         }
 617       }
 618       // collect submeshes of meshed shapes
 619       else if (meshedSM)
 620       {
 621         const int dim = SMESH_Gen::GetShapeDim( shape );
 622         meshedSM[ dim ].push_back( sm );
 623       }
 624     }
 625   }
 626   occgeo.facemeshstatus.SetSize (occgeo.fmap.Extent());
 627   occgeo.facemeshstatus = 0;
 628   occgeo.face_maxh_modified.SetSize(occgeo.fmap.Extent());
 629   occgeo.face_maxh_modified = 0;
 630   occgeo.face_maxh.SetSize(occgeo.fmap.Extent());
 631   occgeo.face_maxh = netgen::mparam.maxh;
 632 }
 633
 634 //================================================================================
 635 /*!
 636  * \brief Return a default min size value suitable for the given geometry.
 637  */
 638 //================================================================================
 639
 640 double NETGENPlugin_Mesher::GetDefaultMinSize(const TopoDS_Shape& geom,
 641                                               const double        maxSize)
 642 {
 643   updateTriangulation( geom );
 644
 645   TopLoc_Location loc;
 646   int i1, i2, i3;
 647   const int* pi[4] = { &i1, &i2, &i3, &i1 };
 648   double minh = 1e100;
 649   Bnd_B3d bb;
 650   TopExp_Explorer fExp( geom, TopAbs_FACE );
 651   for ( ; fExp.More(); fExp.Next() )
 652   {
 653     Handle(Poly_Triangulation) triangulation =
 654       BRep_Tool::Triangulation ( TopoDS::Face( fExp.Current() ), loc);
 655     if ( triangulation.IsNull() ) continue;
 656     const double fTol = BRep_Tool::Tolerance( TopoDS::Face( fExp.Current() ));
 657     const TColgp_Array1OfPnt&   points = triangulation->Nodes();
 658     const Poly_Array1OfTriangle& trias = triangulation->Triangles();
 659     for ( int iT = trias.Lower(); iT <= trias.Upper(); ++iT )
 660     {
 661       trias(iT).Get( i1, i2, i3 );
 662       for ( int j = 0; j < 3; ++j )
 663       {
 664         double dist2 = points(*pi[j]).SquareDistance( points( *pi[j+1] ));
 665         if ( dist2 < minh && fTol*fTol < dist2 )
 666           minh = dist2;
 667         bb.Add( points(*pi[j]));
 668       }
 669     }
 670   }
 671   if ( minh > 0.25 * bb.SquareExtent() ) // simple geometry, rough triangulation
 672   {
 673     minh = 1e-3 * sqrt( bb.SquareExtent());
 674     //cout << "BND BOX minh = " <<minh << endl;
 675   }
 676   else
 677   {
 678     minh = 3 * sqrt( minh ); // triangulation for visualization is rather fine
 679     //cout << "TRIANGULATION minh = " <<minh << endl;
 680   }
 681   if ( minh > 0.5 * maxSize )
 682     minh = maxSize / 3.;
 683
 684   return minh;
 685 }
 686
 687 //================================================================================
 688 /*!
 689  * \brief Restrict size of elements at a given point
 690  */
 691 //================================================================================
 692
 693 void NETGENPlugin_Mesher::RestrictLocalSize(netgen::Mesh& ngMesh, const gp_XYZ& p, const double size)
 694 {
 695   if ( netgen::mparam.minh > size )
 696   {
 697     ngMesh.SetMinimalH( size );
 698     netgen::mparam.minh = size;
 699   }
 700   netgen::Point3d pi(p.X(), p.Y(), p.Z());
 701   ngMesh.RestrictLocalH( pi, size );
 702 }
 703
 704 //================================================================================
 705 /*!
 706  * \brief fill ngMesh with nodes and elements of computed submeshes
 707  */
 708 //================================================================================
 709
 710 bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
 711                                      netgen::Mesh&                  ngMesh,
 712                                      vector<const SMDS_MeshNode*>&  nodeVec,
 713                                      const list< SMESH_subMesh* > & meshedSM,
 714                                      SMESH_MesherHelper*            quadHelper,
 715                                      SMESH_ProxyMesh::Ptr           proxyMesh)
 716 {
 717   TNode2IdMap nodeNgIdMap;
 718   for ( int i = 1; i < nodeVec.size(); ++i )
 719     nodeNgIdMap.insert( make_pair( nodeVec[i], i ));
 720
 721   TopTools_MapOfShape visitedShapes;
 722   map< SMESH_subMesh*, set< int > > visitedEdgeSM2Faces;
 723   set< SMESH_subMesh* > computedSM( meshedSM.begin(), meshedSM.end() );
 724
 725   SMESH_MesherHelper helper (*_mesh);
 726
 727   int faceNgID = ngMesh.GetNFD();
 728
 729   list< SMESH_subMesh* >::const_iterator smIt, smEnd = meshedSM.end();
 730   for ( smIt = meshedSM.begin(); smIt != smEnd; ++smIt )
 731   {
 732     SMESH_subMesh* sm = *smIt;
 733     if ( !visitedShapes.Add( sm->GetSubShape() ))
 734       continue;
 735
 736     const SMESHDS_SubMesh * smDS = sm->GetSubMeshDS();
 737     if ( !smDS ) continue;
 738
 739     switch ( sm->GetSubShape().ShapeType() )
 740     {
 741     case TopAbs_EDGE: { // EDGE
 742       // ----------------------
 743       TopoDS_Edge geomEdge  = TopoDS::Edge( sm->GetSubShape() );
 744       if ( geomEdge.Orientation() >= TopAbs_INTERNAL )
 745         geomEdge.Orientation( TopAbs_FORWARD ); // issue 0020676
 746
 747       // Add ng segments for each not meshed FACE the EDGE bounds
 748       PShapeIteratorPtr fIt = helper.GetAncestors( geomEdge, *sm->GetFather(), TopAbs_FACE );
 749       while ( const TopoDS_Shape * anc = fIt->next() )
 750       {
 751         faceNgID = occgeom.fmap.FindIndex( *anc );
 752         if ( faceNgID < 1 )
 753           continue; // meshed face
 754
 755         int faceSMDSId = helper.GetMeshDS()->ShapeToIndex( *anc );
 756         if ( visitedEdgeSM2Faces[ sm ].count( faceSMDSId ))
 757           continue; // already treated EDGE
 758
 759         TopoDS_Face face = TopoDS::Face( occgeom.fmap( faceNgID ));
 760         if ( face.Orientation() >= TopAbs_INTERNAL )
 761           face.Orientation( TopAbs_FORWARD ); // issue 0020676
 762
 763         // get all meshed EDGEs of the FACE connected to geomEdge (issue 0021140)
 764         helper.SetSubShape( face );
 765         list< TopoDS_Edge > edges = getConnectedEdges( geomEdge, face, computedSM, helper,
 766                                                        visitedEdgeSM2Faces );
 767         if ( edges.empty() )
 768           continue; // wrong ancestor?
 769
 770         // find out orientation of <edges> within <face>
 771         TopoDS_Edge eNotSeam = edges.front();
 772         if ( helper.HasSeam() )
 773         {
 774           list< TopoDS_Edge >::iterator eIt = edges.begin();
 775           while ( helper.IsRealSeam( *eIt )) ++eIt;
 776           if ( eIt != edges.end() )
 777             eNotSeam = *eIt;
 778         }
 779         TopAbs_Orientation fOri = helper.GetSubShapeOri( face, eNotSeam );
 780         bool isForwad = ( fOri == eNotSeam.Orientation() || fOri >= TopAbs_INTERNAL );
 781
 782         // get all nodes from connected <edges>
 783         bool isQuad   = smDS->NbElements() ? smDS->GetElements()->next()->IsQuadratic() : false;
 784         StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad );
 785         const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
 786         int i, nbSeg = fSide.NbSegments();
 787
 788         // remember EDGEs of fSide to treat only once
 789         for ( int iE = 0; iE < fSide.NbEdges(); ++iE )
 790           visitedEdgeSM2Faces[ helper.GetMesh()->GetSubMesh( fSide.Edge(iE )) ].insert(faceSMDSId);
 791
 792         double otherSeamParam = 0;
 793         bool isSeam = false;
 794
 795         // add segments
 796
 797         int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
 798
 799         for ( i = 0; i < nbSeg; ++i )
 800         {
 801           const UVPtStruct& p1 = points[ i ];
 802           const UVPtStruct& p2 = points[ i+1 ];
 803
 804           if ( p1.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX ) //an EDGE begins
 805           {
 806             isSeam = false;
 807             if ( helper.IsRealSeam( p1.node->getshapeId() ))
 808             {
 809               TopoDS_Edge e = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
 810               isSeam = helper.IsRealSeam( e );
 811               if ( isSeam )
 812               {
 813                 otherSeamParam = helper.GetOtherParam( helper.GetPeriodicIndex() & 1 ? p2.u : p2.v );
 814               }
 815             }
 816           }
 817           netgen::Segment seg;
 818           // ng node ids
 819           seg[0] = prevNgId;
 820           seg[1] = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
 821           // node param on curve
 822           seg.epgeominfo[ 0 ].dist = p1.param;
 823           seg.epgeominfo[ 1 ].dist = p2.param;
 824           // uv on face
 825           seg.epgeominfo[ 0 ].u = p1.u;
 826           seg.epgeominfo[ 0 ].v = p1.v;
 827           seg.epgeominfo[ 1 ].u = p2.u;
 828           seg.epgeominfo[ 1 ].v = p2.v;
 829
 830           //geomEdge = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
 831           //seg.epgeominfo[ 0 ].edgenr = seg.epgeominfo[ 1 ].edgenr = occgeom.emap.FindIndex( geomEdge );
 832
 833           //seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
 834           seg.si = faceNgID;                   // = geom.fmap.FindIndex (face);
 835           seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
 836           ngMesh.AddSegment (seg);
 837
 838           SMESH_TNodeXYZ np1( p1.node ), np2( p2.node );
 839           RestrictLocalSize( ngMesh, 0.5*(np1+np2), (np1-np2).Modulus() );
 840
 841 #ifdef DUMP_SEGMENTS
 842           cout << "Segment: " << seg.edgenr << " on SMESH face " << helper.GetMeshDS()->ShapeToIndex( face ) << endl
 843                << "\tface index: " << seg.si << endl
 844                << "\tp1: " << seg[0] << endl
 845                << "\tp2: " << seg[1] << endl
 846                << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
 847                << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
 848             //<< "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
 849                << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
 850                << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl;
 851             //<< "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
 852 #endif
 853           if ( isSeam )
 854           {
 855             if ( helper.GetPeriodicIndex() && 1 ) {
 856               seg.epgeominfo[ 0 ].u = otherSeamParam;
 857               seg.epgeominfo[ 1 ].u = otherSeamParam;
 858               swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
 859             } else {
 860               seg.epgeominfo[ 0 ].v = otherSeamParam;
 861               seg.epgeominfo[ 1 ].v = otherSeamParam;
 862               swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
 863             }
 864             swap (seg[0], seg[1]);
 865             swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
 866             seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
 867             ngMesh.AddSegment (seg);
 868 #ifdef DUMP_SEGMENTS
 869             cout << "Segment: " << seg.edgenr << endl
 870                  << "\t is SEAM (reverse) of the previous. "
 871                  << " Other " << (helper.GetPeriodicIndex() && 1 ? "U" : "V")
 872                  << " = " << otherSeamParam << endl;
 873 #endif
 874           }
 875           else if ( fOri == TopAbs_INTERNAL )
 876           {
 877             swap (seg[0], seg[1]);
 878             swap( seg.epgeominfo[0], seg.epgeominfo[1] );
 879             seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
 880             ngMesh.AddSegment (seg);
 881 #ifdef DUMP_SEGMENTS
 882             cout << "Segment: " << seg.edgenr << endl << "\t is REVERSE of the previous" << endl;
 883 #endif
 884           }
 885         }
 886       } // loop on geomEdge ancestors
 887
 888       if ( quadHelper ) // remember medium nodes of sub-meshes
 889       {
 890         SMDS_ElemIteratorPtr edges = smDS->GetElements();
 891         while ( edges->more() )
 892         {
 893           const SMDS_MeshElement* e = edges->next();
 894           if ( !quadHelper->AddTLinks( static_cast< const SMDS_MeshEdge*>( e )))
 895             break;
 896         }
 897       }
 898
 899       break;
 900     } // case TopAbs_EDGE
 901
 902     case TopAbs_FACE: { // FACE
 903       // ----------------------
 904       const TopoDS_Face& geomFace  = TopoDS::Face( sm->GetSubShape() );
 905       helper.SetSubShape( geomFace );
 906       bool isInternalFace = ( geomFace.Orientation() == TopAbs_INTERNAL );
 907
 908       // Find solids the geomFace bounds
 909       int solidID1 = 0, solidID2 = 0;
 910       StdMeshers_QuadToTriaAdaptor* quadAdaptor =
 911         dynamic_cast<StdMeshers_QuadToTriaAdaptor*>( proxyMesh.get() );
 912       if ( quadAdaptor )
 913       {
 914         solidID1 = occgeom.somap.FindIndex( quadAdaptor->GetShape() );
 915       }
 916       else
 917       {
 918         PShapeIteratorPtr solidIt = helper.GetAncestors( geomFace, *sm->GetFather(), TopAbs_SOLID);
 919         while ( const TopoDS_Shape * solid = solidIt->next() )
 920         {
 921           int id = occgeom.somap.FindIndex ( *solid );
 922           if ( solidID1 && id != solidID1 ) solidID2 = id;
 923           else                              solidID1 = id;
 924         }
 925       }
 926       // Add ng face descriptors of meshed faces
 927       faceNgID++;
 928       ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(faceNgID, solidID1, solidID2, 0));
 929
 930       // if second oreder is required, even already meshed faces must be passed to NETGEN
 931       int fID = occgeom.fmap.Add( geomFace );
 932       while ( fID < faceNgID ) // geomFace is already in occgeom.fmap, add a copy
 933         fID = occgeom.fmap.Add( BRepBuilderAPI_Copy( geomFace, /*copyGeom=*/false ));
 934       // Problem with the second order in a quadrangular mesh remains.
 935       // 1) All quadrangles generated by NETGEN are moved to an inexistent face
 936       //    by FillSMesh() (find "AddFaceDescriptor")
 937       // 2) Temporary triangles generated by StdMeshers_QuadToTriaAdaptor
 938       //    are on faces where quadrangles were.
 939       // Due to these 2 points, wrong geom faces are used while conversion to qudratic
 940       // of the mentioned above quadrangles and triangles
 941
 942       // Orient the face correctly in solidID1 (issue 0020206)
 943       bool reverse = false;
 944       if ( solidID1 ) {
 945         TopoDS_Shape solid = occgeom.somap( solidID1 );
 946         TopAbs_Orientation faceOriInSolid = helper.GetSubShapeOri( solid, geomFace );
 947         if ( faceOriInSolid >= 0 )
 948           reverse =
 949             helper.IsReversedSubMesh( TopoDS::Face( geomFace.Oriented( faceOriInSolid )));
 950       }
 951
 952       // Add surface elements
 953
 954       netgen::Element2d tri(3);
 955       tri.SetIndex ( faceNgID );
 956
 957
 958 #ifdef DUMP_TRIANGLES
 959       cout << "SMESH face " << helper.GetMeshDS()->ShapeToIndex( geomFace )
 960            << " internal="<<isInternalFace << endl;
 961 #endif
 962       if ( proxyMesh )
 963         smDS = proxyMesh->GetSubMesh( geomFace );
 964
 965       SMDS_ElemIteratorPtr faces = smDS->GetElements();
 966       while ( faces->more() )
 967       {
 968         const SMDS_MeshElement* f = faces->next();
 969         if ( f->NbNodes() % 3 != 0 ) // not triangle
 970         {
 971           PShapeIteratorPtr solidIt=helper.GetAncestors(geomFace,*sm->GetFather(),TopAbs_SOLID);
 972           if ( const TopoDS_Shape * solid = solidIt->next() )
 973             sm = _mesh->GetSubMesh( *solid );
 974           SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 975           smError.reset( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"Not triangle submesh"));
 976           smError->myBadElements.push_back( f );
 977           return false;
 978         }
 979
 980         for ( int i = 0; i < 3; ++i )
 981         {
 982           const SMDS_MeshNode* node = f->GetNode( i ), * inFaceNode=0;
 983
 984           // get node UV on face
 985           int shapeID = node->getshapeId();
 986           if ( helper.IsSeamShape( shapeID ))
 987             if ( helper.IsSeamShape( f->GetNodeWrap( i+1 )->getshapeId() ))
 988               inFaceNode = f->GetNodeWrap( i-1 );
 989             else
 990               inFaceNode = f->GetNodeWrap( i+1 );
 991           gp_XY uv = helper.GetNodeUV( geomFace, node, inFaceNode );
 992
 993           int ind = reverse ? 3-i : i+1;
 994           tri.GeomInfoPi(ind).u = uv.X();
 995           tri.GeomInfoPi(ind).v = uv.Y();
 996           tri.PNum      (ind) = ngNodeId( node, ngMesh, nodeNgIdMap );
 997         }
 998
 999         ngMesh.AddSurfaceElement (tri);
1000 #ifdef DUMP_TRIANGLES
1001         cout << tri << endl;
1002 #endif
1003
1004         if ( isInternalFace )
1005         {
1006           swap( tri[1], tri[2] );
1007           ngMesh.AddSurfaceElement (tri);
1008 #ifdef DUMP_TRIANGLES
1009           cout << tri << endl;
1010 #endif
1011         }
1012       }
1013
1014       if ( quadHelper ) // remember medium nodes of sub-meshes
1015       {
1016         SMDS_ElemIteratorPtr faces = smDS->GetElements();
1017         while ( faces->more() )
1018         {
1019           const SMDS_MeshElement* f = faces->next();
1020           if ( !quadHelper->AddTLinks( static_cast< const SMDS_MeshFace*>( f )))
1021             break;
1022         }
1023       }
1024
1025       break;
1026     } // case TopAbs_FACE
1027
1028     case TopAbs_VERTEX: { // VERTEX
1029       // --------------------------
1030       // issue 0021405. Add node only if a VERTEX is shared by a not meshed EDGE,
1031       // else netgen removes a free node and nodeVector becomes invalid
1032       PShapeIteratorPtr ansIt = helper.GetAncestors( sm->GetSubShape(),
1033                                                      *sm->GetFather(),
1034                                                      TopAbs_EDGE );
1035       bool toAdd = false;
1036       while ( const TopoDS_Shape* e = ansIt->next() )
1037       {
1038         SMESH_subMesh* eSub = helper.GetMesh()->GetSubMesh( *e );
1039         if (( toAdd = eSub->IsEmpty() )) break;
1040       }
1041       if ( toAdd )
1042       {
1043         SMDS_NodeIteratorPtr nodeIt = smDS->GetNodes();
1044         if ( nodeIt->more() )
1045           ngNodeId( nodeIt->next(), ngMesh, nodeNgIdMap );
1046       }
1047       break;
1048     }
1049     default:;
1050     } // switch
1051   } // loop on submeshes
1052
1053   // fill nodeVec
1054   nodeVec.resize( ngMesh.GetNP() + 1 );
1055   TNode2IdMap::iterator node_NgId, nodeNgIdEnd = nodeNgIdMap.end();
1056   for ( node_NgId = nodeNgIdMap.begin(); node_NgId != nodeNgIdEnd; ++node_NgId)
1057     nodeVec[ node_NgId->second ] = node_NgId->first;
1058
1059   return true;
1060 }
1061
1062 //================================================================================
1063 /*!
1064  * \brief Duplicate mesh faces on internal geom faces
1065  */
1066 //================================================================================
1067
1068 void NETGENPlugin_Mesher::FixIntFaces(const netgen::OCCGeometry& occgeom,
1069                                       netgen::Mesh&              ngMesh,
1070                                       NETGENPlugin_Internals&    internalShapes)
1071 {
1072   SMESHDS_Mesh* meshDS = internalShapes.getMesh().GetMeshDS();
1073
1074   // find ng indices of internal faces
1075   set<int> ngFaceIds;
1076   for ( int ngFaceID = 1; ngFaceID <= occgeom.fmap.Extent(); ++ngFaceID )
1077   {
1078     int smeshID = meshDS->ShapeToIndex( occgeom.fmap( ngFaceID ));
1079     if ( internalShapes.isInternalShape( smeshID ))
1080       ngFaceIds.insert( ngFaceID );
1081   }
1082   if ( !ngFaceIds.empty() )
1083   {
1084     // duplicate faces
1085     int i, nbFaces = ngMesh.GetNSE();
1086     for (int i = 1; i <= nbFaces; ++i)
1087     {
1088       netgen::Element2d elem = ngMesh.SurfaceElement(i);
1089       if ( ngFaceIds.count( elem.GetIndex() ))
1090       {
1091         swap( elem[1], elem[2] );
1092         ngMesh.AddSurfaceElement (elem);
1093       }
1094     }
1095   }
1096 }
1097
1098 namespace
1099 {
1100   //================================================================================
1101   // define gp_XY_Subtracted pointer to function calling gp_XY::Subtracted(gp_XY)
1102   gp_XY_FunPtr(Subtracted);
1103   //gp_XY_FunPtr(Added);
1104
1105   //================================================================================
1106   /*!
1107    * \brief Evaluate distance between two 2d points along the surface
1108    */
1109   //================================================================================
1110
1111   double evalDist( const gp_XY&                uv1,
1112                    const gp_XY&                uv2,
1113                    const Handle(Geom_Surface)& surf,
1114                    const int                   stopHandler=-1)
1115   {
1116     if ( stopHandler > 0 ) // continue recursion
1117     {
1118       gp_XY mid = SMESH_MesherHelper::GetMiddleUV( surf, uv1, uv2 );
1119       return evalDist( uv1,mid, surf, stopHandler-1 ) + evalDist( mid,uv2, surf, stopHandler-1 );
1120     }
1121     double dist3D = surf->Value( uv1.X(), uv1.Y() ).Distance( surf->Value( uv2.X(), uv2.Y() ));
1122     if ( stopHandler == 0 ) // stop recursion
1123       return dist3D;
1124
1125     // start recursion if necessary
1126     double dist2D = SMESH_MesherHelper::applyIn2D(surf, uv1, uv2, gp_XY_Subtracted, 0).Modulus();
1127     if ( fabs( dist3D - dist2D ) < dist2D * 1e-10 )
1128       return dist3D; // equal parametrization of a planar surface
1129
1130     return evalDist( uv1, uv2, surf, 3 ); // start recursion
1131   }
1132
1133   //================================================================================
1134   /*!
1135    * \brief Data of vertex internal in geom face
1136    */
1137   //================================================================================
1138
1139   struct TIntVData
1140   {
1141     gp_XY uv;        //!< UV in face parametric space
1142     int   ngId;      //!< ng id of corrsponding node
1143     gp_XY uvClose;   //!< UV of closest boundary node
1144     int   ngIdClose; //!< ng id of closest boundary node
1145   };
1146
1147   //================================================================================
1148   /*!
1149    * \brief Data of vertex internal in solid
1150    */
1151   //================================================================================
1152
1153   struct TIntVSoData
1154   {
1155     int   ngId;      //!< ng id of corresponding node
1156     int   ngIdClose; //!< ng id of closest 2d mesh element
1157     int   ngIdCloseN; //!< ng id of closest node of the closest 2d mesh element
1158   };
1159
1160   inline double dist2(const netgen::MeshPoint& p1, const netgen::MeshPoint& p2)
1161   {
1162     return gp_Pnt( NGPOINT_COORDS(p1)).SquareDistance( gp_Pnt( NGPOINT_COORDS(p2)));
1163   }
1164 }
1165
1166 //================================================================================
1167 /*!
1168  * \brief Make netgen take internal vertices in faces into account by adding
1169  *        segments including internal vertices
1170  *
1171  * This function works in supposition that 1D mesh is already computed in ngMesh
1172  */
1173 //================================================================================
1174
1175 void NETGENPlugin_Mesher::AddIntVerticesInFaces(const netgen::OCCGeometry&     occgeom,
1176                                                 netgen::Mesh&                  ngMesh,
1177                                                 vector<const SMDS_MeshNode*>&  nodeVec,
1178                                                 NETGENPlugin_Internals&        internalShapes)
1179 {
1180   if ( nodeVec.size() < ngMesh.GetNP() )
1181     nodeVec.resize( ngMesh.GetNP(), 0 );
1182
1183   SMESHDS_Mesh* meshDS = internalShapes.getMesh().GetMeshDS();
1184   SMESH_MesherHelper helper( internalShapes.getMesh() );
1185
1186   const map<int,list<int> >& face2Vert = internalShapes.getFacesWithVertices();
1187   map<int,list<int> >::const_iterator f2v = face2Vert.begin();
1188   for ( ; f2v != face2Vert.end(); ++f2v )
1189   {
1190     const TopoDS_Face& face = TopoDS::Face( meshDS->IndexToShape( f2v->first ));
1191     if ( face.IsNull() ) continue;
1192     int faceNgID = occgeom.fmap.FindIndex (face);
1193     if ( faceNgID < 0 ) continue;
1194
1195     TopLoc_Location loc;
1196     Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc);
1197
1198     helper.SetSubShape( face );
1199     helper.SetElementsOnShape( true );
1200
1201     // Get data of internal vertices and add them to ngMesh
1202
1203     multimap< double, TIntVData > dist2VData; // sort vertices by distance from boundary nodes
1204
1205     int i, nbSegInit = ngMesh.GetNSeg();
1206
1207     // boundary characteristics
1208     double totSegLen2D = 0;
1209     int totNbSeg = 0;
1210
1211     const list<int>& iVertices = f2v->second;
1212     list<int>::const_iterator iv = iVertices.begin();
1213     for ( int nbV = 0; iv != iVertices.end(); ++iv, nbV++ )
1214     {
1215       TIntVData vData;
1216       // get node on vertex
1217       const TopoDS_Vertex V = TopoDS::Vertex( meshDS->IndexToShape( *iv ));
1218       const SMDS_MeshNode * nV = SMESH_Algo::VertexNode( V, meshDS );
1219       if ( !nV )
1220       {
1221         SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( V );
1222         sm->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1223         nV = SMESH_Algo::VertexNode( V, meshDS );
1224         if ( !nV ) continue;
1225       }
1226       // add ng node
1227       netgen::MeshPoint mp( netgen::Point<3> (nV->X(), nV->Y(), nV->Z()) );
1228       ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
1229       vData.ngId = ngMesh.GetNP();
1230       nodeVec.push_back( nV );
1231
1232       // get node UV
1233       bool uvOK = false;
1234       vData.uv = helper.GetNodeUV( face, nV, 0, &uvOK );
1235       if ( !uvOK ) helper.CheckNodeUV( face, nV, vData.uv, BRep_Tool::Tolerance(V),/*force=*/1);
1236
1237       // loop on all segments of the face to find the node closest to vertex and to count
1238       // average segment 2d length
1239       double closeDist2 = numeric_limits<double>::max(), dist2;
1240       int ngIdLast = 0;
1241       for (i = 1; i <= ngMesh.GetNSeg(); ++i)
1242       {
1243         netgen::Segment & seg = ngMesh.LineSegment(i);
1244         if ( seg.si != faceNgID ) continue;
1245         gp_XY uv[2];
1246         for ( int iEnd = 0; iEnd < 2; ++iEnd)
1247         {
1248           uv[iEnd].SetCoord( seg.epgeominfo[iEnd].u, seg.epgeominfo[iEnd].v );
1249           if ( ngIdLast == seg[ iEnd ] ) continue;
1250           dist2 = helper.applyIn2D(surf, uv[iEnd], vData.uv, gp_XY_Subtracted,0).SquareModulus();
1251           if ( dist2 < closeDist2 )
1252             vData.ngIdClose = seg[ iEnd ], vData.uvClose = uv[iEnd], closeDist2 = dist2;
1253           ngIdLast = seg[ iEnd ];
1254         }
1255         if ( !nbV )
1256         {
1257           totSegLen2D += helper.applyIn2D(surf, uv[0], uv[1], gp_XY_Subtracted, false).Modulus();
1258           totNbSeg++;
1259         }
1260       }
1261       dist2VData.insert( make_pair( closeDist2, vData ));
1262     }
1263
1264     if ( totNbSeg == 0 ) break;
1265     double avgSegLen2d = totSegLen2D / totNbSeg;
1266
1267     // Loop on vertices to add segments
1268
1269     multimap< double, TIntVData >::iterator dist_vData = dist2VData.begin();
1270     for ( ; dist_vData != dist2VData.end(); ++dist_vData )
1271     {
1272       double closeDist2 = dist_vData->first, dist2;
1273       TIntVData & vData = dist_vData->second;
1274
1275       // try to find more close node among segments added for internal vertices
1276       for (i = nbSegInit+1; i <= ngMesh.GetNSeg(); ++i)
1277       {
1278         netgen::Segment & seg = ngMesh.LineSegment(i);
1279         if ( seg.si != faceNgID ) continue;
1280         gp_XY uv[2];
1281         for ( int iEnd = 0; iEnd < 2; ++iEnd)
1282         {
1283           uv[iEnd].SetCoord( seg.epgeominfo[iEnd].u, seg.epgeominfo[iEnd].v );
1284           dist2 = helper.applyIn2D(surf, uv[iEnd], vData.uv, gp_XY_Subtracted,0).SquareModulus();
1285           if ( dist2 < closeDist2 )
1286             vData.ngIdClose = seg[ iEnd ], vData.uvClose = uv[iEnd], closeDist2 = dist2;
1287         }
1288       }
1289       // decide whether to use the closest node as the second end of segment or to
1290       // create a new point
1291       int segEnd1 = vData.ngId;
1292       int segEnd2 = vData.ngIdClose; // to use closest node
1293       gp_XY uvV = vData.uv, uvP = vData.uvClose;
1294       double segLenHint  = ngMesh.GetH( ngMesh.Point( vData.ngId ));
1295       double nodeDist2D  = sqrt( closeDist2 );
1296       double nodeDist3D  = evalDist( vData.uv, vData.uvClose, surf );
1297       bool avgLenOK  = ( avgSegLen2d < 0.75 * nodeDist2D );
1298       bool hintLenOK = ( segLenHint  < 0.75 * nodeDist3D );
1299       //cout << "uvV " << uvV.X() <<","<<uvV.Y() << " ";
1300       if ( hintLenOK || avgLenOK )
1301       {
1302         // create a point between the closest node and V
1303
1304         // how far from V
1305         double r = min( 0.5, ( hintLenOK ? segLenHint/nodeDist3D : avgSegLen2d/nodeDist2D ));
1306         // direction from V to closet node in 2D
1307         gp_Dir2d v2n( helper.applyIn2D(surf, uvP, uvV, gp_XY_Subtracted, false ));
1308         // new point
1309         uvP = vData.uv + r * nodeDist2D * v2n.XY();
1310         gp_Pnt P = surf->Value( uvP.X(), uvP.Y() ).Transformed( loc );
1311
1312         netgen::MeshPoint mp( netgen::Point<3> (P.X(), P.Y(), P.Z()));
1313         ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
1314         segEnd2 = ngMesh.GetNP();
1315         //cout << "Middle " << r << " uv " << uvP.X() << "," << uvP.Y() << "( " << ngMesh.Point(segEnd2).X()<<","<<ngMesh.Point(segEnd2).Y()<<","<<ngMesh.Point(segEnd2).Z()<<" )"<< endl;
1316         SMDS_MeshNode * nP = helper.AddNode(P.X(), P.Y(), P.Z());
1317         nodeVec.push_back( nP );
1318       }
1319       //else cout << "at Node " << " uv " << uvP.X() << "," << uvP.Y() << endl;
1320
1321       // Add the segment
1322       netgen::Segment seg;
1323
1324       if ( segEnd1 > segEnd2 ) swap( segEnd1, segEnd2 ), swap( uvV, uvP );
1325       seg[0] = segEnd1;  // ng node id
1326       seg[1] = segEnd2;  // ng node id
1327       seg.edgenr = ngMesh.GetNSeg() + 1;// segment id
1328       seg.si = faceNgID;
1329
1330       seg.epgeominfo[ 0 ].dist = 0; // param on curve
1331       seg.epgeominfo[ 0 ].u    = uvV.X();
1332       seg.epgeominfo[ 0 ].v    = uvV.Y();
1333       seg.epgeominfo[ 1 ].dist = 1; // param on curve
1334       seg.epgeominfo[ 1 ].u    = uvP.X();
1335       seg.epgeominfo[ 1 ].v    = uvP.Y();
1336
1337 //       seg.epgeominfo[ 0 ].edgenr = 10; //  = geom.emap.FindIndex(edge);
1338 //       seg.epgeominfo[ 1 ].edgenr = 10; //  = geom.emap.FindIndex(edge);
1339
1340       ngMesh.AddSegment (seg);
1341
1342       // add reverse segment
1343       swap (seg[0], seg[1]);
1344       swap( seg.epgeominfo[0], seg.epgeominfo[1] );
1345       seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
1346       ngMesh.AddSegment (seg);
1347     }
1348
1349   }
1350 }
1351
1352 //================================================================================
1353 /*!
1354  * \brief Make netgen take internal vertices in solids into account by adding
1355  *        faces including internal vertices
1356  *
1357  * This function works in supposition that 2D mesh is already computed in ngMesh
1358  */
1359 //================================================================================
1360
1361 void NETGENPlugin_Mesher::AddIntVerticesInSolids(const netgen::OCCGeometry&     occgeom,
1362                                                  netgen::Mesh&                  ngMesh,
1363                                                  vector<const SMDS_MeshNode*>&  nodeVec,
1364                                                  NETGENPlugin_Internals&        internalShapes)
1365 {
1366 #ifdef DUMP_TRIANGLES_SCRIPT
1367   // create a python script making a mesh containing triangles added for internal vertices
1368   ofstream py(DUMP_TRIANGLES_SCRIPT);
1369   py << "import SMESH"<< endl
1370      << "from salome.smesh import smeshBuilder"<<endl
1371      << "smesh = smeshBuilder.New(salome.myStudy)"
1372      << "m = smesh.Mesh(name='triangles')" << endl;
1373 #endif
1374   if ( nodeVec.size() < ngMesh.GetNP() )
1375     nodeVec.resize( ngMesh.GetNP(), 0 );
1376
1377   SMESHDS_Mesh* meshDS = internalShapes.getMesh().GetMeshDS();
1378   SMESH_MesherHelper helper( internalShapes.getMesh() );
1379
1380   const map<int,list<int> >& so2Vert = internalShapes.getSolidsWithVertices();
1381   map<int,list<int> >::const_iterator s2v = so2Vert.begin();
1382   for ( ; s2v != so2Vert.end(); ++s2v )
1383   {
1384     const TopoDS_Shape& solid = meshDS->IndexToShape( s2v->first );
1385     if ( solid.IsNull() ) continue;
1386     int solidNgID = occgeom.somap.FindIndex (solid);
1387     if ( solidNgID < 0 && !occgeom.somap.IsEmpty() ) continue;
1388
1389     helper.SetSubShape( solid );
1390     helper.SetElementsOnShape( true );
1391
1392     // find ng indices of faces within the solid
1393     set<int> ngFaceIds;
1394     for (TopExp_Explorer fExp(solid, TopAbs_FACE); fExp.More(); fExp.Next() )
1395       ngFaceIds.insert( occgeom.fmap.FindIndex( fExp.Current() ));
1396     if ( ngFaceIds.size() == 1 && *ngFaceIds.begin() == 0 )
1397       ngFaceIds.insert( 1 );
1398
1399     // Get data of internal vertices and add them to ngMesh
1400
1401     multimap< double, TIntVSoData > dist2VData; // sort vertices by distance from ng faces
1402
1403     int i, nbFaceInit = ngMesh.GetNSE();
1404
1405     // boundary characteristics
1406     double totSegLen = 0;
1407     int totNbSeg = 0;
1408
1409     const list<int>& iVertices = s2v->second;
1410     list<int>::const_iterator iv = iVertices.begin();
1411     for ( int nbV = 0; iv != iVertices.end(); ++iv, nbV++ )
1412     {
1413       TIntVSoData vData;
1414       const TopoDS_Vertex V = TopoDS::Vertex( meshDS->IndexToShape( *iv ));
1415
1416       // get node on vertex
1417       const SMDS_MeshNode * nV = SMESH_Algo::VertexNode( V, meshDS );
1418       if ( !nV )
1419       {
1420         SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( V );
1421         sm->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1422         nV = SMESH_Algo::VertexNode( V, meshDS );
1423         if ( !nV ) continue;
1424       }
1425       // add ng node
1426       netgen::MeshPoint mpV( netgen::Point<3> (nV->X(), nV->Y(), nV->Z()) );
1427       ngMesh.AddPoint ( mpV, 1, netgen::FIXEDPOINT );
1428       vData.ngId = ngMesh.GetNP();
1429       nodeVec.push_back( nV );
1430
1431       // loop on all 2d elements to find the one closest to vertex and to count
1432       // average segment length
1433       double closeDist2 = numeric_limits<double>::max(), avgDist2;
1434       for (i = 1; i <= ngMesh.GetNSE(); ++i)
1435       {
1436         const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
1437         if ( !ngFaceIds.count( elem.GetIndex() )) continue;
1438         avgDist2 = 0;
1439         multimap< double, int> dist2nID; // sort nodes of element by distance from V
1440         for ( int j = 0; j < elem.GetNP(); ++j)
1441         {
1442           netgen::MeshPoint mp = ngMesh.Point( elem[j] );
1443           double d2 = dist2( mpV, mp );
1444           dist2nID.insert( make_pair( d2, elem[j] ));
1445           avgDist2 += d2 / elem.GetNP();
1446           if ( !nbV )
1447             totNbSeg++, totSegLen+= sqrt( dist2( mp, ngMesh.Point( elem[(j+1)%elem.GetNP()])));
1448         }
1449         double dist = dist2nID.begin()->first; //avgDist2;
1450         if ( dist < closeDist2 )
1451           vData.ngIdClose= i, vData.ngIdCloseN= dist2nID.begin()->second, closeDist2= dist;
1452       }
1453       dist2VData.insert( make_pair( closeDist2, vData ));
1454     }
1455
1456     if ( totNbSeg == 0 ) break;
1457     double avgSegLen = totSegLen / totNbSeg;
1458
1459     // Loop on vertices to add triangles
1460
1461     multimap< double, TIntVSoData >::iterator dist_vData = dist2VData.begin();
1462     for ( ; dist_vData != dist2VData.end(); ++dist_vData )
1463     {
1464       double closeDist2   = dist_vData->first;
1465       TIntVSoData & vData = dist_vData->second;
1466
1467       const netgen::MeshPoint& mpV = ngMesh.Point( vData.ngId );
1468
1469       // try to find more close face among ones added for internal vertices
1470       for (i = nbFaceInit+1; i <= ngMesh.GetNSE(); ++i)
1471       {
1472         double avgDist2 = 0;
1473         multimap< double, int> dist2nID;
1474         const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
1475         for ( int j = 0; j < elem.GetNP(); ++j)
1476         {
1477           double d = dist2( mpV, ngMesh.Point( elem[j] ));
1478           dist2nID.insert( make_pair( d, elem[j] ));
1479           avgDist2 += d / elem.GetNP();
1480           if ( avgDist2 < closeDist2 )
1481             vData.ngIdClose= i, vData.ngIdCloseN= dist2nID.begin()->second, closeDist2= avgDist2;
1482         }
1483       }
1484       // sort nodes of the closest face by angle with vector from V to the closest node
1485       const double tol = numeric_limits<double>::min();
1486       map< double, int > angle2ID;
1487       const netgen::Element2d& closeFace = ngMesh.SurfaceElement( vData.ngIdClose );
1488       netgen::MeshPoint mp[2];
1489       mp[0] = ngMesh.Point( vData.ngIdCloseN );
1490       gp_XYZ p1( NGPOINT_COORDS( mp[0] ));
1491       gp_XYZ pV( NGPOINT_COORDS( mpV ));
1492       gp_Vec v2p1( pV, p1 );
1493       double distN1 = v2p1.Magnitude();
1494       if ( distN1 <= tol ) continue;
1495       v2p1 /= distN1;
1496       for ( int j = 0; j < closeFace.GetNP(); ++j)
1497       {
1498         mp[1] = ngMesh.Point( closeFace[j] );
1499         gp_Vec v2p( pV, gp_Pnt( NGPOINT_COORDS( mp[1] )) );
1500         angle2ID.insert( make_pair( v2p1.Angle( v2p ), closeFace[j]));
1501       }
1502       // get node with angle of 60 degrees or greater
1503       map< double, int >::iterator angle_id = angle2ID.lower_bound( 60. * M_PI / 180. );
1504       if ( angle_id == angle2ID.end() ) angle_id = --angle2ID.end();
1505       const double minAngle = 30. * M_PI / 180.;
1506       const double angle = angle_id->first;
1507       bool angleOK = ( angle > minAngle );
1508
1509       // find points to create a triangle
1510       netgen::Element2d tri(3);
1511       tri.SetIndex ( 1 );
1512       tri[0] = vData.ngId;
1513       tri[1] = vData.ngIdCloseN; // to use the closest nodes
1514       tri[2] = angle_id->second; // to use the node with best angle
1515
1516       // decide whether to use the closest node and the node with best angle or to create new ones
1517       for ( int isBestAngleN = 0; isBestAngleN < 2; ++isBestAngleN )
1518       {
1519         bool createNew = !angleOK, distOK = true;
1520         double distFromV;
1521         int triInd = isBestAngleN ? 2 : 1;
1522         mp[isBestAngleN] = ngMesh.Point( tri[triInd] );
1523         if ( isBestAngleN )
1524         {
1525           if ( angleOK )
1526           {
1527             double distN2 = sqrt( dist2( mpV, mp[isBestAngleN]));
1528             createNew = ( fabs( distN2 - distN1 ) > 0.25 * distN1 );
1529           }
1530           else if ( angle < tol )
1531           {
1532             v2p1.SetX( v2p1.X() + 1e-3 );
1533           }
1534           distFromV = distN1;
1535         }
1536         else
1537         {
1538           double segLenHint = ngMesh.GetH( ngMesh.Point( vData.ngId ));
1539           bool avgLenOK  = ( avgSegLen < 0.75 * distN1 );
1540           bool hintLenOK = ( segLenHint  < 0.75 * distN1 );
1541           createNew = (createNew || avgLenOK || hintLenOK );
1542           // we create a new node not closer than 0.5 to the closest face
1543           // in order not to clash with other close face
1544           double r = min( 0.5, ( hintLenOK ? segLenHint : avgSegLen ) / distN1 );
1545           distFromV = r * distN1;
1546         }
1547         if ( createNew )
1548         {
1549           // create a new point, between the node and the vertex if angleOK
1550           gp_XYZ p( NGPOINT_COORDS( mp[isBestAngleN] ));
1551           gp_Vec v2p( pV, p ); v2p.Normalize();
1552           if ( isBestAngleN && !angleOK )
1553             p = p1 + gp_Dir( v2p.XYZ() - v2p1.XYZ()).XYZ() * distN1 * 0.95;
1554           else
1555             p = pV + v2p.XYZ() * distFromV;
1556
1557           if ( !isBestAngleN ) p1 = p, distN1 = distFromV;
1558
1559           mp[isBestAngleN].SetPoint( netgen::Point<3> (p.X(), p.Y(), p.Z()));
1560           ngMesh.AddPoint ( mp[isBestAngleN], 1, netgen::SURFACEPOINT );
1561           tri[triInd] = ngMesh.GetNP();
1562           nodeVec.push_back( helper.AddNode( p.X(), p.Y(), p.Z()) );
1563         }
1564       }
1565       ngMesh.AddSurfaceElement (tri);
1566       swap( tri[1], tri[2] );
1567       ngMesh.AddSurfaceElement (tri);
1568
1569 #ifdef DUMP_TRIANGLES_SCRIPT
1570       py << "n1 = m.AddNode( "<< mpV.X()<<", "<< mpV.Y()<<", "<< mpV.Z()<<") "<< endl
1571          << "n2 = m.AddNode( "<< mp[0].X()<<", "<< mp[0].Y()<<", "<< mp[0].Z()<<") "<< endl
1572          << "n3 = m.AddNode( "<< mp[1].X()<<", "<< mp[1].Y()<<", "<< mp[1].Z()<<" )" << endl
1573          << "m.AddFace([n1,n2,n3])" << endl;
1574 #endif
1575     } // loop on internal vertices of a solid
1576
1577   } // loop on solids with internal vertices
1578 }
1579
1580 //================================================================================
1581 /*!
1582  * \brief Fill netgen mesh with segments of a FACE
1583  *  \param ngMesh - netgen mesh
1584  *  \param geom - container of OCCT geometry to mesh
1585  *  \param wires - data of nodes on FACE boundary
1586  *  \param helper - mesher helper holding the FACE
1587  *  \param nodeVec - vector of nodes in which node index == netgen ID
1588  *  \retval SMESH_ComputeErrorPtr - error description
1589  */
1590 //================================================================================
1591
1592 SMESH_ComputeErrorPtr
1593 NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
1594                                        netgen::OCCGeometry&             geom,
1595                                        const TSideVector&               wires,
1596                                        SMESH_MesherHelper&              helper,
1597                                        vector< const SMDS_MeshNode* > & nodeVec)
1598 {
1599   // ----------------------------
1600   // Check wires and count nodes
1601   // ----------------------------
1602   int nbNodes = 0;
1603   for ( int iW = 0; iW < wires.size(); ++iW )
1604   {
1605     StdMeshers_FaceSidePtr wire = wires[ iW ];
1606     if ( wire->MissVertexNode() )
1607     {
1608       // Commented for issue 0020960. It worked for the case, let's wait for case where it doesn't.
1609       // It seems that there is no reason for this limitation
1610 //       return TError
1611 //         (new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH, "Missing nodes on vertices"));
1612     }
1613     const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
1614     if ( uvPtVec.size() != wire->NbPoints() )
1615       return SMESH_ComputeError::New(COMPERR_BAD_INPUT_MESH,
1616                                      SMESH_Comment("Unexpected nb of points on wire ") << iW
1617                                      << ": " << uvPtVec.size()<<" != "<<wire->NbPoints());
1618     nbNodes += wire->NbPoints();
1619   }
1620   nodeVec.reserve( nodeVec.size() + nbNodes + 1 );
1621   if ( nodeVec.empty() )
1622     nodeVec.push_back( 0 );
1623
1624   // -----------------
1625   // Fill netgen mesh
1626   // -----------------
1627
1628   const bool wasNgMeshEmpty = ( ngMesh.GetNP() < 1 ); /* true => this method is called by
1629                                                          NETGENPlugin_NETGEN_2D_ONLY */
1630
1631   // map for nodes on vertices since they can be shared between wires
1632   // ( issue 0020676, face_int_box.brep) and nodes built by NETGEN
1633   map<const SMDS_MeshNode*, int > node2ngID;
1634   if ( !wasNgMeshEmpty ) // fill node2ngID with nodes built by NETGEN
1635   {
1636     set< int > subIDs; // ids of sub-shapes of the FACE
1637     for ( int iW = 0; iW < wires.size(); ++iW )
1638     {
1639       StdMeshers_FaceSidePtr wire = wires[ iW ];
1640       for ( int iE = 0, nbE = wire->NbEdges(); iE < nbE; ++iE )
1641       {
1642         subIDs.insert( wire->EdgeID( iE ));
1643         subIDs.insert( helper.GetMeshDS()->ShapeToIndex( wire->FirstVertex( iE )));
1644       }
1645     }
1646     for ( size_t ngID = 1; ngID < nodeVec.size(); ++ngID )
1647       if ( subIDs.count( nodeVec[ngID]->getshapeId() ))
1648         node2ngID.insert( make_pair( nodeVec[ngID], ngID ));
1649   }
1650
1651   const int solidID = 0, faceID = geom.fmap.FindIndex( helper.GetSubShape() );
1652   if ( ngMesh.GetNFD() < 1 )
1653     ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(faceID, solidID, solidID, 0));
1654
1655   for ( int iW = 0; iW < wires.size(); ++iW )
1656   {
1657     StdMeshers_FaceSidePtr wire = wires[ iW ];
1658     const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
1659     const int nbSegments = wire->NbPoints() - 1;
1660
1661     // assure the 1st node to be in node2ngID, which is needed to correctly
1662     // "close chain of segments" (see below) in case if the 1st node is not
1663     // onVertex because it is on a Viscous layer
1664     node2ngID.insert( make_pair( uvPtVec[ 0 ].node, ngMesh.GetNP() + 1 ));
1665
1666     // compute length of every segment
1667     vector<double> segLen( nbSegments );
1668     for ( int i = 0; i < nbSegments; ++i )
1669       segLen[i] = SMESH_TNodeXYZ( uvPtVec[ i ].node ).Distance( uvPtVec[ i+1 ].node );
1670
1671     int edgeID = 1, posID = -2;
1672     bool isInternalWire = false;
1673     double vertexNormPar = 0;
1674     const int prevNbNGSeg = ngMesh.GetNSeg();
1675     for ( int i = 0; i < nbSegments; ++i ) // loop on segments
1676     {
1677       // Add the first point of a segment
1678
1679       const SMDS_MeshNode * n = uvPtVec[ i ].node;
1680       const int posShapeID = n->getshapeId();
1681       bool onVertex = ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX );
1682       bool onEdge   = ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE   );
1683
1684       // skip nodes on degenerated edges
1685       if ( helper.IsDegenShape( posShapeID ) &&
1686            helper.IsDegenShape( uvPtVec[ i+1 ].node->getshapeId() ))
1687         continue;
1688
1689       int ngID1 = ngMesh.GetNP() + 1, ngID2 = ngID1+1;
1690       if ( onVertex || ( !wasNgMeshEmpty && onEdge ))
1691         ngID1 = node2ngID.insert( make_pair( n, ngID1 )).first->second;
1692       if ( ngID1 > ngMesh.GetNP() )
1693       {
1694         netgen::MeshPoint mp( netgen::Point<3> (n->X(), n->Y(), n->Z()) );
1695         ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
1696         nodeVec.push_back( n );
1697       }
1698       else // n is in ngMesh already, and ngID2 in prev segment is wrong
1699       {
1700         ngID2 = ngMesh.GetNP() + 1;
1701         if ( i > 0 ) // prev segment belongs to same wire
1702         {
1703           netgen::Segment& prevSeg = ngMesh.LineSegment( ngMesh.GetNSeg() );
1704           prevSeg[1] = ngID1;
1705         }
1706       }
1707
1708       // Add the segment
1709
1710       netgen::Segment seg;
1711
1712       seg[0]     = ngID1;                // ng node id
1713       seg[1]     = ngID2;                // ng node id
1714       seg.edgenr = ngMesh.GetNSeg() + 1; // ng segment id
1715       seg.si     = faceID;               // = geom.fmap.FindIndex (face);
1716
1717       for ( int iEnd = 0; iEnd < 2; ++iEnd)
1718       {
1719         const UVPtStruct& pnt = uvPtVec[ i + iEnd ];
1720
1721         seg.epgeominfo[ iEnd ].dist = pnt.param; // param on curve
1722         seg.epgeominfo[ iEnd ].u    = pnt.u;
1723         seg.epgeominfo[ iEnd ].v    = pnt.v;
1724
1725         // find out edge id and node parameter on edge
1726         onVertex = ( pnt.normParam + 1e-10 > vertexNormPar );
1727         if ( onVertex || posShapeID != posID )
1728         {
1729           // get edge id
1730           double normParam = pnt.normParam;
1731           if ( onVertex )
1732             normParam = 0.5 * ( uvPtVec[ i ].normParam + uvPtVec[ i+1 ].normParam );
1733           int edgeIndexInWire = wire->EdgeIndex( normParam );
1734           vertexNormPar = wire->LastParameter( edgeIndexInWire );
1735           const TopoDS_Edge& edge = wire->Edge( edgeIndexInWire );
1736           edgeID = geom.emap.FindIndex( edge );
1737           posID  = posShapeID;
1738           isInternalWire = ( edge.Orientation() == TopAbs_INTERNAL );
1739           // if ( onVertex ) // param on curve is different on each of two edges
1740           //   seg.epgeominfo[ iEnd ].dist = helper.GetNodeU( edge, pnt.node );
1741         }
1742         seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
1743       }
1744
1745       ngMesh.AddSegment (seg);
1746       {
1747         // restrict size of elements near the segment
1748         SMESH_TNodeXYZ np1( n ), np2( uvPtVec[ i+1 ].node );
1749         // get an average size of adjacent segments to avoid sharp change of
1750         // element size (regression on issue 0020452, note 0010898)
1751         int   iPrev = SMESH_MesherHelper::WrapIndex( i-1, nbSegments );
1752         int   iNext = SMESH_MesherHelper::WrapIndex( i+1, nbSegments );
1753         double sumH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ];
1754         int   nbSeg = ( int( segLen[ iPrev ] > sumH / 100.)  +
1755                         int( segLen[ i     ] > sumH / 100.)  +
1756                         int( segLen[ iNext ] > sumH / 100.));
1757         if ( nbSeg > 0 )
1758           RestrictLocalSize( ngMesh, 0.5*(np1+np2), sumH / nbSeg );
1759       }
1760       if ( isInternalWire )
1761       {
1762         swap (seg[0], seg[1]);
1763         swap( seg.epgeominfo[0], seg.epgeominfo[1] );
1764         seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
1765         ngMesh.AddSegment (seg);
1766       }
1767     } // loop on segments on a wire
1768
1769     // close chain of segments
1770     if ( nbSegments > 0 )
1771     {
1772       netgen::Segment& lastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() - int( isInternalWire));
1773       const SMDS_MeshNode * lastNode = uvPtVec.back().node;
1774       lastSeg[1] = node2ngID.insert( make_pair( lastNode, lastSeg[1] )).first->second;
1775       if ( lastSeg[1] > ngMesh.GetNP() )
1776       {
1777         netgen::MeshPoint mp( netgen::Point<3> (lastNode->X(), lastNode->Y(), lastNode->Z()) );
1778         ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
1779         nodeVec.push_back( lastNode );
1780       }
1781       if ( isInternalWire )
1782       {
1783         netgen::Segment& realLastSeg = ngMesh.LineSegment( ngMesh.GetNSeg() );
1784         realLastSeg[0] = lastSeg[1];
1785       }
1786     }
1787
1788 #ifdef DUMP_SEGMENTS
1789     cout << "BEGIN WIRE " << iW << endl;
1790     for ( int i = prevNbNGSeg+1; i <= ngMesh.GetNSeg(); ++i )
1791     {
1792       netgen::Segment& seg = ngMesh.LineSegment( i );
1793       if ( i > 1 ) {
1794         netgen::Segment& prevSeg = ngMesh.LineSegment( i-1 );
1795         if ( seg[0] == prevSeg[1] && seg[1] == prevSeg[0] )
1796         {
1797           cout << "Segment: " << seg.edgenr << endl << "\tis REVRESE of the previous one" << endl;
1798           continue;
1799         }
1800       }
1801       cout << "Segment: " << seg.edgenr << endl
1802            << "\tp1: " << seg[0] << endl
1803            << "\tp2: " << seg[1] << endl
1804            << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
1805            << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
1806            << "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
1807            << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
1808            << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl
1809            << "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
1810     }
1811     cout << "--END WIRE " << iW << endl;
1812 #endif
1813
1814   } // loop on WIREs of a FACE
1815
1816   // add a segment instead of an internal vertex
1817   if ( wasNgMeshEmpty )
1818   {
1819     NETGENPlugin_Internals intShapes( *helper.GetMesh(), helper.GetSubShape(), /*is3D=*/false );
1820     AddIntVerticesInFaces( geom, ngMesh, nodeVec, intShapes );
1821   }
1822   ngMesh.CalcSurfacesOfNode();
1823
1824   return TError();
1825 }
1826
1827 //================================================================================
1828 /*!
1829  * \brief Fill SMESH mesh according to contents of netgen mesh
1830  *  \param occgeo - container of OCCT geometry to mesh
1831  *  \param ngMesh - netgen mesh
1832  *  \param initState - bn of entities in netgen mesh before computing
1833  *  \param sMesh - SMESH mesh to fill in
1834  *  \param nodeVec - vector of nodes in which node index == netgen ID
1835  *  \param comment - returns problem description
1836  *  \param quadHelper - holder of medium nodes of sub-meshes
1837  *  \retval int - error
1838  */
1839 //================================================================================
1840
1841 int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
1842                                    netgen::Mesh&                       ngMesh,
1843                                    const NETGENPlugin_ngMeshInfo&      initState,
1844                                    SMESH_Mesh&                         sMesh,
1845                                    std::vector<const SMDS_MeshNode*>&  nodeVec,
1846                                    SMESH_Comment&                      comment,
1847                                    SMESH_MesherHelper*                 quadHelper)
1848 {
1849   int nbNod = ngMesh.GetNP();
1850   int nbSeg = ngMesh.GetNSeg();
1851   int nbFac = ngMesh.GetNSE();
1852   int nbVol = ngMesh.GetNE();
1853
1854   SMESHDS_Mesh* meshDS = sMesh.GetMeshDS();
1855
1856   // quadHelper is used for either
1857   // 1) making quadratic elements when a lower dimention mesh is loaded
1858   //    to SMESH before convertion to quadratic by NETGEN
1859   // 2) sewing of quadratic elements with quadratic elements of sub-meshes
1860   if ( quadHelper && !quadHelper->GetIsQuadratic() && quadHelper->GetTLinkNodeMap().empty() )
1861     quadHelper = 0;
1862
1863   // -------------------------------------
1864   // Create and insert nodes into nodeVec
1865   // -------------------------------------
1866
1867   nodeVec.resize( nbNod + 1 );
1868   int i, nbInitNod = initState._nbNodes;
1869   for (i = nbInitNod+1; i <= nbNod; ++i )
1870   {
1871     const netgen::MeshPoint& ngPoint = ngMesh.Point(i);
1872     SMDS_MeshNode* node = NULL;
1873     TopoDS_Vertex aVert;
1874     // First, netgen creates nodes on vertices in occgeo.vmap,
1875     // so node index corresponds to vertex index
1876     // but (issue 0020776) netgen does not create nodes with equal coordinates
1877     if ( i-nbInitNod <= occgeo.vmap.Extent() )
1878     {
1879       gp_Pnt p ( NGPOINT_COORDS(ngPoint) );
1880       for (int iV = i-nbInitNod; aVert.IsNull() && iV <= occgeo.vmap.Extent(); ++iV)
1881       {
1882         aVert = TopoDS::Vertex( occgeo.vmap( iV ) );
1883         gp_Pnt pV = BRep_Tool::Pnt( aVert );
1884         if ( p.SquareDistance( pV ) > 1e-20 )
1885           aVert.Nullify();
1886         else
1887           node = const_cast<SMDS_MeshNode*>( SMESH_Algo::VertexNode( aVert, meshDS ));
1888       }
1889     }
1890     if (!node) // node not found on vertex
1891     {
1892       node = meshDS->AddNode( NGPOINT_COORDS( ngPoint ));
1893       if (!aVert.IsNull())
1894         meshDS->SetNodeOnVertex(node, aVert);
1895     }
1896     nodeVec[i] = node;
1897   }
1898
1899   // -------------------------------------------
1900   // Create mesh segments along geometric edges
1901   // -------------------------------------------
1902
1903   int nbInitSeg = initState._nbSegments;
1904   for (i = nbInitSeg+1; i <= nbSeg; ++i )
1905   {
1906     const netgen::Segment& seg = ngMesh.LineSegment(i);
1907     TopoDS_Edge aEdge;
1908     int pinds[3] = { seg.pnums[0], seg.pnums[1], seg.pnums[2] };
1909     int nbp = 0;
1910     double param2 = 0;
1911     for (int j=0; j < 3; ++j)
1912     {
1913       int pind = pinds[j];
1914       if (pind <= 0 || !nodeVec_ACCESS(pind))
1915         break;
1916       ++nbp;
1917       double param;
1918       if (j < 2)
1919       {
1920         if (aEdge.IsNull())
1921         {
1922           int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
1923           if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
1924             aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
1925         }
1926         param = seg.epgeominfo[j].dist;
1927         param2 += param;
1928       }
1929       else // middle point
1930       {
1931         param = param2 * 0.5;
1932       }
1933       if (!aEdge.IsNull() && nodeVec_ACCESS(pind)->getshapeId() < 1)
1934       {
1935         meshDS->SetNodeOnEdge(nodeVec_ACCESS(pind), aEdge, param);
1936       }
1937     }
1938     if ( nbp > 1 )
1939     {
1940       SMDS_MeshEdge* edge = 0;
1941       if (nbp == 2) // second order ?
1942       {
1943         if ( meshDS->FindEdge( nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1])))
1944           continue;
1945         if ( quadHelper ) // final mesh must be quadratic
1946           edge = quadHelper->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]));
1947         else
1948           edge = meshDS->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]));
1949       }
1950       else
1951       {
1952         if ( meshDS->FindEdge( nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]),
1953                                nodeVec_ACCESS(pinds[2])))
1954           continue;
1955         edge = meshDS->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]),
1956                                nodeVec_ACCESS(pinds[2]));
1957       }
1958       if (!edge)
1959       {
1960         if ( comment.empty() ) comment << "Cannot create a mesh edge";
1961         MESSAGE("Cannot create a mesh edge");
1962         nbSeg = nbFac = nbVol = 0;
1963         break;
1964       }
1965       if ( !aEdge.IsNull() && edge->getshapeId() < 1 )
1966         meshDS->SetMeshElementOnShape(edge, aEdge);
1967     }
1968     else if ( comment.empty() )
1969     {
1970       comment << "Invalid netgen segment #" << i;
1971     }
1972   }
1973
1974   // ----------------------------------------
1975   // Create mesh faces along geometric faces
1976   // ----------------------------------------
1977
1978   int nbInitFac = initState._nbFaces;
1979   int quadFaceID = ngMesh.GetNFD() + 1;
1980   if ( nbInitFac < nbFac )
1981     // add a faces descriptor to exclude qudrangle elements generated by NETGEN
1982     // from computation of 3D mesh
1983     ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(quadFaceID, /*solid1=*/0, /*solid2=*/0, 0));
1984
1985   vector<const SMDS_MeshNode*> nodes;
1986   for (i = nbInitFac+1; i <= nbFac; ++i )
1987   {
1988     const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
1989     int aGeomFaceInd = elem.GetIndex();
1990     TopoDS_Face aFace;
1991     if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
1992       aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
1993     nodes.clear();
1994     for (int j=1; j <= elem.GetNP(); ++j)
1995     {
1996       int pind = elem.PNum(j);
1997       if ( pind < 1 || pind >= nodeVec.size() )
1998         break;
1999       if ( SMDS_MeshNode* node = nodeVec_ACCESS(pind))
2000       {
2001         nodes.push_back( node );
2002         if (!aFace.IsNull() && node->getshapeId() < 1)
2003         {
2004           const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
2005           meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
2006         }
2007       }
2008     }
2009     if ( nodes.size() != elem.GetNP() )
2010     {
2011       if ( comment.empty() )
2012         comment << "Invalid netgen 2d element #" << i;
2013       continue; // bad node ids
2014     }
2015     SMDS_MeshFace* face = NULL;
2016     switch (elem.GetType())
2017     {
2018     case netgen::TRIG:
2019       if ( quadHelper ) // final mesh must be quadratic
2020         face = quadHelper->AddFace(nodes[0],nodes[1],nodes[2]);
2021       else
2022         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
2023       break;
2024     case netgen::QUAD:
2025       if ( quadHelper ) // final mesh must be quadratic
2026         face = quadHelper->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
2027       else
2028         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
2029       // exclude qudrangle elements from computation of 3D mesh
2030       const_cast< netgen::Element2d& >( elem ).SetIndex( quadFaceID );
2031       break;
2032     case netgen::TRIG6:
2033       nodes[5] = mediumNode( nodes[0],nodes[1],nodes[5], quadHelper );
2034       nodes[3] = mediumNode( nodes[1],nodes[2],nodes[3], quadHelper );
2035       nodes[4] = mediumNode( nodes[2],nodes[0],nodes[4], quadHelper );
2036       face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
2037       break;
2038     case netgen::QUAD8:
2039       nodes[4] = mediumNode( nodes[0],nodes[1],nodes[4], quadHelper );
2040       nodes[7] = mediumNode( nodes[1],nodes[2],nodes[7], quadHelper );
2041       nodes[5] = mediumNode( nodes[2],nodes[3],nodes[5], quadHelper );
2042       nodes[6] = mediumNode( nodes[3],nodes[0],nodes[6], quadHelper );
2043       face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
2044                              nodes[4],nodes[7],nodes[5],nodes[6]);
2045       // exclude qudrangle elements from computation of 3D mesh
2046       const_cast< netgen::Element2d& >( elem ).SetIndex( quadFaceID );
2047       break;
2048     default:
2049       MESSAGE("NETGEN created a face of unexpected type, ignoring");
2050       continue;
2051     }
2052     if (!face)
2053     {
2054       if ( comment.empty() ) comment << "Cannot create a mesh face";
2055       MESSAGE("Cannot create a mesh face");
2056       nbSeg = nbFac = nbVol = 0;
2057       break;
2058     }
2059     if (!aFace.IsNull())
2060       meshDS->SetMeshElementOnShape(face, aFace);
2061   }
2062
2063   // ------------------
2064   // Create tetrahedra
2065   // ------------------
2066
2067   for (i = 1; i <= nbVol; ++i)
2068   {
2069     const netgen::Element& elem = ngMesh.VolumeElement(i);
2070     int aSolidInd = elem.GetIndex();
2071     TopoDS_Solid aSolid;
2072     if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
2073       aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
2074     nodes.clear();
2075     for (int j=1; j <= elem.GetNP(); ++j)
2076     {
2077       int pind = elem.PNum(j);
2078       if ( pind < 1 || pind >= nodeVec.size() )
2079         break;
2080       if ( SMDS_MeshNode* node = nodeVec_ACCESS(pind) )
2081       {
2082         nodes.push_back(node);
2083         if ( !aSolid.IsNull() && node->getshapeId() < 1 )
2084           meshDS->SetNodeInVolume(node, aSolid);
2085       }
2086     }
2087     if ( nodes.size() != elem.GetNP() )
2088     {
2089       if ( comment.empty() )
2090         comment << "Invalid netgen 3d element #" << i;
2091       continue;
2092     }
2093     SMDS_MeshVolume* vol = NULL;
2094     switch (elem.GetType())
2095     {
2096     case netgen::TET:
2097       vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
2098       break;
2099     case netgen::TET10:
2100       nodes[4] = mediumNode( nodes[0],nodes[1],nodes[4], quadHelper );
2101       nodes[7] = mediumNode( nodes[1],nodes[2],nodes[7], quadHelper );
2102       nodes[5] = mediumNode( nodes[2],nodes[0],nodes[5], quadHelper );
2103       nodes[6] = mediumNode( nodes[0],nodes[3],nodes[6], quadHelper );
2104       nodes[8] = mediumNode( nodes[1],nodes[3],nodes[8], quadHelper );
2105       nodes[9] = mediumNode( nodes[2],nodes[3],nodes[9], quadHelper );
2106       vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
2107                               nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
2108       break;
2109     default:
2110       MESSAGE("NETGEN created a volume of unexpected type, ignoring");
2111       continue;
2112     }
2113     if (!vol)
2114     {
2115       if ( comment.empty() ) comment << "Cannot create a mesh volume";
2116       MESSAGE("Cannot create a mesh volume");
2117       nbSeg = nbFac = nbVol = 0;
2118       break;
2119     }
2120     if (!aSolid.IsNull())
2121       meshDS->SetMeshElementOnShape(vol, aSolid);
2122   }
2123   return comment.empty() ? 0 : 1;
2124 }
2125
2126 namespace
2127 {
2128   //================================================================================
2129   /*!
2130    * \brief Restrict size of elements on the given edge
2131    */
2132   //================================================================================
2133
2134   void setLocalSize(const TopoDS_Edge& edge,
2135                     double             size,
2136                     netgen::Mesh&      mesh)
2137   {
2138     const int nb = 1000;
2139     Standard_Real u1, u2;
2140     Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, u1, u2);
2141     if ( curve.IsNull() )
2142     {
2143       TopoDS_Iterator vIt( edge );
2144       if ( !vIt.More() ) return;
2145       gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vIt.Value() ));
2146       NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size );
2147     }
2148     else
2149     {
2150       Standard_Real delta = (u2-u1)/nb;
2151       for(int i=0; i<nb; i++)
2152       {
2153         Standard_Real u = u1 + delta*i;
2154         gp_Pnt p = curve->Value(u);
2155         NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size );
2156         netgen::Point3d pi(p.X(), p.Y(), p.Z());
2157         double resultSize = mesh.GetH(pi);
2158         if ( resultSize - size > 0.1*size )
2159           // netgen does restriction iff oldH/newH > 1.2 (localh.cpp:136)
2160           NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), resultSize/1.201 );
2161       }
2162     }
2163   }
2164
2165   //================================================================================
2166   /*!
2167    * \brief Convert error into text
2168    */
2169   //================================================================================
2170
2171   std::string text(int err)
2172   {
2173     if ( !err )
2174       return string("");
2175     return
2176       SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task;
2177   }
2178
2179   //================================================================================
2180   /*!
2181    * \brief Convert exception into text
2182    */
2183   //================================================================================
2184
2185   std::string text(Standard_Failure& ex)
2186   {
2187     SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
2188     str << " at " << netgen::multithread.task
2189         << ": " << ex.DynamicType()->Name();
2190     if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
2191       str << ": " << ex.GetMessageString();
2192     return str;
2193   }
2194   //================================================================================
2195   /*!
2196    * \brief Convert exception into text
2197    */
2198   //================================================================================
2199
2200   std::string text(netgen::NgException& ex)
2201   {
2202     SMESH_Comment str("NgException");
2203     if ( strlen( netgen::multithread.task ) > 0 )
2204       str << " at " << netgen::multithread.task;
2205     str << ": " << ex.What();
2206     return str;
2207   }
2208
2209   //================================================================================
2210   /*!
2211    * \brief Looks for triangles lying on a SOLID
2212    */
2213   //================================================================================
2214
2215   bool hasBadElemOnSolid( const list<const SMDS_MeshElement*>& elems,
2216                           SMESH_subMesh*                       solidSM )
2217   {
2218     TopTools_IndexedMapOfShape solidSubs;
2219     TopExp::MapShapes( solidSM->GetSubShape(), solidSubs );
2220     SMESHDS_Mesh* mesh = solidSM->GetFather()->GetMeshDS();
2221
2222     list<const SMDS_MeshElement*>::const_iterator e = elems.begin();
2223     for ( ; e != elems.end(); ++e )
2224     {
2225       const SMDS_MeshElement* elem = *e;
2226       if ( elem->GetType() != SMDSAbs_Face )
2227         continue;
2228       int nbNodesOnSolid = 0;
2229       SMDS_NodeIteratorPtr nIt = elem->nodeIterator();
2230       while ( nIt->more() )
2231       {
2232         const SMDS_MeshNode* n = nIt->next();
2233         const TopoDS_Shape&  s = mesh->IndexToShape( n->getshapeId() );
2234         nbNodesOnSolid += ( !s.IsNull() && solidSubs.Contains( s ));
2235         if ( nbNodesOnSolid > 2 )
2236           return true;
2237       }
2238     }
2239     return false;
2240   }
2241
2242   const double edgeMeshingTime = 0.001;
2243   const double faceMeshingTime = 0.019;
2244   const double edgeFaceMeshingTime = edgeMeshingTime + faceMeshingTime;
2245   const double faceOptimizTime = 0.06;
2246   const double voluMeshingTime = 0.15;
2247   const double volOptimizeTime = 0.77;
2248 }
2249
2250 //=============================================================================
2251 /*!
2252  * Here we are going to use the NETGEN mesher
2253  */
2254 //=============================================================================
2255
2256 bool NETGENPlugin_Mesher::Compute()
2257 {
2258   NETGENPlugin_NetgenLibWrapper ngLib;
2259
2260   netgen::MeshingParameters& mparams = netgen::mparam;
2261   MESSAGE("Compute with:\n"
2262           " max size = " << mparams.maxh << "\n"
2263           " segments per edge = " << mparams.segmentsperedge);
2264   MESSAGE("\n"
2265           " growth rate = " << mparams.grading << "\n"
2266           " elements per radius = " << mparams.curvaturesafety << "\n"
2267           " second order = " << mparams.secondorder << "\n"
2268           " quad allowed = " << mparams.quad << "\n"
2269           " surface curvature = " << mparams.uselocalh << "\n"
2270           " fuse edges = " << netgen::merge_solids);
2271
2272   SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
2273   SMESH_MesherHelper quadHelper( *_mesh );
2274   quadHelper.SetIsQuadratic( mparams.secondorder );
2275
2276   static string debugFile = "/tmp/ngMesh.py"; /* to call toPython( ngMesh, debugFile )
2277                                                  while debugging netgen */
2278   // -------------------------
2279   // Prepare OCC geometry
2280   // -------------------------
2281
2282   netgen::OCCGeometry occgeo;
2283   list< SMESH_subMesh* > meshedSM[3]; // for 0-2 dimensions
2284   NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume );
2285   PrepareOCCgeometry( occgeo, _shape, *_mesh, meshedSM, &internals );
2286   _occgeom = &occgeo;
2287
2288   _totalTime = edgeFaceMeshingTime;
2289   if ( _optimize )
2290     _totalTime += faceOptimizTime;
2291   if ( _isVolume )
2292     _totalTime += voluMeshingTime + ( _optimize ? volOptimizeTime : 0 );
2293   double doneTime = 0;
2294   _ticTime = -1;
2295   _progressTic = 1;
2296   _curShapeIndex = -1;
2297
2298   // -------------------------
2299   // Generate the mesh
2300   // -------------------------
2301
2302   _ngMesh = NULL;
2303   NETGENPlugin_ngMeshInfo initState; // it remembers size of ng mesh equal to size of Smesh
2304
2305   SMESH_Comment comment;
2306   int err = 0;
2307
2308   // vector of nodes in which node index == netgen ID
2309   vector< const SMDS_MeshNode* > nodeVec;
2310
2311   {
2312     // ----------------
2313     // compute 1D mesh
2314     // ----------------
2315     if ( _simpleHyp )
2316     {
2317       // not to RestrictLocalH() according to curvature during MESHCONST_ANALYSE
2318       mparams.uselocalh = false;
2319       mparams.grading = 0.8; // not limitited size growth
2320
2321       if ( _simpleHyp->GetNumberOfSegments() )
2322         // nb of segments
2323         mparams.maxh = occgeo.boundingbox.Diam();
2324       else
2325         // segment length
2326         mparams.maxh = _simpleHyp->GetLocalLength();
2327     }
2328
2329     if ( mparams.maxh == 0.0 )
2330       mparams.maxh = occgeo.boundingbox.Diam();
2331     if ( _simpleHyp || ( mparams.minh == 0.0 && _fineness != NETGENPlugin_Hypothesis::UserDefined))
2332       mparams.minh = GetDefaultMinSize( _shape, mparams.maxh );
2333
2334     // Local size on faces
2335     occgeo.face_maxh = mparams.maxh;
2336
2337     // Let netgen create _ngMesh and calculate element size on not meshed shapes
2338 #ifndef NETGEN_V5
2339     char *optstr = 0;
2340 #endif
2341     int startWith = netgen::MESHCONST_ANALYSE;
2342     int endWith   = netgen::MESHCONST_ANALYSE;
2343     try
2344     {
2345       OCC_CATCH_SIGNALS;
2346 #ifdef NETGEN_V5
2347       err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
2348 #else
2349       err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
2350 #endif
2351       if(netgen::multithread.terminate)
2352         return false;
2353
2354       comment << text(err);
2355     }
2356     catch (Standard_Failure& ex)
2357     {
2358       comment << text(ex);
2359     }
2360     err = 0; //- MESHCONST_ANALYSE isn't so important step
2361     if ( !_ngMesh )
2362       return false;
2363     ngLib.setMesh(( Ng_Mesh*) _ngMesh );
2364
2365     _ngMesh->ClearFaceDescriptors(); // we make descriptors our-self
2366
2367     if ( _simpleHyp )
2368     {
2369       // Pass 1D simple parameters to NETGEN
2370       // --------------------------------
2371       int      nbSeg = _simpleHyp->GetNumberOfSegments();
2372       double segSize = _simpleHyp->GetLocalLength();
2373       for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE )
2374       {
2375         const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE));
2376         if ( nbSeg )
2377           segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 );
2378         setLocalSize( e, segSize, *_ngMesh );
2379       }
2380     }
2381     else // if ( ! _simpleHyp )
2382     {
2383       // Local size on vertices and edges
2384       // --------------------------------
2385       for(std::map<int,double>::const_iterator it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++)
2386       {
2387         int key = (*it).first;
2388         double hi = (*it).second;
2389         const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2390         const TopoDS_Edge& e = TopoDS::Edge(shape);
2391         setLocalSize( e, hi, *_ngMesh );
2392       }
2393       for(std::map<int,double>::const_iterator it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++)
2394       {
2395         int key = (*it).first;
2396         double hi = (*it).second;
2397         const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2398         const TopoDS_Vertex& v = TopoDS::Vertex(shape);
2399         gp_Pnt p = BRep_Tool::Pnt(v);
2400         NETGENPlugin_Mesher::RestrictLocalSize( *_ngMesh, p.XYZ(), hi );
2401       }
2402       for(map<int,double>::const_iterator it=FaceId2LocalSize.begin();
2403           it!=FaceId2LocalSize.end(); it++)
2404       {
2405         int key = (*it).first;
2406         double val = (*it).second;
2407         const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2408         int faceNgID = occgeo.fmap.FindIndex(shape);
2409         occgeo.SetFaceMaxH(faceNgID, val);
2410         for ( TopExp_Explorer edgeExp( shape, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
2411           setLocalSize( TopoDS::Edge( edgeExp.Current() ), val, *_ngMesh );
2412       }
2413     }
2414
2415     // Precompute internal edges (issue 0020676) in order to
2416     // add mesh on them correctly (twice) to netgen mesh
2417     if ( !err && internals.hasInternalEdges() )
2418     {
2419       // load internal shapes into OCCGeometry
2420       netgen::OCCGeometry intOccgeo;
2421       internals.getInternalEdges( intOccgeo.fmap, intOccgeo.emap, intOccgeo.vmap, meshedSM );
2422       intOccgeo.boundingbox = occgeo.boundingbox;
2423       intOccgeo.shape = occgeo.shape;
2424       intOccgeo.face_maxh.SetSize(intOccgeo.fmap.Extent());
2425       intOccgeo.face_maxh = netgen::mparam.maxh;
2426       netgen::Mesh *tmpNgMesh = NULL;
2427       try
2428       {
2429         OCC_CATCH_SIGNALS;
2430         // compute local H on internal shapes in the main mesh
2431         //OCCSetLocalMeshSize(intOccgeo, *_ngMesh); it deletes _ngMesh->localH
2432
2433         // let netgen create a temporary mesh
2434 #ifdef NETGEN_V5
2435         netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, mparams, startWith, endWith);
2436 #else
2437         netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, startWith, endWith, optstr);
2438 #endif
2439         if(netgen::multithread.terminate)
2440           return false;
2441
2442         // copy LocalH from the main to temporary mesh
2443         initState.transferLocalH( _ngMesh, tmpNgMesh );
2444
2445         // compute mesh on internal edges
2446         startWith = endWith = netgen::MESHCONST_MESHEDGES;
2447 #ifdef NETGEN_V5
2448         err = netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, mparams, startWith, endWith);
2449 #else
2450         err = netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, startWith, endWith, optstr);
2451 #endif
2452         comment << text(err);
2453       }
2454       catch (Standard_Failure& ex)
2455       {
2456         comment << text(ex);
2457         err = 1;
2458       }
2459       initState.restoreLocalH( tmpNgMesh );
2460
2461       // fill SMESH by netgen mesh
2462       vector< const SMDS_MeshNode* > tmpNodeVec;
2463       FillSMesh( intOccgeo, *tmpNgMesh, initState, *_mesh, tmpNodeVec, comment );
2464       err = ( err || !comment.empty() );
2465
2466       nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)tmpNgMesh);
2467     }
2468
2469     // Fill _ngMesh with nodes and segments of computed submeshes
2470     if ( !err )
2471     {
2472       err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_0D ]) &&
2473                 FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_1D ], &quadHelper));
2474     }
2475     initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2476
2477     // Compute 1d mesh
2478     if (!err)
2479     {
2480       startWith = endWith = netgen::MESHCONST_MESHEDGES;
2481       try
2482       {
2483         OCC_CATCH_SIGNALS;
2484 #ifdef NETGEN_V5
2485         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
2486 #else
2487         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
2488 #endif
2489         if(netgen::multithread.terminate)
2490           return false;
2491
2492         comment << text(err);
2493       }
2494       catch (Standard_Failure& ex)
2495       {
2496         comment << text(ex);
2497         err = 1;
2498       }
2499     }
2500     if ( _isVolume )
2501       _ticTime = ( doneTime += edgeMeshingTime ) / _totalTime / _progressTic;
2502
2503     mparams.uselocalh = true; // restore as it is used at surface optimization
2504
2505     // ---------------------
2506     // compute surface mesh
2507     // ---------------------
2508     if (!err)
2509     {
2510       // Pass 2D simple parameters to NETGEN
2511       if ( _simpleHyp ) {
2512         if ( double area = _simpleHyp->GetMaxElementArea() ) {
2513           // face area
2514           mparams.maxh = sqrt(2. * area/sqrt(3.0));
2515           mparams.grading = 0.4; // moderate size growth
2516         }
2517         else {
2518           // length from edges
2519           if ( _ngMesh->GetNSeg() ) {
2520             double edgeLength = 0;
2521             TopTools_MapOfShape visitedEdges;
2522             for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() )
2523               if( visitedEdges.Add(exp.Current()) )
2524                 edgeLength += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() ));
2525             // we have to multiply length by 2 since for each TopoDS_Edge there
2526             // are double set of NETGEN edges, in other words, we have to
2527             // divide _ngMesh->GetNSeg() by 2.
2528             mparams.maxh = 2*edgeLength / _ngMesh->GetNSeg();
2529           }
2530           else {
2531             mparams.maxh = 1000;
2532           }
2533           mparams.grading = 0.2; // slow size growth
2534         }
2535         mparams.quad = _simpleHyp->GetAllowQuadrangles();
2536         mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
2537         _ngMesh->SetGlobalH (mparams.maxh);
2538         netgen::Box<3> bb = occgeo.GetBoundingBox();
2539         bb.Increase (bb.Diam()/20);
2540         _ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading);
2541       }
2542
2543       // Care of vertices internal in faces (issue 0020676)
2544       if ( internals.hasInternalVertexInFace() )
2545       {
2546         // store computed segments in SMESH in order not to create SMESH
2547         // edges for ng segments added by AddIntVerticesInFaces()
2548         FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
2549         // add segments to faces with internal vertices
2550         AddIntVerticesInFaces( occgeo, *_ngMesh, nodeVec, internals );
2551         initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2552       }
2553
2554       // Build viscous layers
2555       if ( _isViscousLayers2D )
2556       {
2557         if ( !internals.hasInternalVertexInFace() ) {
2558           FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
2559           initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2560         }
2561         SMESH_ProxyMesh::Ptr viscousMesh;
2562         SMESH_MesherHelper   helper( *_mesh );
2563         for ( int faceID = 1; faceID <= occgeo.fmap.Extent(); ++faceID )
2564         {
2565           const TopoDS_Face& F = TopoDS::Face( occgeo.fmap( faceID ));
2566           viscousMesh = StdMeshers_ViscousLayers2D::Compute( *_mesh, F );
2567           if ( !viscousMesh )
2568             return false;
2569           // exclude from computation ng segments built on EDGEs of F
2570           for (int i = 1; i <= _ngMesh->GetNSeg(); i++)
2571           {
2572             netgen::Segment & seg = _ngMesh->LineSegment(i);
2573             if (seg.si == faceID)
2574               seg.si = 0;
2575           }
2576           // add new segments to _ngMesh instead of excluded ones
2577           helper.SetSubShape( F );
2578           TSideVector wires =
2579             StdMeshers_FaceSide::GetFaceWires( F, *_mesh, /*skipMediumNodes=*/true,
2580                                                error, viscousMesh );
2581           error = AddSegmentsToMesh( *_ngMesh, occgeo, wires, helper, nodeVec );
2582
2583           if ( !error ) error = SMESH_ComputeError::New();
2584         }
2585         initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2586       }
2587
2588       // Let netgen compute 2D mesh
2589       startWith = netgen::MESHCONST_MESHSURFACE;
2590       endWith = _optimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
2591       try
2592       {
2593         OCC_CATCH_SIGNALS;
2594 #ifdef NETGEN_V5
2595         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
2596 #else
2597         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
2598 #endif
2599         if(netgen::multithread.terminate)
2600           return false;
2601
2602         comment << text (err);
2603       }
2604       catch (Standard_Failure& ex)
2605       {
2606         comment << text(ex);
2607         //err = 1; -- try to make volumes anyway
2608       }
2609       catch (netgen::NgException exc)
2610       {
2611         comment << text(exc);
2612         //err = 1; -- try to make volumes anyway
2613       }
2614     }
2615     if ( _isVolume )
2616     {
2617       doneTime += faceMeshingTime + ( _optimize ? faceOptimizTime : 0 );
2618       _ticTime = doneTime / _totalTime / _progressTic;
2619     }
2620     // ---------------------
2621     // generate volume mesh
2622     // ---------------------
2623     // Fill _ngMesh with nodes and faces of computed 2D submeshes
2624     if ( !err && _isVolume && ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad ))
2625     {
2626       // load SMESH with computed segments and faces
2627       FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
2628
2629       // compute pyramids on quadrangles
2630       SMESH_ProxyMesh::Ptr proxyMesh;
2631       if ( _mesh->NbQuadrangles() > 0 )
2632         for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS )
2633         {
2634           StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
2635           proxyMesh.reset( Adaptor );
2636
2637           int nbPyrams = _mesh->NbPyramids();
2638           Adaptor->Compute( *_mesh, occgeo.somap(iS) );
2639           if ( nbPyrams != _mesh->NbPyramids() )
2640           {
2641             list< SMESH_subMesh* > quadFaceSM;
2642             for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next())
2643               if ( Adaptor->GetProxySubMesh( face.Current() ))
2644               {
2645                 quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() ));
2646                 meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() );
2647               }
2648             FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, &quadHelper, proxyMesh);
2649           }
2650         }
2651       // fill _ngMesh with faces of sub-meshes
2652       err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ], &quadHelper));
2653       initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2654       //toPython( _ngMesh, "/tmp/ngPython.py");
2655     }
2656     if (!err && _isVolume)
2657     {
2658       // Pass 3D simple parameters to NETGEN
2659       const NETGENPlugin_SimpleHypothesis_3D* simple3d =
2660         dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
2661       if ( simple3d ) {
2662         if ( double vol = simple3d->GetMaxElementVolume() ) {
2663           // max volume
2664           mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
2665           mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
2666         }
2667         else {
2668           // length from faces
2669           mparams.maxh = _ngMesh->AverageH();
2670         }
2671         _ngMesh->SetGlobalH (mparams.maxh);
2672         mparams.grading = 0.4;
2673 #ifdef NETGEN_V5
2674         _ngMesh->CalcLocalH(mparams.grading);
2675 #else
2676         _ngMesh->CalcLocalH();
2677 #endif
2678       }
2679       // Care of vertices internal in solids and internal faces (issue 0020676)
2680       if ( internals.hasInternalVertexInSolid() || internals.hasInternalFaces() )
2681       {
2682         // store computed faces in SMESH in order not to create SMESH
2683         // faces for ng faces added here
2684         FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
2685         // add ng faces to solids with internal vertices
2686         AddIntVerticesInSolids( occgeo, *_ngMesh, nodeVec, internals );
2687         // duplicate mesh faces on internal faces
2688         FixIntFaces( occgeo, *_ngMesh, internals );
2689         initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2690       }
2691       // Let netgen compute 3D mesh
2692       startWith = endWith = netgen::MESHCONST_MESHVOLUME;
2693       try
2694       {
2695         OCC_CATCH_SIGNALS;
2696 #ifdef NETGEN_V5
2697         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
2698 #else
2699         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
2700 #endif
2701         if(netgen::multithread.terminate)
2702           return false;
2703
2704         if ( comment.empty() ) // do not overwrite a previos error
2705           comment << text(err);
2706       }
2707       catch (Standard_Failure& ex)
2708       {
2709         if ( comment.empty() ) // do not overwrite a previos error
2710           comment << text(ex);
2711         err = 1;
2712       }
2713       catch (netgen::NgException exc)
2714       {
2715         if ( comment.empty() ) // do not overwrite a previos error
2716           comment << text(exc);
2717         err = 1;
2718       }
2719       _ticTime = ( doneTime += voluMeshingTime ) / _totalTime / _progressTic;
2720
2721       // Let netgen optimize 3D mesh
2722       if ( !err && _optimize )
2723       {
2724         startWith = endWith = netgen::MESHCONST_OPTVOLUME;
2725         try
2726         {
2727           OCC_CATCH_SIGNALS;
2728 #ifdef NETGEN_V5
2729           err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
2730 #else
2731           err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
2732 #endif
2733           if(netgen::multithread.terminate)
2734             return false;
2735
2736           if ( comment.empty() ) // do not overwrite a previos error
2737             comment << text(err);
2738         }
2739         catch (Standard_Failure& ex)
2740         {
2741           if ( comment.empty() ) // do not overwrite a previos error
2742             comment << text(ex);
2743         }
2744         catch (netgen::NgException exc)
2745         {
2746           if ( comment.empty() ) // do not overwrite a previos error
2747             comment << text(exc);
2748         }
2749       }
2750     }
2751     if (!err && mparams.secondorder > 0)
2752     {
2753       try
2754       {
2755         OCC_CATCH_SIGNALS;
2756         if ( !meshedSM[ MeshDim_1D ].empty() )
2757         {
2758           // remove segments not attached to geometry (IPAL0052479)
2759           for (int i = 1; i <= _ngMesh->GetNSeg(); ++i)
2760           {
2761             const netgen::Segment & seg = _ngMesh->LineSegment (i);
2762             if ( seg.epgeominfo[ 0 ].edgenr == 0 )
2763               _ngMesh->DeleteSegment( i );
2764           }
2765           _ngMesh->Compress();
2766         }
2767         // convert to quadratic
2768         netgen::OCCRefinementSurfaces ref (occgeo);
2769         ref.MakeSecondOrder (*_ngMesh);
2770
2771         // care of elements already loaded to SMESH
2772         // if ( initState._nbSegments > 0 )
2773         //   makeQuadratic( occgeo.emap, _mesh );
2774         // if ( initState._nbFaces > 0 )
2775         //   makeQuadratic( occgeo.fmap, _mesh );
2776       }
2777       catch (Standard_Failure& ex)
2778       {
2779         if ( comment.empty() ) // do not overwrite a previos error
2780           comment << "Exception in netgen at passing to 2nd order ";
2781       }
2782       catch (netgen::NgException exc)
2783       {
2784         if ( comment.empty() ) // do not overwrite a previos error
2785           comment << exc.What();
2786       }
2787     }
2788   }
2789
2790   _ticTime = 0.98 / _progressTic;
2791
2792   int nbNod = _ngMesh->GetNP();
2793   int nbSeg = _ngMesh->GetNSeg();
2794   int nbFac = _ngMesh->GetNSE();
2795   int nbVol = _ngMesh->GetNE();
2796   bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
2797
2798   MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
2799           ", nb nodes: "    << nbNod <<
2800           ", nb segments: " << nbSeg <<
2801           ", nb faces: "    << nbFac <<
2802           ", nb volumes: "  << nbVol);
2803
2804   // Feed back the SMESHDS with the generated Nodes and Elements
2805   if ( true /*isOK*/ ) // get whatever built
2806   {
2807     FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
2808
2809     if ( quadHelper.GetIsQuadratic() ) // remove free nodes
2810       for ( size_t i = 0; i < nodeVec.size(); ++i )
2811         if ( nodeVec[i] && nodeVec[i]->NbInverseElements() == 0 )
2812           _mesh->GetMeshDS()->RemoveFreeNode( nodeVec[i], 0, /*fromGroups=*/false );
2813   }
2814   SMESH_ComputeErrorPtr readErr = ReadErrors(nodeVec);
2815   if ( readErr && !readErr->myBadElements.empty() )
2816     error = readErr;
2817
2818   if ( error->IsOK() && ( !isOK || comment.size() > 0 ))
2819     error->myName = COMPERR_ALGO_FAILED;
2820   if ( !comment.empty() )
2821     error->myComment = comment;
2822
2823   // SetIsAlwaysComputed( true ) to empty sub-meshes, which
2824   // appear if the geometry contains coincident sub-shape due
2825   // to bool merge_solids = 1; in netgen/libsrc/occ/occgenmesh.cpp
2826   const int nbMaps = 2;
2827   const TopTools_IndexedMapOfShape* geoMaps[nbMaps] =
2828     { & occgeo.vmap, & occgeo.emap/*, & occgeo.fmap*/ };
2829   for ( int iMap = 0; iMap < nbMaps; ++iMap )
2830     for (int i = 1; i <= geoMaps[iMap]->Extent(); i++)
2831       if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( geoMaps[iMap]->FindKey(i)))
2832         if ( !sm->IsMeshComputed() )
2833           sm->SetIsAlwaysComputed( true );
2834
2835   // set bad compute error to subshapes of all failed sub-shapes
2836   if ( !error->IsOK() )
2837   {
2838     bool pb2D = false, pb3D = false;
2839     for (int i = 1; i <= occgeo.fmap.Extent(); i++) {
2840       int status = occgeo.facemeshstatus[i-1];
2841       if (status == 1 ) continue;
2842       if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.fmap( i ))) {
2843         SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2844         if ( !smError || smError->IsOK() ) {
2845           if ( status == -1 )
2846             smError.reset( new SMESH_ComputeError( *error ));
2847           else
2848             smError.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, "Ignored" ));
2849           if ( SMESH_Algo::GetMeshError( sm ) == SMESH_Algo::MEr_OK )
2850             smError->myName = COMPERR_WARNING;
2851         }
2852         pb2D = pb2D || smError->IsKO();
2853       }
2854     }
2855     if ( !pb2D ) // all faces are OK
2856       for (int i = 1; i <= occgeo.somap.Extent(); i++)
2857         if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.somap( i )))
2858         {
2859           bool smComputed = nbVol && !sm->IsEmpty();
2860           if ( smComputed && internals.hasInternalVertexInSolid( sm->GetId() ))
2861           {
2862             int nbIntV = internals.getSolidsWithVertices().find( sm->GetId() )->second.size();
2863             SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
2864             smComputed = ( smDS->NbElements() > 0 || smDS->NbNodes() > nbIntV );
2865           }
2866           SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2867           if ( !smComputed && ( !smError || smError->IsOK() ))
2868           {
2869             smError.reset( new SMESH_ComputeError( *error ));
2870             if ( nbVol && SMESH_Algo::GetMeshError( sm ) == SMESH_Algo::MEr_OK )
2871             {
2872               smError->myName = COMPERR_WARNING;
2873             }
2874             else if ( !smError->myBadElements.empty() ) // bad surface mesh
2875             {
2876               if ( !hasBadElemOnSolid( smError->myBadElements, sm ))
2877                 smError.reset();
2878             }
2879           }
2880           pb3D = pb3D || ( smError && smError->IsKO() );
2881         }
2882     if ( !pb2D && !pb3D )
2883       err = 0; // no fatal errors, only warnings
2884   }
2885
2886   ngLib._isComputeOk = !err;
2887
2888   return !err;
2889 }
2890
2891 //=============================================================================
2892 /*!
2893  * Evaluate
2894  */
2895 //=============================================================================
2896 bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
2897 {
2898   netgen::MeshingParameters& mparams = netgen::mparam;
2899
2900
2901   // -------------------------
2902   // Prepare OCC geometry
2903   // -------------------------
2904   netgen::OCCGeometry occgeo;
2905   list< SMESH_subMesh* > meshedSM[4]; // for 0-3 dimensions
2906   NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume );
2907   PrepareOCCgeometry( occgeo, _shape, *_mesh, meshedSM, &internals );
2908
2909   bool tooManyElems = false;
2910   const int hugeNb = std::numeric_limits<int>::max() / 100;
2911
2912   // ----------------
2913   // evaluate 1D
2914   // ----------------
2915   // pass 1D simple parameters to NETGEN
2916   if ( _simpleHyp )
2917   {
2918     // not to RestrictLocalH() according to curvature during MESHCONST_ANALYSE
2919     mparams.uselocalh = false;
2920     mparams.grading = 0.8; // not limitited size growth
2921
2922     if ( _simpleHyp->GetNumberOfSegments() )
2923       // nb of segments
2924       mparams.maxh = occgeo.boundingbox.Diam();
2925     else
2926       // segment length
2927       mparams.maxh = _simpleHyp->GetLocalLength();
2928   }
2929
2930   if ( mparams.maxh == 0.0 )
2931     mparams.maxh = occgeo.boundingbox.Diam();
2932   if ( _simpleHyp || ( mparams.minh == 0.0 && _fineness != NETGENPlugin_Hypothesis::UserDefined))
2933     mparams.minh = GetDefaultMinSize( _shape, mparams.maxh );
2934
2935   // let netgen create _ngMesh and calculate element size on not meshed shapes
2936   NETGENPlugin_NetgenLibWrapper ngLib;
2937   netgen::Mesh *ngMesh = NULL;
2938 #ifndef NETGEN_V5
2939   char *optstr = 0;
2940 #endif
2941   int startWith = netgen::MESHCONST_ANALYSE;
2942   int endWith   = netgen::MESHCONST_MESHEDGES;
2943 #ifdef NETGEN_V5
2944   int err = netgen::OCCGenerateMesh(occgeo, ngMesh, mparams, startWith, endWith);
2945 #else
2946   int err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
2947 #endif
2948
2949   if(netgen::multithread.terminate)
2950     return false;
2951
2952   ngLib.setMesh(( Ng_Mesh*) ngMesh );
2953   if (err) {
2954     if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( _shape ))
2955       sm->GetComputeError().reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED ));
2956     return false;
2957   }
2958   if ( _simpleHyp )
2959   {
2960     // Pass 1D simple parameters to NETGEN
2961     // --------------------------------
2962     int      nbSeg = _simpleHyp->GetNumberOfSegments();
2963     double segSize = _simpleHyp->GetLocalLength();
2964     for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE )
2965     {
2966       const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE));
2967       if ( nbSeg )
2968         segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 );
2969       setLocalSize( e, segSize, *ngMesh );
2970     }
2971   }
2972   else // if ( ! _simpleHyp )
2973   {
2974     // Local size on vertices and edges
2975     // --------------------------------
2976     for(std::map<int,double>::const_iterator it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++)
2977     {
2978       int key = (*it).first;
2979       double hi = (*it).second;
2980       const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2981       const TopoDS_Edge& e = TopoDS::Edge(shape);
2982       setLocalSize( e, hi, *ngMesh );
2983     }
2984     for(std::map<int,double>::const_iterator it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++)
2985     {
2986       int key = (*it).first;
2987       double hi = (*it).second;
2988       const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2989       const TopoDS_Vertex& v = TopoDS::Vertex(shape);
2990       gp_Pnt p = BRep_Tool::Pnt(v);
2991       NETGENPlugin_Mesher::RestrictLocalSize( *ngMesh, p.XYZ(), hi );
2992     }
2993     for(map<int,double>::const_iterator it=FaceId2LocalSize.begin();
2994         it!=FaceId2LocalSize.end(); it++)
2995     {
2996       int key = (*it).first;
2997       double val = (*it).second;
2998       const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2999       int faceNgID = occgeo.fmap.FindIndex(shape);
3000       occgeo.SetFaceMaxH(faceNgID, val);
3001       for ( TopExp_Explorer edgeExp( shape, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
3002         setLocalSize( TopoDS::Edge( edgeExp.Current() ), val, *ngMesh );
3003     }
3004   }
3005   // calculate total nb of segments and length of edges
3006   double fullLen = 0.0;
3007   int fullNbSeg = 0;
3008   int entity = mparams.secondorder > 0 ? SMDSEntity_Quad_Edge : SMDSEntity_Edge;
3009   TopTools_DataMapOfShapeInteger Edge2NbSeg;
3010   for (TopExp_Explorer exp(_shape, TopAbs_EDGE); exp.More(); exp.Next())
3011   {
3012     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
3013     if( !Edge2NbSeg.Bind(E,0) )
3014       continue;
3015
3016     double aLen = SMESH_Algo::EdgeLength(E);
3017     fullLen += aLen;
3018
3019     vector<int>& aVec = aResMap[_mesh->GetSubMesh(E)];
3020     if ( aVec.empty() )
3021       aVec.resize( SMDSEntity_Last, 0);
3022     else
3023       fullNbSeg += aVec[ entity ];
3024   }
3025
3026   // store nb of segments computed by Netgen
3027   NCollection_Map<Link> linkMap;
3028   for (int i = 1; i <= ngMesh->GetNSeg(); ++i )
3029   {
3030     const netgen::Segment& seg = ngMesh->LineSegment(i);
3031     Link link(seg[0], seg[1]);
3032     if ( !linkMap.Add( link )) continue;
3033     int aGeomEdgeInd = seg.epgeominfo[0].edgenr;
3034     if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
3035     {
3036       vector<int>& aVec = aResMap[_mesh->GetSubMesh(occgeo.emap(aGeomEdgeInd))];
3037       aVec[ entity ]++;
3038     }
3039   }
3040   // store nb of nodes on edges computed by Netgen
3041   TopTools_DataMapIteratorOfDataMapOfShapeInteger Edge2NbSegIt(Edge2NbSeg);
3042   for (; Edge2NbSegIt.More(); Edge2NbSegIt.Next())
3043   {
3044     vector<int>& aVec = aResMap[_mesh->GetSubMesh(Edge2NbSegIt.Key())];
3045     if ( aVec[ entity ] > 1 && aVec[ SMDSEntity_Node ] == 0 )
3046       aVec[SMDSEntity_Node] = mparams.secondorder > 0  ? 2*aVec[ entity ]-1 : aVec[ entity ]-1;
3047
3048     fullNbSeg += aVec[ entity ];
3049     Edge2NbSeg( Edge2NbSegIt.Key() ) = aVec[ entity ];
3050   }
3051   if ( fullNbSeg == 0 )
3052     return false;
3053
3054   // ----------------
3055   // evaluate 2D
3056   // ----------------
3057   if ( _simpleHyp ) {
3058     if ( double area = _simpleHyp->GetMaxElementArea() ) {
3059       // face area
3060       mparams.maxh = sqrt(2. * area/sqrt(3.0));
3061       mparams.grading = 0.4; // moderate size growth
3062     }
3063     else {
3064       // length from edges
3065       mparams.maxh = fullLen/fullNbSeg;
3066       mparams.grading = 0.2; // slow size growth
3067     }
3068   }
3069   mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
3070   mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
3071
3072   for (TopExp_Explorer exp(_shape, TopAbs_FACE); exp.More(); exp.Next())
3073   {
3074     TopoDS_Face F = TopoDS::Face( exp.Current() );
3075     SMESH_subMesh *sm = _mesh->GetSubMesh(F);
3076     GProp_GProps G;
3077     BRepGProp::SurfaceProperties(F,G);
3078     double anArea = G.Mass();
3079     tooManyElems = tooManyElems || ( anArea/hugeNb > mparams.maxh*mparams.maxh );
3080     int nb1d = 0;
3081     if ( !tooManyElems )
3082     {
3083       TopTools_MapOfShape egdes;
3084       for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next())
3085         if ( egdes.Add( exp1.Current() ))
3086           nb1d += Edge2NbSeg.Find(exp1.Current());
3087     }
3088     int nbFaces = tooManyElems ? hugeNb : int( 4*anArea / (mparams.maxh*mparams.maxh*sqrt(3.)));
3089     int nbNodes = tooManyElems ? hugeNb : (( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
3090
3091     vector<int> aVec(SMDSEntity_Last, 0);
3092     if( mparams.secondorder > 0 ) {
3093       int nb1d_in = (nbFaces*3 - nb1d) / 2;
3094       aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
3095       aVec[SMDSEntity_Quad_Triangle] = nbFaces;
3096     }
3097     else {
3098       aVec[SMDSEntity_Node] = Max ( nbNodes, 0  );
3099       aVec[SMDSEntity_Triangle] = nbFaces;
3100     }
3101     aResMap[sm].swap(aVec);
3102   }
3103
3104   // ----------------
3105   // evaluate 3D
3106   // ----------------
3107   if(_isVolume) {
3108     // pass 3D simple parameters to NETGEN
3109     const NETGENPlugin_SimpleHypothesis_3D* simple3d =
3110       dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
3111     if ( simple3d ) {
3112       if ( double vol = simple3d->GetMaxElementVolume() ) {
3113         // max volume
3114         mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
3115         mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
3116       }
3117       else {
3118         // using previous length from faces
3119       }
3120       mparams.grading = 0.4;
3121       mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
3122     }
3123     GProp_GProps G;
3124     BRepGProp::VolumeProperties(_shape,G);
3125     double aVolume = G.Mass();
3126     double tetrVol = 0.1179*mparams.maxh*mparams.maxh*mparams.maxh;
3127     tooManyElems = tooManyElems || ( aVolume/hugeNb > tetrVol );
3128     int nbVols = tooManyElems ? hugeNb : int(aVolume/tetrVol);
3129     int nb1d_in = int(( nbVols*6 - fullNbSeg ) / 6 );
3130     vector<int> aVec(SMDSEntity_Last, 0 );
3131     if ( tooManyElems ) // avoid FPE
3132     {
3133       aVec[SMDSEntity_Node] = hugeNb;
3134       aVec[ mparams.secondorder > 0 ? SMDSEntity_Quad_Tetra : SMDSEntity_Tetra] = hugeNb;
3135     }
3136     else
3137     {
3138       if( mparams.secondorder > 0 ) {
3139         aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in;
3140         aVec[SMDSEntity_Quad_Tetra] = nbVols;
3141       }
3142       else {
3143         aVec[SMDSEntity_Node] = nb1d_in/3 + 1;
3144         aVec[SMDSEntity_Tetra] = nbVols;
3145       }
3146     }
3147     SMESH_subMesh *sm = _mesh->GetSubMesh(_shape);
3148     aResMap[sm].swap(aVec);
3149   }
3150
3151   return true;
3152 }
3153
3154 double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder,
3155                                         const int *       algoProgressTic,
3156                                         const double *    algoProgress) const
3157 {
3158   ((int&) _progressTic ) = *algoProgressTic + 1;
3159
3160   if ( !_occgeom ) return 0;
3161
3162   double progress = -1;
3163   if ( !_isVolume )
3164   {
3165     if ( _ticTime < 0 && netgen::multithread.task[0] == 'O'/*Optimizing surface*/ )
3166     {
3167       ((double&) _ticTime ) = edgeFaceMeshingTime / _totalTime / _progressTic;
3168     }
3169     else if ( !_optimize /*&& _occgeom->fmap.Extent() > 1*/ )
3170     {
3171       int doneShapeIndex = -1;
3172       while ( doneShapeIndex+1 < _occgeom->facemeshstatus.Size() &&
3173               _occgeom->facemeshstatus[ doneShapeIndex+1 ])
3174         doneShapeIndex++;
3175       if ( doneShapeIndex+1 != _curShapeIndex )
3176       {
3177         ((int&) _curShapeIndex) = doneShapeIndex+1;
3178         double    doneShapeRate = _curShapeIndex / double( _occgeom->fmap.Extent() );
3179         double         doneTime = edgeMeshingTime + doneShapeRate * faceMeshingTime;
3180         ((double&)    _ticTime) = doneTime / _totalTime / _progressTic;
3181         // cout << "shape " << _curShapeIndex << " _ticTime " << _ticTime
3182         //      << " " << doneTime / _totalTime / _progressTic << endl;
3183       }
3184     }
3185   }
3186   else if ( !_optimize && _occgeom->somap.Extent() > 1 )
3187   {
3188     int curShapeIndex = _curShapeIndex;
3189     if ( _ngMesh->GetNE() > 0 )
3190     {
3191       netgen::Element el = (*_ngMesh)[netgen::ElementIndex( _ngMesh->GetNE()-1 )];
3192       curShapeIndex = el.GetIndex();
3193     }
3194     if ( curShapeIndex != _curShapeIndex )
3195     {
3196       ((int&) _curShapeIndex) = curShapeIndex;
3197       double    doneShapeRate = _curShapeIndex / double( _occgeom->somap.Extent() );
3198       double         doneTime = edgeFaceMeshingTime + doneShapeRate * voluMeshingTime;
3199       ((double&)    _ticTime) = doneTime / _totalTime / _progressTic;
3200       // cout << "shape " << _curShapeIndex << " _ticTime " << _ticTime
3201       //      << " " << doneTime / _totalTime / _progressTic << endl;
3202     }
3203   }
3204   if ( _ticTime > 0 )
3205     progress  = Max( *algoProgressTic * _ticTime, *algoProgress );
3206   if ( progress > 0 )
3207   {
3208     ((int&) *algoProgressTic )++;
3209     ((double&) *algoProgress) = progress;
3210   }
3211   //cout << progress << " "  << *algoProgressTic << " " << netgen::multithread.task << " "<< _ticTime << endl;
3212
3213   return Min( progress, 0.99 );
3214 }
3215
3216 //================================================================================
3217 /*!
3218  * \brief Remove "test.out" and "problemfaces" files in current directory
3219  */
3220 //================================================================================
3221
3222 void NETGENPlugin_Mesher::RemoveTmpFiles()
3223 {
3224   bool rm =  SMESH_File("test.out").remove() ;
3225 #ifndef WIN32
3226   if (rm && netgen::testout)
3227   {
3228     delete netgen::testout;
3229     netgen::testout = 0;
3230   }
3231 #endif
3232   SMESH_File("problemfaces").remove();
3233   SMESH_File("occmesh.rep").remove();
3234 }
3235
3236 //================================================================================
3237 /*!
3238  * \brief Read mesh entities preventing successful computation from "test.out" file
3239  */
3240 //================================================================================
3241
3242 SMESH_ComputeErrorPtr
3243 NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
3244 {
3245   SMESH_ComputeErrorPtr err = SMESH_ComputeError::New
3246     (COMPERR_BAD_INPUT_MESH, "Some edges multiple times in surface mesh");
3247   SMESH_File file("test.out");
3248   vector<int> two(2);
3249   const char* badEdgeStr = " multiple times in surface mesh";
3250   const int   badEdgeStrLen = strlen( badEdgeStr );
3251   while( !file.eof() )
3252   {
3253     if ( strncmp( file, "Edge ", 5 ) == 0 &&
3254          file.getInts( two ) &&
3255          strncmp( file, badEdgeStr, badEdgeStrLen ) == 0 &&
3256          two[0] < nodeVec.size() && two[1] < nodeVec.size())
3257     {
3258       err->myBadElements.push_back( new SMDS_LinearEdge( nodeVec[ two[0]], nodeVec[ two[1]] ));
3259       file += badEdgeStrLen;
3260     }
3261     else if ( strncmp( file, "Intersecting: ", 14 ) == 0 )
3262     {
3263 // Intersecting:
3264 // openelement 18 with open element 126
3265 // 41  36  38
3266 // 69  70  72
3267       vector<int> three1(3), three2(3);
3268       file.getLine();
3269       const char* pos = file;
3270       bool ok = ( strncmp( file, "openelement ", 12 ) == 0 );
3271       ok = ok && file.getInts( two );
3272       ok = ok && file.getInts( three1 );
3273       ok = ok && file.getInts( three2 );
3274       for ( int i = 0; ok && i < 3; ++i )
3275         ok = ( three1[i] < nodeVec.size() && nodeVec[ three1[i]]);
3276       for ( int i = 0; ok && i < 3; ++i )
3277         ok = ( three2[i] < nodeVec.size() && nodeVec[ three2[i]]);
3278       if ( ok )
3279       {
3280         err->myBadElements.push_back( new SMDS_FaceOfNodes( nodeVec[ three1[0]],
3281                                                             nodeVec[ three1[1]],
3282                                                             nodeVec[ three1[2]]));
3283         err->myBadElements.push_back( new SMDS_FaceOfNodes( nodeVec[ three2[0]],
3284                                                             nodeVec[ three2[1]],
3285                                                             nodeVec[ three2[2]]));
3286         err->myComment = "Intersecting triangles";
3287       }
3288       else
3289       {
3290         file.setPos( pos );
3291       }
3292     }
3293     else
3294     {
3295       ++file;
3296     }
3297   }
3298   return err;
3299 }
3300
3301 //================================================================================
3302 /*!
3303  * \brief Write a python script creating an equivalent SALOME mesh.
3304  * This is useful to see what mesh is passed as input for the next step of mesh
3305  * generation (of mesh of higher dimension)
3306  */
3307 //================================================================================
3308
3309 void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh,
3310                                     const std::string&  pyFile)
3311 {
3312   ofstream outfile(pyFile.c_str(), ios::out);
3313   if ( !outfile ) return;
3314
3315   outfile << "import SMESH" << endl
3316           << "from salome.smesh import smeshBuilder" << endl
3317           << "smesh = smeshBuilder.New(salome.myStudy)" << endl
3318           << "mesh = smesh.Mesh()" << endl << endl;
3319
3320   using namespace netgen;
3321   PointIndex pi;
3322   for (pi = PointIndex::BASE;
3323        pi < ngMesh->GetNP()+PointIndex::BASE; pi++)
3324   {
3325     outfile << "mesh.AddNode( ";
3326     outfile << (*ngMesh)[pi](0) << ", ";
3327     outfile << (*ngMesh)[pi](1) << ", ";
3328     outfile << (*ngMesh)[pi](2) << ") ## "<< pi << endl;
3329   }
3330
3331   int nbDom = ngMesh->GetNDomains();
3332   for ( int i = 0; i < nbDom; ++i )
3333     outfile<< "grp" << i+1 << " = mesh.CreateEmptyGroup( SMESH.FACE, 'domain"<< i+1 << "')"<< endl;
3334
3335   SurfaceElementIndex sei;
3336   for (sei = 0; sei < ngMesh->GetNSE(); sei++)
3337   {
3338     outfile << "mesh.AddFace([ ";
3339     Element2d sel = (*ngMesh)[sei];
3340     for (int j = 0; j < sel.GetNP(); j++)
3341       outfile << sel[j] << ( j+1 < sel.GetNP() ? ", " : " ])");
3342     if ( sel.IsDeleted() ) outfile << " ## IsDeleted ";
3343     outfile << endl;
3344
3345     if ((*ngMesh)[sei].GetIndex())
3346     {
3347       if ( int dom1 = ngMesh->GetFaceDescriptor((*ngMesh)[sei].GetIndex ()).DomainIn())
3348         outfile << "grp"<< dom1 <<".Add([ " << (int)sei+1 << " ])" << endl;
3349       if ( int dom2 = ngMesh->GetFaceDescriptor((*ngMesh)[sei].GetIndex ()).DomainOut())
3350         outfile << "grp"<< dom2 <<".Add([ " << (int)sei+1 << " ])" << endl;
3351     }
3352   }
3353
3354   for (ElementIndex ei = 0; ei < ngMesh->GetNE(); ei++)
3355   {
3356     Element el = (*ngMesh)[ei];
3357     outfile << "mesh.AddVolume([ ";
3358     for (int j = 0; j < el.GetNP(); j++)
3359       outfile << el[j] << ( j+1 < el.GetNP() ? ", " : " ])");
3360     outfile << endl;
3361   }
3362
3363   for (int i = 1; i <= ngMesh->GetNSeg(); i++)
3364   {
3365     const Segment & seg = ngMesh->LineSegment (i);
3366     outfile << "mesh.AddEdge([ "
3367             << seg[0] << ", "
3368             << seg[1] << " ])" << endl;
3369   }
3370   cout << "Write " << pyFile << endl;
3371 }
3372
3373 //================================================================================
3374 /*!
3375  * \brief Constructor of NETGENPlugin_ngMeshInfo
3376  */
3377 //================================================================================
3378
3379 NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh):
3380   _copyOfLocalH(0)
3381 {
3382   if ( ngMesh )
3383   {
3384     _nbNodes    = ngMesh->GetNP();
3385     _nbSegments = ngMesh->GetNSeg();
3386     _nbFaces    = ngMesh->GetNSE();
3387     _nbVolumes  = ngMesh->GetNE();
3388   }
3389   else
3390   {
3391     _nbNodes = _nbSegments = _nbFaces = _nbVolumes = 0;
3392   }
3393 }
3394
3395 //================================================================================
3396 /*!
3397  * \brief Copy LocalH member from one netgen mesh to another
3398  */
3399 //================================================================================
3400
3401 void NETGENPlugin_ngMeshInfo::transferLocalH( netgen::Mesh* fromMesh,
3402                                               netgen::Mesh* toMesh )
3403 {
3404   if ( !fromMesh->LocalHFunctionGenerated() ) return;
3405   if ( !toMesh->LocalHFunctionGenerated() )
3406 #ifdef NETGEN_V5
3407     toMesh->CalcLocalH(netgen::mparam.grading);
3408 #else
3409     toMesh->CalcLocalH();
3410 #endif
3411
3412   const size_t size = sizeof( netgen::LocalH );
3413   _copyOfLocalH = new char[ size ];
3414   memcpy( (void*)_copyOfLocalH, (void*)&toMesh->LocalHFunction(), size );
3415   memcpy( (void*)&toMesh->LocalHFunction(), (void*)&fromMesh->LocalHFunction(), size );
3416 }
3417
3418 //================================================================================
3419 /*!
3420  * \brief Restore LocalH member of a netgen mesh
3421  */
3422 //================================================================================
3423
3424 void NETGENPlugin_ngMeshInfo::restoreLocalH( netgen::Mesh* toMesh )
3425 {
3426   if ( _copyOfLocalH )
3427   {
3428     const size_t size = sizeof( netgen::LocalH );
3429     memcpy( (void*)&toMesh->LocalHFunction(), (void*)_copyOfLocalH, size );
3430     delete [] _copyOfLocalH;
3431     _copyOfLocalH = 0;
3432   }
3433 }
3434
3435 //================================================================================
3436 /*!
3437  * \brief Find "internal" sub-shapes
3438  */
3439 //================================================================================
3440
3441 NETGENPlugin_Internals::NETGENPlugin_Internals( SMESH_Mesh&         mesh,
3442                                                 const TopoDS_Shape& shape,
3443                                                 bool                is3D )
3444   : _mesh( mesh ), _is3D( is3D )
3445 {
3446   SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
3447
3448   TopExp_Explorer f,e;
3449   for ( f.Init( shape, TopAbs_FACE ); f.More(); f.Next() )
3450   {
3451     int faceID = meshDS->ShapeToIndex( f.Current() );
3452
3453     // find not computed internal edges
3454
3455     for ( e.Init( f.Current().Oriented(TopAbs_FORWARD), TopAbs_EDGE ); e.More(); e.Next() )
3456       if ( e.Current().Orientation() == TopAbs_INTERNAL )
3457       {
3458         SMESH_subMesh* eSM = mesh.GetSubMesh( e.Current() );
3459         if ( eSM->IsEmpty() )
3460         {
3461           _e2face.insert( make_pair( eSM->GetId(), faceID ));
3462           for ( TopoDS_Iterator v(e.Current()); v.More(); v.Next() )
3463             _e2face.insert( make_pair( meshDS->ShapeToIndex( v.Value() ), faceID ));
3464         }
3465       }
3466
3467     // find internal vertices in a face
3468     set<int> intVV; // issue 0020850 where same vertex is twice in a face
3469     for ( TopoDS_Iterator fSub( f.Current() ); fSub.More(); fSub.Next())
3470       if ( fSub.Value().ShapeType() == TopAbs_VERTEX )
3471       {
3472         int vID = meshDS->ShapeToIndex( fSub.Value() );
3473         if ( intVV.insert( vID ).second )
3474           _f2v[ faceID ].push_back( vID );
3475       }
3476
3477     if ( is3D )
3478     {
3479       // find internal faces and their subshapes where nodes are to be doubled
3480       //  to make a crack with non-sewed borders
3481
3482       if ( f.Current().Orientation() == TopAbs_INTERNAL )
3483       {
3484         _intShapes.insert( meshDS->ShapeToIndex( f.Current() ));
3485
3486         // egdes
3487         list< TopoDS_Shape > edges;
3488         for ( e.Init( f.Current(), TopAbs_EDGE ); e.More(); e.Next())
3489           if ( SMESH_MesherHelper::NbAncestors( e.Current(), mesh, TopAbs_FACE ) > 1 )
3490           {
3491             _intShapes.insert( meshDS->ShapeToIndex( e.Current() ));
3492             edges.push_back( e.Current() );
3493             // find border faces
3494             PShapeIteratorPtr fIt =
3495               SMESH_MesherHelper::GetAncestors( edges.back(),mesh,TopAbs_FACE );
3496             while ( const TopoDS_Shape* pFace = fIt->next() )
3497               if ( !pFace->IsSame( f.Current() ))
3498                 _borderFaces.insert( meshDS->ShapeToIndex( *pFace ));
3499           }
3500         // vertices
3501         // we consider vertex internal if it is shared by more than one internal edge
3502         list< TopoDS_Shape >::iterator edge = edges.begin();
3503         for ( ; edge != edges.end(); ++edge )
3504           for ( TopoDS_Iterator v( *edge ); v.More(); v.Next() )
3505           {
3506             set<int> internalEdges;
3507             PShapeIteratorPtr eIt =
3508               SMESH_MesherHelper::GetAncestors( v.Value(),mesh,TopAbs_EDGE );
3509             while ( const TopoDS_Shape* pEdge = eIt->next() )
3510             {
3511               int edgeID = meshDS->ShapeToIndex( *pEdge );
3512               if ( isInternalShape( edgeID ))
3513                 internalEdges.insert( edgeID );
3514             }
3515             if ( internalEdges.size() > 1 )
3516               _intShapes.insert( meshDS->ShapeToIndex( v.Value() ));
3517           }
3518       }
3519     }
3520   } // loop on geom faces
3521
3522   // find vertices internal in solids
3523   if ( is3D )
3524   {
3525     for ( TopExp_Explorer so(shape, TopAbs_SOLID); so.More(); so.Next())
3526     {
3527       int soID = meshDS->ShapeToIndex( so.Current() );
3528       for ( TopoDS_Iterator soSub( so.Current() ); soSub.More(); soSub.Next())
3529         if ( soSub.Value().ShapeType() == TopAbs_VERTEX )
3530           _s2v[ soID ].push_back( meshDS->ShapeToIndex( soSub.Value() ));
3531     }
3532   }
3533 }
3534
3535 //================================================================================
3536 /*!
3537  * \brief Find mesh faces on non-internal geom faces sharing internal edge
3538  * some nodes of which are to be doubled to make the second border of the "crack"
3539  */
3540 //================================================================================
3541
3542 void NETGENPlugin_Internals::findBorderElements( TIDSortedElemSet & borderElems )
3543 {
3544   if ( _intShapes.empty() ) return;
3545
3546   SMESH_Mesh& mesh = const_cast<SMESH_Mesh&>(_mesh);
3547   SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
3548
3549   // loop on internal geom edges
3550   set<int>::const_iterator intShapeId = _intShapes.begin();
3551   for ( ; intShapeId != _intShapes.end(); ++intShapeId )
3552   {
3553     const TopoDS_Shape& s = meshDS->IndexToShape( *intShapeId );
3554     if ( s.ShapeType() != TopAbs_EDGE ) continue;
3555
3556     // get internal and non-internal geom faces sharing the internal edge <s>
3557     int intFace = 0;
3558     set<int>::iterator bordFace = _borderFaces.end();
3559     PShapeIteratorPtr faces = SMESH_MesherHelper::GetAncestors( s, _mesh, TopAbs_FACE );
3560     while ( const TopoDS_Shape* pFace = faces->next() )
3561     {
3562       int faceID = meshDS->ShapeToIndex( *pFace );
3563       if ( isInternalShape( faceID ))
3564         intFace = faceID;
3565       else
3566         bordFace = _borderFaces.insert( faceID ).first;
3567     }
3568     if ( bordFace == _borderFaces.end() || !intFace ) continue;
3569
3570     // get all links of mesh faces on internal geom face sharing nodes on edge <s>
3571     set< SMESH_OrientedLink > links; //!< links of faces on internal geom face
3572     list<const SMDS_MeshElement*> suspectFaces[2]; //!< mesh faces on border geom faces
3573     int nbSuspectFaces = 0;
3574     SMESHDS_SubMesh* intFaceSM = meshDS->MeshElements( intFace );
3575     if ( !intFaceSM || intFaceSM->NbElements() == 0 ) continue;
3576     SMESH_subMeshIteratorPtr smIt = mesh.GetSubMesh( s )->getDependsOnIterator(true,true);
3577     while ( smIt->more() )
3578     {
3579       SMESHDS_SubMesh* sm = smIt->next()->GetSubMeshDS();
3580       if ( !sm ) continue;
3581       SMDS_NodeIteratorPtr nIt = sm->GetNodes();
3582       while ( nIt->more() )
3583       {
3584         const SMDS_MeshNode* nOnEdge = nIt->next();
3585         SMDS_ElemIteratorPtr fIt = nOnEdge->GetInverseElementIterator(SMDSAbs_Face);
3586         while ( fIt->more() )
3587         {
3588           const SMDS_MeshElement* f = fIt->next();
3589           int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
3590           if ( intFaceSM->Contains( f ))
3591           {
3592             for ( int i = 0; i < nbNodes; ++i )
3593               links.insert( SMESH_OrientedLink( f->GetNode(i), f->GetNode((i+1)%nbNodes)));
3594           }
3595           else
3596           {
3597             int nbDblNodes = 0;
3598             for ( int i = 0; i < nbNodes; ++i )
3599               nbDblNodes += isInternalShape( f->GetNode(i)->getshapeId() );
3600             if ( nbDblNodes )
3601               suspectFaces[ nbDblNodes < 2 ].push_back( f );
3602             nbSuspectFaces++;
3603           }
3604         }
3605       }
3606     }
3607     // suspectFaces[0] having link with same orientation as mesh faces on
3608     // the internal geom face are <borderElems>. suspectFaces[1] have
3609     // only one node on edge <s>, we decide on them later (at the 2nd loop)
3610     // by links of <borderElems> found at the 1st and 2nd loops
3611     set< SMESH_OrientedLink > borderLinks;
3612     for ( int isPostponed = 0; isPostponed < 2; ++isPostponed )
3613     {
3614       list<const SMDS_MeshElement*>::iterator fIt = suspectFaces[isPostponed].begin();
3615       for ( int nbF = 0; fIt != suspectFaces[isPostponed].end(); ++fIt, ++nbF )
3616       {
3617         const SMDS_MeshElement* f = *fIt;
3618         bool isBorder = false, linkFound = false, borderLinkFound = false;
3619         list< SMESH_OrientedLink > faceLinks;
3620         int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
3621         for ( int i = 0; i < nbNodes; ++i )
3622         {
3623           SMESH_OrientedLink link( f->GetNode(i), f->GetNode((i+1)%nbNodes));
3624           faceLinks.push_back( link );
3625           if ( !linkFound )
3626           {
3627             set< SMESH_OrientedLink >::iterator foundLink = links.find( link );
3628             if ( foundLink != links.end() )
3629             {
3630               linkFound= true;
3631               isBorder = ( foundLink->_reversed == link._reversed );
3632               if ( !isBorder && !isPostponed ) break;
3633               faceLinks.pop_back();
3634             }
3635             else if ( isPostponed && !borderLinkFound )
3636             {
3637               foundLink = borderLinks.find( link );
3638               if ( foundLink != borderLinks.end() )
3639               {
3640                 borderLinkFound = true;
3641                 isBorder = ( foundLink->_reversed != link._reversed );
3642               }
3643             }
3644           }
3645         }
3646         if ( isBorder )
3647         {
3648           borderElems.insert( f );
3649           borderLinks.insert( faceLinks.begin(), faceLinks.end() );
3650         }
3651         else if ( !linkFound && !borderLinkFound )
3652         {
3653           suspectFaces[1].push_back( f );
3654           if ( nbF > 2 * nbSuspectFaces )
3655             break; // dead loop protection
3656         }
3657       }
3658     }
3659   }
3660 }
3661
3662 //================================================================================
3663 /*!
3664  * \brief put internal shapes in maps and fill in submeshes to precompute
3665  */
3666 //================================================================================
3667
3668 void NETGENPlugin_Internals::getInternalEdges( TopTools_IndexedMapOfShape& fmap,
3669                                                TopTools_IndexedMapOfShape& emap,
3670                                                TopTools_IndexedMapOfShape& vmap,
3671                                                list< SMESH_subMesh* > smToPrecompute[])
3672 {
3673   if ( !hasInternalEdges() ) return;
3674   map<int,int>::const_iterator ev_face = _e2face.begin();
3675   for ( ; ev_face != _e2face.end(); ++ev_face )
3676   {
3677     const TopoDS_Shape& ev   = _mesh.GetMeshDS()->IndexToShape( ev_face->first );
3678     const TopoDS_Shape& face = _mesh.GetMeshDS()->IndexToShape( ev_face->second );
3679
3680     ( ev.ShapeType() == TopAbs_EDGE ? emap : vmap ).Add( ev );
3681     fmap.Add( face );
3682     //cout<<"INTERNAL EDGE or VERTEX "<<ev_face->first<<" on face "<<ev_face->second<<endl;
3683
3684     smToPrecompute[ MeshDim_1D ].push_back( _mesh.GetSubMeshContaining( ev_face->first ));
3685   }
3686 }
3687
3688 //================================================================================
3689 /*!
3690  * \brief return shapes and submeshes to be meshed and already meshed boundary submeshes
3691  */
3692 //================================================================================
3693
3694 void NETGENPlugin_Internals::getInternalFaces( TopTools_IndexedMapOfShape& fmap,
3695                                                TopTools_IndexedMapOfShape& emap,
3696                                                list< SMESH_subMesh* >&     intFaceSM,
3697                                                list< SMESH_subMesh* >&     boundarySM)
3698 {
3699   if ( !hasInternalFaces() ) return;
3700
3701   // <fmap> and <emap> are for not yet meshed shapes
3702   // <intFaceSM> is for submeshes of faces
3703   // <boundarySM> is for meshed edges and vertices
3704
3705   intFaceSM.clear();
3706   boundarySM.clear();
3707
3708   set<int> shapeIDs ( _intShapes );
3709   if ( !_borderFaces.empty() )
3710     shapeIDs.insert( _borderFaces.begin(), _borderFaces.end() );
3711
3712   set<int>::const_iterator intS = shapeIDs.begin();
3713   for ( ; intS != shapeIDs.end(); ++intS )
3714   {
3715     SMESH_subMesh* sm = _mesh.GetSubMeshContaining( *intS );
3716
3717     if ( sm->GetSubShape().ShapeType() != TopAbs_FACE ) continue;
3718
3719     intFaceSM.push_back( sm );
3720
3721     // add submeshes of not computed internal faces
3722     if ( !sm->IsEmpty() ) continue;
3723
3724     SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(true,true);
3725     while ( smIt->more() )
3726     {
3727       sm = smIt->next();
3728       const TopoDS_Shape& s = sm->GetSubShape();
3729
3730       if ( sm->IsEmpty() )
3731       {
3732         // not yet meshed
3733         switch ( s.ShapeType() ) {
3734         case TopAbs_FACE: fmap.Add ( s ); break;
3735         case TopAbs_EDGE: emap.Add ( s ); break;
3736         default:;
3737         }
3738       }
3739       else
3740       {
3741         if ( s.ShapeType() != TopAbs_FACE )
3742           boundarySM.push_back( sm );
3743       }
3744     }
3745   }
3746 }
3747
3748 //================================================================================
3749 /*!
3750  * \brief Return true if given shape is to be precomputed in order to be correctly
3751  * added to netgen mesh
3752  */
3753 //================================================================================
3754
3755 bool NETGENPlugin_Internals::isShapeToPrecompute(const TopoDS_Shape& s)
3756 {
3757   int shapeID = _mesh.GetMeshDS()->ShapeToIndex( s );
3758   switch ( s.ShapeType() ) {
3759   case TopAbs_FACE  : break; //return isInternalShape( shapeID ) || isBorderFace( shapeID );
3760   case TopAbs_EDGE  : return isInternalEdge( shapeID );
3761   case TopAbs_VERTEX: break;
3762   default:;
3763   }
3764   return false;
3765 }
3766
3767 //================================================================================
3768 /*!
3769  * \brief Return SMESH
3770  */
3771 //================================================================================
3772
3773 SMESH_Mesh& NETGENPlugin_Internals::getMesh() const
3774 {
3775   return const_cast<SMESH_Mesh&>( _mesh );
3776 }
3777
3778 //================================================================================
3779 /*!
3780  * \brief Initialize netgen library
3781  */
3782 //================================================================================
3783
3784 NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper()
3785 {
3786   Ng_Init();
3787
3788   _isComputeOk      = false;
3789   _coutBuffer       = NULL;
3790   if ( !getenv( "KEEP_NETGEN_OUTPUT" ))
3791   {
3792     // redirect all netgen output (mycout,myerr,cout) to _outputFileName
3793     _outputFileName = getOutputFileName();
3794     netgen::mycout  = new ofstream ( _outputFileName.c_str() );
3795     netgen::myerr   = netgen::mycout;
3796     _coutBuffer     = std::cout.rdbuf();
3797 #ifdef _DEBUG_
3798     cout << "NOTE: netgen output is redirected to file " << _outputFileName << endl;
3799 #else
3800     std::cout.rdbuf( netgen::mycout->rdbuf() );
3801 #endif
3802   }
3803
3804   _ngMesh = Ng_NewMesh();
3805 }
3806
3807 //================================================================================
3808 /*!
3809  * \brief Finish using netgen library
3810  */
3811 //================================================================================
3812
3813 NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper()
3814 {
3815   Ng_DeleteMesh( _ngMesh );
3816   Ng_Exit();
3817   NETGENPlugin_Mesher::RemoveTmpFiles();
3818   if ( _coutBuffer )
3819     std::cout.rdbuf( _coutBuffer );
3820 #ifdef _DEBUG_
3821   if( _isComputeOk )
3822 #endif
3823     removeOutputFile();
3824 }
3825
3826 //================================================================================
3827 /*!
3828  * \brief Set netgen mesh to delete at destruction
3829  */
3830 //================================================================================
3831
3832 void NETGENPlugin_NetgenLibWrapper::setMesh( Ng_Mesh* mesh )
3833 {
3834   if ( _ngMesh )
3835     Ng_DeleteMesh( _ngMesh );
3836   _ngMesh = mesh;
3837 }
3838
3839 //================================================================================
3840 /*!
3841  * \brief Return a unique file name
3842  */
3843 //================================================================================
3844
3845 std::string NETGENPlugin_NetgenLibWrapper::getOutputFileName()
3846 {
3847   std::string aTmpDir = SALOMEDS_Tool::GetTmpDir();
3848
3849   TCollection_AsciiString aGenericName = (char*)aTmpDir.c_str();
3850   aGenericName += "NETGEN_";
3851 #ifndef WIN32
3852   aGenericName += getpid();
3853 #else
3854   aGenericName += _getpid();
3855 #endif
3856   aGenericName += "_";
3857   aGenericName += Abs((Standard_Integer)(long) aGenericName.ToCString());
3858   aGenericName += ".out";
3859
3860   return aGenericName.ToCString();
3861 }
3862
3863 //================================================================================
3864 /*!
3865  * \brief Remove file with netgen output
3866  */
3867 //================================================================================
3868
3869 void NETGENPlugin_NetgenLibWrapper::removeOutputFile()
3870 {
3871   if ( !_outputFileName.empty() )
3872   {
3873     if ( netgen::mycout )
3874     {
3875       delete netgen::mycout;
3876       netgen::mycout = 0;
3877       netgen::myerr = 0;
3878     }
3879     string    tmpDir = SALOMEDS_Tool::GetDirFromPath ( _outputFileName );
3880     string aFileName = SALOMEDS_Tool::GetNameFromPath( _outputFileName ) + ".out";
3881     SALOMEDS::ListOfFileNames_var aFiles = new SALOMEDS::ListOfFileNames;
3882     aFiles->length(1);
3883     aFiles[0] = aFileName.c_str();
3884
3885     SALOMEDS_Tool::RemoveTemporaryFiles( tmpDir.c_str(), aFiles.in(), true );
3886   }
3887 }