src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx

   1 //  Copyright (C) 2007-2008  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 // File      : NETGENPlugin_NETGEN_3D.cxx
  24 //             Moved here from SMESH_NETGEN_3D.cxx
  25 // Created   : lundi 27 Janvier 2003
  26 // Author    : Nadir BOUHAMOU (CEA)
  27 // Project   : SALOME
  28 // $Header$
  29 //=============================================================================
  30 //
  31 #include "NETGENPlugin_NETGEN_3D.hxx"
  32
  33 #include "NETGENPlugin_Mesher.hxx"
  34
  35 #include "SMDS_MeshElement.hxx"
  36 #include "SMDS_MeshNode.hxx"
  37 #include "SMESHDS_Mesh.hxx"
  38 #include "SMESH_Comment.hxx"
  39 #include "SMESH_ControlsDef.hxx"
  40 #include "SMESH_Gen.hxx"
  41 #include "SMESH_Mesh.hxx"
  42 #include "SMESH_MesherHelper.hxx"
  43 #include "StdMeshers_QuadToTriaAdaptor.hxx"
  44
  45 #include <BRep_Tool.hxx>
  46 #include <TopExp.hxx>
  47 #include <TopExp_Explorer.hxx>
  48 #include <TopoDS.hxx>
  49
  50 #include <Standard_Failure.hxx>
  51 #include <Standard_ErrorHandler.hxx>
  52
  53 #include "utilities.h"
  54
  55 #include <list>
  56 #include <vector>
  57 #include <map>
  58
  59 /*
  60   Netgen include files
  61 */
  62
  63 namespace nglib {
  64 #include <nglib.h>
  65 }
  66 using namespace nglib;
  67
  68 //=============================================================================
  69 /*!
  70  *
  71  */
  72 //=============================================================================
  73
  74 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
  75                              SMESH_Gen* gen)
  76   : SMESH_3D_Algo(hypId, studyId, gen)
  77 {
  78   MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
  79   _name = "NETGEN_3D";
  80   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  81   _compatibleHypothesis.push_back("MaxElementVolume");
  82
  83   _maxElementVolume = 0.;
  84
  85   _hypMaxElementVolume = NULL;
  86
  87   _requireShape = false; // can work without shape
  88 }
  89
  90 //=============================================================================
  91 /*!
  92  *
  93  */
  94 //=============================================================================
  95
  96 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
  97 {
  98   MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
  99 }
 100
 101 //=============================================================================
 102 /*!
 103  *
 104  */
 105 //=============================================================================
 106
 107 bool NETGENPlugin_NETGEN_3D::CheckHypothesis
 108                          (SMESH_Mesh& aMesh,
 109                           const TopoDS_Shape& aShape,
 110                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 111 {
 112   MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
 113
 114   _hypMaxElementVolume = NULL;
 115   _maxElementVolume = DBL_MAX;
 116
 117   list<const SMESHDS_Hypothesis*>::const_iterator itl;
 118   const SMESHDS_Hypothesis* theHyp;
 119
 120   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
 121   int nbHyp = hyps.size();
 122   if (!nbHyp)
 123   {
 124     aStatus = SMESH_Hypothesis::HYP_OK;
 125     //aStatus = SMESH_Hypothesis::HYP_MISSING;
 126     return true;  // can work with no hypothesis
 127   }
 128
 129   itl = hyps.begin();
 130   theHyp = (*itl); // use only the first hypothesis
 131
 132   string hypName = theHyp->GetName();
 133
 134   bool isOk = false;
 135
 136   if (hypName == "MaxElementVolume")
 137   {
 138     _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
 139     ASSERT(_hypMaxElementVolume);
 140     _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
 141     isOk =true;
 142     aStatus = SMESH_Hypothesis::HYP_OK;
 143   }
 144   else
 145     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 146
 147   return isOk;
 148 }
 149
 150 //=============================================================================
 151 /*!
 152  *Here we are going to use the NETGEN mesher
 153  */
 154 //=============================================================================
 155
 156 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
 157                                      const TopoDS_Shape& aShape)
 158 {
 159   MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
 160
 161   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 162
 163   const int invalid_ID = -1;
 164
 165   SMESH::Controls::Area areaControl;
 166   SMESH::Controls::TSequenceOfXYZ nodesCoords;
 167
 168   // -------------------------------------------------------------------
 169   // get triangles on aShell and make a map of nodes to Netgen node IDs
 170   // -------------------------------------------------------------------
 171
 172   SMESH_MesherHelper helper(aMesh);
 173   SMESH_MesherHelper* myTool = &helper;
 174   bool _quadraticMesh = myTool->IsQuadraticSubMesh(aShape);
 175
 176   typedef map< const SMDS_MeshNode*, int> TNodeToIDMap;
 177   TNodeToIDMap nodeToNetgenID;
 178   list< const SMDS_MeshElement* > triangles;
 179   list< bool >                    isReversed; // orientation of triangles
 180
 181   TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
 182   bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
 183
 184   // for the degeneraged edge: ignore all but one node on it;
 185   // map storing ids of degen edges and vertices and their netgen id:
 186   map< int, int* > degenShapeIdToPtrNgId;
 187   map< int, int* >::iterator shId_ngId;
 188   list< int > degenNgIds;
 189
 190   StdMeshers_QuadToTriaAdaptor Adaptor;
 191   Adaptor.Compute(aMesh,aShape);
 192
 193   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
 194   {
 195     const TopoDS_Shape& aShapeFace = exp.Current();
 196     const SMESHDS_SubMesh * aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
 197     if ( aSubMeshDSFace )
 198     {
 199       bool isRev = false;
 200       if ( checkReverse && helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
 201         // IsReversedSubMesh() can work wrong on strongly curved faces,
 202         // so we use it as less as possible
 203         isRev = SMESH_Algo::IsReversedSubMesh( TopoDS::Face(aShapeFace), meshDS );
 204
 205       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 206       while ( iteratorElem->more() ) // loop on elements on a face
 207       {
 208         // check element
 209         const SMDS_MeshElement* elem = iteratorElem->next();
 210         if ( !elem )
 211           return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 212         bool isTraingle = ( elem->NbNodes()==3 || (_quadraticMesh && elem->NbNodes()==6 ));
 213         if ( !isTraingle ) {
 214           //return error( COMPERR_BAD_INPUT_MESH,
 215           //              SMESH_Comment("Not triangle element ")<<elem->GetID());
 216           // using adaptor
 217           std::list<const SMDS_FaceOfNodes*> faces = Adaptor.GetTriangles(elem);
 218           if(faces.size()==0) {
 219             return error( COMPERR_BAD_INPUT_MESH,
 220                           SMESH_Comment("Not triangles in adaptor for element ")<<elem->GetID());
 221           }
 222           std::list<const SMDS_FaceOfNodes*>::iterator itf = faces.begin();
 223           for(; itf!=faces.end(); itf++ ) {
 224             triangles.push_back( (*itf) );
 225             isReversed.push_back( isRev );
 226             // put triange's nodes to nodeToNetgenID map
 227             SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator();
 228             while ( triangleNodesIt->more() ) {
 229               const SMDS_MeshNode * node =
 230                 static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 231               if(myTool->IsMedium(node))
 232                 continue;
 233               nodeToNetgenID.insert( make_pair( node, invalid_ID ));
 234             }
 235           }
 236         }
 237         else {
 238           // keep a triangle
 239           triangles.push_back( elem );
 240           isReversed.push_back( isRev );
 241           // put elem nodes to nodeToNetgenID map
 242           SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
 243           while ( triangleNodesIt->more() ) {
 244             const SMDS_MeshNode * node =
 245               static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 246             if(myTool->IsMedium(node))
 247               continue;
 248             nodeToNetgenID.insert( make_pair( node, invalid_ID ));
 249           }
 250         }
 251 #ifdef _DEBUG_
 252         // check if a trainge is degenerated
 253         areaControl.GetPoints( elem, nodesCoords );
 254         double area = areaControl.GetValue( nodesCoords );
 255         if ( area <= DBL_MIN ) {
 256           MESSAGE( "Warning: Degenerated " << elem );
 257         }
 258 #endif
 259       }
 260       // look for degeneraged edges and vetices
 261       for (TopExp_Explorer expE(aShapeFace,TopAbs_EDGE);expE.More();expE.Next())
 262       {
 263         TopoDS_Edge aShapeEdge = TopoDS::Edge( expE.Current() );
 264         if ( BRep_Tool::Degenerated( aShapeEdge ))
 265         {
 266           degenNgIds.push_back( invalid_ID );
 267           int* ptrIdOnEdge = & degenNgIds.back();
 268           // remember edge id
 269           int edgeID = meshDS->ShapeToIndex( aShapeEdge );
 270           degenShapeIdToPtrNgId.insert( make_pair( edgeID, ptrIdOnEdge ));
 271           // remember vertex id
 272           int vertexID = meshDS->ShapeToIndex( TopExp::FirstVertex( aShapeEdge ));
 273           degenShapeIdToPtrNgId.insert( make_pair( vertexID, ptrIdOnEdge ));
 274         }
 275       }
 276     }
 277   }
 278   // ---------------------------------
 279   // Feed the Netgen with surface mesh
 280   // ---------------------------------
 281
 282   int Netgen_NbOfNodes = 0;
 283   int Netgen_param2ndOrder = 0;
 284   double Netgen_paramFine = 1.;
 285   double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
 286
 287   double Netgen_point[3];
 288   int Netgen_triangle[3];
 289   int Netgen_tetrahedron[4];
 290
 291   Ng_Init();
 292
 293   Ng_Mesh * Netgen_mesh = Ng_NewMesh();
 294
 295   // set nodes and remember thier netgen IDs
 296   bool isDegen = false, hasDegen = !degenShapeIdToPtrNgId.empty();
 297   TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 298   for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 299   {
 300     const SMDS_MeshNode* node = n_id->first;
 301
 302     // ignore nodes on degenerated edge
 303     if ( hasDegen ) {
 304       int shapeId = node->GetPosition()->GetShapeId();
 305       shId_ngId = degenShapeIdToPtrNgId.find( shapeId );
 306       isDegen = ( shId_ngId != degenShapeIdToPtrNgId.end() );
 307       if ( isDegen && *(shId_ngId->second) != invalid_ID ) {
 308         n_id->second = *(shId_ngId->second);
 309         continue;
 310       }
 311     }
 312     Netgen_point [ 0 ] = node->X();
 313     Netgen_point [ 1 ] = node->Y();
 314     Netgen_point [ 2 ] = node->Z();
 315     Ng_AddPoint(Netgen_mesh, Netgen_point);
 316     n_id->second = ++Netgen_NbOfNodes; // set netgen ID
 317
 318     if ( isDegen ) // all nodes on a degen edge get one netgen ID
 319       *(shId_ngId->second) = n_id->second;
 320   }
 321
 322   // set triangles
 323   list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
 324   list< bool >::iterator                 reverse = isReversed.begin();
 325   for ( ; tria != triangles.end(); ++tria, ++reverse )
 326   {
 327     int i = 0;
 328     SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
 329     while ( triangleNodesIt->more() ) {
 330       const SMDS_MeshNode * node =
 331         static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 332       if(myTool->IsMedium(node))
 333         continue;
 334       Netgen_triangle[ *reverse ? 2 - i : i ] = nodeToNetgenID[ node ];
 335       ++i;
 336     }
 337     if ( !hasDegen ||
 338          // ignore degenerated triangles, they have 2 or 3 same ids
 339          (Netgen_triangle[0] != Netgen_triangle[1] &&
 340           Netgen_triangle[0] != Netgen_triangle[2] &&
 341           Netgen_triangle[2] != Netgen_triangle[1] ))
 342     {
 343       Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 344     }
 345   }
 346
 347   // -------------------------
 348   // Generate the volume mesh
 349   // -------------------------
 350
 351   Ng_Meshing_Parameters Netgen_param;
 352
 353   Netgen_param.secondorder = Netgen_param2ndOrder;
 354   Netgen_param.fineness = Netgen_paramFine;
 355   Netgen_param.maxh = Netgen_paramSize;
 356
 357   Ng_Result status;
 358
 359   try {
 360 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 361     OCC_CATCH_SIGNALS;
 362 #endif
 363     status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
 364   }
 365   catch (Standard_Failure& exc) {
 366     error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
 367     status = NG_VOLUME_FAILURE;
 368   }
 369   catch (...) {
 370     error("Exception in Ng_GenerateVolumeMesh()");
 371     status = NG_VOLUME_FAILURE;
 372   }
 373   if ( GetComputeError()->IsOK() ) {
 374     switch ( status ) {
 375     case NG_SURFACE_INPUT_ERROR:error( status, "NG_SURFACE_INPUT_ERROR");
 376     case NG_VOLUME_FAILURE:     error( status, "NG_VOLUME_FAILURE");
 377     case NG_STL_INPUT_ERROR:    error( status, "NG_STL_INPUT_ERROR");
 378     case NG_SURFACE_FAILURE:    error( status, "NG_SURFACE_FAILURE");
 379     case NG_FILE_NOT_FOUND:     error( status, "NG_FILE_NOT_FOUND");
 380     };
 381   }
 382
 383   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 384
 385   int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
 386
 387   MESSAGE("End of Volume Mesh Generation. status=" << status <<
 388           ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
 389           ", nb tetra: " << Netgen_NbOfTetra);
 390
 391   // -------------------------------------------------------------------
 392   // Feed back the SMESHDS with the generated Nodes and Volume Elements
 393   // -------------------------------------------------------------------
 394
 395   bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
 396   if ( isOK )
 397   {
 398     // vector of nodes in which node index == netgen ID
 399     vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
 400     // insert old nodes into nodeVec
 401     for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
 402       nodeVec.at( n_id->second ) = n_id->first;
 403     }
 404     // create and insert new nodes into nodeVec
 405     int nodeIndex = Netgen_NbOfNodes + 1;
 406     int shapeID = meshDS->ShapeToIndex( aShape );
 407     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 408     {
 409       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 410       SMDS_MeshNode * node = meshDS->AddNode(Netgen_point[0],
 411                                              Netgen_point[1],
 412                                              Netgen_point[2]);
 413       meshDS->SetNodeInVolume(node, shapeID);
 414       nodeVec.at(nodeIndex) = node;
 415     }
 416
 417     // create tetrahedrons
 418     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 419     {
 420       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 421       SMDS_MeshVolume * elt = myTool->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
 422                                                  nodeVec.at( Netgen_tetrahedron[1] ),
 423                                                  nodeVec.at( Netgen_tetrahedron[2] ),
 424                                                  nodeVec.at( Netgen_tetrahedron[3] ));
 425       meshDS->SetMeshElementOnShape(elt, shapeID );
 426     }
 427   }
 428
 429   Ng_DeleteMesh(Netgen_mesh);
 430   Ng_Exit();
 431
 432   NETGENPlugin_Mesher::RemoveTmpFiles();
 433
 434   return (status == NG_OK);
 435 }
 436
 437 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
 438                                      SMESH_MesherHelper* aHelper)
 439 {
 440   MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
 441   const int invalid_ID = -1;
 442   bool _quadraticMesh = false;
 443   typedef map< const SMDS_MeshNode*, int> TNodeToIDMap;
 444   TNodeToIDMap nodeToNetgenID;
 445   list< const SMDS_MeshElement* > triangles;
 446   SMESHDS_Mesh* MeshDS = aHelper->GetMeshDS();
 447
 448   SMESH_MesherHelper::MType MeshType = aHelper->IsQuadraticMesh();
 449
 450   if(MeshType == SMESH_MesherHelper::COMP)
 451     return error( COMPERR_BAD_INPUT_MESH,
 452                   SMESH_Comment("Mesh with linear and quadratic elements given."));
 453   else if (MeshType == SMESH_MesherHelper::QUADRATIC)
 454     _quadraticMesh = true;
 455
 456   StdMeshers_QuadToTriaAdaptor Adaptor;
 457   Adaptor.Compute(aMesh);
 458
 459   SMDS_FaceIteratorPtr iteratorFace = MeshDS->facesIterator();
 460   while(iteratorFace->more()) {
 461     // check element
 462     const SMDS_MeshElement* elem = iteratorFace->next();
 463     if ( !elem )
 464       return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 465     bool isTraingle = ( elem->NbNodes()==3 || (_quadraticMesh && elem->NbNodes()==6 ));
 466     if ( !isTraingle ) {
 467       //return error( COMPERR_BAD_INPUT_MESH,
 468       //              SMESH_Comment("Not triangle element ")<<elem->GetID());
 469       // using adaptor
 470       std::list<const SMDS_FaceOfNodes*> faces = Adaptor.GetTriangles(elem);
 471       if(faces.size()==0) {
 472         return error( COMPERR_BAD_INPUT_MESH,
 473                       SMESH_Comment("Not triangles in adaptor for element ")<<elem->GetID());
 474       }
 475       std::list<const SMDS_FaceOfNodes*>::iterator itf = faces.begin();
 476       for(; itf!=faces.end(); itf++ ) {
 477         triangles.push_back( (*itf) );
 478         // put triange's nodes to nodeToNetgenID map
 479         SMDS_ElemIteratorPtr triangleNodesIt = (*itf)->nodesIterator();
 480         while ( triangleNodesIt->more() ) {
 481           const SMDS_MeshNode * node =
 482             static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 483           if(aHelper->IsMedium(node))
 484             continue;
 485           nodeToNetgenID.insert( make_pair( node, invalid_ID ));
 486         }
 487       }
 488     }
 489     else {
 490       // keep a triangle
 491       triangles.push_back( elem );
 492       // put elem nodes to nodeToNetgenID map
 493       SMDS_ElemIteratorPtr triangleNodesIt = elem->nodesIterator();
 494       while ( triangleNodesIt->more() ) {
 495         const SMDS_MeshNode * node =
 496           static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 497         if(aHelper->IsMedium(node))
 498           continue;
 499         nodeToNetgenID.insert( make_pair( node, invalid_ID ));
 500       }
 501     }
 502   }
 503
 504   // ---------------------------------
 505   // Feed the Netgen with surface mesh
 506   // ---------------------------------
 507
 508   int Netgen_NbOfNodes = 0;
 509   int Netgen_param2ndOrder = 0;
 510   double Netgen_paramFine = 1.;
 511   double Netgen_paramSize = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
 512
 513   double Netgen_point[3];
 514   int Netgen_triangle[3];
 515   int Netgen_tetrahedron[4];
 516
 517   Ng_Init();
 518
 519   Ng_Mesh * Netgen_mesh = Ng_NewMesh();
 520
 521     // set nodes and remember thier netgen IDs
 522
 523   TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 524   for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 525   {
 526     const SMDS_MeshNode* node = n_id->first;
 527
 528     Netgen_point [ 0 ] = node->X();
 529     Netgen_point [ 1 ] = node->Y();
 530     Netgen_point [ 2 ] = node->Z();
 531     Ng_AddPoint(Netgen_mesh, Netgen_point);
 532     n_id->second = ++Netgen_NbOfNodes; // set netgen ID
 533
 534   }
 535
 536   // set triangles
 537   list< const SMDS_MeshElement* >::iterator tria = triangles.begin();
 538   for ( ; tria != triangles.end(); ++tria)
 539   {
 540     int i = 0;
 541     SMDS_ElemIteratorPtr triangleNodesIt = (*tria)->nodesIterator();
 542     while ( triangleNodesIt->more() ) {
 543       const SMDS_MeshNode * node =
 544         static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 545       if(aHelper->IsMedium(node))
 546         continue;
 547       Netgen_triangle[ i ] = nodeToNetgenID[ node ];
 548       ++i;
 549     }
 550
 551     Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 552   }
 553
 554   // -------------------------
 555   // Generate the volume mesh
 556   // -------------------------
 557
 558   Ng_Meshing_Parameters Netgen_param;
 559
 560   Netgen_param.secondorder = Netgen_param2ndOrder;
 561   Netgen_param.fineness = Netgen_paramFine;
 562   Netgen_param.maxh = Netgen_paramSize;
 563
 564   Ng_Result status;
 565
 566   try {
 567 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 568     OCC_CATCH_SIGNALS;
 569 #endif
 570     status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
 571   }
 572   catch (Standard_Failure& exc) {
 573     error(COMPERR_OCC_EXCEPTION, exc.GetMessageString());
 574     status = NG_VOLUME_FAILURE;
 575   }
 576   catch (...) {
 577     error("Bad mesh input!!!");
 578     status = NG_VOLUME_FAILURE;
 579   }
 580   if ( GetComputeError()->IsOK() ) {
 581     error( status, "Bad mesh input!!!");
 582   }
 583
 584   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 585
 586   int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
 587
 588   MESSAGE("End of Volume Mesh Generation. status=" << status <<
 589           ", nb new nodes: " << Netgen_NbOfNodesNew - Netgen_NbOfNodes <<
 590           ", nb tetra: " << Netgen_NbOfTetra);
 591
 592   // -------------------------------------------------------------------
 593   // Feed back the SMESHDS with the generated Nodes and Volume Elements
 594   // -------------------------------------------------------------------
 595
 596   bool isOK = ( Netgen_NbOfTetra > 0 );// get whatever built
 597   if ( isOK )
 598   {
 599     // vector of nodes in which node index == netgen ID
 600     vector< const SMDS_MeshNode* > nodeVec ( Netgen_NbOfNodesNew + 1 );
 601     // insert old nodes into nodeVec
 602     for ( n_id = nodeToNetgenID.begin(); n_id != nodeToNetgenID.end(); ++n_id ) {
 603       nodeVec.at( n_id->second ) = n_id->first;
 604     }
 605     // create and insert new nodes into nodeVec
 606     int nodeIndex = Netgen_NbOfNodes + 1;
 607
 608     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 609     {
 610       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 611       SMDS_MeshNode * node = aHelper->AddNode(Netgen_point[0],
 612                                               Netgen_point[1],
 613                                               Netgen_point[2]);
 614       nodeVec.at(nodeIndex) = node;
 615     }
 616
 617     // create tetrahedrons
 618     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 619     {
 620       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 621       aHelper->AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
 622                           nodeVec.at( Netgen_tetrahedron[1] ),
 623                           nodeVec.at( Netgen_tetrahedron[2] ),
 624                           nodeVec.at( Netgen_tetrahedron[3] ));
 625     }
 626   }
 627
 628   Ng_DeleteMesh(Netgen_mesh);
 629   Ng_Exit();
 630
 631   NETGENPlugin_Mesher::RemoveTmpFiles();
 632
 633   return (status == NG_OK);
 634 }