src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx

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