0020682: EDF 1222 SMESH: 3D mesh from a skin mesh and with volumic cells

* Redesign to solve pb with internal face (StudyFiss_bugNetgen3D.hdf)

src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx

@@ -42,11 +42,12 @@
 #include "SMESH_MeshEditor.hxx"
 #include "StdMeshers_QuadToTriaAdaptor.hxx"
 
+#include <BRepGProp.hxx>
 #include <BRep_Tool.hxx>
 #include <GProp_GProps.hxx>
-#include <BRepGProp.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
+#include <TopTools_ListIteratorOfListOfShape.hxx>
 #include <TopoDS.hxx>
 
 #include <Standard_Failure.hxx>
@@ -150,6 +151,158 @@ bool NETGENPlugin_NETGEN_3D::CheckHypothesis
   return isOk;
 }
 
+namespace
+{
+  //================================================================================
+  /*!
+   * \brief It correctly initializes netgen library at constructor and
+   *        correctly finishes using netgen library at destructor
+   */
+  //================================================================================
+
+  struct TNetgenLibWrapper
+  {
+    Ng_Mesh* _ngMesh;
+    TNetgenLibWrapper()
+    {
+      Ng_Init();
+      _ngMesh = Ng_NewMesh();
+    }
+    ~TNetgenLibWrapper()
+    {
+      Ng_DeleteMesh( _ngMesh );
+      Ng_Exit();
+      NETGENPlugin_Mesher::RemoveTmpFiles();
+    }
+  };
+
+  //================================================================================
+  /*!
+   * \brief Find mesh faces on non-internal geom faces sharing internal edge
+   *        some nodes of which are to be doubled to make the second border of the "crack"
+   */
+  //================================================================================
+
+  void findBorders( const set<int>&     internalShapeIds,
+                    SMESH_MesherHelper& helper,
+                    TIDSortedElemSet &  borderElems,
+                    set<int> &          borderFaceIds )
+  {
+    SMESH_Mesh*     mesh = helper.GetMesh();
+    SMESHDS_Mesh* meshDS = helper.GetMeshDS();
+
+    // loop on internal geom edges
+    set<int>::const_iterator intShapeId = internalShapeIds.begin();
+    for ( ; intShapeId != internalShapeIds.end(); ++intShapeId )
+    {
+      const TopoDS_Shape& s = meshDS->IndexToShape( *intShapeId );
+      if ( s.ShapeType() != TopAbs_EDGE ) continue;
+
+      // get internal and non-internal geom faces sharing the internal edge <s>
+      int intFace = 0;
+      set<int>::iterator bordFace = borderFaceIds.end();
+      TopTools_ListIteratorOfListOfShape ancIt( mesh->GetAncestors( s ));
+      for ( ; ancIt.More(); ancIt.Next() )
+      {
+        if ( ancIt.Value().ShapeType() != TopAbs_FACE ) continue;
+        int faceID = meshDS->ShapeToIndex( ancIt.Value() );
+        if ( internalShapeIds.count( faceID ))
+          intFace = faceID;
+        else
+          bordFace = borderFaceIds.insert( faceID ).first;
+      }
+      if ( bordFace == borderFaceIds.end() || !intFace ) continue;
+
+      // get all links of mesh faces on internal geom face sharing nodes on edge <s>
+      set< SMESH_OrientedLink > links; //!< links of faces on internal geom face
+      TIDSortedElemSet suspectFaces;   //!< mesh faces on border geom faces
+      SMESHDS_SubMesh* intFaceSM = meshDS->MeshElements( intFace );
+      if ( !intFaceSM || intFaceSM->NbElements() == 0 ) continue;
+      SMESH_subMeshIteratorPtr smIt = mesh->GetSubMesh( s )->getDependsOnIterator(true,true);
+      while ( smIt->more() )
+      {
+        SMESHDS_SubMesh* sm = smIt->next()->GetSubMeshDS();
+        if ( !sm ) continue;
+        SMDS_NodeIteratorPtr nIt = sm->GetNodes();
+        while ( nIt->more() )
+        {
+          const SMDS_MeshNode* nOnEdge = nIt->next();
+          SMDS_ElemIteratorPtr fIt = nOnEdge->GetInverseElementIterator(SMDSAbs_Face);
+          while ( fIt->more() )
+          {
+            const SMDS_MeshElement* f = fIt->next();
+            if ( intFaceSM->Contains( f ))
+            {
+              int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
+              for ( int i = 0; i < nbNodes; ++i )
+                links.insert( SMESH_OrientedLink( f->GetNode(i), f->GetNode((i+1)%nbNodes)));
+            }
+            else
+            {
+              suspectFaces.insert( f );
+            }
+          }
+        }
+      }
+      // <suspectFaces> having link with same orientation as mesh faces on
+      // the internal geom face are <borderElems>.
+      // Some of them have only one node on edge s, we collect them to decide
+      // later by links of found <borderElems>
+      TIDSortedElemSet posponedFaces;
+      set< SMESH_OrientedLink > borderLinks;
+      TIDSortedElemSet::iterator fIt = suspectFaces.begin();
+      for ( ; fIt != suspectFaces.end(); ++fIt )
+      {
+        const SMDS_MeshElement* f = *fIt;
+        bool linkFound = false, isBorder = false;
+        list< SMESH_OrientedLink > faceLinks;
+        int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
+        for ( int i = 0; i < nbNodes; ++i )
+        {
+          SMESH_OrientedLink link( f->GetNode(i), f->GetNode((i+1)%nbNodes));
+          faceLinks.push_back( link );
+          if ( !linkFound )
+          {
+            set< SMESH_OrientedLink >::iterator foundLink = links.find( link );
+            if ( foundLink != links.end() )
+            {
+              linkFound= true;
+              isBorder = ( foundLink->_reversed == link._reversed );
+              if ( !isBorder ) break;
+            }
+          }
+        }
+        if ( isBorder )
+        {
+          borderElems.insert( f );
+          borderLinks.insert( faceLinks.begin(), faceLinks.end() );
+        }
+        else if ( !linkFound )
+        {
+          posponedFaces.insert( f );
+        }
+      }
+      // decide on posponedFaces
+      for ( fIt = posponedFaces.begin(); fIt != posponedFaces.end(); ++fIt )
+      {
+        const SMDS_MeshElement* f = *fIt;
+        int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
+        for ( int i = 0; i < nbNodes; ++i )
+        {
+          SMESH_OrientedLink link( f->GetNode(i), f->GetNode((i+1)%nbNodes));
+          set< SMESH_OrientedLink >::iterator foundLink = borderLinks.find( link );
+          if ( foundLink != borderLinks.end() )
+          {
+            if ( foundLink->_reversed != link._reversed )
+              borderElems.insert( f );
+            break;
+          }
+        }
+      }
+    }
+  }
+}
+
 //=============================================================================
 /*!
  *Here we are going to use the NETGEN mesher
@@ -163,122 +316,8 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
 
   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 
-  const int invalid_ID = -1;
-
-  SMESH::Controls::Area areaControl;
-  SMESH::Controls::TSequenceOfXYZ nodesCoords;
-
-  // -------------------------------------------------------------------
-  // get triangles on aShell and make a map of nodes to Netgen node IDs
-  // -------------------------------------------------------------------
-
   SMESH_MesherHelper helper(aMesh);
-  SMESH_MesherHelper* myTool = &helper;
+  bool _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
-  bool _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
-
-  typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
-  TNodeToIDMap nodeToNetgenID;
-  list< const SMDS_MeshElement* > triangles;
-  list< bool >                    isReversed; // orientation of triangles
-
-  TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
-  bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
-
-  // for the degeneraged edge: ignore all but one node on it;
-  // map storing ids of degen edges and vertices and their netgen id:
-  map< int, int* > degenShapeIdToPtrNgId;
-  map< int, int* >::iterator shId_ngId;
-  list< int > degenNgIds;
-
-  StdMeshers_QuadToTriaAdaptor Adaptor;
-  Adaptor.Compute(aMesh,aShape);
-
-  for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
-  {
-    const TopoDS_Shape& aShapeFace = exp.Current();
-    const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
-    if ( aSubMeshDSFace )
-    {
-      bool isRev = false;
-      if ( checkReverse && helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
-        // IsReversedSubMesh() can work wrong on strongly curved faces,
-        // so we use it as less as possible
-        isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
-
-      SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
-      while ( iteratorElem->more() ) // loop on elements on a face
-      {
-        // check element
-        const SMDS_MeshElement* elem = iteratorElem->next();
-        if ( !elem )
-          return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
-        bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 ));
-        if ( !isTraingle ) {
-          // using adaptor
-          const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
-          if(faces==0) {
-            return error( COMPERR_BAD_INPUT_MESH,
-                          SMESH_Comment("No triangles in adaptor for element ")<<elem->GetID());
-          }
-          list<const SMDS_FaceOfNodes*>::const_iterator itf = faces->begin();
-          for(; itf!=faces->end(); itf++ ) {
-            triangles.push_back( (*itf) );
-            isReversed.push_back( isRev );
-            // put triange's nodes to nodeToNetgenID map
-            SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator();
-            while ( triangleNodesIt->more() ) {
-              const SMDS_MeshNode * node =
-                static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
-              if(myTool->IsMedium(node))
-                continue;
-              nodeToNetgenID.insert( make_pair( node, invalid_ID ));
-            }
-          }
-        }
-        else {
-          // keep a triangle
-          triangles.push_back( elem );
-          isReversed.push_back( isRev );
-          // put elem nodes to nodeToNetgenID map
-          SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
-          while ( triangleNodesIt->more() ) {
-            const SMDS_MeshNode * node =
-              static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
-            if(myTool->IsMedium(node))
-              continue;
-            nodeToNetgenID.insert( make_pair( node, invalid_ID ));
-          }
-        }
-#ifdef _DEBUG_
-        // check if a trainge is degenerated
-        areaControl.GetPoints( elem, nodesCoords );
-        double area = areaControl.GetValue( nodesCoords );
-        if ( area <= DBL_MIN ) {
-          MESSAGE( "Warning: Degenerated " << elem );
-        }
-#endif
-      }
-      // look for degeneraged edges and vetices
-      for (TopExp_Explorer expE(aShapeFace,TopAbs_EDGE);expE.More();expE.Next())
-      {
-        TopoDS_Edge aShapeEdge = TopoDS::Edge( expE.Current() );
-        if ( BRep_Tool::Degenerated( aShapeEdge ))
-        {
-          degenNgIds.push_back( invalid_ID );
-          int* ptrIdOnEdge = & degenNgIds.back();
-          // remember edge id
-          int edgeID = meshDS->ShapeToIndex( aShapeEdge );
-          degenShapeIdToPtrNgId.insert( make_pair( edgeID, ptrIdOnEdge ));
-          // remember vertex id
-          int vertexID = meshDS->ShapeToIndex( TopExp::FirstVertex( aShapeEdge ));
-          degenShapeIdToPtrNgId.insert( make_pair( vertexID, ptrIdOnEdge ));
-        }
-      }
-    }
-  }
-  // ---------------------------------
-  // Feed the Netgen with surface mesh
-  // ---------------------------------
 
   int Netgen_NbOfNodes     = 0;
   int Netgen_param2ndOrder = 0;
@@ -289,61 +328,184 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
   int Netgen_triangle[3];
   int Netgen_tetrahedron[4];
 
-  Ng_Init();
+  TNetgenLibWrapper ngLib;
+  Ng_Mesh * Netgen_mesh = ngLib._ngMesh;
 
-  Ng_Mesh * Netgen_mesh = Ng_NewMesh();
+  // maps of 1) ordinary nodes and 2) doubled nodes on internal shapes
+  typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
+  typedef TNodeToIDMap::value_type                     TN2ID;
+  TNodeToIDMap nodeToNetgenID[2];
 
-  // set nodes and remember thier netgen IDs
-  bool isDegen = false, hasDegen = !degenShapeIdToPtrNgId.empty();
-  TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
-  for ( ; n_id != nodeToNetgenID.end(); ++n_id )
   {
-    const SMDS_MeshNode* node = n_id->first;
+    const int invalid_ID = -1;
 
-    // ignore nodes on degenerated edge
+    SMESH::Controls::Area areaControl;
-    if ( hasDegen ) {
+    SMESH::Controls::TSequenceOfXYZ nodesCoords;
-      int shapeId = node->GetPosition()->GetShapeId();
+
-      shId_ngId = degenShapeIdToPtrNgId.find( shapeId );
+    // --------------------------------------------------------------------
-      isDegen = ( shId_ngId != degenShapeIdToPtrNgId.end() );
+    // Issue 0020676 (StudyFiss_bugNetgen3D.hdf). Pb with internal face.
-      if ( isDegen && *(shId_ngId->second) != invalid_ID ) {
+    // Find internal geom faces, edges and vertices. 
-        n_id->second = *(shId_ngId->second);
+    // Nodes and faces built on the found internal shapes
-        continue;
+    // will be doubled in Netgen input to make two borders of the "crack".
+    // --------------------------------------------------------------------
+
+    set<int> internalShapeIds;
+    set<int> borderFaceIds; //!< non-internal geom faces sharing internal edge
+
+    // mesh faces on <borderFaceIds>, having nodes on internal edge that are
+    // to be replaced by doubled nodes
+    TIDSortedElemSet borderElems;
+
+    // find "internal" faces and edges
+    TopExp_Explorer exFa(aShape,TopAbs_FACE), exEd, exVe;
+    for ( ; exFa.More(); exFa.Next())
+    {
+      if ( exFa.Current().Orientation() == TopAbs_INTERNAL )
+      {
+        internalShapeIds.insert( meshDS->ShapeToIndex( exFa.Current() ));
+        for ( exEd.Init( exFa.Current(), TopAbs_EDGE ); exEd.More(); exEd.Next())
+          if ( helper.NbAncestors( exEd.Current(), aMesh, TopAbs_FACE ) > 1 )
+            internalShapeIds.insert( meshDS->ShapeToIndex( exEd.Current() ));
       }
     }
+    if ( !internalShapeIds.empty() )
+    {
+      // find internal vertices,
+      // we consider vertex internal if it is shared by more than one internal edge
+      TopTools_ListIteratorOfListOfShape ancIt;
+      set<int>::iterator intShapeId = internalShapeIds.begin();
+      for ( ; intShapeId != internalShapeIds.end(); ++intShapeId )
+      {
+        const TopoDS_Shape& s = meshDS->IndexToShape( *intShapeId );
+        if ( s.ShapeType() != TopAbs_EDGE ) continue;
+        for ( exVe.Init( s, TopAbs_VERTEX ); exVe.More(); exVe.Next())
+        {
+          set<int> internalEdges;
+          for ( ancIt.Initialize( aMesh.GetAncestors( exVe.Current() ));
+                ancIt.More(); ancIt.Next() )
+          {
+            if ( ancIt.Value().ShapeType() != TopAbs_EDGE ) continue;
+            int edgeID = meshDS->ShapeToIndex( ancIt.Value() );
+            if ( internalShapeIds.count( edgeID ))
+              internalEdges.insert( edgeID );
+          }
+          if ( internalEdges.size() > 1 )
+            internalShapeIds.insert( meshDS->ShapeToIndex( exVe.Current() ));
+        }
+      }
+      // find border shapes and mesh elements
+      findBorders( internalShapeIds, helper, borderElems, borderFaceIds );
+    }
+
+    // ---------------------------------
+    // Feed the Netgen with surface mesh
+    // ---------------------------------
+
+    TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
+    bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
+
+    StdMeshers_QuadToTriaAdaptor Adaptor;
+    if ( aMesh.NbQuadrangles() > 0 )
+      Adaptor.Compute(aMesh,aShape);
+
+    for ( exFa.ReInit(); exFa.More(); exFa.Next())
+    {
+      const TopoDS_Shape& aShapeFace = exFa.Current();
+      int faceID = meshDS->ShapeToIndex( aShapeFace );
+      bool isInternalFace = internalShapeIds.count( faceID );
+      bool isBorderFace   = borderFaceIds.count( faceID );
+      bool isRev = false;
+      if ( checkReverse && !isInternalFace &&
+           helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
+        // IsReversedSubMesh() can work wrong on strongly curved faces,
+        // so we use it as less as possible
+        isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
+
+      const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
+      if ( !aSubMeshDSFace ) continue;
+      SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
+      while ( iteratorElem->more() ) // loop on elements on a geom face
+      {
+        // check mesh face
+        const SMDS_MeshElement* elem = iteratorElem->next();
+        if ( !elem )
+          return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
+        vector< const SMDS_MeshElement* > trias;
+        bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 ));
+        if ( !isTraingle )
+        {
+          // use adaptor to convert quadrangle face into triangles
+          const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
+          if(faces==0)
+            return error( COMPERR_BAD_INPUT_MESH,
+                          SMESH_Comment("No triangles in adaptor for element ")<<elem->GetID());
+          trias.assign( faces->begin(), faces->end() );
+        }
+        else
+        {
+          trias.push_back( elem );
+        }
+        // Add nodes of triangles and triangles them-selves to netgen mesh
+
+        // a triangle on internal face is added twice,
+        // on border face, once but with doubled nodes
+        bool isBorder = ( isBorderFace && borderElems.count( elem ));
+        int nbDblLoops = ( isInternalFace && isTraingle || isBorder ) ? 2 : 1;
+
+        for ( int i = 0; i < trias.size(); ++i )
+        {
+          bool reverse = isRev;
+          for ( int isDblF = isBorder; isDblF < nbDblLoops; ++isDblF, reverse = !reverse )
+          {
+            // add three nodes of triangle
+            bool hasDegen = false;
+            for ( int iN = 0; iN < 3; ++iN )
+            {
+              const SMDS_MeshNode* node = trias[i]->GetNode( iN );
+              int shapeID = node->GetPosition()->GetShapeId();
+              if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
+                   helper.IsDegenShape( shapeID ))
+              {
+                // ignore all nodes on degeneraged edge and use node on its vertex instead
+                TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
+                node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
+                hasDegen = true;
+              }
+              bool isDblN = isDblF && internalShapeIds.count( shapeID );
+              int& ngID = nodeToNetgenID[isDblN].insert(TN2ID( node, invalid_ID )).first->second;
+              if ( ngID == invalid_ID )
+              {
+                ngID = ++Netgen_NbOfNodes;
                 Netgen_point [ 0 ] = node->X();
                 Netgen_point [ 1 ] = node->Y();
                 Netgen_point [ 2 ] = node->Z();
                 Ng_AddPoint(Netgen_mesh, Netgen_point);
-    n_id->second = ++Netgen_NbOfNodes; // set netgen ID
-
-    if ( isDegen ) // all nodes on a degen edge get one netgen ID
-      *(shId_ngId->second) = n_id->second;
               }
-
+              Netgen_triangle[ iN ] = ngID;
-  // set triangles
-  list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
-  list< bool >::iterator                 reverse = isReversed.begin();
-  for ( ; tria != triangles.end(); ++tria, ++reverse )
-  {
-    int i = 0;
-    SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
-    while ( triangleNodesIt->more() ) {
-      const SMDS_MeshNode * node =
-        static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
-      if(myTool->IsMedium(node))
-        continue;
-      Netgen_triangle[ *reverse ? 2 - i : i ] = nodeToNetgenID[ node ];
-      ++i;
             }
-    if ( !hasDegen ||
+            // add triangle
-         // ignore degenerated triangles, they have 2 or 3 same ids
+            if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
-         (Netgen_triangle[0] != Netgen_triangle[1] &&
+                              Netgen_triangle[0] == Netgen_triangle[2] ||
-          Netgen_triangle[0] != Netgen_triangle[2] &&
+                              Netgen_triangle[2] == Netgen_triangle[1] ))
-          Netgen_triangle[2] != Netgen_triangle[1] ))
+              break;
-    {
+            if ( reverse )
+              swap( Netgen_triangle[1], Netgen_triangle[2] );
+
             Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
           }
         }
+#ifdef _DEBUG_
+        // check if a trainge is degenerated
+        areaControl.GetPoints( elem, nodesCoords );
+        double area = areaControl.GetValue( nodesCoords );
+        if ( area <= DBL_MIN ) {
+          MESSAGE( "Warning: Degenerated " << elem );
+        }
+#endif
+      } // loop on elements on a face
+    } // loop on faces of a SOLID or SHELL
+
+  }
 
   // -------------------------
   // Generate the volume mesh
@@ -393,15 +555,26 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
   // Feed back the SMESHDS with the generated Nodes and Volume Elements
   // -------------------------------------------------------------------
 
-  bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
-  if ( isOK )
-  {
   // vector of nodes in which node index == netgen ID
   vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
   // insert old nodes into nodeVec
-    for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
+  for ( int isDbl = 0; isDbl < 2; ++isDbl )
-      nodeVec.at( n_id->second ) = n_id->first;
+  {
+    TNodeToIDMap::iterator n_id = nodeToNetgenID[isDbl].begin();
+    for ( ; n_id != nodeToNetgenID[isDbl].end(); ++n_id )
+      nodeVec[ n_id->second ] = n_id->first;
+    nodeToNetgenID[isDbl].clear();
   }
+  if ( status == NG_VOLUME_FAILURE )
+  {
+    SMESH_ComputeErrorPtr err = NETGENPlugin_Mesher::readErrors(nodeVec);
+    if ( err && !err->myBadElements.empty() )
+      error( err );
+  }
+
+  bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
+  if ( isOK )
+  {
     // create and insert new nodes into nodeVec
     int nodeIndex = Netgen_NbOfNodes + 1;
     int shapeID = meshDS->ShapeToIndex( aShape );
@@ -419,7 +592,7 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
     {
       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
-      SMDS_MeshVolume * elt = myTool->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
+      SMDS_MeshVolume * elt = helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
                                                 nodeVec.at( Netgen_tetrahedron[1] ),
                                                 nodeVec.at( Netgen_tetrahedron[2] ),
                                                 nodeVec.at( Netgen_tetrahedron[3] ));
@@ -427,14 +600,15 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
     }
   }
 
-  Ng_DeleteMesh(Netgen_mesh);
-  Ng_Exit();
-
-  NETGENPlugin_Mesher::RemoveTmpFiles();
-
   return (status == NG_OK);
 }
 
+//================================================================================
+/*!
+ * \brief Compute tetrahedral mesh from 2D mesh without geometry
+ */
+//================================================================================
+
 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
                                      SMESH_MesherHelper* aHelper)
 {
@@ -462,42 +636,33 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
   while( fIt->more()) sortedFaces.insert( fIt->next() );
 
   TIDSortedElemSet::iterator itFace = sortedFaces.begin(), fEnd = sortedFaces.end();
-  for ( ; itFace != fEnd; ++itFace ) {
+  for ( ; itFace != fEnd; ++itFace )
+  {
     // check element
     const SMDS_MeshElement* elem = *itFace;
     if ( !elem )
       return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
+
+    vector< const SMDS_MeshElement* > trias;
     bool isTraingle = ( elem->NbNodes() == ( elem->IsQuadratic() ? 6 : 3 ));
     if ( !isTraingle ) {
       // using adaptor
       const list<const SMDS_FaceOfNodes*>* faces = Adaptor.GetTriangles(elem);
-      if(faces==0) {
+      if(faces==0)
-        continue; // Issue 0020682. There already can be 3d mesh
+        return error( COMPERR_BAD_INPUT_MESH,
-      }
+                      SMESH_Comment("No triangles in adaptor for element ")<<elem->GetID());
-      list<const SMDS_FaceOfNodes*>::const_iterator itf = faces->begin();
+      trias.assign( faces->begin(), faces->end() );
-      for(; itf!=faces->end(); itf++ ) {
-        triangles.push_back( (*itf) );
-        // put triange's nodes to nodeToNetgenID map
-        SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator();
-        while ( triangleNodesIt->more() ) {
-          const SMDS_MeshNode * node =
-            static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
-          if(aHelper->IsMedium(node))
-            continue;
-          nodeToNetgenID.insert( make_pair( node, invalid_ID ));
-        }
-      }
     }
     else {
-      // keep a triangle
+      trias.push_back( elem );
-      triangles.push_back( elem );
+    }
+    for ( int i = 0; i < trias.size(); ++i )
+    {
+      triangles.push_back( trias[i] );
+      for ( int iN = 0; iN < 3; ++iN )
+      {
+        const SMDS_MeshNode* node = trias[i]->GetNode( iN );
         // put elem nodes to nodeToNetgenID map
-      SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
-      while ( triangleNodesIt->more() ) {
-        const SMDS_MeshNode * node =
-          static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
-        if(aHelper->IsMedium(node))
-          continue;
         nodeToNetgenID.insert( make_pair( node, invalid_ID ));
       }
     }
@@ -516,9 +681,8 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
   int Netgen_triangle[3];
   int Netgen_tetrahedron[4];
 
-  Ng_Init();
+  TNetgenLibWrapper ngLib;
-
+  Ng_Mesh * Netgen_mesh = ngLib._ngMesh;
-  Ng_Mesh * Netgen_mesh = Ng_NewMesh();
 
   // set nodes and remember thier netgen IDs
   
@@ -625,11 +789,6 @@ bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
     }
   }
 
-  Ng_DeleteMesh(Netgen_mesh);
-  Ng_Exit();
-  
-  NETGENPlugin_Mesher::RemoveTmpFiles();
-  
   return (status == NG_OK);
 }
 
@@ -646,8 +805,8 @@ bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh,
 {
   int nbtri = 0, nbqua = 0;
   double fullArea = 0.0;
-  for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
+  for (TopExp_Explorer expF(aShape, TopAbs_FACE); expF.More(); expF.Next()) {
-    TopoDS_Face F = TopoDS::Face( exp.Current() );
+    TopoDS_Face F = TopoDS::Face( expF.Current() );
     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
     MapShapeNbElemsItr anIt = aResMap.find(sm);
     if( anIt==aResMap.end() ) {
@@ -669,12 +828,12 @@ bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh,
   bool IsQuadratic = false;
   bool IsFirst = true;
   TopTools_MapOfShape tmpMap;
-  for (TopExp_Explorer exp(aShape, TopAbs_EDGE); exp.More(); exp.Next()) {
+  for (TopExp_Explorer expF(aShape, TopAbs_EDGE); expF.More(); expF.Next()) {
-    TopoDS_Edge E = TopoDS::Edge(exp.Current());
+    TopoDS_Edge E = TopoDS::Edge(expF.Current());
     if( tmpMap.Contains(E) )
       continue;
     tmpMap.Add(E);
-    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
+    SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(expF.Current());
     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
     if( anIt==aResMap.end() ) {
       SMESH_ComputeErrorPtr& smError = aSubMesh->GetComputeError();