Salome HOME
Revert "Merge branch 'yan/parallel_mesh2'"
[plugins/netgenplugin.git] / src / NETGENPlugin / NETGENPlugin_Mesher.cxx
index 1d6ad11616c6d3cb2907b3efb76fe89688595439..732a2441bc6c4d28bded942e0565bacdde988d20 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2020  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2022  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -599,6 +599,18 @@ void NETGENPlugin_Mesher::SetDefaultParameters()
   _fineness               = NETGENPlugin_Hypothesis::GetDefaultFineness();
   mparams.uselocalh       = NETGENPlugin_Hypothesis::GetDefaultSurfaceCurvature();
   netgen::merge_solids    = NETGENPlugin_Hypothesis::GetDefaultFuseEdges();
+
+#ifdef NETGEN_V6
+
+  mparams.nthreads = std::thread::hardware_concurrency();
+
+  if ( getenv( "SALOME_NETGEN_DISABLE_MULTITHREADING" ))
+  {
+    mparams.nthreads = 1;
+    mparams.parallel_meshing = false;
+  }
+
+#endif
 }
 
 //=============================================================================
@@ -669,6 +681,12 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
       }
     }
   }
+
+#ifdef NETGEN_V6
+
+  netgen::mparam.closeedgefac = 2;
+
+#endif
 }
 
 //=============================================================================
