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 sunH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ];
1676         int   nbSeg = ( int( segLen[ iPrev ] > sunH / 100.)  +
1677                         int( segLen[ i     ] > sunH / 100.)  +
1678                         int( segLen[ iNext ] > sunH / 100.));
1679         if ( nbSeg > 0 )
1680           RestrictLocalSize( ngMesh, 0.5*(np1+np2), sunH / 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   const double edgeMeshingTime = 0.001;
2099   const double faceMeshingTime = 0.019;
2100   const double edgeFaceMeshingTime = edgeMeshingTime + faceMeshingTime;
2101   const double faceOptimizTime = 0.06;
2102   const double voluMeshingTime = 0.15;
2103   const double volOptimizeTime = 0.77;
2104 }
2105
2106 //=============================================================================
2107 /*!
2108  * Here we are going to use the NETGEN mesher
2109  */
2110 //=============================================================================
2111
2112 bool NETGENPlugin_Mesher::Compute()
2113 {
2114   NETGENPlugin_NetgenLibWrapper ngLib;
2115
2116   netgen::MeshingParameters& mparams = netgen::mparam;
2117   MESSAGE("Compute with:\n"
2118           " max size = " << mparams.maxh << "\n"
2119           " segments per edge = " << mparams.segmentsperedge);
2120   MESSAGE("\n"
2121           " growth rate = " << mparams.grading << "\n"
2122           " elements per radius = " << mparams.curvaturesafety << "\n"
2123           " second order = " << mparams.secondorder << "\n"
2124           " quad allowed = " << mparams.quad << "\n"
2125           " surface curvature = " << mparams.uselocalh << "\n"
2126           " fuse edges = " << netgen::merge_solids);
2127
2128   SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
2129
2130   static string debugFile = "/tmp/ngMesh.py"; /* to call toPython( ngMesh, debugFile )
2131                                                  while debugging netgen */
2132   // -------------------------
2133   // Prepare OCC geometry
2134   // -------------------------
2135
2136   netgen::OCCGeometry occgeo;
2137   list< SMESH_subMesh* > meshedSM[3]; // for 0-2 dimensions
2138   NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume );
2139   PrepareOCCgeometry( occgeo, _shape, *_mesh, meshedSM, &internals );
2140   _occgeom = &occgeo;
2141
2142   _totalTime = edgeFaceMeshingTime;
2143   if ( _optimize )
2144     _totalTime += faceOptimizTime;
2145   if ( _isVolume )
2146     _totalTime += voluMeshingTime + ( _optimize ? volOptimizeTime : 0 );
2147   double doneTime = 0;
2148   _ticTime = -1;
2149   _progressTic = 1;
2150   _curShapeIndex = -1;
2151
2152   // -------------------------
2153   // Generate the mesh
2154   // -------------------------
2155
2156   _ngMesh = NULL;
2157   NETGENPlugin_ngMeshInfo initState; // it remembers size of ng mesh equal to size of Smesh
2158
2159   SMESH_Comment comment;
2160   int err = 0;
2161
2162   // vector of nodes in which node index == netgen ID
2163   vector< const SMDS_MeshNode* > nodeVec;
2164
2165   {
2166     // ----------------
2167     // compute 1D mesh
2168     // ----------------
2169     if ( _simpleHyp )
2170     {
2171       // not to RestrictLocalH() according to curvature during MESHCONST_ANALYSE
2172       mparams.uselocalh = false;
2173       mparams.grading = 0.8; // not limitited size growth
2174
2175       if ( _simpleHyp->GetNumberOfSegments() )
2176         // nb of segments
2177         mparams.maxh = occgeo.boundingbox.Diam();
2178       else
2179         // segment length
2180         mparams.maxh = _simpleHyp->GetLocalLength();
2181     }
2182
2183     if ( mparams.maxh == 0.0 )
2184       mparams.maxh = occgeo.boundingbox.Diam();
2185     if ( _simpleHyp || ( mparams.minh == 0.0 && _fineness != NETGENPlugin_Hypothesis::UserDefined))
2186       mparams.minh = GetDefaultMinSize( _shape, mparams.maxh );
2187
2188     // Local size on faces
2189     occgeo.face_maxh = mparams.maxh;
2190
2191     // Let netgen create _ngMesh and calculate element size on not meshed shapes
2192 #ifndef NETGEN_V5
2193     char *optstr = 0;
2194 #endif
2195     int startWith = netgen::MESHCONST_ANALYSE;
2196     int endWith   = netgen::MESHCONST_ANALYSE;
2197     try
2198     {
2199       OCC_CATCH_SIGNALS;
2200 #ifdef NETGEN_V5
2201       err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
2202 #else
2203       err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
2204 #endif
2205       if(netgen::multithread.terminate)
2206         return false;
2207
2208       comment << text(err);
2209     }
2210     catch (Standard_Failure& ex)
2211     {
2212       comment << text(ex);
2213     }
2214     err = 0; //- MESHCONST_ANALYSE isn't so important step
2215     if ( !_ngMesh )
2216       return false;
2217     ngLib.setMesh(( Ng_Mesh*) _ngMesh );
2218
2219     _ngMesh->ClearFaceDescriptors(); // we make descriptors our-self
2220
2221     if ( _simpleHyp )
2222     {
2223       // Pass 1D simple parameters to NETGEN
2224       // --------------------------------
2225       int      nbSeg = _simpleHyp->GetNumberOfSegments();
2226       double segSize = _simpleHyp->GetLocalLength();
2227       for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE )
2228       {
2229         const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE));
2230         if ( nbSeg )
2231           segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 );
2232         setLocalSize( e, segSize, *_ngMesh );
2233       }
2234     }
2235     else // if ( ! _simpleHyp )
2236     {
2237       // Local size on vertices and edges
2238       // --------------------------------
2239       for(std::map<int,double>::const_iterator it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++)
2240       {
2241         int key = (*it).first;
2242         double hi = (*it).second;
2243         const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2244         const TopoDS_Edge& e = TopoDS::Edge(shape);
2245         setLocalSize( e, hi, *_ngMesh );
2246       }
2247       for(std::map<int,double>::const_iterator it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++)
2248       {
2249         int key = (*it).first;
2250         double hi = (*it).second;
2251         const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2252         const TopoDS_Vertex& v = TopoDS::Vertex(shape);
2253         gp_Pnt p = BRep_Tool::Pnt(v);
2254         NETGENPlugin_Mesher::RestrictLocalSize( *_ngMesh, p.XYZ(), hi );
2255       }
2256       for(map<int,double>::const_iterator it=FaceId2LocalSize.begin();
2257           it!=FaceId2LocalSize.end(); it++)
2258       {
2259         int key = (*it).first;
2260         double val = (*it).second;
2261         const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2262         int faceNgID = occgeo.fmap.FindIndex(shape);
2263         occgeo.SetFaceMaxH(faceNgID, val);
2264         for ( TopExp_Explorer edgeExp( shape, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
2265           setLocalSize( TopoDS::Edge( edgeExp.Current() ), val, *_ngMesh );
2266       }
2267     }
2268
2269     // Precompute internal edges (issue 0020676) in order to
2270     // add mesh on them correctly (twice) to netgen mesh
2271     if ( !err && internals.hasInternalEdges() )
2272     {
2273       // load internal shapes into OCCGeometry
2274       netgen::OCCGeometry intOccgeo;
2275       internals.getInternalEdges( intOccgeo.fmap, intOccgeo.emap, intOccgeo.vmap, meshedSM );
2276       intOccgeo.boundingbox = occgeo.boundingbox;
2277       intOccgeo.shape = occgeo.shape;
2278       intOccgeo.face_maxh.SetSize(intOccgeo.fmap.Extent());
2279       intOccgeo.face_maxh = netgen::mparam.maxh;
2280       netgen::Mesh *tmpNgMesh = NULL;
2281       try
2282       {
2283         OCC_CATCH_SIGNALS;
2284         // compute local H on internal shapes in the main mesh
2285         //OCCSetLocalMeshSize(intOccgeo, *_ngMesh); it deletes _ngMesh->localH
2286
2287         // let netgen create a temporary mesh
2288 #ifdef NETGEN_V5
2289         netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, mparams, startWith, endWith);
2290 #else
2291         netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, startWith, endWith, optstr);
2292 #endif
2293         if(netgen::multithread.terminate)
2294           return false;
2295
2296         // copy LocalH from the main to temporary mesh
2297         initState.transferLocalH( _ngMesh, tmpNgMesh );
2298
2299         // compute mesh on internal edges
2300         startWith = endWith = netgen::MESHCONST_MESHEDGES;
2301 #ifdef NETGEN_V5
2302         err = netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, mparams, startWith, endWith);
2303 #else
2304         err = netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, startWith, endWith, optstr);
2305 #endif
2306         comment << text(err);
2307       }
2308       catch (Standard_Failure& ex)
2309       {
2310         comment << text(ex);
2311         err = 1;
2312       }
2313       initState.restoreLocalH( tmpNgMesh );
2314
2315       // fill SMESH by netgen mesh
2316       vector< const SMDS_MeshNode* > tmpNodeVec;
2317       FillSMesh( intOccgeo, *tmpNgMesh, initState, *_mesh, tmpNodeVec, comment );
2318       err = ( err || !comment.empty() );
2319
2320       nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)tmpNgMesh);
2321     }
2322
2323     // Fill _ngMesh with nodes and segments of computed submeshes
2324     if ( !err )
2325     {
2326       err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_0D ]) &&
2327                 FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_1D ]));
2328     }
2329     initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2330
2331     // Compute 1d mesh
2332     if (!err)
2333     {
2334       startWith = endWith = netgen::MESHCONST_MESHEDGES;
2335       try
2336       {
2337         OCC_CATCH_SIGNALS;
2338 #ifdef NETGEN_V5
2339         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
2340 #else
2341         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
2342 #endif
2343         if(netgen::multithread.terminate)
2344           return false;
2345
2346         comment << text(err);
2347       }
2348       catch (Standard_Failure& ex)
2349       {
2350         comment << text(ex);
2351         err = 1;
2352       }
2353     }
2354     if ( _isVolume )
2355       _ticTime = ( doneTime += edgeMeshingTime ) / _totalTime / _progressTic;
2356
2357     mparams.uselocalh = true; // restore as it is used at surface optimization
2358
2359     // ---------------------
2360     // compute surface mesh
2361     // ---------------------
2362     if (!err)
2363     {
2364       // Pass 2D simple parameters to NETGEN
2365       if ( _simpleHyp ) {
2366         if ( double area = _simpleHyp->GetMaxElementArea() ) {
2367           // face area
2368           mparams.maxh = sqrt(2. * area/sqrt(3.0));
2369           mparams.grading = 0.4; // moderate size growth
2370         }
2371         else {
2372           // length from edges
2373           if ( _ngMesh->GetNSeg() ) {
2374             double edgeLength = 0;
2375             TopTools_MapOfShape visitedEdges;
2376             for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() )
2377               if( visitedEdges.Add(exp.Current()) )
2378                 edgeLength += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() ));
2379             // we have to multiply length by 2 since for each TopoDS_Edge there
2380             // are double set of NETGEN edges, in other words, we have to
2381             // divide _ngMesh->GetNSeg() by 2.
2382             mparams.maxh = 2*edgeLength / _ngMesh->GetNSeg();
2383           }
2384           else {
2385             mparams.maxh = 1000;
2386           }
2387           mparams.grading = 0.2; // slow size growth
2388         }
2389         mparams.quad = _simpleHyp->GetAllowQuadrangles();
2390         mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
2391         _ngMesh->SetGlobalH (mparams.maxh);
2392         netgen::Box<3> bb = occgeo.GetBoundingBox();
2393         bb.Increase (bb.Diam()/20);
2394         _ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading);
2395       }
2396
2397       // Care of vertices internal in faces (issue 0020676)
2398       if ( internals.hasInternalVertexInFace() )
2399       {
2400         // store computed segments in SMESH in order not to create SMESH
2401         // edges for ng segments added by AddIntVerticesInFaces()
2402         FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
2403         // add segments to faces with internal vertices
2404         AddIntVerticesInFaces( occgeo, *_ngMesh, nodeVec, internals );
2405         initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2406       }
2407
2408       // Build viscous layers
2409       if ( _isViscousLayers2D )
2410       {
2411         if ( !internals.hasInternalVertexInFace() ) {
2412           FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
2413           initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2414         }
2415         SMESH_ProxyMesh::Ptr viscousMesh;
2416         SMESH_MesherHelper   helper( *_mesh );
2417         for ( int faceID = 1; faceID <= occgeo.fmap.Extent(); ++faceID )
2418         {
2419           const TopoDS_Face& F = TopoDS::Face( occgeo.fmap( faceID ));
2420           viscousMesh = StdMeshers_ViscousLayers2D::Compute( *_mesh, F );
2421           if ( !viscousMesh )
2422             return false;
2423           // exclude from computation ng segments built on EDGEs of F
2424           for (int i = 1; i <= _ngMesh->GetNSeg(); i++)
2425           {
2426             netgen::Segment & seg = _ngMesh->LineSegment(i);
2427             if (seg.si == faceID)
2428               seg.si = 0;
2429           }
2430           // add new segments to _ngMesh instead of excluded ones
2431           helper.SetSubShape( F );
2432           TSideVector wires =
2433             StdMeshers_FaceSide::GetFaceWires( F, *_mesh, /*skipMediumNodes=*/true,
2434                                                error, viscousMesh );
2435           error = AddSegmentsToMesh( *_ngMesh, occgeo, wires, helper, nodeVec );
2436
2437           if ( !error ) error = SMESH_ComputeError::New();
2438         }
2439         initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2440       }
2441
2442       // Let netgen compute 2D mesh
2443       startWith = netgen::MESHCONST_MESHSURFACE;
2444       endWith = _optimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
2445       try
2446       {
2447         OCC_CATCH_SIGNALS;
2448 #ifdef NETGEN_V5
2449         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
2450 #else
2451         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
2452 #endif
2453         if(netgen::multithread.terminate)
2454           return false;
2455
2456         comment << text (err);
2457       }
2458       catch (Standard_Failure& ex)
2459       {
2460         comment << text(ex);
2461         //err = 1; -- try to make volumes anyway
2462       }
2463       catch (netgen::NgException exc)
2464       {
2465         comment << text(exc);
2466         //err = 1; -- try to make volumes anyway
2467       }
2468     }
2469     if ( _isVolume )
2470     {
2471       doneTime += faceMeshingTime + ( _optimize ? faceOptimizTime : 0 );
2472       _ticTime = doneTime / _totalTime / _progressTic;
2473     }
2474     // ---------------------
2475     // generate volume mesh
2476     // ---------------------
2477     // Fill _ngMesh with nodes and faces of computed 2D submeshes
2478     if ( !err && _isVolume && ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad ))
2479     {
2480       // load SMESH with computed segments and faces
2481       FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
2482
2483       // compute pyramids on quadrangles
2484       SMESH_ProxyMesh::Ptr proxyMesh;
2485       if ( _mesh->NbQuadrangles() > 0 )
2486         for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS )
2487         {
2488           StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
2489           proxyMesh.reset( Adaptor );
2490
2491           int nbPyrams = _mesh->NbPyramids();
2492           Adaptor->Compute( *_mesh, occgeo.somap(iS) );
2493           if ( nbPyrams != _mesh->NbPyramids() )
2494           {
2495             list< SMESH_subMesh* > quadFaceSM;
2496             for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next())
2497               if ( Adaptor->GetProxySubMesh( face.Current() ))
2498               {
2499                 quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() ));
2500                 meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() );
2501               }
2502             FillNgMesh(occgeo, *_ngMesh, nodeVec, quadFaceSM, proxyMesh);
2503           }
2504         }
2505       // fill _ngMesh with faces of sub-meshes
2506       err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ]));
2507       initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2508       //toPython( _ngMesh, "/tmp/ngPython.py");
2509     }
2510     if (!err && _isVolume)
2511     {
2512       // Pass 3D simple parameters to NETGEN
2513       const NETGENPlugin_SimpleHypothesis_3D* simple3d =
2514         dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
2515       if ( simple3d ) {
2516         if ( double vol = simple3d->GetMaxElementVolume() ) {
2517           // max volume
2518           mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
2519           mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
2520         }
2521         else {
2522           // length from faces
2523           mparams.maxh = _ngMesh->AverageH();
2524         }
2525         _ngMesh->SetGlobalH (mparams.maxh);
2526         mparams.grading = 0.4;
2527 #ifdef NETGEN_V5
2528         _ngMesh->CalcLocalH(mparams.grading);
2529 #else
2530         _ngMesh->CalcLocalH();
2531 #endif
2532       }
2533       // Care of vertices internal in solids and internal faces (issue 0020676)
2534       if ( internals.hasInternalVertexInSolid() || internals.hasInternalFaces() )
2535       {
2536         // store computed faces in SMESH in order not to create SMESH
2537         // faces for ng faces added here
2538         FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment );
2539         // add ng faces to solids with internal vertices
2540         AddIntVerticesInSolids( occgeo, *_ngMesh, nodeVec, internals );
2541         // duplicate mesh faces on internal faces
2542         FixIntFaces( occgeo, *_ngMesh, internals );
2543         initState = NETGENPlugin_ngMeshInfo(_ngMesh);
2544       }
2545       // Let netgen compute 3D mesh
2546       startWith = endWith = netgen::MESHCONST_MESHVOLUME;
2547       try
2548       {
2549         OCC_CATCH_SIGNALS;
2550 #ifdef NETGEN_V5
2551         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
2552 #else
2553         err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
2554 #endif
2555         if(netgen::multithread.terminate)
2556           return false;
2557
2558         if ( comment.empty() ) // do not overwrite a previos error
2559           comment << text(err);
2560       }
2561       catch (Standard_Failure& ex)
2562       {
2563         if ( comment.empty() ) // do not overwrite a previos error
2564           comment << text(ex);
2565         err = 1;
2566       }
2567       catch (netgen::NgException exc)
2568       {
2569         if ( comment.empty() ) // do not overwrite a previos error
2570           comment << text(exc);
2571         err = 1;
2572       }
2573       _ticTime = ( doneTime += voluMeshingTime ) / _totalTime / _progressTic;
2574
2575       // Let netgen optimize 3D mesh
2576       if ( !err && _optimize )
2577       {
2578         startWith = endWith = netgen::MESHCONST_OPTVOLUME;
2579         try
2580         {
2581           OCC_CATCH_SIGNALS;
2582 #ifdef NETGEN_V5
2583           err = netgen::OCCGenerateMesh(occgeo, _ngMesh, mparams, startWith, endWith);
2584 #else
2585           err = netgen::OCCGenerateMesh(occgeo, _ngMesh, startWith, endWith, optstr);
2586 #endif
2587           if(netgen::multithread.terminate)
2588             return false;
2589
2590           if ( comment.empty() ) // do not overwrite a previos error
2591             comment << text(err);
2592         }
2593         catch (Standard_Failure& ex)
2594         {
2595           if ( comment.empty() ) // do not overwrite a previos error
2596             comment << text(ex);
2597         }
2598         catch (netgen::NgException exc)
2599         {
2600           if ( comment.empty() ) // do not overwrite a previos error
2601             comment << text(exc);
2602         }
2603       }
2604     }
2605     if (!err && mparams.secondorder > 0)
2606     {
2607       try
2608       {
2609         OCC_CATCH_SIGNALS;
2610         netgen::OCCRefinementSurfaces ref (occgeo);
2611         ref.MakeSecondOrder (*_ngMesh);
2612       }
2613       catch (Standard_Failure& ex)
2614       {
2615         if ( comment.empty() ) // do not overwrite a previos error
2616           comment << "Exception in netgen at passing to 2nd order ";
2617       }
2618       catch (netgen::NgException exc)
2619       {
2620         if ( comment.empty() ) // do not overwrite a previos error
2621           comment << exc.What();
2622       }
2623     }
2624   }
2625
2626   _ticTime = 0.98 / _progressTic;
2627
2628   int nbNod = _ngMesh->GetNP();
2629   int nbSeg = _ngMesh->GetNSeg();
2630   int nbFac = _ngMesh->GetNSE();
2631   int nbVol = _ngMesh->GetNE();
2632   bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
2633
2634   MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
2635           ", nb nodes: "    << nbNod <<
2636           ", nb segments: " << nbSeg <<
2637           ", nb faces: "    << nbFac <<
2638           ", nb volumes: "  << nbVol);
2639
2640   // Feed back the SMESHDS with the generated Nodes and Elements
2641   if ( true /*isOK*/ ) // get whatever built
2642     FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment ); //!<
2643
2644   SMESH_ComputeErrorPtr readErr = ReadErrors(nodeVec);
2645   if ( readErr && !readErr->myBadElements.empty() )
2646     error = readErr;
2647
2648   if ( error->IsOK() && ( !isOK || comment.size() > 0 ))
2649     error->myName = COMPERR_ALGO_FAILED;
2650   if ( !comment.empty() )
2651     error->myComment = comment;
2652
2653   // SetIsAlwaysComputed( true ) to empty sub-meshes, which
2654   // appear if the geometry contains coincident sub-shape due
2655   // to bool merge_solids = 1; in netgen/libsrc/occ/occgenmesh.cpp
2656   const int nbMaps = 2;
2657   const TopTools_IndexedMapOfShape* geoMaps[nbMaps] =
2658     { & occgeo.vmap, & occgeo.emap/*, & occgeo.fmap*/ };
2659   for ( int iMap = 0; iMap < nbMaps; ++iMap )
2660     for (int i = 1; i <= geoMaps[iMap]->Extent(); i++)
2661       if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( geoMaps[iMap]->FindKey(i)))
2662         if ( !sm->IsMeshComputed() )
2663           sm->SetIsAlwaysComputed( true );
2664
2665   // set bad compute error to subshapes of all failed sub-shapes
2666   if ( !error->IsOK() )
2667   {
2668     bool pb2D = false, pb3D = false;
2669     for (int i = 1; i <= occgeo.fmap.Extent(); i++) {
2670       int status = occgeo.facemeshstatus[i-1];
2671       if (status == 1 ) continue;
2672       if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.fmap( i ))) {
2673         SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2674         if ( !smError || smError->IsOK() ) {
2675           if ( status == -1 )
2676             smError.reset( new SMESH_ComputeError( *error ));
2677           else
2678             smError.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, "Ignored" ));
2679           if ( SMESH_Algo::GetMeshError( sm ) == SMESH_Algo::MEr_OK )
2680             smError->myName = COMPERR_WARNING;
2681         }
2682         pb2D = pb2D || smError->IsKO();
2683       }
2684     }
2685     if ( !pb2D ) // all faces are OK
2686       for (int i = 1; i <= occgeo.somap.Extent(); i++)
2687         if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.somap( i )))
2688         {
2689           bool smComputed = nbVol && !sm->IsEmpty();
2690           if ( smComputed && internals.hasInternalVertexInSolid( sm->GetId() ))
2691           {
2692             int nbIntV = internals.getSolidsWithVertices().find( sm->GetId() )->second.size();
2693             SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
2694             smComputed = ( smDS->NbElements() > 0 || smDS->NbNodes() > nbIntV );
2695           }
2696           SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2697           if ( !smComputed && ( !smError || smError->IsOK() ))
2698           {
2699             smError.reset( new SMESH_ComputeError( *error ));
2700             if ( nbVol && SMESH_Algo::GetMeshError( sm ) == SMESH_Algo::MEr_OK )
2701               smError->myName = COMPERR_WARNING;
2702           }
2703           pb3D = pb3D || ( smError && smError->IsKO() );
2704         }
2705     if ( !pb2D && !pb3D )
2706       err = 0; // no fatal errors, only warnings
2707   }
2708
2709   ngLib._isComputeOk = !err;
2710
2711   return !err;
2712 }
2713
2714 //=============================================================================
2715 /*!
2716  * Evaluate
2717  */
2718 //=============================================================================
2719 bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
2720 {
2721   netgen::MeshingParameters& mparams = netgen::mparam;
2722
2723
2724   // -------------------------
2725   // Prepare OCC geometry
2726   // -------------------------
2727   netgen::OCCGeometry occgeo;
2728   list< SMESH_subMesh* > meshedSM[4]; // for 0-3 dimensions
2729   NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume );
2730   PrepareOCCgeometry( occgeo, _shape, *_mesh, meshedSM, &internals );
2731
2732   bool tooManyElems = false;
2733   const int hugeNb = std::numeric_limits<int>::max() / 100;
2734
2735   // ----------------
2736   // evaluate 1D
2737   // ----------------
2738   // pass 1D simple parameters to NETGEN
2739   if ( _simpleHyp )
2740   {
2741     // not to RestrictLocalH() according to curvature during MESHCONST_ANALYSE
2742     mparams.uselocalh = false;
2743     mparams.grading = 0.8; // not limitited size growth
2744
2745     if ( _simpleHyp->GetNumberOfSegments() )
2746       // nb of segments
2747       mparams.maxh = occgeo.boundingbox.Diam();
2748     else
2749       // segment length
2750       mparams.maxh = _simpleHyp->GetLocalLength();
2751   }
2752
2753   if ( mparams.maxh == 0.0 )
2754     mparams.maxh = occgeo.boundingbox.Diam();
2755   if ( _simpleHyp || ( mparams.minh == 0.0 && _fineness != NETGENPlugin_Hypothesis::UserDefined))
2756     mparams.minh = GetDefaultMinSize( _shape, mparams.maxh );
2757
2758   // let netgen create _ngMesh and calculate element size on not meshed shapes
2759   NETGENPlugin_NetgenLibWrapper ngLib;
2760   netgen::Mesh *ngMesh = NULL;
2761 #ifndef NETGEN_V5
2762   char *optstr = 0;
2763 #endif
2764   int startWith = netgen::MESHCONST_ANALYSE;
2765   int endWith   = netgen::MESHCONST_MESHEDGES;
2766 #ifdef NETGEN_V5
2767   int err = netgen::OCCGenerateMesh(occgeo, ngMesh, mparams, startWith, endWith);
2768 #else
2769   int err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
2770 #endif
2771
2772   if(netgen::multithread.terminate)
2773     return false;
2774
2775   ngLib.setMesh(( Ng_Mesh*) ngMesh );
2776   if (err) {
2777     if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( _shape ))
2778       sm->GetComputeError().reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED ));
2779     return false;
2780   }
2781   if ( _simpleHyp )
2782   {
2783     // Pass 1D simple parameters to NETGEN
2784     // --------------------------------
2785     int      nbSeg = _simpleHyp->GetNumberOfSegments();
2786     double segSize = _simpleHyp->GetLocalLength();
2787     for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE )
2788     {
2789       const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE));
2790       if ( nbSeg )
2791         segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 );
2792       setLocalSize( e, segSize, *ngMesh );
2793     }
2794   }
2795   else // if ( ! _simpleHyp )
2796   {
2797     // Local size on vertices and edges
2798     // --------------------------------
2799     for(std::map<int,double>::const_iterator it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++)
2800     {
2801       int key = (*it).first;
2802       double hi = (*it).second;
2803       const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2804       const TopoDS_Edge& e = TopoDS::Edge(shape);
2805       setLocalSize( e, hi, *ngMesh );
2806     }
2807     for(std::map<int,double>::const_iterator it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++)
2808     {
2809       int key = (*it).first;
2810       double hi = (*it).second;
2811       const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2812       const TopoDS_Vertex& v = TopoDS::Vertex(shape);
2813       gp_Pnt p = BRep_Tool::Pnt(v);
2814       NETGENPlugin_Mesher::RestrictLocalSize( *ngMesh, p.XYZ(), hi );
2815     }
2816     for(map<int,double>::const_iterator it=FaceId2LocalSize.begin();
2817         it!=FaceId2LocalSize.end(); it++)
2818     {
2819       int key = (*it).first;
2820       double val = (*it).second;
2821       const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
2822       int faceNgID = occgeo.fmap.FindIndex(shape);
2823       occgeo.SetFaceMaxH(faceNgID, val);
2824       for ( TopExp_Explorer edgeExp( shape, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
2825         setLocalSize( TopoDS::Edge( edgeExp.Current() ), val, *ngMesh );
2826     }
2827   }
2828   // calculate total nb of segments and length of edges
2829   double fullLen = 0.0;
2830   int fullNbSeg = 0;
2831   int entity = mparams.secondorder > 0 ? SMDSEntity_Quad_Edge : SMDSEntity_Edge;
2832   TopTools_DataMapOfShapeInteger Edge2NbSeg;
2833   for (TopExp_Explorer exp(_shape, TopAbs_EDGE); exp.More(); exp.Next())
2834   {
2835     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
2836     if( !Edge2NbSeg.Bind(E,0) )
2837       continue;
2838
2839     double aLen = SMESH_Algo::EdgeLength(E);
2840     fullLen += aLen;
2841
2842     vector<int>& aVec = aResMap[_mesh->GetSubMesh(E)];
2843     if ( aVec.empty() )
2844       aVec.resize( SMDSEntity_Last, 0);
2845     else
2846       fullNbSeg += aVec[ entity ];
2847   }
2848
2849   // store nb of segments computed by Netgen
2850   NCollection_Map<Link> linkMap;
2851   for (int i = 1; i <= ngMesh->GetNSeg(); ++i )
2852   {
2853     const netgen::Segment& seg = ngMesh->LineSegment(i);
2854     Link link(seg[0], seg[1]);
2855     if ( !linkMap.Add( link )) continue;
2856     int aGeomEdgeInd = seg.epgeominfo[0].edgenr;
2857     if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
2858     {
2859       vector<int>& aVec = aResMap[_mesh->GetSubMesh(occgeo.emap(aGeomEdgeInd))];
2860       aVec[ entity ]++;
2861     }
2862   }
2863   // store nb of nodes on edges computed by Netgen
2864   TopTools_DataMapIteratorOfDataMapOfShapeInteger Edge2NbSegIt(Edge2NbSeg);
2865   for (; Edge2NbSegIt.More(); Edge2NbSegIt.Next())
2866   {
2867     vector<int>& aVec = aResMap[_mesh->GetSubMesh(Edge2NbSegIt.Key())];
2868     if ( aVec[ entity ] > 1 && aVec[ SMDSEntity_Node ] == 0 )
2869       aVec[SMDSEntity_Node] = mparams.secondorder > 0  ? 2*aVec[ entity ]-1 : aVec[ entity ]-1;
2870
2871     fullNbSeg += aVec[ entity ];
2872     Edge2NbSeg( Edge2NbSegIt.Key() ) = aVec[ entity ];
2873   }
2874   if ( fullNbSeg == 0 )
2875     return false;
2876
2877   // ----------------
2878   // evaluate 2D
2879   // ----------------
2880   if ( _simpleHyp ) {
2881     if ( double area = _simpleHyp->GetMaxElementArea() ) {
2882       // face area
2883       mparams.maxh = sqrt(2. * area/sqrt(3.0));
2884       mparams.grading = 0.4; // moderate size growth
2885     }
2886     else {
2887       // length from edges
2888       mparams.maxh = fullLen/fullNbSeg;
2889       mparams.grading = 0.2; // slow size growth
2890     }
2891   }
2892   mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
2893   mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
2894
2895   for (TopExp_Explorer exp(_shape, TopAbs_FACE); exp.More(); exp.Next())
2896   {
2897     TopoDS_Face F = TopoDS::Face( exp.Current() );
2898     SMESH_subMesh *sm = _mesh->GetSubMesh(F);
2899     GProp_GProps G;
2900     BRepGProp::SurfaceProperties(F,G);
2901     double anArea = G.Mass();
2902     tooManyElems = tooManyElems || ( anArea/hugeNb > mparams.maxh*mparams.maxh );
2903     int nb1d = 0;
2904     if ( !tooManyElems )
2905     {
2906       TopTools_MapOfShape egdes;
2907       for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next())
2908         if ( egdes.Add( exp1.Current() ))
2909           nb1d += Edge2NbSeg.Find(exp1.Current());
2910     }
2911     int nbFaces = tooManyElems ? hugeNb : int( 4*anArea / (mparams.maxh*mparams.maxh*sqrt(3.)));
2912     int nbNodes = tooManyElems ? hugeNb : (( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
2913
2914     vector<int> aVec(SMDSEntity_Last, 0);
2915     if( mparams.secondorder > 0 ) {
2916       int nb1d_in = (nbFaces*3 - nb1d) / 2;
2917       aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
2918       aVec[SMDSEntity_Quad_Triangle] = nbFaces;
2919     }
2920     else {
2921       aVec[SMDSEntity_Node] = Max ( nbNodes, 0  );
2922       aVec[SMDSEntity_Triangle] = nbFaces;
2923     }
2924     aResMap[sm].swap(aVec);
2925   }
2926
2927   // ----------------
2928   // evaluate 3D
2929   // ----------------
2930   if(_isVolume) {
2931     // pass 3D simple parameters to NETGEN
2932     const NETGENPlugin_SimpleHypothesis_3D* simple3d =
2933       dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
2934     if ( simple3d ) {
2935       if ( double vol = simple3d->GetMaxElementVolume() ) {
2936         // max volume
2937         mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
2938         mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
2939       }
2940       else {
2941         // using previous length from faces
2942       }
2943       mparams.grading = 0.4;
2944       mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
2945     }
2946     GProp_GProps G;
2947     BRepGProp::VolumeProperties(_shape,G);
2948     double aVolume = G.Mass();
2949     double tetrVol = 0.1179*mparams.maxh*mparams.maxh*mparams.maxh;
2950     tooManyElems = tooManyElems || ( aVolume/hugeNb > tetrVol );
2951     int nbVols = tooManyElems ? hugeNb : int(aVolume/tetrVol);
2952     int nb1d_in = int(( nbVols*6 - fullNbSeg ) / 6 );
2953     vector<int> aVec(SMDSEntity_Last, 0 );
2954     if ( tooManyElems ) // avoid FPE
2955     {
2956       aVec[SMDSEntity_Node] = hugeNb;
2957       aVec[ mparams.secondorder > 0 ? SMDSEntity_Quad_Tetra : SMDSEntity_Tetra] = hugeNb;
2958     }
2959     else
2960     {
2961       if( mparams.secondorder > 0 ) {
2962         aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in;
2963         aVec[SMDSEntity_Quad_Tetra] = nbVols;
2964       }
2965       else {
2966         aVec[SMDSEntity_Node] = nb1d_in/3 + 1;
2967         aVec[SMDSEntity_Tetra] = nbVols;
2968       }
2969     }
2970     SMESH_subMesh *sm = _mesh->GetSubMesh(_shape);
2971     aResMap[sm].swap(aVec);
2972   }
2973
2974   return true;
2975 }
2976
2977 double NETGENPlugin_Mesher::GetProgress(const SMESH_Algo* holder,
2978                                         const int *       algoProgressTic,
2979                                         const double *    algoProgress) const
2980 {
2981   ((int&) _progressTic ) = *algoProgressTic + 1;
2982
2983   if ( !_occgeom ) return 0;
2984
2985   double progress = -1;
2986   if ( !_isVolume )
2987   {
2988     if ( _ticTime < 0 && netgen::multithread.task[0] == 'O'/*Optimizing surface*/ )
2989     {
2990       ((double&) _ticTime ) = edgeFaceMeshingTime / _totalTime / _progressTic;
2991     }
2992     else if ( !_optimize /*&& _occgeom->fmap.Extent() > 1*/ )
2993     {
2994       int doneShapeIndex = -1;
2995       while ( doneShapeIndex+1 < _occgeom->facemeshstatus.Size() &&
2996               _occgeom->facemeshstatus[ doneShapeIndex+1 ])
2997         doneShapeIndex++;
2998       if ( doneShapeIndex+1 != _curShapeIndex )
2999       {
3000         ((int&) _curShapeIndex) = doneShapeIndex+1;
3001         double    doneShapeRate = _curShapeIndex / double( _occgeom->fmap.Extent() );
3002         double         doneTime = edgeMeshingTime + doneShapeRate * faceMeshingTime;
3003         ((double&)    _ticTime) = doneTime / _totalTime / _progressTic;
3004         // cout << "shape " << _curShapeIndex << " _ticTime " << _ticTime
3005         //      << " " << doneTime / _totalTime / _progressTic << endl;
3006       }
3007     }
3008   }
3009   else if ( !_optimize && _occgeom->somap.Extent() > 1 )
3010   {
3011     int curShapeIndex = _curShapeIndex;
3012     if ( _ngMesh->GetNE() > 0 )
3013     {
3014       netgen::Element el = (*_ngMesh)[netgen::ElementIndex( _ngMesh->GetNE()-1 )];
3015       curShapeIndex = el.GetIndex();
3016     }
3017     if ( curShapeIndex != _curShapeIndex )
3018     {
3019       ((int&) _curShapeIndex) = curShapeIndex;
3020       double    doneShapeRate = _curShapeIndex / double( _occgeom->somap.Extent() );
3021       double         doneTime = edgeFaceMeshingTime + doneShapeRate * voluMeshingTime;
3022       ((double&)    _ticTime) = doneTime / _totalTime / _progressTic;
3023       // cout << "shape " << _curShapeIndex << " _ticTime " << _ticTime
3024       //      << " " << doneTime / _totalTime / _progressTic << endl;
3025     }
3026   }
3027   if ( _ticTime > 0 )
3028     progress  = Max( *algoProgressTic * _ticTime, *algoProgress );
3029   if ( progress > 0 )
3030   {
3031     ((int&) *algoProgressTic )++;
3032     ((double&) *algoProgress) = progress;
3033   }
3034   //cout << progress << " "  << *algoProgressTic << " " << netgen::multithread.task << " "<< _ticTime << endl;
3035
3036   return Min( progress, 0.99 );
3037 }
3038
3039 //================================================================================
3040 /*!
3041  * \brief Remove "test.out" and "problemfaces" files in current directory
3042  */
3043 //================================================================================
3044
3045 void NETGENPlugin_Mesher::RemoveTmpFiles()
3046 {
3047   bool rm =  SMESH_File("test.out").remove() ;
3048 #ifndef WIN32
3049   if (rm && netgen::testout)
3050   {
3051     delete netgen::testout;
3052     netgen::testout = 0;
3053   }
3054 #endif
3055   SMESH_File("problemfaces").remove();
3056   SMESH_File("occmesh.rep").remove();
3057 }
3058
3059 //================================================================================
3060 /*!
3061  * \brief Read mesh entities preventing successful computation from "test.out" file
3062  */
3063 //================================================================================
3064
3065 SMESH_ComputeErrorPtr
3066 NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
3067 {
3068   SMESH_ComputeErrorPtr err = SMESH_ComputeError::New
3069     (COMPERR_BAD_INPUT_MESH, "Some edges multiple times in surface mesh");
3070   SMESH_File file("test.out");
3071   vector<int> two(2);
3072   const char* badEdgeStr = " multiple times in surface mesh";
3073   const int   badEdgeStrLen = strlen( badEdgeStr );
3074   while( !file.eof() )
3075   {
3076     if ( strncmp( file, "Edge ", 5 ) == 0 &&
3077          file.getInts( two ) &&
3078          strncmp( file, badEdgeStr, badEdgeStrLen ) == 0 &&
3079          two[0] < nodeVec.size() && two[1] < nodeVec.size())
3080     {
3081       err->myBadElements.push_back( new SMDS_LinearEdge( nodeVec[ two[0]], nodeVec[ two[1]] ));
3082       file += badEdgeStrLen;
3083     }
3084     else if ( strncmp( file, "Intersecting: ", 14 ) == 0 )
3085     {
3086 // Intersecting:
3087 // openelement 18 with open element 126
3088 // 41  36  38
3089 // 69  70  72
3090       vector<int> three1(3), three2(3);
3091       file.getLine();
3092       const char* pos = file;
3093       bool ok = ( strncmp( file, "openelement ", 12 ) == 0 );
3094       ok = ok && file.getInts( two );
3095       ok = ok && file.getInts( three1 );
3096       ok = ok && file.getInts( three2 );
3097       for ( int i = 0; ok && i < 3; ++i )
3098         ok = ( three1[i] < nodeVec.size() && nodeVec[ three1[i]]);
3099       for ( int i = 0; ok && i < 3; ++i )
3100         ok = ( three2[i] < nodeVec.size() && nodeVec[ three2[i]]);
3101       if ( ok )
3102       {
3103         err->myBadElements.push_back( new SMDS_FaceOfNodes( nodeVec[ three1[0]],
3104                                                             nodeVec[ three1[1]],
3105                                                             nodeVec[ three1[2]]));
3106         err->myBadElements.push_back( new SMDS_FaceOfNodes( nodeVec[ three2[0]],
3107                                                             nodeVec[ three2[1]],
3108                                                             nodeVec[ three2[2]]));
3109         err->myComment = "Intersecting triangles";
3110       }
3111       else
3112       {
3113         file.setPos( pos );
3114       }
3115     }
3116     else
3117     {
3118       ++file;
3119     }
3120   }
3121   return err;
3122 }
3123
3124 //================================================================================
3125 /*!
3126  * \brief Write a python script creating an equivalent SALOME mesh.
3127  * This is useful to see what mesh is passed as input for the next step of mesh
3128  * generation (of mesh of higher dimension)
3129  */
3130 //================================================================================
3131
3132 void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh,
3133                                     const std::string&  pyFile)
3134 {
3135   ofstream outfile(pyFile.c_str(), ios::out);
3136   if ( !outfile ) return;
3137
3138   outfile << "import SMESH" << endl
3139           << "from salome.smesh import smeshBuilder" << endl
3140           << "smesh = smeshBuilder.New(salome.myStudy)" << endl
3141           << "mesh = smesh.Mesh()" << endl << endl;
3142
3143   using namespace netgen;
3144   PointIndex pi;
3145   for (pi = PointIndex::BASE;
3146        pi < ngMesh->GetNP()+PointIndex::BASE; pi++)
3147   {
3148     outfile << "mesh.AddNode( ";
3149     outfile << (*ngMesh)[pi](0) << ", ";
3150     outfile << (*ngMesh)[pi](1) << ", ";
3151     outfile << (*ngMesh)[pi](2) << ") ## "<< pi << endl;
3152   }
3153
3154   int nbDom = ngMesh->GetNDomains();
3155   for ( int i = 0; i < nbDom; ++i )
3156     outfile<< "grp" << i+1 << " = mesh.CreateEmptyGroup( SMESH.FACE, 'domain"<< i+1 << "')"<< endl;
3157
3158   SurfaceElementIndex sei;
3159   for (sei = 0; sei < ngMesh->GetNSE(); sei++)
3160   {
3161     outfile << "mesh.AddFace([ ";
3162     Element2d sel = (*ngMesh)[sei];
3163     for (int j = 0; j < sel.GetNP(); j++)
3164       outfile << sel[j] << ( j+1 < sel.GetNP() ? ", " : " ])");
3165     if ( sel.IsDeleted() ) outfile << " ## IsDeleted ";
3166     outfile << endl;
3167
3168     if ((*ngMesh)[sei].GetIndex())
3169     {
3170       if ( int dom1 = ngMesh->GetFaceDescriptor((*ngMesh)[sei].GetIndex ()).DomainIn())
3171         outfile << "grp"<< dom1 <<".Add([ " << (int)sei+1 << " ])" << endl;
3172       if ( int dom2 = ngMesh->GetFaceDescriptor((*ngMesh)[sei].GetIndex ()).DomainOut())
3173         outfile << "grp"<< dom2 <<".Add([ " << (int)sei+1 << " ])" << endl;
3174     }
3175   }
3176
3177   for (ElementIndex ei = 0; ei < ngMesh->GetNE(); ei++)
3178   {
3179     Element el = (*ngMesh)[ei];
3180     outfile << "mesh.AddVolume([ ";
3181     for (int j = 0; j < el.GetNP(); j++)
3182       outfile << el[j] << ( j+1 < el.GetNP() ? ", " : " ])");
3183     outfile << endl;
3184   }
3185
3186   for (int i = 1; i <= ngMesh->GetNSeg(); i++)
3187   {
3188     const Segment & seg = ngMesh->LineSegment (i);
3189     outfile << "mesh.AddEdge([ "
3190             << seg[0] << ", "
3191             << seg[1] << " ])" << endl;
3192   }
3193   cout << "Write " << pyFile << endl;
3194 }
3195
3196 //================================================================================
3197 /*!
3198  * \brief Constructor of NETGENPlugin_ngMeshInfo
3199  */
3200 //================================================================================
3201
3202 NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh):
3203   _copyOfLocalH(0)
3204 {
3205   if ( ngMesh )
3206   {
3207     _nbNodes    = ngMesh->GetNP();
3208     _nbSegments = ngMesh->GetNSeg();
3209     _nbFaces    = ngMesh->GetNSE();
3210     _nbVolumes  = ngMesh->GetNE();
3211   }
3212   else
3213   {
3214     _nbNodes = _nbSegments = _nbFaces = _nbVolumes = 0;
3215   }
3216 }
3217
3218 //================================================================================
3219 /*!
3220  * \brief Copy LocalH member from one netgen mesh to another
3221  */
3222 //================================================================================
3223
3224 void NETGENPlugin_ngMeshInfo::transferLocalH( netgen::Mesh* fromMesh,
3225                                               netgen::Mesh* toMesh )
3226 {
3227   if ( !fromMesh->LocalHFunctionGenerated() ) return;
3228   if ( !toMesh->LocalHFunctionGenerated() )
3229 #ifdef NETGEN_V5
3230     toMesh->CalcLocalH(netgen::mparam.grading);
3231 #else
3232     toMesh->CalcLocalH();
3233 #endif
3234
3235   const size_t size = sizeof( netgen::LocalH );
3236   _copyOfLocalH = new char[ size ];
3237   memcpy( (void*)_copyOfLocalH, (void*)&toMesh->LocalHFunction(), size );
3238   memcpy( (void*)&toMesh->LocalHFunction(), (void*)&fromMesh->LocalHFunction(), size );
3239 }
3240
3241 //================================================================================
3242 /*!
3243  * \brief Restore LocalH member of a netgen mesh
3244  */
3245 //================================================================================
3246
3247 void NETGENPlugin_ngMeshInfo::restoreLocalH( netgen::Mesh* toMesh )
3248 {
3249   if ( _copyOfLocalH )
3250   {
3251     const size_t size = sizeof( netgen::LocalH );
3252     memcpy( (void*)&toMesh->LocalHFunction(), (void*)_copyOfLocalH, size );
3253     delete [] _copyOfLocalH;
3254     _copyOfLocalH = 0;
3255   }
3256 }
3257
3258 //================================================================================
3259 /*!
3260  * \brief Find "internal" sub-shapes
3261  */
3262 //================================================================================
3263
3264 NETGENPlugin_Internals::NETGENPlugin_Internals( SMESH_Mesh&         mesh,
3265                                                 const TopoDS_Shape& shape,
3266                                                 bool                is3D )
3267   : _mesh( mesh ), _is3D( is3D )
3268 {
3269   SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
3270
3271   TopExp_Explorer f,e;
3272   for ( f.Init( shape, TopAbs_FACE ); f.More(); f.Next() )
3273   {
3274     int faceID = meshDS->ShapeToIndex( f.Current() );
3275
3276     // find not computed internal edges
3277
3278     for ( e.Init( f.Current().Oriented(TopAbs_FORWARD), TopAbs_EDGE ); e.More(); e.Next() )
3279       if ( e.Current().Orientation() == TopAbs_INTERNAL )
3280       {
3281         SMESH_subMesh* eSM = mesh.GetSubMesh( e.Current() );
3282         if ( eSM->IsEmpty() )
3283         {
3284           _e2face.insert( make_pair( eSM->GetId(), faceID ));
3285           for ( TopoDS_Iterator v(e.Current()); v.More(); v.Next() )
3286             _e2face.insert( make_pair( meshDS->ShapeToIndex( v.Value() ), faceID ));
3287         }
3288       }
3289
3290     // find internal vertices in a face
3291     set<int> intVV; // issue 0020850 where same vertex is twice in a face
3292     for ( TopoDS_Iterator fSub( f.Current() ); fSub.More(); fSub.Next())
3293       if ( fSub.Value().ShapeType() == TopAbs_VERTEX )
3294       {
3295         int vID = meshDS->ShapeToIndex( fSub.Value() );
3296         if ( intVV.insert( vID ).second )
3297           _f2v[ faceID ].push_back( vID );
3298       }
3299
3300     if ( is3D )
3301     {
3302       // find internal faces and their subshapes where nodes are to be doubled
3303       //  to make a crack with non-sewed borders
3304
3305       if ( f.Current().Orientation() == TopAbs_INTERNAL )
3306       {
3307         _intShapes.insert( meshDS->ShapeToIndex( f.Current() ));
3308
3309         // egdes
3310         list< TopoDS_Shape > edges;
3311         for ( e.Init( f.Current(), TopAbs_EDGE ); e.More(); e.Next())
3312           if ( SMESH_MesherHelper::NbAncestors( e.Current(), mesh, TopAbs_FACE ) > 1 )
3313           {
3314             _intShapes.insert( meshDS->ShapeToIndex( e.Current() ));
3315             edges.push_back( e.Current() );
3316             // find border faces
3317             PShapeIteratorPtr fIt =
3318               SMESH_MesherHelper::GetAncestors( edges.back(),mesh,TopAbs_FACE );
3319             while ( const TopoDS_Shape* pFace = fIt->next() )
3320               if ( !pFace->IsSame( f.Current() ))
3321                 _borderFaces.insert( meshDS->ShapeToIndex( *pFace ));
3322           }
3323         // vertices
3324         // we consider vertex internal if it is shared by more than one internal edge
3325         list< TopoDS_Shape >::iterator edge = edges.begin();
3326         for ( ; edge != edges.end(); ++edge )
3327           for ( TopoDS_Iterator v( *edge ); v.More(); v.Next() )
3328           {
3329             set<int> internalEdges;
3330             PShapeIteratorPtr eIt =
3331               SMESH_MesherHelper::GetAncestors( v.Value(),mesh,TopAbs_EDGE );
3332             while ( const TopoDS_Shape* pEdge = eIt->next() )
3333             {
3334               int edgeID = meshDS->ShapeToIndex( *pEdge );
3335               if ( isInternalShape( edgeID ))
3336                 internalEdges.insert( edgeID );
3337             }
3338             if ( internalEdges.size() > 1 )
3339               _intShapes.insert( meshDS->ShapeToIndex( v.Value() ));
3340           }
3341       }
3342     }
3343   } // loop on geom faces
3344
3345   // find vertices internal in solids
3346   if ( is3D )
3347   {
3348     for ( TopExp_Explorer so(shape, TopAbs_SOLID); so.More(); so.Next())
3349     {
3350       int soID = meshDS->ShapeToIndex( so.Current() );
3351       for ( TopoDS_Iterator soSub( so.Current() ); soSub.More(); soSub.Next())
3352         if ( soSub.Value().ShapeType() == TopAbs_VERTEX )
3353           _s2v[ soID ].push_back( meshDS->ShapeToIndex( soSub.Value() ));
3354     }
3355   }
3356 }
3357
3358 //================================================================================
3359 /*!
3360  * \brief Find mesh faces on non-internal geom faces sharing internal edge
3361  * some nodes of which are to be doubled to make the second border of the "crack"
3362  */
3363 //================================================================================
3364
3365 void NETGENPlugin_Internals::findBorderElements( TIDSortedElemSet & borderElems )
3366 {
3367   if ( _intShapes.empty() ) return;
3368
3369   SMESH_Mesh& mesh = const_cast<SMESH_Mesh&>(_mesh);
3370   SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
3371
3372   // loop on internal geom edges
3373   set<int>::const_iterator intShapeId = _intShapes.begin();
3374   for ( ; intShapeId != _intShapes.end(); ++intShapeId )
3375   {
3376     const TopoDS_Shape& s = meshDS->IndexToShape( *intShapeId );
3377     if ( s.ShapeType() != TopAbs_EDGE ) continue;
3378
3379     // get internal and non-internal geom faces sharing the internal edge <s>
3380     int intFace = 0;
3381     set<int>::iterator bordFace = _borderFaces.end();
3382     PShapeIteratorPtr faces = SMESH_MesherHelper::GetAncestors( s, _mesh, TopAbs_FACE );
3383     while ( const TopoDS_Shape* pFace = faces->next() )
3384     {
3385       int faceID = meshDS->ShapeToIndex( *pFace );
3386       if ( isInternalShape( faceID ))
3387         intFace = faceID;
3388       else
3389         bordFace = _borderFaces.insert( faceID ).first;
3390     }
3391     if ( bordFace == _borderFaces.end() || !intFace ) continue;
3392
3393     // get all links of mesh faces on internal geom face sharing nodes on edge <s>
3394     set< SMESH_OrientedLink > links; //!< links of faces on internal geom face
3395     list<const SMDS_MeshElement*> suspectFaces[2]; //!< mesh faces on border geom faces
3396     int nbSuspectFaces = 0;
3397     SMESHDS_SubMesh* intFaceSM = meshDS->MeshElements( intFace );
3398     if ( !intFaceSM || intFaceSM->NbElements() == 0 ) continue;
3399     SMESH_subMeshIteratorPtr smIt = mesh.GetSubMesh( s )->getDependsOnIterator(true,true);
3400     while ( smIt->more() )
3401     {
3402       SMESHDS_SubMesh* sm = smIt->next()->GetSubMeshDS();
3403       if ( !sm ) continue;
3404       SMDS_NodeIteratorPtr nIt = sm->GetNodes();
3405       while ( nIt->more() )
3406       {
3407         const SMDS_MeshNode* nOnEdge = nIt->next();
3408         SMDS_ElemIteratorPtr fIt = nOnEdge->GetInverseElementIterator(SMDSAbs_Face);
3409         while ( fIt->more() )
3410         {
3411           const SMDS_MeshElement* f = fIt->next();
3412           int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
3413           if ( intFaceSM->Contains( f ))
3414           {
3415             for ( int i = 0; i < nbNodes; ++i )
3416               links.insert( SMESH_OrientedLink( f->GetNode(i), f->GetNode((i+1)%nbNodes)));
3417           }
3418           else
3419           {
3420             int nbDblNodes = 0;
3421             for ( int i = 0; i < nbNodes; ++i )
3422               nbDblNodes += isInternalShape( f->GetNode(i)->getshapeId() );
3423             if ( nbDblNodes )
3424               suspectFaces[ nbDblNodes < 2 ].push_back( f );
3425             nbSuspectFaces++;
3426           }
3427         }
3428       }
3429     }
3430     // suspectFaces[0] having link with same orientation as mesh faces on
3431     // the internal geom face are <borderElems>. suspectFaces[1] have
3432     // only one node on edge <s>, we decide on them later (at the 2nd loop)
3433     // by links of <borderElems> found at the 1st and 2nd loops
3434     set< SMESH_OrientedLink > borderLinks;
3435     for ( int isPostponed = 0; isPostponed < 2; ++isPostponed )
3436     {
3437       list<const SMDS_MeshElement*>::iterator fIt = suspectFaces[isPostponed].begin();
3438       for ( int nbF = 0; fIt != suspectFaces[isPostponed].end(); ++fIt, ++nbF )
3439       {
3440         const SMDS_MeshElement* f = *fIt;
3441         bool isBorder = false, linkFound = false, borderLinkFound = false;
3442         list< SMESH_OrientedLink > faceLinks;
3443         int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
3444         for ( int i = 0; i < nbNodes; ++i )
3445         {
3446           SMESH_OrientedLink link( f->GetNode(i), f->GetNode((i+1)%nbNodes));
3447           faceLinks.push_back( link );
3448           if ( !linkFound )
3449           {
3450             set< SMESH_OrientedLink >::iterator foundLink = links.find( link );
3451             if ( foundLink != links.end() )
3452             {
3453               linkFound= true;
3454               isBorder = ( foundLink->_reversed == link._reversed );
3455               if ( !isBorder && !isPostponed ) break;
3456               faceLinks.pop_back();
3457             }
3458             else if ( isPostponed && !borderLinkFound )
3459             {
3460               foundLink = borderLinks.find( link );
3461               if ( foundLink != borderLinks.end() )
3462               {
3463                 borderLinkFound = true;
3464                 isBorder = ( foundLink->_reversed != link._reversed );
3465               }
3466             }
3467           }
3468         }
3469         if ( isBorder )
3470         {
3471           borderElems.insert( f );
3472           borderLinks.insert( faceLinks.begin(), faceLinks.end() );
3473         }
3474         else if ( !linkFound && !borderLinkFound )
3475         {
3476           suspectFaces[1].push_back( f );
3477           if ( nbF > 2 * nbSuspectFaces )
3478             break; // dead loop protection
3479         }
3480       }
3481     }
3482   }
3483 }
3484
3485 //================================================================================
3486 /*!
3487  * \brief put internal shapes in maps and fill in submeshes to precompute
3488  */
3489 //================================================================================
3490
3491 void NETGENPlugin_Internals::getInternalEdges( TopTools_IndexedMapOfShape& fmap,
3492                                                TopTools_IndexedMapOfShape& emap,
3493                                                TopTools_IndexedMapOfShape& vmap,
3494                                                list< SMESH_subMesh* > smToPrecompute[])
3495 {
3496   if ( !hasInternalEdges() ) return;
3497   map<int,int>::const_iterator ev_face = _e2face.begin();
3498   for ( ; ev_face != _e2face.end(); ++ev_face )
3499   {
3500     const TopoDS_Shape& ev   = _mesh.GetMeshDS()->IndexToShape( ev_face->first );
3501     const TopoDS_Shape& face = _mesh.GetMeshDS()->IndexToShape( ev_face->second );
3502
3503     ( ev.ShapeType() == TopAbs_EDGE ? emap : vmap ).Add( ev );
3504     fmap.Add( face );
3505     //cout<<"INTERNAL EDGE or VERTEX "<<ev_face->first<<" on face "<<ev_face->second<<endl;
3506
3507     smToPrecompute[ MeshDim_1D ].push_back( _mesh.GetSubMeshContaining( ev_face->first ));
3508   }
3509 }
3510
3511 //================================================================================
3512 /*!
3513  * \brief return shapes and submeshes to be meshed and already meshed boundary submeshes
3514  */
3515 //================================================================================
3516
3517 void NETGENPlugin_Internals::getInternalFaces( TopTools_IndexedMapOfShape& fmap,
3518                                                TopTools_IndexedMapOfShape& emap,
3519                                                list< SMESH_subMesh* >&     intFaceSM,
3520                                                list< SMESH_subMesh* >&     boundarySM)
3521 {
3522   if ( !hasInternalFaces() ) return;
3523
3524   // <fmap> and <emap> are for not yet meshed shapes
3525   // <intFaceSM> is for submeshes of faces
3526   // <boundarySM> is for meshed edges and vertices
3527
3528   intFaceSM.clear();
3529   boundarySM.clear();
3530
3531   set<int> shapeIDs ( _intShapes );
3532   if ( !_borderFaces.empty() )
3533     shapeIDs.insert( _borderFaces.begin(), _borderFaces.end() );
3534
3535   set<int>::const_iterator intS = shapeIDs.begin();
3536   for ( ; intS != shapeIDs.end(); ++intS )
3537   {
3538     SMESH_subMesh* sm = _mesh.GetSubMeshContaining( *intS );
3539
3540     if ( sm->GetSubShape().ShapeType() != TopAbs_FACE ) continue;
3541
3542     intFaceSM.push_back( sm );
3543
3544     // add submeshes of not computed internal faces
3545     if ( !sm->IsEmpty() ) continue;
3546
3547     SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(true,true);
3548     while ( smIt->more() )
3549     {
3550       sm = smIt->next();
3551       const TopoDS_Shape& s = sm->GetSubShape();
3552
3553       if ( sm->IsEmpty() )
3554       {
3555         // not yet meshed
3556         switch ( s.ShapeType() ) {
3557         case TopAbs_FACE: fmap.Add ( s ); break;
3558         case TopAbs_EDGE: emap.Add ( s ); break;
3559         default:;
3560         }
3561       }
3562       else
3563       {
3564         if ( s.ShapeType() != TopAbs_FACE )
3565           boundarySM.push_back( sm );
3566       }
3567     }
3568   }
3569 }
3570
3571 //================================================================================
3572 /*!
3573  * \brief Return true if given shape is to be precomputed in order to be correctly
3574  * added to netgen mesh
3575  */
3576 //================================================================================
3577
3578 bool NETGENPlugin_Internals::isShapeToPrecompute(const TopoDS_Shape& s)
3579 {
3580   int shapeID = _mesh.GetMeshDS()->ShapeToIndex( s );
3581   switch ( s.ShapeType() ) {
3582   case TopAbs_FACE  : break; //return isInternalShape( shapeID ) || isBorderFace( shapeID );
3583   case TopAbs_EDGE  : return isInternalEdge( shapeID );
3584   case TopAbs_VERTEX: break;
3585   default:;
3586   }
3587   return false;
3588 }
3589
3590 //================================================================================
3591 /*!
3592  * \brief Return SMESH
3593  */
3594 //================================================================================
3595
3596 SMESH_Mesh& NETGENPlugin_Internals::getMesh() const
3597 {
3598   return const_cast<SMESH_Mesh&>( _mesh );
3599 }
3600
3601 //================================================================================
3602 /*!
3603  * \brief Initialize netgen library
3604  */
3605 //================================================================================
3606
3607 NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper()
3608 {
3609   Ng_Init();
3610
3611   // redirect all netgen output (mycout,myerr,cout) to _outputFileName
3612   _isComputeOk    = false;
3613   _outputFileName = getOutputFileName();
3614   netgen::mycout  = new ofstream ( _outputFileName.c_str() );
3615   netgen::myerr   = netgen::mycout;
3616   _coutBuffer     = std::cout.rdbuf();
3617 #ifdef _DEBUG_
3618   cout << "NOTE: netgen output is redirected to file " << _outputFileName << endl;
3619 #else
3620   std::cout.rdbuf( netgen::mycout->rdbuf() );
3621 #endif
3622
3623   _ngMesh = Ng_NewMesh();
3624 }
3625
3626 //================================================================================
3627 /*!
3628  * \brief Finish using netgen library
3629  */
3630 //================================================================================
3631
3632 NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper()
3633 {
3634   Ng_DeleteMesh( _ngMesh );
3635   Ng_Exit();
3636   NETGENPlugin_Mesher::RemoveTmpFiles();
3637   std::cout.rdbuf( _coutBuffer );
3638 #ifdef _DEBUG_
3639   if( _isComputeOk )
3640 #endif
3641     removeOutputFile();
3642 }
3643
3644 //================================================================================
3645 /*!
3646  * \brief Set netgen mesh to delete at destruction
3647  */
3648 //================================================================================
3649
3650 void NETGENPlugin_NetgenLibWrapper::setMesh( Ng_Mesh* mesh )
3651 {
3652   if ( _ngMesh )
3653     Ng_DeleteMesh( _ngMesh );
3654   _ngMesh = mesh;
3655 }
3656
3657 //================================================================================
3658 /*!
3659  * \brief Return a unique file name
3660  */
3661 //================================================================================
3662
3663 std::string NETGENPlugin_NetgenLibWrapper::getOutputFileName()
3664 {
3665   std::string aTmpDir = SALOMEDS_Tool::GetTmpDir();
3666
3667   TCollection_AsciiString aGenericName = (char*)aTmpDir.c_str();
3668   aGenericName += "NETGEN_";
3669 #ifndef WIN32
3670   aGenericName += getpid();
3671 #else
3672   aGenericName += _getpid();
3673 #endif
3674   aGenericName += "_";
3675   aGenericName += Abs((Standard_Integer)(long) aGenericName.ToCString());
3676   aGenericName += ".out";
3677
3678   return aGenericName.ToCString();
3679 }
3680
3681 //================================================================================
3682 /*!
3683  * \brief Remove file with netgen output
3684  */
3685 //================================================================================
3686
3687 void NETGENPlugin_NetgenLibWrapper::removeOutputFile()
3688 {
3689   string tmpDir = SALOMEDS_Tool::GetDirFromPath( _outputFileName );
3690   SALOMEDS::ListOfFileNames_var aFiles = new SALOMEDS::ListOfFileNames;
3691   aFiles->length(1);
3692   std::string aFileName = SALOMEDS_Tool::GetNameFromPath( _outputFileName ) + ".out";
3693   aFiles[0] = aFileName.c_str();
3694   if ( netgen::mycout)
3695   {
3696     delete netgen::mycout;
3697     netgen::mycout = 0;
3698     netgen::myerr = 0;
3699   }
3700
3701   SALOMEDS_Tool::RemoveTemporaryFiles( tmpDir.c_str(), aFiles.in(), true );
3702 #ifdef _DEBUG_
3703   //cout << "NOTE: netgen output log was REMOVED       " << _outputFileName << endl;
3704 #endif
3705 }