src/NETGENPlugin/NETGENPlugin_Mesher.cxx

   1 // Copyright (C) 2007-2012  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.
  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 <utilities.h>
  47
  48 #include <vector>
  49 #include <limits>
  50
  51 #include <BRep_Tool.hxx>
  52 #include <Bnd_B3d.hxx>
  53 #include <GCPnts_AbscissaPoint.hxx>
  54 #include <GeomAdaptor_Curve.hxx>
  55 #include <NCollection_Map.hxx>
  56 #include <OSD_File.hxx>
  57 #include <OSD_Path.hxx>
  58 #include <Standard_ErrorHandler.hxx>
  59 #include <Standard_ProgramError.hxx>
  60 #include <TCollection_AsciiString.hxx>
  61 #include <TopExp.hxx>
  62 #include <TopExp_Explorer.hxx>
  63 #include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
  64 #include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
  65 #include <TopTools_DataMapOfShapeInteger.hxx>
  66 #include <TopTools_DataMapOfShapeShape.hxx>
  67 #include <TopTools_ListIteratorOfListOfShape.hxx>
  68 #include <TopTools_MapOfShape.hxx>
  69 #include <TopoDS.hxx>
  70
  71 // Netgen include files
  72 #ifndef OCCGEOMETRY
  73 #define OCCGEOMETRY
  74 #endif
  75 #include <occgeom.hpp>
  76 #include <meshing.hpp>
  77 //#include <ngexception.hpp>
  78 namespace netgen {
  79   extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
  80   extern MeshingParameters mparam;
  81   extern volatile multithreadt multithread;
  82 }
  83
  84 using namespace nglib;
  85 using namespace std;
  86
  87 #ifdef _DEBUG_
  88 #define nodeVec_ACCESS(index) ((SMDS_MeshNode*) nodeVec.at((index)))
  89 #else
  90 #define nodeVec_ACCESS(index) ((SMDS_MeshNode*) nodeVec[index])
  91 #endif
  92
  93 #ifdef NETGEN_NEW
  94 #define NGPOINT_COORDS(p) p(0),p(1),p(2)
  95 #else
  96 #define NGPOINT_COORDS(p) p.X(),p.Y(),p.Z()
  97 #endif
  98
  99 // dump elements added to ng mesh
 100 //#define DUMP_SEGMENTS
 101 //#define DUMP_TRIANGLES
 102 //#define DUMP_TRIANGLES_SCRIPT "/tmp/trias.py" //!< debug addIntVerticesInSolids()
 103
 104 TopTools_IndexedMapOfShape ShapesWithLocalSize;
 105 std::map<int,double> VertexId2LocalSize;
 106 std::map<int,double> EdgeId2LocalSize;
 107 std::map<int,double> FaceId2LocalSize;
 108
 109 //=============================================================================
 110 /*!
 111  *
 112  */
 113 //=============================================================================
 114
 115 NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
 116                                           const TopoDS_Shape& aShape,
 117                                           const bool isVolume)
 118   : _mesh    (mesh),
 119     _shape   (aShape),
 120     _isVolume(isVolume),
 121     _optimize(true),
 122     _fineness(NETGENPlugin_Hypothesis::GetDefaultFineness()),
 123     _simpleHyp(NULL)
 124 {
 125   defaultParameters();
 126   ShapesWithLocalSize.Clear();
 127   VertexId2LocalSize.clear();
 128   EdgeId2LocalSize.clear();
 129   FaceId2LocalSize.clear();
 130 }
 131
 132 //================================================================================
 133 /*!
 134  * \brief Initialize global NETGEN parameters with default values
 135  */
 136 //================================================================================
 137
 138 void NETGENPlugin_Mesher::defaultParameters()
 139 {
 140   netgen::MeshingParameters& mparams = netgen::mparam;
 141   // maximal mesh edge size
 142   mparams.maxh = 0;//NETGENPlugin_Hypothesis::GetDefaultMaxSize();
 143   mparams.minh = 0;
 144   // minimal number of segments per edge
 145   mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
 146   // rate of growth of size between elements
 147   mparams.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
 148   // safety factor for curvatures (elements per radius)
 149   mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
 150   // create elements of second order
 151   mparams.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder() ? 1 : 0;
 152   // quad-dominated surface meshing
 153   if (_isVolume)
 154     mparams.quad = 0;
 155   else
 156     mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0;
 157   _fineness = NETGENPlugin_Hypothesis::GetDefaultFineness();
 158 }
 159
 160 //=============================================================================
 161 /*!
 162  *
 163  */
 164 //=============================================================================
 165 void SetLocalSize(TopoDS_Shape GeomShape, double LocalSize)
 166 {
 167   TopAbs_ShapeEnum GeomType = GeomShape.ShapeType();
 168   if (GeomType == TopAbs_COMPOUND) {
 169     for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()) {
 170       SetLocalSize(it.Value(), LocalSize);
 171     }
 172     return;
 173   }
 174   int key;
 175   if (! ShapesWithLocalSize.Contains(GeomShape))
 176     key = ShapesWithLocalSize.Add(GeomShape);
 177   else
 178     key = ShapesWithLocalSize.FindIndex(GeomShape);
 179   if (GeomType == TopAbs_VERTEX) {
 180     VertexId2LocalSize[key] = LocalSize;
 181   } else if (GeomType == TopAbs_EDGE) {
 182     EdgeId2LocalSize[key] = LocalSize;
 183   } else if (GeomType == TopAbs_FACE) {
 184     FaceId2LocalSize[key] = LocalSize;
 185   }
 186 }
 187
 188 //=============================================================================
 189 /*!
 190  * Pass parameters to NETGEN
 191  */
 192 //=============================================================================
 193 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
 194 {
 195   if (hyp)
 196   {
 197     netgen::MeshingParameters& mparams = netgen::mparam;
 198     // Initialize global NETGEN parameters:
 199     // maximal mesh segment size
 200     mparams.maxh = hyp->GetMaxSize();
 201     // maximal mesh element linear size
 202     mparams.minh = hyp->GetMinSize();
 203     // minimal number of segments per edge
 204     mparams.segmentsperedge = hyp->GetNbSegPerEdge();
 205     // rate of growth of size between elements
 206     mparams.grading = hyp->GetGrowthRate();
 207     // safety factor for curvatures (elements per radius)
 208     mparams.curvaturesafety = hyp->GetNbSegPerRadius();
 209     // create elements of second order
 210     mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0;
 211     // quad-dominated surface meshing
 212     // only triangles are allowed for volumic mesh
 213     if (!_isVolume)
 214       mparams.quad = static_cast<const NETGENPlugin_Hypothesis_2D*>
 215         (hyp)->GetQuadAllowed() ? 1 : 0;
 216     _optimize = hyp->GetOptimize();
 217     _fineness = hyp->GetFineness();
 218     _simpleHyp = NULL;
 219
 220     SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen();
 221     CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager");
 222     SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject);
 223     SALOMEDS::Study_var myStudy = aStudyMgr->GetStudyByID(hyp->GetStudyId());
 224
 225     const NETGENPlugin_Hypothesis::TLocalSize localSizes = hyp->GetLocalSizesAndEntries();
 226     NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin();
 227     for (it ; it != localSizes.end() ; it++)
 228       {
 229         std::string entry = (*it).first;
 230         double val = (*it).second;
 231         // --
 232         GEOM::GEOM_Object_var aGeomObj;
 233         TopoDS_Shape S = TopoDS_Shape();
 234         SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
 235         SALOMEDS::GenericAttribute_var anAttr;
 236         if (!aSObj->_is_nil() && aSObj->FindAttribute(anAttr, "AttributeIOR")) {
 237           SALOMEDS::AttributeIOR_var anIOR = SALOMEDS::AttributeIOR::_narrow(anAttr);
 238           CORBA::String_var aVal = anIOR->Value();
 239           CORBA::Object_var obj = myStudy->ConvertIORToObject(aVal);
 240           aGeomObj = GEOM::GEOM_Object::_narrow(obj);
 241         }
 242         if ( !aGeomObj->_is_nil() )
 243           S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
 244         // --
 245         SetLocalSize(S, val);
 246       }
 247   }
 248 }
 249
 250 //=============================================================================
 251 /*!
 252  * Pass simple parameters to NETGEN
 253  */
 254 //=============================================================================
 255
 256 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp)
 257 {
 258   _simpleHyp = hyp;
 259   if ( _simpleHyp )
 260     defaultParameters();
 261 }
 262
 263 //=============================================================================
 264 /*!
 265  *  Link - a pair of integer numbers
 266  */
 267 //=============================================================================
 268 struct Link
 269 {
 270   int n1, n2;
 271   Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
 272   Link() : n1(0), n2(0) {}
 273 };
 274
 275 int HashCode(const Link& aLink, int aLimit)
 276 {
 277   return HashCode(aLink.n1 + aLink.n2, aLimit);
 278 }
 279
 280 Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
 281 {
 282   return (aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ||
 283           aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1);
 284 }
 285
 286 namespace
 287 {
 288   //================================================================================
 289   /*!
 290    * \brief return id of netgen point corresponding to SMDS node
 291    */
 292   //================================================================================
 293   typedef map< const SMDS_MeshNode*, int > TNode2IdMap;
 294
 295   int ngNodeId( const SMDS_MeshNode* node,
 296                 netgen::Mesh&        ngMesh,
 297                 TNode2IdMap&         nodeNgIdMap)
 298   {
 299     int newNgId = ngMesh.GetNP() + 1;
 300
 301     TNode2IdMap::iterator node_id = nodeNgIdMap.insert( make_pair( node, newNgId )).first;
 302
 303     if ( node_id->second == newNgId)
 304     {
 305 #if defined(DUMP_SEGMENTS) || defined(DUMP_TRIANGLES)
 306       cout << "Ng " << newNgId << " - " << node;
 307 #endif
 308       netgen::MeshPoint p( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
 309       ngMesh.AddPoint( p );
 310     }
 311     return node_id->second;
 312   }
 313
 314   //================================================================================
 315   /*!
 316    * \brief Return computed EDGEs connected to the given one
 317    */
 318   //================================================================================
 319
 320   list< TopoDS_Edge > getConnectedEdges( const TopoDS_Edge&                 edge,
 321                                          const TopoDS_Face&                 face,
 322                                          const set< SMESH_subMesh* > &      computedSM,
 323                                          const SMESH_MesherHelper&          helper,
 324                                          map< SMESH_subMesh*, set< int > >& addedEdgeSM2Faces)
 325   {
 326     // get ordered EDGEs
 327     TopoDS_Vertex v1;
 328     list< TopoDS_Edge > edges;
 329     list< int > nbEdgesInWire;
 330     int nbWires = SMESH_Block::GetOrderedEdges( face, v1, edges, nbEdgesInWire);
 331
 332     // find <edge> within <edges>
 333     list< TopoDS_Edge >::iterator eItFwd = edges.begin();
 334     for ( ; eItFwd != edges.end(); ++eItFwd )
 335       if ( edge.IsSame( *eItFwd ))
 336         break;
 337     if ( eItFwd == edges.end()) return list< TopoDS_Edge>();
 338
 339     if ( eItFwd->Orientation() >= TopAbs_INTERNAL )
 340     {
 341       // connected INTERNAL edges returned from GetOrderedEdges() are wrongly oriented
 342       // so treat each INTERNAL edge separately
 343       TopoDS_Edge e = *eItFwd;
 344       edges.clear();
 345       edges.push_back( e );
 346       return edges;
 347     }
 348
 349     // get all computed EDGEs connected to <edge>
 350
 351     list< TopoDS_Edge >::iterator eItBack = eItFwd, ePrev;
 352     TopoDS_Vertex vCommon;
 353     TopTools_MapOfShape eAdded; // map used not to add a seam edge twice to <edges>
 354     eAdded.Add( edge );
 355
 356     // put edges before <edge> to <edges> back
 357     while ( edges.begin() != eItFwd )
 358       edges.splice( edges.end(), edges, edges.begin() );
 359
 360     // search forward
 361     ePrev = eItFwd;
 362     while ( ++eItFwd != edges.end() )
 363     {
 364       SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( *eItFwd );
 365
 366       bool connected = TopExp::CommonVertex( *ePrev, *eItFwd, vCommon );
 367       bool computed  = sm->IsMeshComputed();
 368       bool added     = addedEdgeSM2Faces[ sm ].count( helper.GetSubShapeID() );
 369       bool doubled   = !eAdded.Add( *eItFwd );
 370       bool orientOK  = (( ePrev ->Orientation() < TopAbs_INTERNAL ) ==
 371                         ( eItFwd->Orientation() < TopAbs_INTERNAL )    );
 372       if ( !connected || !computed || !orientOK || added || doubled )
 373       {
 374         // stop advancement; move edges from tail to head
 375         while ( edges.back() != *ePrev )
 376           edges.splice( edges.begin(), edges, --edges.end() );
 377         break;
 378       }
 379       ePrev = eItFwd;
 380     }
 381     // search backward
 382     while ( eItBack != edges.begin() )
 383     {
 384       ePrev = eItBack;
 385       --eItBack;
 386       SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( *eItBack );
 387
 388       bool connected = TopExp::CommonVertex( *ePrev, *eItBack, vCommon );
 389       bool computed  = sm->IsMeshComputed();
 390       bool added     = addedEdgeSM2Faces[ sm ].count( helper.GetSubShapeID() );
 391       bool doubled   = !eAdded.Add( *eItBack );
 392       bool orientOK  = (( ePrev  ->Orientation() < TopAbs_INTERNAL ) ==
 393                         ( eItBack->Orientation() < TopAbs_INTERNAL )    );
 394       if ( !connected || !computed || !orientOK || added || doubled)
 395       {
 396         // stop advancement
 397         edges.erase( edges.begin(), ePrev );
 398         break;
 399       }
 400     }
 401     if ( edges.front() != edges.back() )
 402     {
 403       // assure that the 1st vertex is meshed
 404       TopoDS_Edge eLast = edges.back();
 405       while ( !SMESH_Algo::VertexNode( SMESH_MesherHelper::IthVertex( 0, edges.front()), helper.GetMeshDS())
 406               &&
 407               edges.front() != eLast )
 408         edges.splice( edges.end(), edges, edges.begin() );
 409     }
 410     return edges;
 411   }
 412
 413   //================================================================================
 414   /*!
 415    * \brief Make triangulation of a shape precise enough
 416    */
 417   //================================================================================
 418
 419   void updateTriangulation( const TopoDS_Shape& shape )
 420   {
 421     // static set< Poly_Triangulation* > updated;
 422
 423     // TopLoc_Location loc;
 424     // TopExp_Explorer fExp( shape, TopAbs_FACE );
 425     // for ( ; fExp.More(); fExp.Next() )
 426     // {
 427     //   Handle(Poly_Triangulation) triangulation =
 428     //     BRep_Tool::Triangulation ( TopoDS::Face( fExp.Current() ), loc);
 429     //   if ( triangulation.IsNull() ||
 430     //        updated.insert( triangulation.operator->() ).second )
 431     //   {
 432     //     BRepTools::Clean (shape);
 433         try {
 434 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 435           OCC_CATCH_SIGNALS;
 436 #endif
 437           BRepMesh_IncrementalMesh e(shape, 0.01, true);
 438
 439         }
 440         catch (Standard_Failure)
 441         {
 442         }
 443   //       updated.erase( triangulation.operator->() );
 444   //       triangulation = BRep_Tool::Triangulation ( TopoDS::Face( fExp.Current() ), loc);
 445   //       updated.insert( triangulation.operator->() );
 446   //     }
 447   //   }
 448   }
 449 }
 450
 451 //================================================================================
 452 /*!
 453  * \brief Initialize netgen::OCCGeometry with OCCT shape
 454  */
 455 //================================================================================
 456
 457 void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry&     occgeo,
 458                                              const TopoDS_Shape&      shape,
 459                                              SMESH_Mesh&              mesh,
 460                                              list< SMESH_subMesh* > * meshedSM,
 461                                              NETGENPlugin_Internals*  intern)
 462 {
 463   updateTriangulation( shape );
 464
 465   Bnd_Box bb;
 466   BRepBndLib::Add (shape, bb);
 467   double x1,y1,z1,x2,y2,z2;
 468   bb.Get (x1,y1,z1,x2,y2,z2);
 469   MESSAGE("shape bounding box:\n" <<
 470           "(" << x1 << " " << y1 << " " << z1 << ") " <<
 471           "(" << x2 << " " << y2 << " " << z2 << ")");
 472   netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
 473   netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
 474   occgeo.boundingbox = netgen::Box<3> (p1,p2);
 475
 476   occgeo.shape = shape;
 477   occgeo.changed = 1;
 478
 479   // fill maps of shapes of occgeo with not yet meshed subshapes
 480
 481   // get root submeshes
 482   list< SMESH_subMesh* > rootSM;
 483   if ( SMESH_subMesh* sm = mesh.GetSubMeshContaining( shape )) {
 484     rootSM.push_back( sm );
 485   }
 486   else {
 487     for ( TopoDS_Iterator it( shape ); it.More(); it.Next() )
 488       rootSM.push_back( mesh.GetSubMesh( it.Value() ));
 489   }
 490
 491   // add subshapes of empty submeshes
 492   list< SMESH_subMesh* >::iterator rootIt = rootSM.begin(), rootEnd = rootSM.end();
 493   for ( ; rootIt != rootEnd; ++rootIt ) {
 494     SMESH_subMesh * root = *rootIt;
 495     SMESH_subMeshIteratorPtr smIt = root->getDependsOnIterator(/*includeSelf=*/true,
 496                                                                /*complexShapeFirst=*/true);
 497     // to find a right orientation of subshapes (PAL20462)
 498     TopTools_IndexedMapOfShape subShapes;
 499     TopExp::MapShapes(root->GetSubShape(), subShapes);
 500     while ( smIt->more() )
 501     {
 502       SMESH_subMesh* sm = smIt->next();
 503       TopoDS_Shape shape = sm->GetSubShape();
 504       if ( intern && intern->isShapeToPrecompute( shape ))
 505         continue;
 506       if ( !meshedSM || sm->IsEmpty() )
 507       {
 508         if ( shape.ShapeType() != TopAbs_VERTEX )
 509           shape = subShapes( subShapes.FindIndex( shape ));// shape -> index -> oriented shape
 510         if ( shape.Orientation() >= TopAbs_INTERNAL )
 511           shape.Orientation( TopAbs_FORWARD ); // isuue 0020676
 512         switch ( shape.ShapeType() ) {
 513         case TopAbs_FACE  : occgeo.fmap.Add( shape ); break;
 514         case TopAbs_EDGE  : occgeo.emap.Add( shape ); break;
 515         case TopAbs_VERTEX: occgeo.vmap.Add( shape ); break;
 516         case TopAbs_SOLID :occgeo.somap.Add( shape ); break;
 517         default:;
 518         }
 519       }
 520       // collect submeshes of meshed shapes
 521       else if (meshedSM)
 522       {
 523         const int dim = SMESH_Gen::GetShapeDim( shape );
 524         meshedSM[ dim ].push_back( sm );
 525       }
 526     }
 527   }
 528   occgeo.facemeshstatus.SetSize (occgeo.fmap.Extent());
 529   occgeo.facemeshstatus = 0;
 530 #ifdef NETGEN_NEW
 531   occgeo.face_maxh_modified.SetSize(occgeo.fmap.Extent());
 532   occgeo.face_maxh_modified = 0;
 533   occgeo.face_maxh.SetSize(occgeo.fmap.Extent());
 534   occgeo.face_maxh = netgen::mparam.maxh;
 535 #endif
 536 }
 537
 538 //================================================================================
 539 /*!
 540  * \brief Return a default min size value suitable for the given geometry.
 541  */
 542 //================================================================================
 543
 544 double NETGENPlugin_Mesher::GetDefaultMinSize(const TopoDS_Shape& geom,
 545                                               const double        maxSize)
 546 {
 547   updateTriangulation( geom );
 548
 549   TopLoc_Location loc;
 550   int i1, i2, i3;
 551   const int* pi[4] = { &i1, &i2, &i3, &i1 };
 552   double minh = 1e100;
 553   Bnd_B3d bb;
 554   TopExp_Explorer fExp( geom, TopAbs_FACE );
 555   for ( ; fExp.More(); fExp.Next() )
 556   {
 557     Handle(Poly_Triangulation) triangulation =
 558       BRep_Tool::Triangulation ( TopoDS::Face( fExp.Current() ), loc);
 559     if ( triangulation.IsNull() ) continue;
 560     const double fTol = BRep_Tool::Tolerance( TopoDS::Face( fExp.Current() ));
 561     const TColgp_Array1OfPnt&   points = triangulation->Nodes();
 562     const Poly_Array1OfTriangle& trias = triangulation->Triangles();
 563     for ( int iT = trias.Lower(); iT <= trias.Upper(); ++iT )
 564     {
 565       trias(iT).Get( i1, i2, i3 );
 566       for ( int j = 0; j < 3; ++j )
 567       {
 568         double dist2 = points(*pi[j]).SquareDistance( points( *pi[j+1] ));
 569         if ( dist2 < minh && fTol*fTol < dist2 )
 570           minh = dist2;
 571         bb.Add( points(*pi[j]));
 572       }
 573     }
 574   }
 575   if ( minh > 0.25 * bb.SquareExtent() ) // simple geometry, rough triangulation
 576   {
 577     minh = 1e-3 * sqrt( bb.SquareExtent());
 578     //cout << "BND BOX minh = " <<minh << endl;
 579   }
 580   else
 581   {
 582     minh = 3 * sqrt( minh ); // triangulation for visualization is rather fine
 583     //cout << "TRIANGULATION minh = " <<minh << endl;
 584   }
 585   if ( minh > 0.5 * maxSize )
 586     minh = maxSize / 3.;
 587
 588   return minh;
 589 }
 590
 591 //================================================================================
 592 /*!
 593  * \brief Restrict size of elements at a given point
 594  */
 595 //================================================================================
 596
 597 void NETGENPlugin_Mesher::RestrictLocalSize(netgen::Mesh& ngMesh, const gp_XYZ& p, const double size)
 598 {
 599   if ( netgen::mparam.minh > size )
 600   {
 601     ngMesh.SetMinimalH( size );
 602     netgen::mparam.minh = size;
 603   }
 604   netgen::Point3d pi(p.X(), p.Y(), p.Z());
 605   ngMesh.RestrictLocalH( pi, size );
 606 }
 607
 608 //================================================================================
 609 /*!
 610  * \brief fill ngMesh with nodes and elements of computed submeshes
 611  */
 612 //================================================================================
 613
 614 bool NETGENPlugin_Mesher::fillNgMesh(const netgen::OCCGeometry&     occgeom,
 615                                      netgen::Mesh&                  ngMesh,
 616                                      vector<const SMDS_MeshNode*>&  nodeVec,
 617                                      const list< SMESH_subMesh* > & meshedSM)
 618 {
 619   TNode2IdMap nodeNgIdMap;
 620   for ( int i = 1; i < nodeVec.size(); ++i )
 621     nodeNgIdMap.insert( make_pair( nodeVec[i], i ));
 622
 623   TopTools_MapOfShape visitedShapes;
 624   map< SMESH_subMesh*, set< int > > visitedEdgeSM2Faces;
 625   set< SMESH_subMesh* > computedSM( meshedSM.begin(), meshedSM.end() );
 626
 627   SMESH_MesherHelper helper (*_mesh);
 628
 629   int faceNgID = occgeom.fmap.Extent();
 630
 631   list< SMESH_subMesh* >::const_iterator smIt, smEnd = meshedSM.end();
 632   for ( smIt = meshedSM.begin(); smIt != smEnd; ++smIt )
 633   {
 634     SMESH_subMesh* sm = *smIt;
 635     if ( !visitedShapes.Add( sm->GetSubShape() ))
 636       continue;
 637
 638     SMESHDS_SubMesh * smDS = sm->GetSubMeshDS();
 639     if ( !smDS ) continue;
 640
 641     switch ( sm->GetSubShape().ShapeType() )
 642     {
 643     case TopAbs_EDGE: { // EDGE
 644       // ----------------------
 645       TopoDS_Edge geomEdge  = TopoDS::Edge( sm->GetSubShape() );
 646       if ( geomEdge.Orientation() >= TopAbs_INTERNAL )
 647         geomEdge.Orientation( TopAbs_FORWARD ); // issue 0020676
 648
 649       // Add ng segments for each not meshed FACE the EDGE bounds
 650       PShapeIteratorPtr fIt = helper.GetAncestors( geomEdge, *sm->GetFather(), TopAbs_FACE );
 651       while ( const TopoDS_Shape * anc = fIt->next() )
 652       {
 653         faceNgID = occgeom.fmap.FindIndex( *anc );
 654         if ( faceNgID < 1 )
 655           continue; // meshed face
 656
 657         int faceSMDSId = helper.GetMeshDS()->ShapeToIndex( *anc );
 658         if ( visitedEdgeSM2Faces[ sm ].count( faceSMDSId ))
 659           continue; // already treated EDGE
 660
 661         TopoDS_Face face = TopoDS::Face( occgeom.fmap( faceNgID ));
 662         if ( face.Orientation() >= TopAbs_INTERNAL )
 663           face.Orientation( TopAbs_FORWARD ); // issue 0020676
 664
 665         // get all meshed EDGEs of the FACE connected to geomEdge (issue 0021140)
 666         helper.SetSubShape( face );
 667         list< TopoDS_Edge > edges = getConnectedEdges( geomEdge, face, computedSM, helper,
 668                                                        visitedEdgeSM2Faces );
 669         if ( edges.empty() )
 670           continue; // wrong ancestor?
 671
 672         // find out orientation of <edges> within <face>
 673         TopoDS_Edge eNotSeam = edges.front();
 674         if ( helper.HasSeam() )
 675         {
 676           list< TopoDS_Edge >::iterator eIt = edges.begin();
 677           while ( helper.IsRealSeam( *eIt )) ++eIt;
 678           if ( eIt != edges.end() )
 679             eNotSeam = *eIt;
 680         }
 681         TopAbs_Orientation fOri = helper.GetSubShapeOri( face, eNotSeam );
 682         bool isForwad = ( fOri == eNotSeam.Orientation() || fOri >= TopAbs_INTERNAL );
 683
 684         // get all nodes from connected <edges>
 685         bool isQuad   = smDS->NbElements() ? smDS->GetElements()->next()->IsQuadratic() : false;
 686         StdMeshers_FaceSide fSide( face, edges, _mesh, isForwad, isQuad );
 687         const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
 688         int i, nbSeg = fSide.NbSegments();
 689
 690         // remember EDGEs of fSide to treat only once
 691         for ( int iE = 0; iE < fSide.NbEdges(); ++iE )
 692           visitedEdgeSM2Faces[ helper.GetMesh()->GetSubMesh( fSide.Edge(iE )) ].insert(faceSMDSId);
 693
 694         double otherSeamParam = 0;
 695         bool isSeam = false;
 696
 697         // add segments
 698
 699         int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
 700
 701         for ( i = 0; i < nbSeg; ++i )
 702         {
 703           const UVPtStruct& p1 = points[ i ];
 704           const UVPtStruct& p2 = points[ i+1 ];
 705
 706           if ( p1.node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX ) //an EDGE begins
 707           {
 708             isSeam = false;
 709             if ( helper.IsRealSeam( p1.node->getshapeId() ))
 710             {
 711               TopoDS_Edge e = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
 712               isSeam = helper.IsRealSeam( e );
 713               if ( isSeam )
 714               {
 715                 otherSeamParam = helper.GetOtherParam( helper.GetPeriodicIndex() & 1 ? p2.u : p2.v );
 716               }
 717             }
 718           }
 719           netgen::Segment seg;
 720           // ng node ids
 721           seg[0] = prevNgId;
 722           seg[1] = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
 723           // node param on curve
 724           seg.epgeominfo[ 0 ].dist = p1.param;
 725           seg.epgeominfo[ 1 ].dist = p2.param;
 726           // uv on face
 727           seg.epgeominfo[ 0 ].u = p1.u;
 728           seg.epgeominfo[ 0 ].v = p1.v;
 729           seg.epgeominfo[ 1 ].u = p2.u;
 730           seg.epgeominfo[ 1 ].v = p2.v;
 731
 732           //geomEdge = fSide.Edge( fSide.EdgeIndex( 0.5 * ( p1.normParam + p2.normParam )));
 733           //seg.epgeominfo[ 0 ].edgenr = seg.epgeominfo[ 1 ].edgenr = occgeom.emap.FindIndex( geomEdge );
 734
 735           //seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
 736           seg.si = faceNgID;                   // = geom.fmap.FindIndex (face);
 737           seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
 738           ngMesh.AddSegment (seg);
 739
 740           SMESH_TNodeXYZ np1( p1.node ), np2( p2.node );
 741           RestrictLocalSize( ngMesh, 0.5*(np1+np2), (np1-np2).Modulus() );
 742
 743 #ifdef DUMP_SEGMENTS
 744           cout << "Segment: " << seg.edgenr << " on SMESH face " << helper.GetMeshDS()->ShapeToIndex( face ) << endl
 745                << "\tface index: " << seg.si << endl
 746                << "\tp1: " << seg[0] << endl
 747                << "\tp2: " << seg[1] << endl
 748                << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
 749                << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
 750             //<< "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
 751                << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
 752                << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl;
 753             //<< "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
 754 #endif
 755           if ( isSeam )
 756           {
 757             if ( helper.GetPeriodicIndex() && 1 ) {
 758               seg.epgeominfo[ 0 ].u = otherSeamParam;
 759               seg.epgeominfo[ 1 ].u = otherSeamParam;
 760               swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
 761             } else {
 762               seg.epgeominfo[ 0 ].v = otherSeamParam;
 763               seg.epgeominfo[ 1 ].v = otherSeamParam;
 764               swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
 765             }
 766             swap (seg[0], seg[1]);
 767             swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
 768             seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
 769             ngMesh.AddSegment (seg);
 770 #ifdef DUMP_SEGMENTS
 771             cout << "Segment: " << seg.edgenr << endl
 772                  << "\t is SEAM (reverse) of the previous. "
 773                  << " Other " << (helper.GetPeriodicIndex() && 1 ? "U" : "V")
 774                  << " = " << otherSeamParam << endl;
 775 #endif
 776           }
 777           else if ( fOri == TopAbs_INTERNAL )
 778           {
 779             swap (seg[0], seg[1]);
 780             swap( seg.epgeominfo[0], seg.epgeominfo[1] );
 781             seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
 782             ngMesh.AddSegment (seg);
 783 #ifdef DUMP_SEGMENTS
 784             cout << "Segment: " << seg.edgenr << endl << "\t is REVERSE of the previous" << endl;
 785 #endif
 786           }
 787         }
 788       } // loop on geomEdge ancestors
 789
 790       break;
 791     } // case TopAbs_EDGE
 792
 793     case TopAbs_FACE: { // FACE
 794       // ----------------------
 795       const TopoDS_Face& geomFace  = TopoDS::Face( sm->GetSubShape() );
 796       helper.SetSubShape( geomFace );
 797       bool isInternalFace = ( geomFace.Orientation() == TopAbs_INTERNAL );
 798
 799       // Find solids the geomFace bounds
 800       int solidID1 = 0, solidID2 = 0;
 801       PShapeIteratorPtr solidIt = helper.GetAncestors( geomFace, *sm->GetFather(), TopAbs_SOLID);
 802       while ( const TopoDS_Shape * solid = solidIt->next() )
 803       {
 804         int id = occgeom.somap.FindIndex ( *solid );
 805         if ( solidID1 && id != solidID1 ) solidID2 = id;
 806         else                              solidID1 = id;
 807       }
 808       faceNgID++;
 809       //_faceDescriptors[ faceNgID ].first  = solidID1;
 810       //_faceDescriptors[ faceNgID ].second = solidID2;
 811       // Add ng face descriptors of meshed faces
 812       ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(faceNgID, solidID1, solidID2, 0));
 813
 814       // Orient the face correctly in solidID1 (issue 0020206)
 815       bool reverse = false;
 816       if ( solidID1 ) {
 817         TopoDS_Shape solid = occgeom.somap( solidID1 );
 818         TopAbs_Orientation faceOriInSolid = helper.GetSubShapeOri( solid, geomFace );
 819         if ( faceOriInSolid >= 0 )
 820           reverse = SMESH_Algo::IsReversedSubMesh
 821             ( TopoDS::Face( geomFace.Oriented( faceOriInSolid )), helper.GetMeshDS() );
 822       }
 823
 824       // Add surface elements
 825
 826       netgen::Element2d tri(3);
 827       tri.SetIndex ( faceNgID );
 828
 829
 830 #ifdef DUMP_TRIANGLES
 831       cout << "SMESH face " << helper.GetMeshDS()->ShapeToIndex( geomFace )
 832            << " internal="<<isInternalFace<< " border="<<isBorderFace << endl;
 833 #endif
 834       SMDS_ElemIteratorPtr faces = smDS->GetElements();
 835       while ( faces->more() )
 836       {
 837         const SMDS_MeshElement* f = faces->next();
 838         if ( f->NbNodes() % 3 != 0 ) // not triangle
 839         {
 840           PShapeIteratorPtr solidIt=helper.GetAncestors(geomFace,*sm->GetFather(),TopAbs_SOLID);
 841           if ( const TopoDS_Shape * solid = solidIt->next() )
 842             sm = _mesh->GetSubMesh( *solid );
 843           SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 844           smError.reset( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"Not triangle submesh"));
 845           smError->myBadElements.push_back( f );
 846           return false;
 847         }
 848
 849         for ( int i = 0; i < 3; ++i )
 850         {
 851           const SMDS_MeshNode* node = f->GetNode( i ), * inFaceNode=0;
 852
 853           // get node UV on face
 854           int shapeID = node->getshapeId();
 855           if ( helper.IsSeamShape( shapeID ))
 856             if ( helper.IsSeamShape( f->GetNodeWrap( i+1 )->getshapeId() ))
 857               inFaceNode = f->GetNodeWrap( i-1 );
 858             else
 859               inFaceNode = f->GetNodeWrap( i+1 );
 860           gp_XY uv = helper.GetNodeUV( geomFace, node, inFaceNode );
 861
 862           int ind = reverse ? 3-i : i+1;
 863           tri.GeomInfoPi(ind).u = uv.X();
 864           tri.GeomInfoPi(ind).v = uv.Y();
 865           tri.PNum      (ind) = ngNodeId( node, ngMesh, nodeNgIdMap );
 866         }
 867
 868         ngMesh.AddSurfaceElement (tri);
 869 #ifdef DUMP_TRIANGLES
 870         cout << tri << endl;
 871 #endif
 872
 873         if ( isInternalFace )
 874         {
 875           swap( tri[1], tri[2] );
 876           ngMesh.AddSurfaceElement (tri);
 877 #ifdef DUMP_TRIANGLES
 878         cout << tri << endl;
 879 #endif
 880         }
 881       }
 882       break;
 883     } // case TopAbs_FACE
 884
 885     case TopAbs_VERTEX: { // VERTEX
 886       // --------------------------
 887       // issue 0021405. Add node only if a VERTEX is shared by a not meshed EDGE,
 888       // else netgen removes a free node and nodeVector becomes invalid
 889       PShapeIteratorPtr ansIt = helper.GetAncestors( sm->GetSubShape(),
 890                                                      *sm->GetFather(),
 891                                                      TopAbs_EDGE );
 892       bool toAdd = false;
 893       while ( const TopoDS_Shape* e = ansIt->next() )
 894       {
 895         SMESH_subMesh* eSub = helper.GetMesh()->GetSubMesh( *e );
 896         if (( toAdd = eSub->IsEmpty() )) break;
 897       }
 898       if ( toAdd )
 899       {
 900         SMDS_NodeIteratorPtr nodeIt = smDS->GetNodes();
 901         if ( nodeIt->more() )
 902           ngNodeId( nodeIt->next(), ngMesh, nodeNgIdMap );
 903       }
 904       break;
 905     }
 906     default:;
 907     } // switch
 908   } // loop on submeshes
 909
 910   // fill nodeVec
 911   nodeVec.resize( ngMesh.GetNP() + 1 );
 912   TNode2IdMap::iterator node_NgId, nodeNgIdEnd = nodeNgIdMap.end();
 913   for ( node_NgId = nodeNgIdMap.begin(); node_NgId != nodeNgIdEnd; ++node_NgId)
 914     nodeVec[ node_NgId->second ] = node_NgId->first;
 915
 916   return true;
 917 }
 918
 919 //================================================================================
 920 /*!
 921  * \brief Duplicate mesh faces on internal geom faces
 922  */
 923 //================================================================================
 924
 925 void NETGENPlugin_Mesher::fixIntFaces(const netgen::OCCGeometry& occgeom,
 926                                       netgen::Mesh&              ngMesh,
 927                                       NETGENPlugin_Internals&    internalShapes)
 928 {
 929   SMESHDS_Mesh* meshDS = internalShapes.getMesh().GetMeshDS();
 930
 931   // find ng indices of internal faces
 932   set<int> ngFaceIds;
 933   for ( int ngFaceID = 1; ngFaceID <= occgeom.fmap.Extent(); ++ngFaceID )
 934   {
 935     int smeshID = meshDS->ShapeToIndex( occgeom.fmap( ngFaceID ));
 936     if ( internalShapes.isInternalShape( smeshID ))
 937       ngFaceIds.insert( ngFaceID );
 938   }
 939   if ( !ngFaceIds.empty() )
 940   {
 941     // duplicate faces
 942     int i, nbFaces = ngMesh.GetNSE();
 943     for (int i = 1; i <= nbFaces; ++i)
 944     {
 945       netgen::Element2d elem = ngMesh.SurfaceElement(i);
 946       if ( ngFaceIds.count( elem.GetIndex() ))
 947       {
 948         swap( elem[1], elem[2] );
 949         ngMesh.AddSurfaceElement (elem);
 950       }
 951     }
 952   }
 953 }
 954
 955 namespace
 956 {
 957   //================================================================================
 958   // define gp_XY_Subtracted pointer to function calling gp_XY::Subtracted(gp_XY)
 959   gp_XY_FunPtr(Subtracted);
 960   //gp_XY_FunPtr(Added);
 961
 962   //================================================================================
 963   /*!
 964    * \brief Evaluate distance between two 2d points along the surface
 965    */
 966   //================================================================================
 967
 968   double evalDist( const gp_XY&                uv1,
 969                    const gp_XY&                uv2,
 970                    const Handle(Geom_Surface)& surf,
 971                    const int                   stopHandler=-1)
 972   {
 973     if ( stopHandler > 0 ) // continue recursion
 974     {
 975       gp_XY mid = SMESH_MesherHelper::GetMiddleUV( surf, uv1, uv2 );
 976       return evalDist( uv1,mid, surf, stopHandler-1 ) + evalDist( mid,uv2, surf, stopHandler-1 );
 977     }
 978     double dist3D = surf->Value( uv1.X(), uv1.Y() ).Distance( surf->Value( uv2.X(), uv2.Y() ));
 979     if ( stopHandler == 0 ) // stop recursion
 980       return dist3D;
 981
 982     // start recursion if necessary
 983     double dist2D = SMESH_MesherHelper::applyIn2D(surf, uv1, uv2, gp_XY_Subtracted, 0).Modulus();
 984     if ( fabs( dist3D - dist2D ) < dist2D * 1e-10 )
 985       return dist3D; // equal parametrization of a planar surface
 986
 987     return evalDist( uv1, uv2, surf, 3 ); // start recursion
 988   }
 989
 990   //================================================================================
 991   /*!
 992    * \brief Data of vertex internal in geom face
 993    */
 994   //================================================================================
 995
 996   struct TIntVData
 997   {
 998     gp_XY uv;        //!< UV in face parametric space
 999     int   ngId;      //!< ng id of corrsponding node
1000     gp_XY uvClose;   //!< UV of closest boundary node
1001     int   ngIdClose; //!< ng id of closest boundary node
1002   };
1003
1004   //================================================================================
1005   /*!
1006    * \brief Data of vertex internal in solid
1007    */
1008   //================================================================================
1009
1010   struct TIntVSoData
1011   {
1012     int   ngId;      //!< ng id of corresponding node
1013     int   ngIdClose; //!< ng id of closest 2d mesh element
1014     int   ngIdCloseN; //!< ng id of closest node of the closest 2d mesh element
1015   };
1016
1017   inline double dist2(const netgen::MeshPoint& p1, const netgen::MeshPoint& p2)
1018   {
1019     return gp_Pnt( NGPOINT_COORDS(p1)).SquareDistance( gp_Pnt( NGPOINT_COORDS(p2)));
1020   }
1021 }
1022
1023 //================================================================================
1024 /*!
1025  * \brief Make netgen take internal vertices in faces into account by adding
1026  *        segments including internal vertices
1027  *
1028  * This function works in supposition that 1D mesh is already computed in ngMesh
1029  */
1030 //================================================================================
1031
1032 void NETGENPlugin_Mesher::addIntVerticesInFaces(const netgen::OCCGeometry&     occgeom,
1033                                                 netgen::Mesh&                  ngMesh,
1034                                                 vector<const SMDS_MeshNode*>&  nodeVec,
1035                                                 NETGENPlugin_Internals&        internalShapes)
1036 {
1037   if ( nodeVec.size() < ngMesh.GetNP() )
1038     nodeVec.resize( ngMesh.GetNP(), 0 );
1039
1040   SMESHDS_Mesh* meshDS = internalShapes.getMesh().GetMeshDS();
1041   SMESH_MesherHelper helper( internalShapes.getMesh() );
1042
1043   const map<int,list<int> >& face2Vert = internalShapes.getFacesWithVertices();
1044   map<int,list<int> >::const_iterator f2v = face2Vert.begin();
1045   for ( ; f2v != face2Vert.end(); ++f2v )
1046   {
1047     const TopoDS_Face& face = TopoDS::Face( meshDS->IndexToShape( f2v->first ));
1048     if ( face.IsNull() ) continue;
1049     int faceNgID = occgeom.fmap.FindIndex (face);
1050     if ( faceNgID < 0 ) continue;
1051
1052     TopLoc_Location loc;
1053     Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc);
1054
1055     helper.SetSubShape( face );
1056     helper.SetElementsOnShape( true );
1057
1058     // Get data of internal vertices and add them to ngMesh
1059
1060     multimap< double, TIntVData > dist2VData; // sort vertices by distance from boundary nodes
1061
1062     int i, nbSegInit = ngMesh.GetNSeg();
1063
1064     // boundary characteristics
1065     double totSegLen2D = 0;
1066     int totNbSeg = 0;
1067
1068     const list<int>& iVertices = f2v->second;
1069     list<int>::const_iterator iv = iVertices.begin();
1070     for ( int nbV = 0; iv != iVertices.end(); ++iv, nbV++ )
1071     {
1072       TIntVData vData;
1073       // get node on vertex
1074       const TopoDS_Vertex V = TopoDS::Vertex( meshDS->IndexToShape( *iv ));
1075       const SMDS_MeshNode * nV = SMESH_Algo::VertexNode( V, meshDS );
1076       if ( !nV )
1077       {
1078         SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( V );
1079         sm->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1080         nV = SMESH_Algo::VertexNode( V, meshDS );
1081         if ( !nV ) continue;
1082       }
1083       // add ng node
1084       netgen::MeshPoint mp( netgen::Point<3> (nV->X(), nV->Y(), nV->Z()) );
1085       ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
1086       vData.ngId = ngMesh.GetNP();
1087       nodeVec.push_back( nV );
1088
1089       // get node UV
1090       bool uvOK = false;
1091       vData.uv = helper.GetNodeUV( face, nV, 0, &uvOK );
1092       if ( !uvOK ) helper.CheckNodeUV( face, nV, vData.uv, BRep_Tool::Tolerance(V),/*force=*/1);
1093
1094       // loop on all segments of the face to find the node closest to vertex and to count
1095       // average segment 2d length
1096       double closeDist2 = numeric_limits<double>::max(), dist2;
1097       int ngIdLast = 0;
1098       for (i = 1; i <= ngMesh.GetNSeg(); ++i)
1099       {
1100         netgen::Segment & seg = ngMesh.LineSegment(i);
1101         if ( seg.si != faceNgID ) continue;
1102         gp_XY uv[2];
1103         for ( int iEnd = 0; iEnd < 2; ++iEnd)
1104         {
1105           uv[iEnd].SetCoord( seg.epgeominfo[iEnd].u, seg.epgeominfo[iEnd].v );
1106           if ( ngIdLast == seg[ iEnd ] ) continue;
1107           dist2 = helper.applyIn2D(surf, uv[iEnd], vData.uv, gp_XY_Subtracted,0).SquareModulus();
1108           if ( dist2 < closeDist2 )
1109             vData.ngIdClose = seg[ iEnd ], vData.uvClose = uv[iEnd], closeDist2 = dist2;
1110           ngIdLast = seg[ iEnd ];
1111         }
1112         if ( !nbV )
1113         {
1114           totSegLen2D += helper.applyIn2D(surf, uv[0], uv[1], gp_XY_Subtracted, false).Modulus();
1115           totNbSeg++;
1116         }
1117       }
1118       dist2VData.insert( make_pair( closeDist2, vData ));
1119     }
1120
1121     if ( totNbSeg == 0 ) break;
1122     double avgSegLen2d = totSegLen2D / totNbSeg;
1123
1124     // Loop on vertices to add segments
1125
1126     multimap< double, TIntVData >::iterator dist_vData = dist2VData.begin();
1127     for ( ; dist_vData != dist2VData.end(); ++dist_vData )
1128     {
1129       double closeDist2 = dist_vData->first, dist2;
1130       TIntVData & vData = dist_vData->second;
1131
1132       // try to find more close node among segments added for internal vertices
1133       for (i = nbSegInit+1; i <= ngMesh.GetNSeg(); ++i)
1134       {
1135         netgen::Segment & seg = ngMesh.LineSegment(i);
1136         if ( seg.si != faceNgID ) continue;
1137         gp_XY uv[2];
1138         for ( int iEnd = 0; iEnd < 2; ++iEnd)
1139         {
1140           uv[iEnd].SetCoord( seg.epgeominfo[iEnd].u, seg.epgeominfo[iEnd].v );
1141           dist2 = helper.applyIn2D(surf, uv[iEnd], vData.uv, gp_XY_Subtracted,0).SquareModulus();
1142           if ( dist2 < closeDist2 )
1143             vData.ngIdClose = seg[ iEnd ], vData.uvClose = uv[iEnd], closeDist2 = dist2;
1144         }
1145       }
1146       // decide whether to use the closest node as the second end of segment or to
1147       // create a new point
1148       int segEnd1 = vData.ngId;
1149       int segEnd2 = vData.ngIdClose; // to use closest node
1150       gp_XY uvV = vData.uv, uvP = vData.uvClose;
1151       double segLenHint  = ngMesh.GetH( ngMesh.Point( vData.ngId ));
1152       double nodeDist2D  = sqrt( closeDist2 );
1153       double nodeDist3D  = evalDist( vData.uv, vData.uvClose, surf );
1154       bool avgLenOK  = ( avgSegLen2d < 0.75 * nodeDist2D );
1155       bool hintLenOK = ( segLenHint  < 0.75 * nodeDist3D );
1156       //cout << "uvV " << uvV.X() <<","<<uvV.Y() << " ";
1157       if ( hintLenOK || avgLenOK )
1158       {
1159         // create a point between the closest node and V
1160
1161         // how far from V
1162         double r = min( 0.5, ( hintLenOK ? segLenHint/nodeDist3D : avgSegLen2d/nodeDist2D ));
1163         // direction from V to closet node in 2D
1164         gp_Dir2d v2n( helper.applyIn2D(surf, uvP, uvV, gp_XY_Subtracted, false ));
1165         // new point
1166         uvP = vData.uv + r * nodeDist2D * v2n.XY();
1167         gp_Pnt P = surf->Value( uvP.X(), uvP.Y() ).Transformed( loc );
1168
1169         netgen::MeshPoint mp( netgen::Point<3> (P.X(), P.Y(), P.Z()));
1170         ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
1171         segEnd2 = ngMesh.GetNP();
1172         //cout << "Middle " << r << " uv " << uvP.X() << "," << uvP.Y() << "( " << ngMesh.Point(segEnd2).X()<<","<<ngMesh.Point(segEnd2).Y()<<","<<ngMesh.Point(segEnd2).Z()<<" )"<< endl;
1173         SMDS_MeshNode * nP = helper.AddNode(P.X(), P.Y(), P.Z());
1174         nodeVec.push_back( nP );
1175       }
1176       //else cout << "at Node " << " uv " << uvP.X() << "," << uvP.Y() << endl;
1177
1178       // Add the segment
1179       netgen::Segment seg;
1180
1181       if ( segEnd1 > segEnd2 ) swap( segEnd1, segEnd2 ), swap( uvV, uvP );
1182       seg[0] = segEnd1;  // ng node id
1183       seg[1] = segEnd2;  // ng node id
1184       seg.edgenr = ngMesh.GetNSeg() + 1;// segment id
1185       seg.si = faceNgID;
1186
1187       seg.epgeominfo[ 0 ].dist = 0; // param on curve
1188       seg.epgeominfo[ 0 ].u    = uvV.X();
1189       seg.epgeominfo[ 0 ].v    = uvV.Y();
1190       seg.epgeominfo[ 1 ].dist = 1; // param on curve
1191       seg.epgeominfo[ 1 ].u    = uvP.X();
1192       seg.epgeominfo[ 1 ].v    = uvP.Y();
1193
1194 //       seg.epgeominfo[ 0 ].edgenr = 10; //  = geom.emap.FindIndex(edge);
1195 //       seg.epgeominfo[ 1 ].edgenr = 10; //  = geom.emap.FindIndex(edge);
1196
1197       ngMesh.AddSegment (seg);
1198
1199       // add reverse segment
1200       swap (seg[0], seg[1]);
1201       swap( seg.epgeominfo[0], seg.epgeominfo[1] );
1202       seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
1203       ngMesh.AddSegment (seg);
1204     }
1205
1206   }
1207 }
1208
1209 //================================================================================
1210 /*!
1211  * \brief Make netgen take internal vertices in solids into account by adding
1212  *        faces including internal vertices
1213  *
1214  * This function works in supposition that 2D mesh is already computed in ngMesh
1215  */
1216 //================================================================================
1217
1218 void NETGENPlugin_Mesher::addIntVerticesInSolids(const netgen::OCCGeometry&     occgeom,
1219                                                  netgen::Mesh&                  ngMesh,
1220                                                  vector<const SMDS_MeshNode*>&  nodeVec,
1221                                                  NETGENPlugin_Internals&        internalShapes)
1222 {
1223 #ifdef DUMP_TRIANGLES_SCRIPT
1224   // create a python script making a mesh containing triangles added for internal vertices
1225   ofstream py(DUMP_TRIANGLES_SCRIPT);
1226   py << "from smesh import * "<< endl
1227      << "m = Mesh(name='triangles')" << endl;
1228 #endif
1229   if ( nodeVec.size() < ngMesh.GetNP() )
1230     nodeVec.resize( ngMesh.GetNP(), 0 );
1231
1232   SMESHDS_Mesh* meshDS = internalShapes.getMesh().GetMeshDS();
1233   SMESH_MesherHelper helper( internalShapes.getMesh() );
1234
1235   const map<int,list<int> >& so2Vert = internalShapes.getSolidsWithVertices();
1236   map<int,list<int> >::const_iterator s2v = so2Vert.begin();
1237   for ( ; s2v != so2Vert.end(); ++s2v )
1238   {
1239     const TopoDS_Shape& solid = meshDS->IndexToShape( s2v->first );
1240     if ( solid.IsNull() ) continue;
1241     int solidNgID = occgeom.somap.FindIndex (solid);
1242     if ( solidNgID < 0 && !occgeom.somap.IsEmpty() ) continue;
1243
1244     helper.SetSubShape( solid );
1245     helper.SetElementsOnShape( true );
1246
1247     // find ng indices of faces within the solid
1248     set<int> ngFaceIds;
1249     for (TopExp_Explorer fExp(solid, TopAbs_FACE); fExp.More(); fExp.Next() )
1250       ngFaceIds.insert( occgeom.fmap.FindIndex( fExp.Current() ));
1251     if ( ngFaceIds.size() == 1 && *ngFaceIds.begin() == 0 )
1252       ngFaceIds.insert( 1 );
1253
1254     // Get data of internal vertices and add them to ngMesh
1255
1256     multimap< double, TIntVSoData > dist2VData; // sort vertices by distance from ng faces
1257
1258     int i, nbFaceInit = ngMesh.GetNSE();
1259
1260     // boundary characteristics
1261     double totSegLen = 0;
1262     int totNbSeg = 0;
1263
1264     const list<int>& iVertices = s2v->second;
1265     list<int>::const_iterator iv = iVertices.begin();
1266     for ( int nbV = 0; iv != iVertices.end(); ++iv, nbV++ )
1267     {
1268       TIntVSoData vData;
1269       const TopoDS_Vertex V = TopoDS::Vertex( meshDS->IndexToShape( *iv ));
1270
1271       // get node on vertex
1272       const SMDS_MeshNode * nV = SMESH_Algo::VertexNode( V, meshDS );
1273       if ( !nV )
1274       {
1275         SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( V );
1276         sm->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1277         nV = SMESH_Algo::VertexNode( V, meshDS );
1278         if ( !nV ) continue;
1279       }
1280       // add ng node
1281       netgen::MeshPoint mpV( netgen::Point<3> (nV->X(), nV->Y(), nV->Z()) );
1282       ngMesh.AddPoint ( mpV, 1, netgen::FIXEDPOINT );
1283       vData.ngId = ngMesh.GetNP();
1284       nodeVec.push_back( nV );
1285
1286       // loop on all 2d elements to find the one closest to vertex and to count
1287       // average segment length
1288       double closeDist2 = numeric_limits<double>::max(), avgDist2;
1289       for (i = 1; i <= ngMesh.GetNSE(); ++i)
1290       {
1291         const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
1292         if ( !ngFaceIds.count( elem.GetIndex() )) continue;
1293         avgDist2 = 0;
1294         multimap< double, int> dist2nID; // sort nodes of element by distance from V
1295         for ( int j = 0; j < elem.GetNP(); ++j)
1296         {
1297           netgen::MeshPoint mp = ngMesh.Point( elem[j] );
1298           double d2 = dist2( mpV, mp );
1299           dist2nID.insert( make_pair( d2, elem[j] ));
1300           avgDist2 += d2 / elem.GetNP();
1301           if ( !nbV )
1302             totNbSeg++, totSegLen+= sqrt( dist2( mp, ngMesh.Point( elem[(j+1)%elem.GetNP()])));
1303         }
1304         double dist = dist2nID.begin()->first; //avgDist2;
1305         if ( dist < closeDist2 )
1306           vData.ngIdClose= i, vData.ngIdCloseN= dist2nID.begin()->second, closeDist2= dist;
1307       }
1308       dist2VData.insert( make_pair( closeDist2, vData ));
1309     }
1310
1311     if ( totNbSeg == 0 ) break;
1312     double avgSegLen = totSegLen / totNbSeg;
1313
1314     // Loop on vertices to add triangles
1315
1316     multimap< double, TIntVSoData >::iterator dist_vData = dist2VData.begin();
1317     for ( ; dist_vData != dist2VData.end(); ++dist_vData )
1318     {
1319       double closeDist2   = dist_vData->first;
1320       TIntVSoData & vData = dist_vData->second;
1321
1322       const netgen::MeshPoint& mpV = ngMesh.Point( vData.ngId );
1323
1324       // try to find more close face among ones added for internal vertices
1325       for (i = nbFaceInit+1; i <= ngMesh.GetNSE(); ++i)
1326       {
1327         double avgDist2 = 0;
1328         multimap< double, int> dist2nID;
1329         const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
1330         for ( int j = 0; j < elem.GetNP(); ++j)
1331         {
1332           double d = dist2( mpV, ngMesh.Point( elem[j] ));
1333           dist2nID.insert( make_pair( d, elem[j] ));
1334           avgDist2 += d / elem.GetNP();
1335           if ( avgDist2 < closeDist2 )
1336             vData.ngIdClose= i, vData.ngIdCloseN= dist2nID.begin()->second, closeDist2= avgDist2;
1337         }
1338       }
1339       // sort nodes of the closest face by angle with vector from V to the closest node
1340       const double tol = numeric_limits<double>::min();
1341       map< double, int > angle2ID;
1342       const netgen::Element2d& closeFace = ngMesh.SurfaceElement( vData.ngIdClose );
1343       netgen::MeshPoint mp[2];
1344       mp[0] = ngMesh.Point( vData.ngIdCloseN );
1345       gp_XYZ p1( NGPOINT_COORDS( mp[0] ));
1346       gp_XYZ pV( NGPOINT_COORDS( mpV ));
1347       gp_Vec v2p1( pV, p1 );
1348       double distN1 = v2p1.Magnitude();
1349       if ( distN1 <= tol ) continue;
1350       v2p1 /= distN1;
1351       for ( int j = 0; j < closeFace.GetNP(); ++j)
1352       {
1353         mp[1] = ngMesh.Point( closeFace[j] );
1354         gp_Vec v2p( pV, gp_Pnt( NGPOINT_COORDS( mp[1] )) );
1355         angle2ID.insert( make_pair( v2p1.Angle( v2p ), closeFace[j]));
1356       }
1357       // get node with angle of 60 degrees or greater
1358       map< double, int >::iterator angle_id = angle2ID.lower_bound( 60. * M_PI / 180. );
1359       if ( angle_id == angle2ID.end() ) angle_id = --angle2ID.end();
1360       const double minAngle = 30. * M_PI / 180.;
1361       const double angle = angle_id->first;
1362       bool angleOK = ( angle > minAngle );
1363
1364       // find points to create a triangle
1365       netgen::Element2d tri(3);
1366       tri.SetIndex ( 1 );
1367       tri[0] = vData.ngId;
1368       tri[1] = vData.ngIdCloseN; // to use the closest nodes
1369       tri[2] = angle_id->second; // to use the node with best angle
1370
1371       // decide whether to use the closest node and the node with best angle or to create new ones
1372       for ( int isBestAngleN = 0; isBestAngleN < 2; ++isBestAngleN )
1373       {
1374         bool createNew = !angleOK, distOK = true;
1375         double distFromV;
1376         int triInd = isBestAngleN ? 2 : 1;
1377         mp[isBestAngleN] = ngMesh.Point( tri[triInd] );
1378         if ( isBestAngleN )
1379         {
1380           if ( angleOK )
1381           {
1382             double distN2 = sqrt( dist2( mpV, mp[isBestAngleN]));
1383             createNew = ( fabs( distN2 - distN1 ) > 0.25 * distN1 );
1384           }
1385           else if ( angle < tol )
1386           {
1387             v2p1.SetX( v2p1.X() + 1e-3 );
1388           }
1389           distFromV = distN1;
1390         }
1391         else
1392         {
1393           double segLenHint = ngMesh.GetH( ngMesh.Point( vData.ngId ));
1394           bool avgLenOK  = ( avgSegLen < 0.75 * distN1 );
1395           bool hintLenOK = ( segLenHint  < 0.75 * distN1 );
1396           createNew = (createNew || avgLenOK || hintLenOK );
1397           // we create a new node not closer than 0.5 to the closest face
1398           // in order not to clash with other close face
1399           double r = min( 0.5, ( hintLenOK ? segLenHint : avgSegLen ) / distN1 );
1400           distFromV = r * distN1;
1401         }
1402         if ( createNew )
1403         {
1404           // create a new point, between the node and the vertex if angleOK
1405           gp_XYZ p( NGPOINT_COORDS( mp[isBestAngleN] ));
1406           gp_Vec v2p( pV, p ); v2p.Normalize();
1407           if ( isBestAngleN && !angleOK )
1408             p = p1 + gp_Dir( v2p.XYZ() - v2p1.XYZ()).XYZ() * distN1 * 0.95;
1409           else
1410             p = pV + v2p.XYZ() * distFromV;
1411
1412           if ( !isBestAngleN ) p1 = p, distN1 = distFromV;
1413
1414           mp[isBestAngleN].SetPoint( netgen::Point<3> (p.X(), p.Y(), p.Z()));
1415           ngMesh.AddPoint ( mp[isBestAngleN], 1, netgen::SURFACEPOINT );
1416           tri[triInd] = ngMesh.GetNP();
1417           nodeVec.push_back( helper.AddNode( p.X(), p.Y(), p.Z()) );
1418         }
1419       }
1420       ngMesh.AddSurfaceElement (tri);
1421       swap( tri[1], tri[2] );
1422       ngMesh.AddSurfaceElement (tri);
1423
1424 #ifdef DUMP_TRIANGLES_SCRIPT
1425       py << "n1 = m.AddNode( "<< mpV.X()<<", "<< mpV.Y()<<", "<< mpV.Z()<<") "<< endl
1426          << "n2 = m.AddNode( "<< mp[0].X()<<", "<< mp[0].Y()<<", "<< mp[0].Z()<<") "<< endl
1427          << "n3 = m.AddNode( "<< mp[1].X()<<", "<< mp[1].Y()<<", "<< mp[1].Z()<<" )" << endl
1428          << "m.AddFace([n1,n2,n3])" << endl;
1429 #endif
1430     } // loop on internal vertices of a solid
1431
1432   } // loop on solids with internal vertices
1433 }
1434
1435 //================================================================================
1436 /*!
1437  * \brief Fill SMESH mesh according to contents of netgen mesh
1438  *  \param occgeo - container of OCCT geometry to mesh
1439  *  \param ngMesh - netgen mesh
1440  *  \param initState - bn of entities in netgen mesh before computing
1441  *  \param sMesh - SMESH mesh to fill in
1442  *  \param nodeVec - vector of nodes in which node index == netgen ID
1443  *  \retval int - error
1444  */
1445 //================================================================================
1446
1447 int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
1448                                    const netgen::Mesh&                 ngMesh,
1449                                    const NETGENPlugin_ngMeshInfo&      initState,
1450                                    SMESH_Mesh&                         sMesh,
1451                                    std::vector<const SMDS_MeshNode*>&  nodeVec,
1452                                    SMESH_Comment&                      comment)
1453 {
1454   int nbNod = ngMesh.GetNP();
1455   int nbSeg = ngMesh.GetNSeg();
1456   int nbFac = ngMesh.GetNSE();
1457   int nbVol = ngMesh.GetNE();
1458
1459   SMESHDS_Mesh* meshDS = sMesh.GetMeshDS();
1460
1461   // create and insert nodes into nodeVec
1462   nodeVec.resize( nbNod + 1 );
1463   int i, nbInitNod = initState._nbNodes;
1464   for (i = nbInitNod+1; i <= nbNod; ++i )
1465   {
1466     const netgen::MeshPoint& ngPoint = ngMesh.Point(i);
1467     SMDS_MeshNode* node = NULL;
1468     TopoDS_Vertex aVert;
1469     // First, netgen creates nodes on vertices in occgeo.vmap,
1470     // so node index corresponds to vertex index
1471     // but (issue 0020776) netgen does not create nodes with equal coordinates
1472     if ( i-nbInitNod <= occgeo.vmap.Extent() )
1473     {
1474       gp_Pnt p ( NGPOINT_COORDS(ngPoint) );
1475       for (int iV = i-nbInitNod; aVert.IsNull() && iV <= occgeo.vmap.Extent(); ++iV)
1476       {
1477         aVert = TopoDS::Vertex( occgeo.vmap( iV ) );
1478         gp_Pnt pV = BRep_Tool::Pnt( aVert );
1479         if ( p.SquareDistance( pV ) > 1e-20 )
1480           aVert.Nullify();
1481         else
1482           node = const_cast<SMDS_MeshNode*>( SMESH_Algo::VertexNode( aVert, meshDS ));
1483       }
1484     }
1485     if (!node) // node not found on vertex
1486     {
1487       node = meshDS->AddNode( NGPOINT_COORDS( ngPoint ));
1488       if (!aVert.IsNull())
1489         meshDS->SetNodeOnVertex(node, aVert);
1490     }
1491     nodeVec[i] = node;
1492   }
1493
1494   // create mesh segments along geometric edges
1495   int nbInitSeg = initState._nbSegments;
1496   for (i = nbInitSeg+1; i <= nbSeg; ++i )
1497   {
1498     const netgen::Segment& seg = ngMesh.LineSegment(i);
1499     TopoDS_Edge aEdge;
1500 #ifdef NETGEN_NEW
1501     int pinds[3] = { seg.pnums[0], seg.pnums[1], seg.pnums[2] };
1502 #else
1503     int pinds[3] = { seg.p1, seg.p2, seg.pmid };
1504 #endif
1505     int nbp = 0;
1506     double param2 = 0;
1507     for (int j=0; j < 3; ++j)
1508     {
1509       int pind = pinds[j];
1510       if (pind <= 0 || !nodeVec_ACCESS(pind))
1511         break;
1512       ++nbp;
1513       double param;
1514       if (j < 2)
1515       {
1516         if (aEdge.IsNull())
1517         {
1518           int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
1519           if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
1520             aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
1521         }
1522         param = seg.epgeominfo[j].dist;
1523         param2 += param;
1524       }
1525       else // middle point
1526       {
1527         param = param2 * 0.5;
1528       }
1529       if (!aEdge.IsNull() && nodeVec_ACCESS(pind)->getshapeId() < 1)
1530       {
1531         meshDS->SetNodeOnEdge(nodeVec_ACCESS(pind), aEdge, param);
1532       }
1533     }
1534     if ( nbp > 1 )
1535     {
1536       SMDS_MeshEdge* edge = 0;
1537       if (nbp == 2) // second order ?
1538       {
1539         if ( meshDS->FindEdge( nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1])))
1540           continue;
1541         edge = meshDS->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]));
1542       }
1543       else
1544       {
1545         if ( meshDS->FindEdge( nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]),
1546                                nodeVec_ACCESS(pinds[2])))
1547           continue;
1548         edge = meshDS->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]),
1549                                nodeVec_ACCESS(pinds[2]));
1550       }
1551       if (!edge)
1552       {
1553         if ( comment.empty() ) comment << "Cannot create a mesh edge";
1554         MESSAGE("Cannot create a mesh edge");
1555         nbSeg = nbFac = nbVol = 0;
1556         break;
1557       }
1558       if ( !aEdge.IsNull() && edge->getshapeId() < 1 )
1559         meshDS->SetMeshElementOnShape(edge, aEdge);
1560     }
1561     else if ( comment.empty() )
1562     {
1563       comment << "Invalid netgen segment #" << i;
1564     }
1565   }
1566
1567   // create mesh faces along geometric faces
1568   int nbInitFac = initState._nbFaces;
1569   for (i = nbInitFac+1; i <= nbFac; ++i )
1570   {
1571     const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
1572     int aGeomFaceInd = elem.GetIndex();
1573     TopoDS_Face aFace;
1574     if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
1575       aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
1576     vector<SMDS_MeshNode*> nodes;
1577     for (int j=1; j <= elem.GetNP(); ++j)
1578     {
1579       int pind = elem.PNum(j);
1580       if ( pind < 1 || pind >= nodeVec.size() )
1581         break;
1582       if ( SMDS_MeshNode* node = nodeVec_ACCESS(pind))
1583       {
1584         nodes.push_back(node);
1585         if (!aFace.IsNull() && node->getshapeId() < 1)
1586         {
1587           const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
1588           meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
1589         }
1590       }
1591     }
1592     if ( nodes.size() != elem.GetNP() )
1593     {
1594       if ( comment.empty() )
1595         comment << "Invalid netgen 2d element #" << i;
1596       continue; // bad node ids
1597     }
1598     SMDS_MeshFace* face = NULL;
1599     switch (elem.GetType())
1600     {
1601     case netgen::TRIG:
1602       face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
1603       break;
1604     case netgen::QUAD:
1605       face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
1606       break;
1607     case netgen::TRIG6:
1608       face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
1609       break;
1610     case netgen::QUAD8:
1611       face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
1612                              nodes[4],nodes[7],nodes[5],nodes[6]);
1613       break;
1614     default:
1615       MESSAGE("NETGEN created a face of unexpected type, ignoring");
1616       continue;
1617     }
1618     if (!face)
1619     {
1620       if ( comment.empty() ) comment << "Cannot create a mesh face";
1621       MESSAGE("Cannot create a mesh face");
1622       nbSeg = nbFac = nbVol = 0;
1623       break;
1624     }
1625     if (!aFace.IsNull())
1626       meshDS->SetMeshElementOnShape(face, aFace);
1627   }
1628
1629   // create tetrahedra
1630   for (i = 1; i <= nbVol; ++i)
1631   {
1632     const netgen::Element& elem = ngMesh.VolumeElement(i);
1633     int aSolidInd = elem.GetIndex();
1634     TopoDS_Solid aSolid;
1635     if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
1636       aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
1637     vector<SMDS_MeshNode*> nodes;
1638     for (int j=1; j <= elem.GetNP(); ++j)
1639     {
1640       int pind = elem.PNum(j);
1641       if ( pind < 1 || pind >= nodeVec.size() )
1642         break;
1643       if ( SMDS_MeshNode* node = nodeVec_ACCESS(pind) )
1644       {
1645         nodes.push_back(node);
1646         if ( !aSolid.IsNull() && node->getshapeId() < 1 )
1647           meshDS->SetNodeInVolume(node, aSolid);
1648       }
1649     }
1650     if ( nodes.size() != elem.GetNP() )
1651     {
1652       if ( comment.empty() )
1653         comment << "Invalid netgen 3d element #" << i;
1654       continue;
1655     }
1656     SMDS_MeshVolume* vol = NULL;
1657     switch (elem.GetType())
1658     {
1659     case netgen::TET:
1660       vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
1661       break;
1662     case netgen::TET10:
1663       vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
1664                               nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
1665       break;
1666     default:
1667       MESSAGE("NETGEN created a volume of unexpected type, ignoring");
1668       continue;
1669     }
1670     if (!vol)
1671     {
1672       if ( comment.empty() ) comment << "Cannot create a mesh volume";
1673       MESSAGE("Cannot create a mesh volume");
1674       nbSeg = nbFac = nbVol = 0;
1675       break;
1676     }
1677     if (!aSolid.IsNull())
1678       meshDS->SetMeshElementOnShape(vol, aSolid);
1679   }
1680   return comment.empty() ? 0 : 1;
1681 }
1682
1683 namespace
1684 {
1685   //================================================================================
1686   /*!
1687    * \brief Restrict size of elements on the given edge
1688    */
1689   //================================================================================
1690
1691   void setLocalSize(const TopoDS_Edge& edge,
1692                     double             size,
1693                     netgen::Mesh&      mesh)
1694   {
1695     const int nb = 1000;
1696     Standard_Real u1, u2;
1697     Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, u1, u2);
1698     if ( curve.IsNull() )
1699     {
1700       TopoDS_Iterator vIt( edge );
1701       if ( !vIt.More() ) return;
1702       gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( vIt.Value() ));
1703       NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size );
1704     }
1705     else
1706     {
1707       Standard_Real delta = (u2-u1)/nb;
1708       for(int i=0; i<nb; i++)
1709       {
1710         Standard_Real u = u1 + delta*i;
1711         gp_Pnt p = curve->Value(u);
1712         NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), size );
1713         netgen::Point3d pi(p.X(), p.Y(), p.Z());
1714         double resultSize = mesh.GetH(pi);
1715         if ( resultSize - size > 0.1*size )
1716           // netgen does restriction iff oldH/newH > 1.2 (localh.cpp:136)
1717           NETGENPlugin_Mesher::RestrictLocalSize( mesh, p.XYZ(), resultSize/1.201 );
1718       }
1719     }
1720   }
1721
1722   //================================================================================
1723   /*!
1724    * \brief Convert error into text
1725    */
1726   //================================================================================
1727
1728   std::string text(int err)
1729   {
1730     if ( !err )
1731       return string("");
1732     return
1733       SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task;
1734   }
1735
1736   //================================================================================
1737   /*!
1738    * \brief Convert exception into text
1739    */
1740   //================================================================================
1741
1742   std::string text(Standard_Failure& ex)
1743   {
1744     SMESH_Comment str("Exception in netgen::OCCGenerateMesh()");
1745     str << " at " << netgen::multithread.task
1746         << ": " << ex.DynamicType()->Name();
1747     if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
1748       str << ": " << ex.GetMessageString();
1749     return str;
1750   }
1751   //================================================================================
1752   /*!
1753    * \brief Convert exception into text
1754    */
1755   //================================================================================
1756
1757   std::string text(netgen::NgException& ex)
1758   {
1759     SMESH_Comment str("NgException");
1760     str << " at " << netgen::multithread.task << ": " << ex.What();
1761     return str;
1762   }
1763 }
1764
1765 //=============================================================================
1766 /*!
1767  * Here we are going to use the NETGEN mesher
1768  */
1769 //=============================================================================
1770
1771 bool NETGENPlugin_Mesher::Compute()
1772 {
1773   NETGENPlugin_NetgenLibWrapper ngLib;
1774
1775   netgen::MeshingParameters& mparams = netgen::mparam;
1776   MESSAGE("Compute with:\n"
1777           " max size = " << mparams.maxh << "\n"
1778           " segments per edge = " << mparams.segmentsperedge);
1779   MESSAGE("\n"
1780           " growth rate = " << mparams.grading << "\n"
1781           " elements per radius = " << mparams.curvaturesafety << "\n"
1782           " second order = " << mparams.secondorder << "\n"
1783           " quad allowed = " << mparams.quad);
1784
1785   SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
1786
1787
1788   // -------------------------
1789   // Prepare OCC geometry
1790   // -------------------------
1791
1792   netgen::OCCGeometry occgeo;
1793   list< SMESH_subMesh* > meshedSM[3]; // for 0-2 dimensions
1794   NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume );
1795   PrepareOCCgeometry( occgeo, _shape, *_mesh, meshedSM, &internals );
1796
1797   // -------------------------
1798   // Generate the mesh
1799   // -------------------------
1800
1801   netgen::Mesh *ngMesh = NULL;
1802   NETGENPlugin_ngMeshInfo initState; // it remembers size of ng mesh equal to size of Smesh
1803
1804   SMESH_Comment comment;
1805   int err = 0;
1806
1807   // vector of nodes in which node index == netgen ID
1808   vector< const SMDS_MeshNode* > nodeVec;
1809
1810   {
1811     // ----------------
1812     // compute 1D mesh
1813     // ----------------
1814     if ( _simpleHyp )
1815     {
1816       // not to RestrictLocalH() according to curvature during MESHCONST_ANALYSE
1817       mparams.uselocalh = false;
1818       mparams.grading = 0.8; // not limitited size growth
1819
1820       if ( _simpleHyp->GetNumberOfSegments() )
1821         // nb of segments
1822         mparams.maxh = occgeo.boundingbox.Diam();
1823       else
1824         // segment length
1825         mparams.maxh = _simpleHyp->GetLocalLength();
1826     }
1827
1828     if ( mparams.maxh == 0.0 )
1829       mparams.maxh = occgeo.boundingbox.Diam();
1830     if ( _simpleHyp || ( mparams.minh == 0.0 && _fineness != NETGENPlugin_Hypothesis::UserDefined))
1831       mparams.minh = GetDefaultMinSize( _shape, mparams.maxh );
1832
1833 #ifdef NETGEN_NEW
1834     // Local size on faces
1835     occgeo.face_maxh = mparams.maxh;
1836     for(map<int,double>::const_iterator it=FaceId2LocalSize.begin();
1837         it!=FaceId2LocalSize.end(); it++)
1838     {
1839       int key = (*it).first;
1840       double val = (*it).second;
1841       const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
1842       int faceNgID = occgeo.fmap.FindIndex(shape);
1843       occgeo.SetFaceMaxH(faceNgID, val);
1844     }
1845 #endif
1846
1847     // Let netgen create ngMesh and calculate element size on not meshed shapes
1848     char *optstr = 0;
1849     int startWith = netgen::MESHCONST_ANALYSE;
1850     int endWith   = netgen::MESHCONST_ANALYSE;
1851     try
1852     {
1853       OCC_CATCH_SIGNALS;
1854       err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
1855 #ifdef WITH_SMESH_CANCEL_COMPUTE
1856       if(netgen::multithread.terminate)
1857         return false;
1858 #endif
1859       comment << text(err);
1860     }
1861     catch (Standard_Failure& ex)
1862     {
1863       comment << text(ex);
1864     }
1865     err = 0; //- MESHCONST_ANALYSE isn't so important step
1866     if ( !ngMesh )
1867       return false;
1868     ngLib.setMesh(( Ng_Mesh*) ngMesh );
1869
1870     if ( _simpleHyp )
1871     {
1872       // Pass 1D simple parameters to NETGEN
1873       // --------------------------------
1874       int      nbSeg = _simpleHyp->GetNumberOfSegments();
1875       double segSize = _simpleHyp->GetLocalLength();
1876       for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE )
1877       {
1878         const TopoDS_Edge& e = TopoDS::Edge( occgeo.emap(iE));
1879         if ( nbSeg )
1880           segSize = SMESH_Algo::EdgeLength( e ) / ( nbSeg - 0.4 );
1881         setLocalSize( e, segSize, *ngMesh );
1882       }
1883     }
1884     else // if ( ! _simpleHyp )
1885     {
1886       // Local size on vertices and edges
1887       // --------------------------------
1888       for(std::map<int,double>::const_iterator it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++)
1889       {
1890         int key = (*it).first;
1891         double hi = (*it).second;
1892         const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
1893         const TopoDS_Edge& e = TopoDS::Edge(shape);
1894         setLocalSize( e, hi, *ngMesh );
1895       }
1896       for(std::map<int,double>::const_iterator it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++)
1897       {
1898         int key = (*it).first;
1899         double hi = (*it).second;
1900         const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
1901         const TopoDS_Vertex& v = TopoDS::Vertex(shape);
1902         gp_Pnt p = BRep_Tool::Pnt(v);
1903         NETGENPlugin_Mesher::RestrictLocalSize( *ngMesh, p.XYZ(), hi );
1904       }
1905     }
1906
1907     // Precompute internal edges (issue 0020676) in order to
1908     // add mesh on them correctly (twice) to netgen mesh
1909     if ( !err && internals.hasInternalEdges() )
1910     {
1911       // load internal shapes into OCCGeometry
1912       netgen::OCCGeometry intOccgeo;
1913       internals.getInternalEdges( intOccgeo.fmap, intOccgeo.emap, intOccgeo.vmap, meshedSM );
1914       intOccgeo.boundingbox = occgeo.boundingbox;
1915       intOccgeo.shape = occgeo.shape;
1916 #ifdef NETGEN_NEW
1917       intOccgeo.face_maxh.SetSize(intOccgeo.fmap.Extent());
1918       intOccgeo.face_maxh = netgen::mparam.maxh;
1919 #endif
1920
1921       // let netgen compute element size by the main geometry in temporary mesh
1922       netgen::Mesh *tmpNgMesh = NULL;
1923       try
1924       {
1925         OCC_CATCH_SIGNALS;
1926         netgen::OCCGenerateMesh(occgeo, tmpNgMesh, startWith, endWith, optstr);
1927 #ifdef WITH_SMESH_CANCEL_COMPUTE
1928         if(netgen::multithread.terminate)
1929           return false;
1930 #endif
1931         // compute mesh on internal edges
1932         endWith = netgen::MESHCONST_MESHEDGES;
1933         err = netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, startWith, endWith, optstr);
1934         comment << text(err);
1935       }
1936       catch (Standard_Failure& ex)
1937       {
1938         comment << text(ex);
1939         err = 1;
1940       }
1941       // fill SMESH by netgen mesh
1942       vector< const SMDS_MeshNode* > tmpNodeVec;
1943       FillSMesh( intOccgeo, *tmpNgMesh, initState, *_mesh, tmpNodeVec, comment );
1944       err = ( err || !comment.empty() );
1945
1946       nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)tmpNgMesh);
1947     }
1948
1949     // Fill ngMesh with nodes and segments of computed submeshes
1950     if ( !err )
1951     {
1952       _faceDescriptors.clear();
1953       err = ! ( fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_0D ]) &&
1954                 fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_1D ]));
1955     }
1956     initState = NETGENPlugin_ngMeshInfo(ngMesh);
1957
1958     // Compute 1d mesh
1959     if (!err)
1960     {
1961       startWith = endWith = netgen::MESHCONST_MESHEDGES;
1962       try
1963       {
1964         OCC_CATCH_SIGNALS;
1965         err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
1966 #ifdef WITH_SMESH_CANCEL_COMPUTE
1967         if(netgen::multithread.terminate)
1968           return false;
1969 #endif
1970         comment << text(err);
1971       }
1972       catch (Standard_Failure& ex)
1973       {
1974         comment << text(ex);
1975         err = 1;
1976       }
1977     }
1978     mparams.uselocalh = true; // restore as it is used at surface optimization
1979
1980     // ---------------------
1981     // compute surface mesh
1982     // ---------------------
1983     if (!err)
1984     {
1985       // Pass 2D simple parameters to NETGEN
1986       if ( _simpleHyp ) {
1987         if ( double area = _simpleHyp->GetMaxElementArea() ) {
1988           // face area
1989           mparams.maxh = sqrt(2. * area/sqrt(3.0));
1990           mparams.grading = 0.4; // moderate size growth
1991         }
1992         else {
1993           // length from edges
1994           if ( ngMesh->GetNSeg() ) {
1995             double edgeLength = 0;
1996             TopTools_MapOfShape visitedEdges;
1997             for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() )
1998               if( visitedEdges.Add(exp.Current()) )
1999                 edgeLength += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() ));
2000             // we have to multiply length by 2 since for each TopoDS_Edge there
2001             // are double set of NETGEN edges, in other words, we have to
2002             // divide ngMesh->GetNSeg() by 2.
2003             mparams.maxh = 2*edgeLength / ngMesh->GetNSeg();
2004           }
2005           else {
2006             mparams.maxh = 1000;
2007           }
2008           mparams.grading = 0.2; // slow size growth
2009         }
2010         mparams.quad = _simpleHyp->GetAllowQuadrangles();
2011         mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
2012         ngMesh->SetGlobalH (mparams.maxh);
2013         netgen::Box<3> bb = occgeo.GetBoundingBox();
2014         bb.Increase (bb.Diam()/20);
2015         ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading);
2016       }
2017
2018       // Care of vertices internal in faces (issue 0020676)
2019       if ( internals.hasInternalVertexInFace() )
2020       {
2021         // store computed segments in SMESH in order not to create SMESH
2022         // edges for ng segments added by addIntVerticesInFaces()
2023         FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
2024         // add segments to faces with internal vertices
2025         addIntVerticesInFaces( occgeo, *ngMesh, nodeVec, internals );
2026         initState = NETGENPlugin_ngMeshInfo(ngMesh);
2027       }
2028
2029       // Let netgen compute 2D mesh
2030       startWith = netgen::MESHCONST_MESHSURFACE;
2031       endWith = _optimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
2032       try
2033       {
2034         OCC_CATCH_SIGNALS;
2035         err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
2036 #ifdef WITH_SMESH_CANCEL_COMPUTE
2037         if(netgen::multithread.terminate)
2038           return false;
2039 #endif
2040         comment << text (err);
2041       }
2042       catch (Standard_Failure& ex)
2043       {
2044         comment << text(ex);
2045         //err = 1; -- try to make volumes anyway
2046       }
2047       catch (netgen::NgException exc)
2048       {
2049         comment << text(exc);
2050         //err = 1; -- try to make volumes anyway
2051       }
2052     }
2053     // ---------------------
2054     // generate volume mesh
2055     // ---------------------
2056       // Fill ngMesh with nodes and faces of computed 2D submeshes
2057     if ( !err && _isVolume && !meshedSM[ MeshDim_2D ].empty() )
2058     {
2059       // load SMESH with computed segments and faces
2060       FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
2061       // fill ng mesh
2062       err = ! ( fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_2D ]));
2063       initState = NETGENPlugin_ngMeshInfo(ngMesh);
2064     }
2065     if (!err && _isVolume)
2066     {
2067       // Pass 3D simple parameters to NETGEN
2068       const NETGENPlugin_SimpleHypothesis_3D* simple3d =
2069         dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
2070       if ( simple3d ) {
2071         if ( double vol = simple3d->GetMaxElementVolume() ) {
2072           // max volume
2073           mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
2074           mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
2075         }
2076         else {
2077           // length from faces
2078           mparams.maxh = ngMesh->AverageH();
2079         }
2080         ngMesh->SetGlobalH (mparams.maxh);
2081         mparams.grading = 0.4;
2082         ngMesh->CalcLocalH();
2083       }
2084       // Care of vertices internal in solids and internal faces (issue 0020676)
2085       if ( internals.hasInternalVertexInSolid() || internals.hasInternalFaces() )
2086       {
2087         // store computed faces in SMESH in order not to create SMESH
2088         // faces for ng faces added here
2089         FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
2090         // add ng faces to solids with internal vertices
2091         addIntVerticesInSolids( occgeo, *ngMesh, nodeVec, internals );
2092         // duplicate mesh faces on internal faces
2093         fixIntFaces( occgeo, *ngMesh, internals );
2094         initState = NETGENPlugin_ngMeshInfo(ngMesh);
2095       }
2096       // Let netgen compute 3D mesh
2097       startWith = endWith = netgen::MESHCONST_MESHVOLUME;
2098       try
2099       {
2100         OCC_CATCH_SIGNALS;
2101         err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
2102 #ifdef WITH_SMESH_CANCEL_COMPUTE
2103         if(netgen::multithread.terminate)
2104           return false;
2105 #endif
2106         if ( comment.empty() ) // do not overwrite a previos error
2107           comment << text(err);
2108       }
2109       catch (Standard_Failure& ex)
2110       {
2111         if ( comment.empty() ) // do not overwrite a previos error
2112           comment << text(ex);
2113         err = 1;
2114       }
2115       catch (netgen::NgException exc)
2116       {
2117         if ( comment.empty() ) // do not overwrite a previos error
2118           comment << text(exc);
2119         err = 1;
2120       }
2121       // Let netgen optimize 3D mesh
2122       if ( !err && _optimize )
2123       {
2124         startWith = endWith = netgen::MESHCONST_OPTVOLUME;
2125         try
2126         {
2127           OCC_CATCH_SIGNALS;
2128           err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
2129 #ifdef WITH_SMESH_CANCEL_COMPUTE
2130           if(netgen::multithread.terminate)
2131             return false;
2132 #endif
2133           if ( comment.empty() ) // do not overwrite a previos error
2134             comment << text(err);
2135         }
2136         catch (Standard_Failure& ex)
2137         {
2138           if ( comment.empty() ) // do not overwrite a previos error
2139             comment << text(ex);
2140         }
2141         catch (netgen::NgException exc)
2142         {
2143           if ( comment.empty() ) // do not overwrite a previos error
2144             comment << text(exc);
2145         }
2146       }
2147     }
2148     if (!err && mparams.secondorder > 0)
2149     {
2150       try
2151       {
2152         OCC_CATCH_SIGNALS;
2153         netgen::OCCRefinementSurfaces ref (occgeo);
2154         ref.MakeSecondOrder (*ngMesh);
2155       }
2156       catch (Standard_Failure& ex)
2157       {
2158         if ( comment.empty() ) // do not overwrite a previos error
2159           comment << "Exception in netgen at passing to 2nd order ";
2160       }
2161       catch (netgen::NgException exc)
2162       {
2163         if ( comment.empty() ) // do not overwrite a previos error
2164           comment << exc.What();
2165       }
2166     }
2167   }
2168   int nbNod = ngMesh->GetNP();
2169   int nbSeg = ngMesh->GetNSeg();
2170   int nbFac = ngMesh->GetNSE();
2171   int nbVol = ngMesh->GetNE();
2172   bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
2173
2174   MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
2175           ", nb nodes: "    << nbNod <<
2176           ", nb segments: " << nbSeg <<
2177           ", nb faces: "    << nbFac <<
2178           ", nb volumes: "  << nbVol);
2179
2180   // Feed back the SMESHDS with the generated Nodes and Elements
2181   if ( true /*isOK*/ ) // get whatever built
2182     FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment ); //!<
2183
2184   SMESH_ComputeErrorPtr readErr = readErrors(nodeVec);
2185   if ( readErr && !readErr->myBadElements.empty() )
2186     error = readErr;
2187
2188   if ( error->IsOK() && ( !isOK || comment.size() > 0 ))
2189     error->myName = COMPERR_ALGO_FAILED;
2190   if ( !comment.empty() )
2191     error->myComment = comment;
2192
2193   // SetIsAlwaysComputed( true ) to empty sub-meshes, which
2194   // appear if the geometry contains coincident sub-shape due
2195   // to bool merge_solids = 1; in netgen/libsrc/occ/occgenmesh.cpp
2196   const int nbMaps = 2;
2197   const TopTools_IndexedMapOfShape* geoMaps[nbMaps] =
2198     { & occgeo.vmap, & occgeo.emap/*, & occgeo.fmap*/ };
2199   for ( int iMap = 0; iMap < nbMaps; ++iMap )
2200     for (int i = 1; i <= geoMaps[iMap]->Extent(); i++)
2201       if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( geoMaps[iMap]->FindKey(i)))
2202         if ( !sm->IsMeshComputed() )
2203           sm->SetIsAlwaysComputed( true );
2204
2205   // set bad compute error to subshapes of all failed sub-shapes
2206   if ( !error->IsOK() )
2207   {
2208     bool pb2D = false, pb3D = false;
2209     for (int i = 1; i <= occgeo.fmap.Extent(); i++) {
2210       int status = occgeo.facemeshstatus[i-1];
2211       if (status == 1 ) continue;
2212       if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.fmap( i ))) {
2213         SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2214         if ( !smError || smError->IsOK() ) {
2215           if ( status == -1 )
2216             smError.reset( new SMESH_ComputeError( *error ));
2217           else
2218             smError.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, "Ignored" ));
2219           if ( SMESH_Algo::GetMeshError( sm ) == SMESH_Algo::MEr_OK )
2220             smError->myName = COMPERR_WARNING;
2221         }
2222         pb2D = pb2D || smError->IsKO();
2223       }
2224     }
2225     if ( !pb2D ) // all faces are OK
2226       for (int i = 1; i <= occgeo.somap.Extent(); i++)
2227         if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.somap( i )))
2228         {
2229           bool smComputed = !sm->IsEmpty();
2230           if ( smComputed && internals.hasInternalVertexInSolid( sm->GetId() ))
2231           {
2232             int nbIntV = internals.getSolidsWithVertices().find( sm->GetId() )->second.size();
2233             SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
2234             smComputed = ( smDS->NbElements() > 0 || smDS->NbNodes() > nbIntV );
2235           }
2236           SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
2237           if ( !smComputed && ( !smError || smError->IsOK() ))
2238           {
2239             smError.reset( new SMESH_ComputeError( *error ));
2240             if ( SMESH_Algo::GetMeshError( sm ) == SMESH_Algo::MEr_OK )
2241               smError->myName = COMPERR_WARNING;
2242           }
2243           pb3D = pb3D || ( smError && smError->IsKO() );
2244         }
2245     if ( !pb2D && !pb3D )
2246       err = 0; // no fatal errors, only warnings
2247   }
2248
2249   return !err;
2250 }
2251
2252 //=============================================================================
2253 /*!
2254  * Evaluate
2255  */
2256 //=============================================================================
2257 bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
2258 {
2259   netgen::MeshingParameters& mparams = netgen::mparam;
2260
2261
2262   // -------------------------
2263   // Prepare OCC geometry
2264   // -------------------------
2265   netgen::OCCGeometry occgeo;
2266   PrepareOCCgeometry( occgeo, _shape, *_mesh );
2267
2268   bool tooManyElems = false;
2269   const int hugeNb = std::numeric_limits<int>::max() / 100;
2270
2271   // ----------------
2272   // evaluate 1D
2273   // ----------------
2274   // pass 1D simple parameters to NETGEN
2275   if ( _simpleHyp ) {
2276     if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
2277       // nb of segments
2278       mparams.segmentsperedge = nbSeg + 0.1;
2279       mparams.maxh = occgeo.boundingbox.Diam();
2280       mparams.minh = GetDefaultMinSize( _shape, mparams.maxh );
2281       mparams.grading = 0.01;
2282     }
2283     else {
2284       // segment length
2285       mparams.segmentsperedge = 1;
2286       mparams.maxh = _simpleHyp->GetLocalLength();
2287     }
2288   }
2289   // let netgen create ngMesh and calculate element size on not meshed shapes
2290   NETGENPlugin_NetgenLibWrapper ngLib;
2291   netgen::Mesh *ngMesh = NULL;
2292   char *optstr = 0;
2293   int startWith = netgen::MESHCONST_ANALYSE;
2294   int endWith   = netgen::MESHCONST_MESHEDGES;
2295   int err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
2296 #ifdef WITH_SMESH_CANCEL_COMPUTE
2297   if(netgen::multithread.terminate)
2298     return false;
2299 #endif
2300   ngLib.setMesh(( Ng_Mesh*) ngMesh );
2301   if (err) {
2302     if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( _shape ))
2303       sm->GetComputeError().reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED ));
2304     return false;
2305   }
2306
2307   // calculate total nb of segments and length of edges
2308   double fullLen = 0.0;
2309   int fullNbSeg = 0;
2310   int entity = mparams.secondorder > 0 ? SMDSEntity_Quad_Edge : SMDSEntity_Edge;
2311   TopTools_DataMapOfShapeInteger Edge2NbSeg;
2312   for (TopExp_Explorer exp(_shape, TopAbs_EDGE); exp.More(); exp.Next())
2313   {
2314     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
2315     if( !Edge2NbSeg.Bind(E,0) )
2316       continue;
2317
2318     double aLen = SMESH_Algo::EdgeLength(E);
2319     fullLen += aLen;
2320
2321     vector<int>& aVec = aResMap[_mesh->GetSubMesh(E)];
2322     if ( aVec.empty() )
2323       aVec.resize( SMDSEntity_Last, 0);
2324     else
2325       fullNbSeg += aVec[ entity ];
2326   }
2327
2328   // store nb of segments computed by Netgen
2329   NCollection_Map<Link> linkMap;
2330   for (int i = 1; i <= ngMesh->GetNSeg(); ++i )
2331   {
2332     const netgen::Segment& seg = ngMesh->LineSegment(i);
2333     Link link(seg[0], seg[1]);
2334     if ( !linkMap.Add( link )) continue;
2335     int aGeomEdgeInd = seg.epgeominfo[0].edgenr;
2336     if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
2337     {
2338       vector<int>& aVec = aResMap[_mesh->GetSubMesh(occgeo.emap(aGeomEdgeInd))];
2339       aVec[ entity ]++;
2340     }
2341   }
2342   // store nb of nodes on edges computed by Netgen
2343   TopTools_DataMapIteratorOfDataMapOfShapeInteger Edge2NbSegIt(Edge2NbSeg);
2344   for (; Edge2NbSegIt.More(); Edge2NbSegIt.Next())
2345   {
2346     vector<int>& aVec = aResMap[_mesh->GetSubMesh(Edge2NbSegIt.Key())];
2347     if ( aVec[ entity ] > 1 && aVec[ SMDSEntity_Node ] == 0 )
2348       aVec[SMDSEntity_Node] = mparams.secondorder > 0  ? 2*aVec[ entity ]-1 : aVec[ entity ]-1;
2349
2350     fullNbSeg += aVec[ entity ];
2351     Edge2NbSeg( Edge2NbSegIt.Key() ) = aVec[ entity ];
2352   }
2353
2354   // ----------------
2355   // evaluate 2D
2356   // ----------------
2357   if ( _simpleHyp ) {
2358     if ( double area = _simpleHyp->GetMaxElementArea() ) {
2359       // face area
2360       mparams.maxh = sqrt(2. * area/sqrt(3.0));
2361       mparams.grading = 0.4; // moderate size growth
2362     }
2363     else {
2364       // length from edges
2365       mparams.maxh = fullLen/fullNbSeg;
2366       mparams.grading = 0.2; // slow size growth
2367     }
2368   }
2369   mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
2370   mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
2371
2372   for (TopExp_Explorer exp(_shape, TopAbs_FACE); exp.More(); exp.Next())
2373   {
2374     TopoDS_Face F = TopoDS::Face( exp.Current() );
2375     SMESH_subMesh *sm = _mesh->GetSubMesh(F);
2376     GProp_GProps G;
2377     BRepGProp::SurfaceProperties(F,G);
2378     double anArea = G.Mass();
2379     tooManyElems = tooManyElems || ( anArea/hugeNb > mparams.maxh*mparams.maxh );
2380     int nb1d = 0;
2381     if ( !tooManyElems )
2382     {
2383       TopTools_MapOfShape egdes;
2384       for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next())
2385         if ( egdes.Add( exp1.Current() ))
2386           nb1d += Edge2NbSeg.Find(exp1.Current());
2387     }
2388     int nbFaces = tooManyElems ? hugeNb : int( 4*anArea / (mparams.maxh*mparams.maxh*sqrt(3.)));
2389     int nbNodes = tooManyElems ? hugeNb : (( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
2390
2391     vector<int> aVec(SMDSEntity_Last, 0);
2392     if( mparams.secondorder > 0 ) {
2393       int nb1d_in = (nbFaces*3 - nb1d) / 2;
2394       aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
2395       aVec[SMDSEntity_Quad_Triangle] = nbFaces;
2396     }
2397     else {
2398       aVec[SMDSEntity_Node] = Max ( nbNodes, 0  );
2399       aVec[SMDSEntity_Triangle] = nbFaces;
2400     }
2401     aResMap[sm].swap(aVec);
2402   }
2403
2404   // ----------------
2405   // evaluate 3D
2406   // ----------------
2407   if(_isVolume) {
2408     // pass 3D simple parameters to NETGEN
2409     const NETGENPlugin_SimpleHypothesis_3D* simple3d =
2410       dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
2411     if ( simple3d ) {
2412       if ( double vol = simple3d->GetMaxElementVolume() ) {
2413         // max volume
2414         mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
2415         mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
2416       }
2417       else {
2418         // using previous length from faces
2419       }
2420       mparams.grading = 0.4;
2421       mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
2422     }
2423     GProp_GProps G;
2424     BRepGProp::VolumeProperties(_shape,G);
2425     double aVolume = G.Mass();
2426     double tetrVol = 0.1179*mparams.maxh*mparams.maxh*mparams.maxh;
2427     tooManyElems = tooManyElems || ( aVolume/hugeNb > tetrVol );
2428     int nbVols = tooManyElems ? hugeNb : int(aVolume/tetrVol);
2429     int nb1d_in = int(( nbVols*6 - fullNbSeg ) / 6 );
2430     vector<int> aVec(SMDSEntity_Last, 0 );
2431     if ( tooManyElems ) // avoid FPE
2432     {
2433       aVec[SMDSEntity_Node] = hugeNb;
2434       aVec[ mparams.secondorder > 0 ? SMDSEntity_Quad_Tetra : SMDSEntity_Tetra] = hugeNb;
2435     }
2436     else
2437     {
2438       if( mparams.secondorder > 0 ) {
2439         aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in;
2440         aVec[SMDSEntity_Quad_Tetra] = nbVols;
2441       }
2442       else {
2443         aVec[SMDSEntity_Node] = nb1d_in/3 + 1;
2444         aVec[SMDSEntity_Tetra] = nbVols;
2445       }
2446     }
2447     SMESH_subMesh *sm = _mesh->GetSubMesh(_shape);
2448     aResMap[sm].swap(aVec);
2449   }
2450
2451   return true;
2452 }
2453
2454 //================================================================================
2455 /*!
2456  * \brief Remove "test.out" and "problemfaces" files in current directory
2457  */
2458 //================================================================================
2459
2460 void NETGENPlugin_Mesher::RemoveTmpFiles()
2461 {
2462   if ( SMESH_File("test.out").remove() && netgen::testout)
2463   {
2464     delete netgen::testout;
2465     netgen::testout = 0;
2466   }
2467   SMESH_File("problemfaces").remove();
2468   SMESH_File("occmesh.rep").remove();
2469 }
2470
2471 //================================================================================
2472 /*!
2473  * \brief Read mesh entities preventing successful computation from "test.out" file
2474  */
2475 //================================================================================
2476
2477 SMESH_ComputeErrorPtr
2478 NETGENPlugin_Mesher::readErrors(const vector<const SMDS_MeshNode* >& nodeVec)
2479 {
2480   SMESH_ComputeErrorPtr err = SMESH_ComputeError::New
2481     (COMPERR_BAD_INPUT_MESH, "Some edges multiple times in surface mesh");
2482   SMESH_File file("test.out");
2483   vector<int> two(2);
2484   const char* badEdgeStr = " multiple times in surface mesh";
2485   const int   badEdgeStrLen = strlen( badEdgeStr );
2486   while( !file.eof() )
2487   {
2488     if ( strncmp( file, "Edge ", 5 ) == 0 &&
2489          file.getInts( two ) &&
2490          strncmp( file, badEdgeStr, badEdgeStrLen ) == 0 &&
2491          two[0] < nodeVec.size() && two[1] < nodeVec.size())
2492     {
2493       err->myBadElements.push_back( new SMDS_LinearEdge( nodeVec[ two[0]], nodeVec[ two[1]] ));
2494       file += badEdgeStrLen;
2495     }
2496     else if ( strncmp( file, "Intersecting: ", 14 ) == 0 )
2497     {
2498 // Intersecting:
2499 // openelement 18 with open element 126
2500 // 41  36  38
2501 // 69  70  72
2502       vector<int> three1(3), three2(3);
2503       file.getLine();
2504       const char* pos = file;
2505       bool ok = ( strncmp( file, "openelement ", 12 ) == 0 );
2506       ok = ok && file.getInts( two );
2507       ok = ok && file.getInts( three1 );
2508       ok = ok && file.getInts( three2 );
2509       for ( int i = 0; ok && i < 3; ++i )
2510         ok = ( three1[i] < nodeVec.size() && nodeVec[ three1[i]]);
2511       for ( int i = 0; ok && i < 3; ++i )
2512         ok = ( three2[i] < nodeVec.size() && nodeVec[ three2[i]]);
2513       if ( ok )
2514       {
2515         err->myBadElements.push_back( new SMDS_FaceOfNodes( nodeVec[ three1[0]],
2516                                                             nodeVec[ three1[1]],
2517                                                             nodeVec[ three1[2]]));
2518         err->myBadElements.push_back( new SMDS_FaceOfNodes( nodeVec[ three2[0]],
2519                                                             nodeVec[ three2[1]],
2520                                                             nodeVec[ three2[2]]));
2521         err->myComment = "Intersecting triangles";
2522       }
2523       else
2524       {
2525         file.setPos( pos );
2526       }
2527     }
2528     else
2529     {
2530       ++file;
2531     }
2532   }
2533   return err;
2534 }
2535
2536 //================================================================================
2537 /*!
2538  * \brief Constructor of NETGENPlugin_ngMeshInfo
2539  */
2540 //================================================================================
2541
2542 NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh)
2543 {
2544   if ( ngMesh )
2545   {
2546     _nbNodes    = ngMesh->GetNP();
2547     _nbSegments = ngMesh->GetNSeg();
2548     _nbFaces    = ngMesh->GetNSE();
2549     _nbVolumes  = ngMesh->GetNE();
2550   }
2551   else
2552   {
2553     _nbNodes = _nbSegments = _nbFaces = _nbVolumes = 0;
2554   }
2555 }
2556
2557 //================================================================================
2558 /*!
2559  * \brief Find "internal" sub-shapes
2560  */
2561 //================================================================================
2562
2563 NETGENPlugin_Internals::NETGENPlugin_Internals( SMESH_Mesh&         mesh,
2564                                                 const TopoDS_Shape& shape,
2565                                                 bool                is3D )
2566   : _mesh( mesh ), _is3D( is3D )
2567 {
2568   SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
2569
2570   TopExp_Explorer f,e;
2571   for ( f.Init( shape, TopAbs_FACE ); f.More(); f.Next() )
2572   {
2573     int faceID = meshDS->ShapeToIndex( f.Current() );
2574
2575     // find not computed internal edges
2576
2577     for ( e.Init( f.Current().Oriented(TopAbs_FORWARD), TopAbs_EDGE ); e.More(); e.Next() )
2578       if ( e.Current().Orientation() == TopAbs_INTERNAL )
2579       {
2580         SMESH_subMesh* eSM = mesh.GetSubMesh( e.Current() );
2581         if ( eSM->IsEmpty() )
2582         {
2583           _e2face.insert( make_pair( eSM->GetId(), faceID ));
2584           for ( TopoDS_Iterator v(e.Current()); v.More(); v.Next() )
2585             _e2face.insert( make_pair( meshDS->ShapeToIndex( v.Value() ), faceID ));
2586         }
2587       }
2588
2589     // find internal vertices in a face
2590     set<int> intVV; // issue 0020850 where same vertex is twice in a face
2591     for ( TopoDS_Iterator fSub( f.Current() ); fSub.More(); fSub.Next())
2592       if ( fSub.Value().ShapeType() == TopAbs_VERTEX )
2593       {
2594         int vID = meshDS->ShapeToIndex( fSub.Value() );
2595         if ( intVV.insert( vID ).second )
2596           _f2v[ faceID ].push_back( vID );
2597       }
2598
2599     if ( is3D )
2600     {
2601       // find internal faces and their subshapes where nodes are to be doubled
2602       //  to make a crack with non-sewed borders
2603
2604       if ( f.Current().Orientation() == TopAbs_INTERNAL )
2605       {
2606         _intShapes.insert( meshDS->ShapeToIndex( f.Current() ));
2607
2608         // egdes
2609         list< TopoDS_Shape > edges;
2610         for ( e.Init( f.Current(), TopAbs_EDGE ); e.More(); e.Next())
2611           if ( SMESH_MesherHelper::NbAncestors( e.Current(), mesh, TopAbs_FACE ) > 1 )
2612           {
2613             _intShapes.insert( meshDS->ShapeToIndex( e.Current() ));
2614             edges.push_back( e.Current() );
2615             // find border faces
2616             PShapeIteratorPtr fIt =
2617               SMESH_MesherHelper::GetAncestors( edges.back(),mesh,TopAbs_FACE );
2618             while ( const TopoDS_Shape* pFace = fIt->next() )
2619               if ( !pFace->IsSame( f.Current() ))
2620                 _borderFaces.insert( meshDS->ShapeToIndex( *pFace ));
2621           }
2622         // vertices
2623         // we consider vertex internal if it is shared by more than one internal edge
2624         list< TopoDS_Shape >::iterator edge = edges.begin();
2625         for ( ; edge != edges.end(); ++edge )
2626           for ( TopoDS_Iterator v( *edge ); v.More(); v.Next() )
2627           {
2628             set<int> internalEdges;
2629             PShapeIteratorPtr eIt =
2630               SMESH_MesherHelper::GetAncestors( v.Value(),mesh,TopAbs_EDGE );
2631             while ( const TopoDS_Shape* pEdge = eIt->next() )
2632             {
2633               int edgeID = meshDS->ShapeToIndex( *pEdge );
2634               if ( isInternalShape( edgeID ))
2635                 internalEdges.insert( edgeID );
2636             }
2637             if ( internalEdges.size() > 1 )
2638               _intShapes.insert( meshDS->ShapeToIndex( v.Value() ));
2639           }
2640       }
2641     }
2642   } // loop on geom faces
2643
2644   // find vertices internal in solids
2645   if ( is3D )
2646   {
2647     for ( TopExp_Explorer so(shape, TopAbs_SOLID); so.More(); so.Next())
2648     {
2649       int soID = meshDS->ShapeToIndex( so.Current() );
2650       for ( TopoDS_Iterator soSub( so.Current() ); soSub.More(); soSub.Next())
2651         if ( soSub.Value().ShapeType() == TopAbs_VERTEX )
2652           _s2v[ soID ].push_back( meshDS->ShapeToIndex( soSub.Value() ));
2653     }
2654   }
2655 }
2656
2657 //================================================================================
2658 /*!
2659  * \brief Find mesh faces on non-internal geom faces sharing internal edge
2660  * some nodes of which are to be doubled to make the second border of the "crack"
2661  */
2662 //================================================================================
2663
2664 void NETGENPlugin_Internals::findBorderElements( TIDSortedElemSet & borderElems )
2665 {
2666   if ( _intShapes.empty() ) return;
2667
2668   SMESH_Mesh& mesh = const_cast<SMESH_Mesh&>(_mesh);
2669   SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
2670
2671   // loop on internal geom edges
2672   set<int>::const_iterator intShapeId = _intShapes.begin();
2673   for ( ; intShapeId != _intShapes.end(); ++intShapeId )
2674   {
2675     const TopoDS_Shape& s = meshDS->IndexToShape( *intShapeId );
2676     if ( s.ShapeType() != TopAbs_EDGE ) continue;
2677
2678     // get internal and non-internal geom faces sharing the internal edge <s>
2679     int intFace = 0;
2680     set<int>::iterator bordFace = _borderFaces.end();
2681     PShapeIteratorPtr faces = SMESH_MesherHelper::GetAncestors( s, _mesh, TopAbs_FACE );
2682     while ( const TopoDS_Shape* pFace = faces->next() )
2683     {
2684       int faceID = meshDS->ShapeToIndex( *pFace );
2685       if ( isInternalShape( faceID ))
2686         intFace = faceID;
2687       else
2688         bordFace = _borderFaces.insert( faceID ).first;
2689     }
2690     if ( bordFace == _borderFaces.end() || !intFace ) continue;
2691
2692     // get all links of mesh faces on internal geom face sharing nodes on edge <s>
2693     set< SMESH_OrientedLink > links; //!< links of faces on internal geom face
2694     list<const SMDS_MeshElement*> suspectFaces[2]; //!< mesh faces on border geom faces
2695     int nbSuspectFaces = 0;
2696     SMESHDS_SubMesh* intFaceSM = meshDS->MeshElements( intFace );
2697     if ( !intFaceSM || intFaceSM->NbElements() == 0 ) continue;
2698     SMESH_subMeshIteratorPtr smIt = mesh.GetSubMesh( s )->getDependsOnIterator(true,true);
2699     while ( smIt->more() )
2700     {
2701       SMESHDS_SubMesh* sm = smIt->next()->GetSubMeshDS();
2702       if ( !sm ) continue;
2703       SMDS_NodeIteratorPtr nIt = sm->GetNodes();
2704       while ( nIt->more() )
2705       {
2706         const SMDS_MeshNode* nOnEdge = nIt->next();
2707         SMDS_ElemIteratorPtr fIt = nOnEdge->GetInverseElementIterator(SMDSAbs_Face);
2708         while ( fIt->more() )
2709         {
2710           const SMDS_MeshElement* f = fIt->next();
2711           int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
2712           if ( intFaceSM->Contains( f ))
2713           {
2714             for ( int i = 0; i < nbNodes; ++i )
2715               links.insert( SMESH_OrientedLink( f->GetNode(i), f->GetNode((i+1)%nbNodes)));
2716           }
2717           else
2718           {
2719             int nbDblNodes = 0;
2720             for ( int i = 0; i < nbNodes; ++i )
2721               nbDblNodes += isInternalShape( f->GetNode(i)->getshapeId() );
2722             if ( nbDblNodes )
2723               suspectFaces[ nbDblNodes < 2 ].push_back( f );
2724             nbSuspectFaces++;
2725           }
2726         }
2727       }
2728     }
2729     // suspectFaces[0] having link with same orientation as mesh faces on
2730     // the internal geom face are <borderElems>. suspectFaces[1] have
2731     // only one node on edge <s>, we decide on them later (at the 2nd loop)
2732     // by links of <borderElems> found at the 1st and 2nd loops
2733     set< SMESH_OrientedLink > borderLinks;
2734     for ( int isPostponed = 0; isPostponed < 2; ++isPostponed )
2735     {
2736       list<const SMDS_MeshElement*>::iterator fIt = suspectFaces[isPostponed].begin();
2737       for ( int nbF = 0; fIt != suspectFaces[isPostponed].end(); ++fIt, ++nbF )
2738       {
2739         const SMDS_MeshElement* f = *fIt;
2740         bool isBorder = false, linkFound = false, borderLinkFound = false;
2741         list< SMESH_OrientedLink > faceLinks;
2742         int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
2743         for ( int i = 0; i < nbNodes; ++i )
2744         {
2745           SMESH_OrientedLink link( f->GetNode(i), f->GetNode((i+1)%nbNodes));
2746           faceLinks.push_back( link );
2747           if ( !linkFound )
2748           {
2749             set< SMESH_OrientedLink >::iterator foundLink = links.find( link );
2750             if ( foundLink != links.end() )
2751             {
2752               linkFound= true;
2753               isBorder = ( foundLink->_reversed == link._reversed );
2754               if ( !isBorder && !isPostponed ) break;
2755               faceLinks.pop_back();
2756             }
2757             else if ( isPostponed && !borderLinkFound )
2758             {
2759               foundLink = borderLinks.find( link );
2760               if ( foundLink != borderLinks.end() )
2761               {
2762                 borderLinkFound = true;
2763                 isBorder = ( foundLink->_reversed != link._reversed );
2764               }
2765             }
2766           }
2767         }
2768         if ( isBorder )
2769         {
2770           borderElems.insert( f );
2771           borderLinks.insert( faceLinks.begin(), faceLinks.end() );
2772         }
2773         else if ( !linkFound && !borderLinkFound )
2774         {
2775           suspectFaces[1].push_back( f );
2776           if ( nbF > 2 * nbSuspectFaces )
2777             break; // dead loop protection
2778         }
2779       }
2780     }
2781   }
2782 }
2783
2784 //================================================================================
2785 /*!
2786  * \brief put internal shapes in maps and fill in submeshes to precompute
2787  */
2788 //================================================================================
2789
2790 void NETGENPlugin_Internals::getInternalEdges( TopTools_IndexedMapOfShape& fmap,
2791                                                TopTools_IndexedMapOfShape& emap,
2792                                                TopTools_IndexedMapOfShape& vmap,
2793                                                list< SMESH_subMesh* > smToPrecompute[])
2794 {
2795   if ( !hasInternalEdges() ) return;
2796   map<int,int>::const_iterator ev_face = _e2face.begin();
2797   for ( ; ev_face != _e2face.end(); ++ev_face )
2798   {
2799     const TopoDS_Shape& ev   = _mesh.GetMeshDS()->IndexToShape( ev_face->first );
2800     const TopoDS_Shape& face = _mesh.GetMeshDS()->IndexToShape( ev_face->second );
2801
2802     ( ev.ShapeType() == TopAbs_EDGE ? emap : vmap ).Add( ev );
2803     fmap.Add( face );
2804     //cout<<"INTERNAL EDGE or VERTEX "<<ev_face->first<<" on face "<<ev_face->second<<endl;
2805
2806     smToPrecompute[ MeshDim_1D ].push_back( _mesh.GetSubMeshContaining( ev_face->first ));
2807   }
2808 }
2809
2810 //================================================================================
2811 /*!
2812  * \brief return shapes and submeshes to be meshed and already meshed boundary submeshes
2813  */
2814 //================================================================================
2815
2816 void NETGENPlugin_Internals::getInternalFaces( TopTools_IndexedMapOfShape& fmap,
2817                                                TopTools_IndexedMapOfShape& emap,
2818                                                list< SMESH_subMesh* >&     intFaceSM,
2819                                                list< SMESH_subMesh* >&     boundarySM)
2820 {
2821   if ( !hasInternalFaces() ) return;
2822
2823   // <fmap> and <emap> are for not yet meshed shapes
2824   // <intFaceSM> is for submeshes of faces
2825   // <boundarySM> is for meshed edges and vertices
2826
2827   intFaceSM.clear();
2828   boundarySM.clear();
2829
2830   set<int> shapeIDs ( _intShapes );
2831   if ( !_borderFaces.empty() )
2832     shapeIDs.insert( _borderFaces.begin(), _borderFaces.end() );
2833
2834   set<int>::const_iterator intS = shapeIDs.begin();
2835   for ( ; intS != shapeIDs.end(); ++intS )
2836   {
2837     SMESH_subMesh* sm = _mesh.GetSubMeshContaining( *intS );
2838
2839     if ( sm->GetSubShape().ShapeType() != TopAbs_FACE ) continue;
2840
2841     intFaceSM.push_back( sm );
2842
2843     // add submeshes of not computed internal faces
2844     if ( !sm->IsEmpty() ) continue;
2845
2846     SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(true,true);
2847     while ( smIt->more() )
2848     {
2849       sm = smIt->next();
2850       const TopoDS_Shape& s = sm->GetSubShape();
2851
2852       if ( sm->IsEmpty() )
2853       {
2854         // not yet meshed
2855         switch ( s.ShapeType() ) {
2856         case TopAbs_FACE: fmap.Add ( s ); break;
2857         case TopAbs_EDGE: emap.Add ( s ); break;
2858         default:;
2859         }
2860       }
2861       else
2862       {
2863         if ( s.ShapeType() != TopAbs_FACE )
2864           boundarySM.push_back( sm );
2865       }
2866     }
2867   }
2868 }
2869
2870 //================================================================================
2871 /*!
2872  * \brief Return true if given shape is to be precomputed in order to be correctly
2873  * added to netgen mesh
2874  */
2875 //================================================================================
2876
2877 bool NETGENPlugin_Internals::isShapeToPrecompute(const TopoDS_Shape& s)
2878 {
2879   int shapeID = _mesh.GetMeshDS()->ShapeToIndex( s );
2880   switch ( s.ShapeType() ) {
2881   case TopAbs_FACE  : break; //return isInternalShape( shapeID ) || isBorderFace( shapeID );
2882   case TopAbs_EDGE  : return isInternalEdge( shapeID );
2883   case TopAbs_VERTEX: break;
2884   default:;
2885   }
2886   return false;
2887 }
2888
2889 //================================================================================
2890 /*!
2891  * \brief Return SMESH
2892  */
2893 //================================================================================
2894
2895 SMESH_Mesh& NETGENPlugin_Internals::getMesh() const
2896 {
2897   return const_cast<SMESH_Mesh&>( _mesh );
2898 }
2899
2900 //================================================================================
2901 /*!
2902  * \brief Initialize netgen library
2903  */
2904 //================================================================================
2905
2906 NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper()
2907 {
2908   Ng_Init();
2909   _ngMesh = Ng_NewMesh();
2910 }
2911
2912 //================================================================================
2913 /*!
2914  * \brief Finish using netgen library
2915  */
2916 //================================================================================
2917
2918 NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper()
2919 {
2920   Ng_DeleteMesh( _ngMesh );
2921   Ng_Exit();
2922   NETGENPlugin_Mesher::RemoveTmpFiles();
2923 }
2924
2925 //================================================================================
2926 /*!
2927  * \brief Set netgen mesh to delete at destruction
2928  */
2929 //================================================================================
2930
2931 void NETGENPlugin_NetgenLibWrapper::setMesh( Ng_Mesh* mesh )
2932 {
2933   if ( _ngMesh )
2934     Ng_DeleteMesh( _ngMesh );
2935   _ngMesh = mesh;
2936 }