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