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