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