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