@@ -732,7 +750,7 @@ void NETGENPlugin_Mesher::SetLocalSize( netgen::OCCGeometry& occgeo,
     if ( faceNgID >= 1 )
     {
 #ifdef NETGEN_V6
-      occgeo.SetFaceMaxH(faceNgID-1, val);
+      occgeo.SetFaceMaxH(faceNgID-1, val, netgen::mparam);
 #else
       occgeo.SetFaceMaxH(faceNgID, val);
 #endif
@@ -801,7 +819,7 @@ void NETGENPlugin_Mesher::SetLocalSizeForChordalError( netgen::OCCGeometry& occg
       double maxCurv = Max( Abs( surfProp.MaxCurvature()), Abs( surfProp.MinCurvature() ));
       double    size = elemSizeForChordalError( _chordalError, 1 / maxCurv );
 #ifdef NETGEN_V6
-      occgeo.SetFaceMaxH( i-1, size * sizeCoef );
+      occgeo.SetFaceMaxH( i-1, size * sizeCoef, netgen::mparam );
 #else
       occgeo.SetFaceMaxH( i, size * sizeCoef );
 #endif
@@ -1162,7 +1180,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
         const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
         if ( points.empty() )
           return false; // invalid node params?
-        int i, nbSeg = fSide.NbSegments();
+        smIdType i, nbSeg = fSide.NbSegments();
 
         // remember EDGEs of fSide to treat only once
         for ( int iE = 0; iE < fSide.NbEdges(); ++iE )
@@ -1297,11 +1315,11 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
           if ( id ) solidSMDSIDs[ bool( solidSMDSIDs[0] )] = meshDS->ShapeToIndex( *solid );
         }
       }
+      bool isShrunk = true;
       if ( proxyMesh && proxyMesh->GetProxySubMesh( geomFace ))
       {
         // if a proxy sub-mesh contains temporary faces, then these faces
         // should be used to mesh only one SOLID
-        bool hasTmp = false;
         smDS = proxyMesh->GetSubMesh( geomFace );
         SMDS_ElemIteratorPtr faces = smDS->GetElements();
         while ( faces->more() )
@@ -1309,7 +1327,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
           const SMDS_MeshElement* f = faces->next();
           if ( proxyMesh->IsTemporary( f ))
           {
-            hasTmp = true;
+            isShrunk = false;
             if ( solidSMDSIDs[1] && proxyMesh->HasPrismsOnTwoSides( meshDS->MeshElements( geomFace )))
               break;
             else
@@ -1328,7 +1346,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
           }
         }
         const int fID = occgeom.fmap.FindIndex( geomFace );
-        if ( !hasTmp ) // shrunk mesh
+        if ( isShrunk ) // shrunk mesh
         {
           // move netgen points according to moved nodes
           SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
@@ -1380,7 +1398,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
       {
         solidSMDSIDs[1] = 0;
       }
-      const bool hasVLOn2Sides = ( solidSMDSIDs[1] > 0 );
+      const bool hasVLOn2Sides = ( solidSMDSIDs[1] > 0 && !isShrunk );
 
       // Add ng face descriptors of meshed faces
       faceNgID++;
@@ -1440,7 +1458,7 @@ bool NETGENPlugin_Mesher::FillNgMesh(netgen::OCCGeometry&           occgeom,
         const SMDS_MeshElement* f = faces->next();
         if ( f->NbNodes() % 3 != 0 ) // not triangle
         {
-          PShapeIteratorPtr solidIt=helper.GetAncestors(geomFace,*sm->GetFather(),TopAbs_SOLID);
+          PShapeIteratorPtr solidIt = helper.GetAncestors( geomFace,*sm->GetFather(),TopAbs_SOLID);
           if ( const TopoDS_Shape * solid = solidIt->next() )
             sm = _mesh->GetSubMesh( *solid );
           SMESH_BadInputElements* badElems =
@@ -2202,7 +2220,7 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
   // ----------------------------
   // Check wires and count nodes
   // ----------------------------
-  int nbNodes = 0;
+  smIdType nbNodes = 0;
   for ( size_t iW = 0; iW < wires.size(); ++iW )
   {
     StdMeshers_FaceSidePtr wire = wires[ iW ];
@@ -2259,7 +2277,7 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
   {
     StdMeshers_FaceSidePtr       wire = wires[ iW ];
     const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
-    const int              nbSegments = wire->NbPoints() - 1;
+    const smIdType         nbSegments = wire->NbPoints() - 1;
 
     // assure the 1st node to be in node2ngID, which is needed to correctly
     // "close chain of segments" (see below) in case if the 1st node is not
@@ -2351,8 +2369,8 @@ NETGENPlugin_Mesher::AddSegmentsToMesh(netgen::Mesh&                    ngMesh,
         SMESH_TNodeXYZ np1( n ), np2( uvPtVec[ i+1 ].node );
         // get an average size of adjacent segments to avoid sharp change of
         // element size (regression on issue 0020452, note 0010898)
-        int   iPrev = SMESH_MesherHelper::WrapIndex( i-1, nbSegments );
-        int   iNext = SMESH_MesherHelper::WrapIndex( i+1, nbSegments );
+        int   iPrev = SMESH_MesherHelper::WrapIndex( i-1, (int) nbSegments );
+        int   iNext = SMESH_MesherHelper::WrapIndex( i+1, (int) nbSegments );
         double sumH = segLen[ iPrev ] + segLen[ i ] + segLen[ iNext ];
         int   nbSeg = ( int( segLen[ iPrev ] > sumH / 100.)  +
                         int( segLen[ i     ] > sumH / 100.)  +
@@ -2465,14 +2483,14 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
   if ( quadHelper && !quadHelper->GetIsQuadratic() && quadHelper->GetTLinkNodeMap().empty() )
     quadHelper = 0;
 
-  int i, nbInitNod = initState._nbNodes;
+  int ngID, nbInitNod = initState._nbNodes;
   if ( initState._elementsRemoved )
   {
     // PAL23427. Update nodeVec to track removal of netgen free points as a result
     // of removal of faces in FillNgMesh() in the case of a shrunk sub-mesh
-    int ngID, nodeVecSize = nodeVec.size();
+    size_t i, nodeVecSize = nodeVec.size();
     const double eps = std::numeric_limits<double>::min();
-    for ( ngID = i = 1; i < nodeVecSize; ++ngID, ++i )
+    for ( i = ngID = 1; i < nodeVecSize; ++ngID, ++i )
     {
       gp_Pnt ngPnt( NGPOINT_COORDS( ngMesh.Point( ngID )));
       gp_Pnt node ( SMESH_NodeXYZ (nodeVec_ACCESS(i) ));
@@ -2494,7 +2512,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
 
   if ( nbNod > nbInitNod )
     nodeVec.resize( nbNod + 1 );
-  for ( i = nbInitNod+1; i <= nbNod; ++i )
+  for ( int i = nbInitNod+1; i <= nbNod; ++i )
   {
     const netgen::MeshPoint& ngPoint = ngMesh.Point(i);
     SMDS_MeshNode* node = NULL;
@@ -2529,7 +2547,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
   // -------------------------------------------
 
   int nbInitSeg = initState._nbSegments;
-  for (i = nbInitSeg+1; i <= nbSeg; ++i )
+  for ( int i = nbInitSeg+1; i <= nbSeg; ++i )
   {
     const netgen::Segment& seg = ngMesh.LineSegment(i);
     TopoDS_Edge aEdge;
@@ -2611,7 +2629,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
     ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(quadFaceID, /*solid1=*/0, /*solid2=*/0, 0));
 
   vector<const SMDS_MeshNode*> nodes;
-  for (i = nbInitFac+1; i <= nbFac; ++i )
+  for ( int i = nbInitFac+1; i <= nbFac; ++i )
   {
     const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
     const int        aGeomFaceInd = elem.GetIndex();
@@ -2692,7 +2710,7 @@ int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry&          occgeo,
   // Create tetrahedra
   // ------------------
 
-  for ( i = 1; i <= nbVol; ++i )
+  for ( int i = 1; i <= nbVol; ++i )
   {
     const netgen::Element& elem = ngMesh.VolumeElement(i);      
     int aSolidInd = elem.GetIndex();
@@ -2957,7 +2975,7 @@ bool NETGENPlugin_Mesher::Compute()
     {
       // Pass 1D simple parameters to NETGEN
       // --------------------------------
-      int      nbSeg = _simpleHyp->GetNumberOfSegments();
+      double nbSeg   = (double) _simpleHyp->GetNumberOfSegments();
       double segSize = _simpleHyp->GetLocalLength();
       for ( int iE = 1; iE <= occgeo.emap.Extent(); ++iE )
       {
@@ -3184,7 +3202,7 @@ bool NETGENPlugin_Mesher::Compute()
       FillSMesh( occgeo, *_ngMesh, initState, *_mesh, nodeVec, comment, &quadHelper );
 
       // compute prismatic boundary volumes
-      int nbQuad = _mesh->NbQuadrangles();
+      smIdType nbQuad = _mesh->NbQuadrangles();
       SMESH_ProxyMesh::Ptr viscousMesh;
       if ( _viscousLayersHyp )
       {
@@ -3228,7 +3246,7 @@ bool NETGENPlugin_Mesher::Compute()
       // fill _ngMesh with faces of sub-meshes
       err = ! ( FillNgMesh(occgeo, *_ngMesh, nodeVec, meshedSM[ MeshDim_2D ], &quadHelper));
       initState = NETGENPlugin_ngMeshInfo(_ngMesh, /*checkRemovedElems=*/true);
-      // toPython( _ngMesh );
+      // toPython( _ngMesh )
     }
     if (!err && _isVolume)
     {
@@ -3442,9 +3460,9 @@ bool NETGENPlugin_Mesher::Compute()
           bool smComputed = nbVol && !sm->IsEmpty();
           if ( smComputed && internals.hasInternalVertexInSolid( sm->GetId() ))
           {
-            int nbIntV = internals.getSolidsWithVertices().find( sm->GetId() )->second.size();
+            size_t nbIntV = internals.getSolidsWithVertices().find( sm->GetId() )->second.size();
             SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
-            smComputed = ( smDS->NbElements() > 0 || smDS->NbNodes() > nbIntV );
+            smComputed = ( smDS->NbElements() > 0 || smDS->NbNodes() > (smIdType) nbIntV );
           }
           SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
           if ( !smComputed && ( !smError || smError->IsOK() ))
@@ -3552,7 +3570,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
   // }
   // calculate total nb of segments and length of edges
   double fullLen = 0.0;
-  int fullNbSeg = 0;
+  smIdType fullNbSeg = 0;
   int entity = mparams.secondorder > 0 ? SMDSEntity_Quad_Edge : SMDSEntity_Edge;
   TopTools_DataMapOfShapeInteger Edge2NbSeg;
   for (TopExp_Explorer exp(_shape, TopAbs_EDGE); exp.More(); exp.Next())
@@ -3564,7 +3582,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
     double aLen = SMESH_Algo::EdgeLength(E);
     fullLen += aLen;
 
-    vector<int>& aVec = aResMap[_mesh->GetSubMesh(E)];
+    vector<smIdType>& aVec = aResMap[_mesh->GetSubMesh(E)];
     if ( aVec.empty() )
       aVec.resize( SMDSEntity_Last, 0);
     else
@@ -3581,7 +3599,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
     int aGeomEdgeInd = seg.epgeominfo[0].edgenr;
     if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
     {
-      vector<int>& aVec = aResMap[_mesh->GetSubMesh(occgeo.emap(aGeomEdgeInd))];
+      vector<smIdType>& aVec = aResMap[_mesh->GetSubMesh(occgeo.emap(aGeomEdgeInd))];
       aVec[ entity ]++;
     }
   }
@@ -3589,12 +3607,12 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
   TopTools_DataMapIteratorOfDataMapOfShapeInteger Edge2NbSegIt(Edge2NbSeg);
   for (; Edge2NbSegIt.More(); Edge2NbSegIt.Next())
   {
-    vector<int>& aVec = aResMap[_mesh->GetSubMesh(Edge2NbSegIt.Key())];
+    vector<smIdType>& aVec = aResMap[_mesh->GetSubMesh(Edge2NbSegIt.Key())];
     if ( aVec[ entity ] > 1 && aVec[ SMDSEntity_Node ] == 0 )
       aVec[SMDSEntity_Node] = mparams.secondorder > 0  ? 2*aVec[ entity ]-1 : aVec[ entity ]-1;
 
     fullNbSeg += aVec[ entity ];
-    Edge2NbSeg( Edge2NbSegIt.Key() ) = aVec[ entity ];
+    Edge2NbSeg( Edge2NbSegIt.Key() ) = (int) aVec[ entity ];
   }
   if ( fullNbSeg == 0 )
     return false;
@@ -3610,12 +3628,12 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
     }
     else {
       // length from edges
-      mparams.maxh = fullLen/fullNbSeg;
+      mparams.maxh = fullLen / double( fullNbSeg );
       mparams.grading = 0.2; // slow size growth
     }
   }
   mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
-  mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
+  mparams.maxh = min( mparams.maxh, fullLen / double( fullNbSeg ) * (1. + mparams.grading));
 
   for (TopExp_Explorer exp(_shape, TopAbs_FACE); exp.More(); exp.Next())
   {
@@ -3636,7 +3654,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
     int nbFaces = tooManyElems ? hugeNb : int( 4*anArea / (mparams.maxh*mparams.maxh*sqrt(3.)));
     int nbNodes = tooManyElems ? hugeNb : (( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
 
-    vector<int> aVec(SMDSEntity_Last, 0);
+    vector<smIdType> aVec(SMDSEntity_Last, 0);
     if( mparams.secondorder > 0 ) {
       int nb1d_in = (nbFaces*3 - nb1d) / 2;
       aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
@@ -3666,7 +3684,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
         // using previous length from faces
       }
       mparams.grading = 0.4;
-      mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
+      mparams.maxh = min( mparams.maxh, fullLen / double( fullNbSeg ) * (1. + mparams.grading));
     }
     GProp_GProps G;
     BRepGProp::VolumeProperties(_shape,G);
@@ -3675,7 +3693,7 @@ bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
     tooManyElems = tooManyElems || ( aVolume/hugeNb > tetrVol );
     int nbVols = tooManyElems ? hugeNb : int(aVolume/tetrVol);
     int nb1d_in = int(( nbVols*6 - fullNbSeg ) / 6 );
-    vector<int> aVec(SMDSEntity_Last, 0 );
+    vector<smIdType> aVec(SMDSEntity_Last, 0 );
     if ( tooManyElems ) // avoid FPE
     {
       aVec[SMDSEntity_Node] = hugeNb;
@@ -3787,8 +3805,8 @@ NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
   vector<int> two(2);
   vector<int> three1(3), three2(3);
   const char* badEdgeStr = " multiple times in surface mesh";
-  const int   badEdgeStrLen = strlen( badEdgeStr );
-  const int   nbNodes = nodeVec.size();
+  const int   badEdgeStrLen = (int) strlen( badEdgeStr );
+  const int   nbNodes = (int) nodeVec.size();
 
   while( !file.eof() )
   {
@@ -3798,7 +3816,7 @@ NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
          two[0] < nbNodes  &&  two[1] < nbNodes )
     {
       err->myBadElements.push_back( new SMDS_LinearEdge( nodeVec[ two[0]], nodeVec[ two[1]] ));
-      file += badEdgeStrLen;
+      file += (int) badEdgeStrLen;
     }
     else if ( strncmp( file, "Intersecting: ", 14 ) == 0 )
     {
@@ -3814,7 +3832,7 @@ NETGENPlugin_Mesher::ReadErrors(const vector<const SMDS_MeshNode* >& nodeVec)
       ok = ok && file.getInts( three2 );
       for ( int i = 0; ok && i < 3; ++i )
         ok = ( three1[i] < nbNodes && nodeVec[ three1[i]]);
-      for ( int i = 0; ok && i < 3; ++i ) 
+      for ( int i = 0; ok && i < 3; ++i )
         ok = ( three2[i] < nbNodes && nodeVec[ three2[i]]);
       if ( ok )
       {
@@ -3865,6 +3883,61 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
           << "mesh = smesh.Mesh()" << std::endl << std::endl;
 
   using namespace netgen;
+
+#ifdef NETGEN_V6
+
+  for ( int i = 1; i <= ngMesh->GetNP(); i++)
+  {
+    const Point3d & p = ngMesh->Point(i);
+    outfile << "mesh.AddNode( ";
+    outfile << p.X() << ", ";
+    outfile << p.Y() << ", ";
+    outfile << p.Z() << ") ## "<< i << std::endl;
+  }
+
+  int nbDom = ngMesh->GetNDomains();
+  for ( int i = 0; i < nbDom; ++i )
+    outfile<< "grp" << i+1 << " = mesh.CreateEmptyGroup( SMESH.FACE, 'domain"<< i+1 << "')"<< std::endl;
+
+  int nbDel = 0;
+  for (int i = 1; i <= ngMesh->GetNSE(); i++)
+  {
+    outfile << "mesh.AddFace([ ";
+    Element2d sel = ngMesh->SurfaceElement(i);
+    for (int j = 1; j <= sel.GetNP(); j++)
+      outfile << sel.PNum(j) << ( j < sel.GetNP() ? ", " : " ])");
+    if ( sel.IsDeleted() ) outfile << " ## IsDeleted ";
+    outfile << std::endl;
+    nbDel += sel.IsDeleted();
+
+    if (sel.GetIndex())
+    {
+      if ( int dom1 = ngMesh->GetFaceDescriptor(sel.GetIndex ()).DomainIn())
+        outfile << "grp"<< dom1 <<".Add([ " << i - nbDel << " ])" << std::endl;
+      if ( int dom2 = ngMesh->GetFaceDescriptor(sel.GetIndex ()).DomainOut())
+        outfile << "grp"<< dom2 <<".Add([ " << i - nbDel << " ])" << std::endl;
+    }
+  }
+
+  for (int i = 1; i <= ngMesh->GetNE(); i++)
+  {
+    Element el = ngMesh->VolumeElement(i);
+    outfile << "mesh.AddVolume([ ";
+    for (int j = 1; j <= el.GetNP(); j++)
+      outfile << el.PNum(j) << ( j < el.GetNP() ? ", " : " ])");
+    outfile << std::endl;
+  }
+
+  for (int i = 1; i <= ngMesh->GetNSeg(); i++)
+  {
+    const Segment & seg = ngMesh->LineSegment (i);
+    outfile << "mesh.AddEdge([ "
+            << seg[0]+1 << ", "
+            << seg[1]+1 << " ])" << std::endl;
+  }
+
+#else  //////// V 5
+
   PointIndex pi;
   for (pi = PointIndex::BASE; 
        pi < ngMesh->GetNP()+PointIndex::BASE; pi++)
@@ -3879,6 +3952,7 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
   for ( int i = 0; i < nbDom; ++i )
     outfile<< "grp" << i+1 << " = mesh.CreateEmptyGroup( SMESH.FACE, 'domain"<< i+1 << "')"<< std::endl;
 
+  int nbDel = 0;
   SurfaceElementIndex sei;
   for (sei = 0; sei < ngMesh->GetNSE(); sei++)
   {
@@ -3888,13 +3962,14 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
       outfile << sel[j] << ( j+1 < sel.GetNP() ? ", " : " ])");
     if ( sel.IsDeleted() ) outfile << " ## IsDeleted ";
     outfile << std::endl;
+    nbDel += sel.IsDeleted();
 
     if ((*ngMesh)[sei].GetIndex())
     {
       if ( int dom1 = ngMesh->GetFaceDescriptor((*ngMesh)[sei].GetIndex ()).DomainIn())
-        outfile << "grp"<< dom1 <<".Add([ " << (int)sei+1 << " ])" << std::endl;
+        outfile << "grp"<< dom1 <<".Add([ " << (int)sei+1 - nbDel << " ])" << std::endl;
       if ( int dom2 = ngMesh->GetFaceDescriptor((*ngMesh)[sei].GetIndex ()).DomainOut())
-        outfile << "grp"<< dom2 <<".Add([ " << (int)sei+1 << " ])" << std::endl;
+        outfile << "grp"<< dom2 <<".Add([ " << (int)sei+1 - nbDel  << " ])" << std::endl;
     }
   }
 
@@ -3914,6 +3989,9 @@ void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh )
             << seg[0] << ", "
             << seg[1] << " ])" << std::endl;
   }
+
+#endif
+
   std::cout << "Write " << pyFile << std::endl;
 }
 
@@ -4426,10 +4504,11 @@ int NETGENPlugin_NetgenLibWrapper::GenerateMesh( netgen::OCCGeometry& occgeo,
                                                  netgen::Mesh* & ngMesh )
 {
   int err = 0;
-#ifdef NETGEN_V6
-
   if ( !ngMesh )
     ngMesh = new netgen::Mesh;
+
+#ifdef NETGEN_V6
+
   ngMesh->SetGeometry( shared_ptr<netgen::NetgenGeometry>( &occgeo, &NOOP_Deleter ));
 
   netgen::mparam.perfstepsstart = startWith;