src/NETGENPlugin/NETGENPlugin_Mesher.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 //  NETGENPlugin : C++ implementation
  23 // File      : NETGENPlugin_Mesher.cxx
  24 // Author    : Michael Sazonov (OCN)
  25 // Date      : 31/03/2006
  26 // Project   : SALOME
  27 //=============================================================================
  28 //
  29 #include "NETGENPlugin_Mesher.hxx"
  30 #include "NETGENPlugin_Hypothesis_2D.hxx"
  31 #include "NETGENPlugin_SimpleHypothesis_3D.hxx"
  32
  33 #include <SMESH_Mesh.hxx>
  34 #include <SMESH_Comment.hxx>
  35 #include <SMESH_ComputeError.hxx>
  36 #include <SMESH_subMesh.hxx>
  37 #include <SMESH_MesherHelper.hxx>
  38 #include <SMESHDS_Mesh.hxx>
  39 #include <SMDS_MeshElement.hxx>
  40 #include <SMDS_MeshNode.hxx>
  41 #include <utilities.h>
  42
  43 #include <vector>
  44 #include <limits>
  45
  46 #include <BRep_Tool.hxx>
  47 #include <TopExp.hxx>
  48 #include <TopExp_Explorer.hxx>
  49 #include <TopoDS.hxx>
  50 #include <NCollection_Map.hxx>
  51 #include <OSD_Path.hxx>
  52 #include <OSD_File.hxx>
  53 #include <TCollection_AsciiString.hxx>
  54 #include <TopTools_ListIteratorOfListOfShape.hxx>
  55 #include <TopTools_DataMapOfShapeInteger.hxx>
  56 #include <Standard_ErrorHandler.hxx>
  57 #include <Standard_ProgramError.hxx>
  58
  59 // Netgen include files
  60 namespace nglib {
  61 #include <nglib.h>
  62 }
  63 #define OCCGEOMETRY
  64 #include <occgeom.hpp>
  65 #include <meshing.hpp>
  66 //#include <ngexception.hpp>
  67 namespace netgen {
  68   extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
  69   extern MeshingParameters mparam;
  70 }
  71
  72 using namespace std;
  73
  74 static void removeFile( const TCollection_AsciiString& fileName )
  75 {
  76   try {
  77     OSD_File( fileName ).Remove();
  78   }
  79   catch ( Standard_ProgramError ) {
  80     MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
  81   }
  82 }
  83
  84 //=============================================================================
  85 /*!
  86  *
  87  */
  88 //=============================================================================
  89
  90 NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
  91                                           const TopoDS_Shape& aShape,
  92                                           const bool isVolume)
  93   : _mesh    (mesh),
  94     _shape   (aShape),
  95     _isVolume(isVolume),
  96     _optimize(true),
  97     _simpleHyp(NULL)
  98 {
  99   defaultParameters();
 100 }
 101
 102 //================================================================================
 103 /*!
 104  * \brief Initialize global NETGEN parameters with default values
 105  */
 106 //================================================================================
 107
 108 void NETGENPlugin_Mesher::defaultParameters()
 109 {
 110 #ifdef WNT
 111   netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
 112 #else
 113   netgen::MeshingParameters& mparams = netgen::mparam;
 114 #endif
 115   // maximal mesh edge size
 116   mparams.maxh = NETGENPlugin_Hypothesis::GetDefaultMaxSize();
 117   // minimal number of segments per edge
 118   mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
 119   // rate of growth of size between elements
 120   mparams.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
 121   // safety factor for curvatures (elements per radius)
 122   mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
 123   // create elements of second order
 124   mparams.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder() ? 1 : 0;
 125   // quad-dominated surface meshing
 126   if (_isVolume)
 127     mparams.quad = 0;
 128   else
 129     mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0;
 130 }
 131
 132 //=============================================================================
 133 /*!
 134  * Pass parameters to NETGEN
 135  */
 136 //=============================================================================
 137 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
 138 {
 139   if (hyp)
 140   {
 141 #ifdef WNT
 142     netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
 143 #else
 144     netgen::MeshingParameters& mparams = netgen::mparam;
 145 #endif
 146     // Initialize global NETGEN parameters:
 147     // maximal mesh segment size
 148     mparams.maxh = hyp->GetMaxSize();
 149     // minimal number of segments per edge
 150     mparams.segmentsperedge = hyp->GetNbSegPerEdge();
 151     // rate of growth of size between elements
 152     mparams.grading = hyp->GetGrowthRate();
 153     // safety factor for curvatures (elements per radius)
 154     mparams.curvaturesafety = hyp->GetNbSegPerRadius();
 155     // create elements of second order
 156     mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0;
 157     // quad-dominated surface meshing
 158     // only triangles are allowed for volumic mesh
 159     if (!_isVolume)
 160       mparams.quad = static_cast<const NETGENPlugin_Hypothesis_2D*>
 161         (hyp)->GetQuadAllowed() ? 1 : 0;
 162     _optimize = hyp->GetOptimize();
 163     _simpleHyp = NULL;
 164   }
 165 }
 166
 167 //=============================================================================
 168 /*!
 169  * Pass simple parameters to NETGEN
 170  */
 171 //=============================================================================
 172
 173 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp)
 174 {
 175   _simpleHyp = hyp;
 176   if ( _simpleHyp )
 177     defaultParameters();
 178 }
 179
 180 //=============================================================================
 181 /*!
 182  *  Link - a pair of integer numbers
 183  */
 184 //=============================================================================
 185 struct Link
 186 {
 187   int n1, n2;
 188   Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
 189   Link() : n1(0), n2(0) {}
 190 };
 191
 192 int HashCode(const Link& aLink, int aLimit)
 193 {
 194   return HashCode(aLink.n1 + aLink.n2, aLimit);
 195 }
 196
 197 Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
 198 {
 199   return (aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ||
 200           aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1);
 201 }
 202
 203 //================================================================================
 204 /*!
 205  * \brief Initialize netgen::OCCGeometry with OCCT shape
 206  */
 207 //================================================================================
 208
 209 void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry&     occgeo,
 210                                              const TopoDS_Shape&      shape,
 211                                              SMESH_Mesh&              mesh,
 212                                              list< SMESH_subMesh* > * meshedSM)
 213 {
 214   BRepTools::Clean (shape);
 215   try {
 216 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
 217     OCC_CATCH_SIGNALS;
 218 #endif
 219     BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, 0.01, true);
 220   } catch (Standard_Failure) {
 221   }
 222   Bnd_Box bb;
 223   BRepBndLib::Add (shape, bb);
 224   double x1,y1,z1,x2,y2,z2;
 225   bb.Get (x1,y1,z1,x2,y2,z2);
 226   MESSAGE("shape bounding box:\n" <<
 227           "(" << x1 << " " << y1 << " " << z1 << ") " <<
 228           "(" << x2 << " " << y2 << " " << z2 << ")");
 229   netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
 230   netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
 231   occgeo.boundingbox = netgen::Box<3> (p1,p2);
 232
 233   occgeo.shape = shape;
 234   occgeo.changed = 1;
 235   //occgeo.BuildFMap();
 236
 237   // fill maps of shapes of occgeo with not yet meshed subshapes
 238
 239   // get root submeshes
 240   list< SMESH_subMesh* > rootSM;
 241   if ( SMESH_subMesh* sm = mesh.GetSubMeshContaining( shape )) {
 242     rootSM.push_back( sm );
 243   }
 244   else {
 245     for ( TopoDS_Iterator it( shape ); it.More(); it.Next() )
 246       rootSM.push_back( mesh.GetSubMesh( it.Value() ));
 247   }
 248
 249   // add subshapes of empty submeshes
 250   list< SMESH_subMesh* >::iterator rootIt = rootSM.begin(), rootEnd = rootSM.end();
 251   for ( ; rootIt != rootEnd; ++rootIt ) {
 252     SMESH_subMesh * root = *rootIt;
 253     SMESH_subMeshIteratorPtr smIt = root->getDependsOnIterator(/*includeSelf=*/true,
 254                                                                /*complexShapeFirst=*/true);
 255     // to find a right orientation of subshapes (PAL20462)
 256     TopTools_IndexedMapOfShape subShapes;
 257     TopExp::MapShapes(root->GetSubShape(), subShapes);
 258     while ( smIt->more() ) {
 259       SMESH_subMesh* sm = smIt->next();
 260       if ( !meshedSM || sm->IsEmpty() ) {
 261         TopoDS_Shape shape = sm->GetSubShape();
 262         if ( shape.ShapeType() != TopAbs_VERTEX )
 263           shape = subShapes( subShapes.FindIndex( shape ));// - shape->index->oriented shape
 264         switch ( shape.ShapeType() ) {
 265         case TopAbs_FACE  : occgeo.fmap.Add( shape ); break;
 266         case TopAbs_EDGE  : occgeo.emap.Add( shape ); break;
 267         case TopAbs_VERTEX: occgeo.vmap.Add( shape ); break;
 268         case TopAbs_SOLID :occgeo.somap.Add( shape ); break;
 269         default:;
 270         }
 271       }
 272       // collect submeshes of meshed shapes
 273       else if (meshedSM) {
 274         meshedSM->push_back( sm );
 275       }
 276     }
 277   }
 278   occgeo.facemeshstatus.SetSize (occgeo.fmap.Extent());
 279   occgeo.facemeshstatus = 0;
 280
 281 }
 282
 283 //================================================================================
 284 /*!
 285  * \brief return id of netgen point corresponding to SMDS node
 286  */
 287 //================================================================================
 288 typedef map< const SMDS_MeshNode*, int > TNode2IdMap;
 289
 290 static int ngNodeId( const SMDS_MeshNode* node,
 291                      netgen::Mesh&        ngMesh,
 292                      TNode2IdMap&         nodeNgIdMap)
 293 {
 294   int newNgId = ngMesh.GetNP() + 1;
 295
 296   pair< TNode2IdMap::iterator, bool > it_isNew = nodeNgIdMap.insert( make_pair( node, newNgId ));
 297
 298   if ( it_isNew.second ) {
 299     netgen::MeshPoint p( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
 300     ngMesh.AddPoint( p );
 301   }
 302   return it_isNew.first->second;
 303 }
 304
 305 //================================================================================
 306 /*!
 307  * \brief fill ngMesh with nodes and elements of computed submeshes
 308  */
 309 //================================================================================
 310
 311 bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry&           occgeom,
 312                                      netgen::Mesh&                  ngMesh,
 313                                      vector<SMDS_MeshNode*>&        nodeVec,
 314                                      const list< SMESH_subMesh* > & meshedSM)
 315 {
 316   TNode2IdMap nodeNgIdMap;
 317
 318   TopTools_MapOfShape visitedShapes;
 319
 320   SMESH_MesherHelper helper (*_mesh);
 321
 322   int faceID = occgeom.fmap.Extent();
 323   _faceDescriptors.clear();
 324
 325   list< SMESH_subMesh* >::const_iterator smIt, smEnd = meshedSM.end();
 326   for ( smIt = meshedSM.begin(); smIt != smEnd; ++smIt )
 327   {
 328     SMESH_subMesh* sm = *smIt;
 329     if ( !visitedShapes.Add( sm->GetSubShape() ))
 330       continue;
 331
 332     SMESHDS_SubMesh * smDS = sm->GetSubMeshDS();
 333
 334     switch ( sm->GetSubShape().ShapeType() )
 335     {
 336     case TopAbs_EDGE: { // EDGE
 337       // ----------------------
 338       const TopoDS_Edge& geomEdge  = TopoDS::Edge( sm->GetSubShape() );
 339
 340       // Add ng segments for each not meshed face the edge bounds
 341       TopTools_MapOfShape visitedAncestors;
 342       const TopTools_ListOfShape& ancestors = _mesh->GetAncestors( geomEdge );
 343       TopTools_ListIteratorOfListOfShape ancestorIt ( ancestors );
 344       for ( ; ancestorIt.More(); ancestorIt.Next() )
 345       {
 346         const TopoDS_Shape & ans = ancestorIt.Value();
 347         if ( ans.ShapeType() != TopAbs_FACE || !visitedAncestors.Add( ans ))
 348           continue;
 349         const TopoDS_Face& face = TopoDS::Face( ans );
 350
 351         int faceID = occgeom.fmap.FindIndex( face );
 352         if ( faceID < 1 )
 353           continue; // meshed face
 354
 355         // find out orientation of geomEdge within face
 356         bool isForwad = false;
 357         for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() ) {
 358           if ( geomEdge.IsSame( exp.Current() )) {
 359             isForwad = ( exp.Current().Orientation() == geomEdge.Orientation() );
 360             break;
 361           }
 362         }
 363         bool isQuad = smDS->GetElements()->next()->IsQuadratic();
 364
 365         // get all nodes from geomEdge
 366         StdMeshers_FaceSide fSide( face, geomEdge, _mesh, isForwad, isQuad );
 367         const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
 368         int i, nbSeg = fSide.NbSegments();
 369
 370         double otherSeamParam = 0;
 371         helper.SetSubShape( face );
 372         bool isSeam = helper.IsRealSeam( geomEdge );
 373         if ( isSeam )
 374           otherSeamParam =
 375             helper.GetOtherParam( helper.GetPeriodicIndex() == 1 ? points[0].u : points[0].v );
 376
 377         // add segments
 378
 379         int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
 380
 381         for ( i = 0; i < nbSeg; ++i )
 382         {
 383           const UVPtStruct& p1 = points[ i ];
 384           const UVPtStruct& p2 = points[ i+1 ];
 385
 386           netgen::Segment seg;
 387           // ng node ids
 388           seg.p1 = prevNgId;
 389           seg.p2 = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
 390           // node param on curve
 391           seg.epgeominfo[ 0 ].dist = p1.param;
 392           seg.epgeominfo[ 1 ].dist = p2.param;
 393           // uv on face
 394           seg.epgeominfo[ 0 ].u = p1.u;
 395           seg.epgeominfo[ 0 ].v = p1.v;
 396           seg.epgeominfo[ 1 ].u = p2.u;
 397           seg.epgeominfo[ 1 ].v = p2.v;
 398
 399           //seg.epgeominfo[ iEnd ].edgenr = edgeID; //  = geom.emap.FindIndex(edge);
 400           seg.si = faceID;                   // = geom.fmap.FindIndex (face);
 401           seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
 402           ngMesh.AddSegment (seg);
 403
 404           if ( isSeam )
 405           {
 406             if ( helper.GetPeriodicIndex() == 1 ) {
 407               seg.epgeominfo[ 0 ].u = otherSeamParam;
 408               seg.epgeominfo[ 1 ].u = otherSeamParam;
 409               swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
 410             } else {
 411               seg.epgeominfo[ 0 ].v = otherSeamParam;
 412               seg.epgeominfo[ 1 ].v = otherSeamParam;
 413               swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
 414             }
 415             swap (seg.p1, seg.p2);
 416             swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
 417             seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
 418             ngMesh.AddSegment (seg);
 419           }
 420         }
 421       } // loop on geomEdge ancestors
 422
 423       break;
 424     } // case TopAbs_EDGE
 425
 426     case TopAbs_FACE: { // FACE
 427       // ----------------------
 428       const TopoDS_Face& geomFace  = TopoDS::Face( sm->GetSubShape() );
 429       helper.SetSubShape( geomFace );
 430
 431       // Find solids the geomFace bounds
 432       int solidID1 = 0, solidID2 = 0;
 433       const TopTools_ListOfShape& ancestors = _mesh->GetAncestors( geomFace );
 434       TopTools_ListIteratorOfListOfShape ancestorIt ( ancestors );
 435       for ( ; ancestorIt.More(); ancestorIt.Next() )
 436       {
 437         const TopoDS_Shape & solid = ancestorIt.Value();
 438         if ( solid.ShapeType() == TopAbs_SOLID  ) {
 439           int id = occgeom.somap.FindIndex ( solid );
 440           if ( solidID1 && id != solidID1 ) solidID2 = id;
 441           else                              solidID1 = id;
 442         }
 443       }
 444       faceID++;
 445       _faceDescriptors[ faceID ].first  = solidID1;
 446       _faceDescriptors[ faceID ].second = solidID2;
 447
 448       // Orient the face correctly in solidID1 (issue 0020206)
 449       bool reverse = false;
 450       if ( solidID1 ) {
 451         TopoDS_Shape solid = occgeom.somap( solidID1 );
 452         for ( TopExp_Explorer f( solid, TopAbs_FACE ); f.More(); f.Next() ) {
 453           if ( geomFace.IsSame( f.Current() )) {
 454             reverse = SMESH_Algo::IsReversedSubMesh( TopoDS::Face( f.Current()), helper.GetMeshDS() );
 455             break;
 456           }
 457         }
 458       }
 459
 460       // Add surface elements
 461       SMDS_ElemIteratorPtr faces = smDS->GetElements();
 462       while ( faces->more() ) {
 463
 464         const SMDS_MeshElement* f = faces->next();
 465         if ( f->NbNodes() % 3 != 0 ) { // not triangle
 466           for ( ancestorIt.Initialize(ancestors); ancestorIt.More(); ancestorIt.Next() )
 467             if ( ancestorIt.Value().ShapeType() == TopAbs_SOLID  ) {
 468               sm = _mesh->GetSubMesh( ancestorIt.Value() );
 469               break;
 470             }
 471           SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 472           smError.reset( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"Not triangle submesh"));
 473           smError->myBadElements.push_back( f );
 474           return false;
 475         }
 476
 477         netgen::Element2d tri(3);
 478         tri.SetIndex ( faceID );
 479
 480         for ( int i = 0; i < 3; ++i ) {
 481           const SMDS_MeshNode* node = f->GetNode( i ), * inFaceNode=0;
 482           if ( helper.IsSeamShape( node->GetPosition()->GetShapeId() ))
 483             if ( helper.IsSeamShape( f->GetNodeWrap( i+1 )->GetPosition()->GetShapeId() ))
 484               inFaceNode = f->GetNodeWrap( i-1 );
 485             else
 486               inFaceNode = f->GetNodeWrap( i+1 );
 487
 488           gp_XY uv = helper.GetNodeUV( geomFace, node, inFaceNode );
 489           if ( reverse ) {
 490             tri.GeomInfoPi(3-i).u = uv.X();
 491             tri.GeomInfoPi(3-i).v = uv.Y();
 492             tri.PNum      (3-i) = ngNodeId( node, ngMesh, nodeNgIdMap );
 493           } else {
 494             tri.GeomInfoPi(i+1).u = uv.X();
 495             tri.GeomInfoPi(i+1).v = uv.Y();
 496             tri.PNum      (i+1) = ngNodeId( node, ngMesh, nodeNgIdMap );
 497           }
 498         }
 499
 500         ngMesh.AddSurfaceElement (tri);
 501
 502       }
 503       break;
 504     } //
 505
 506     case TopAbs_VERTEX: { // VERTEX
 507       // --------------------------
 508       SMDS_NodeIteratorPtr nodeIt = smDS->GetNodes();
 509       if ( nodeIt->more() )
 510         ngNodeId( nodeIt->next(), ngMesh, nodeNgIdMap );
 511       break;
 512     }
 513     default:;
 514     } // switch
 515   } // loop on submeshes
 516
 517   // fill nodeVec
 518   nodeVec.resize( ngMesh.GetNP() + 1 );
 519   TNode2IdMap::iterator node_NgId, nodeNgIdEnd = nodeNgIdMap.end();
 520   for ( node_NgId = nodeNgIdMap.begin(); node_NgId != nodeNgIdEnd; ++node_NgId)
 521     nodeVec[ node_NgId->second ] = (SMDS_MeshNode*) node_NgId->first;
 522
 523   return true;
 524 }
 525
 526 //=============================================================================
 527 /*!
 528  * Here we are going to use the NETGEN mesher
 529  */
 530 //=============================================================================
 531 bool NETGENPlugin_Mesher::Compute()
 532 {
 533 #ifdef WNT
 534   netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
 535 #else
 536   netgen::MeshingParameters& mparams = netgen::mparam;
 537 #endif
 538   MESSAGE("Compute with:\n"
 539           " max size = " << mparams.maxh << "\n"
 540           " segments per edge = " << mparams.segmentsperedge);
 541   MESSAGE("\n"
 542           " growth rate = " << mparams.grading << "\n"
 543           " elements per radius = " << mparams.curvaturesafety << "\n"
 544           " second order = " << mparams.secondorder << "\n"
 545           " quad allowed = " << mparams.quad);
 546
 547   SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
 548   nglib::Ng_Init();
 549
 550   // -------------------------
 551   // Prepare OCC geometry
 552   // -------------------------
 553
 554   netgen::OCCGeometry occgeo;
 555   list< SMESH_subMesh* > meshedSM;
 556   PrepareOCCgeometry( occgeo, _shape, *_mesh, &meshedSM );
 557
 558   // -------------------------
 559   // Generate the mesh
 560   // -------------------------
 561
 562   netgen::Mesh *ngMesh = NULL;
 563
 564   SMESH_Comment comment;
 565   int err = 0;
 566   int nbInitNod = 0;
 567   int nbInitSeg = 0;
 568   int nbInitFac = 0;
 569   // vector of nodes in which node index == netgen ID
 570   vector< SMDS_MeshNode* > nodeVec;
 571   try
 572   {
 573     // ----------------
 574     // compute 1D mesh
 575     // ----------------
 576     // pass 1D simple parameters to NETGEN
 577     if ( _simpleHyp ) {
 578       if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
 579         // nb of segments
 580         mparams.segmentsperedge = nbSeg + 0.1;
 581         mparams.maxh = occgeo.boundingbox.Diam();
 582         mparams.grading = 0.01;
 583       }
 584       else {
 585         // segment length
 586         mparams.segmentsperedge = 1;
 587         mparams.maxh = _simpleHyp->GetLocalLength();
 588       }
 589     }
 590     // let netgen create ngMesh and calculate element size on not meshed shapes
 591     char *optstr = 0;
 592     int startWith = netgen::MESHCONST_ANALYSE;
 593     int endWith   = netgen::MESHCONST_ANALYSE;
 594     err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
 595     if (err) comment << "Error in netgen::OCCGenerateMesh() at MESHCONST_ANALYSE step";
 596
 597     // fill ngMesh with nodes and elements of computed submeshes
 598     err = ! fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM);
 599     nbInitNod = ngMesh->GetNP();
 600     nbInitSeg = ngMesh->GetNSeg();
 601     nbInitFac = ngMesh->GetNSE();
 602
 603     // compute mesh
 604     if (!err)
 605     {
 606       startWith = endWith = netgen::MESHCONST_MESHEDGES;
 607       err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
 608       if (err) comment << "Error in netgen::OCCGenerateMesh() at 1D mesh generation";
 609     }
 610     // ---------------------
 611     // compute surface mesh
 612     // ---------------------
 613     if (!err)
 614     {
 615       // pass 2D simple parameters to NETGEN
 616       if ( _simpleHyp ) {
 617         if ( double area = _simpleHyp->GetMaxElementArea() ) {
 618           // face area
 619           mparams.maxh = sqrt(2. * area/sqrt(3.0));
 620           mparams.grading = 0.4; // moderate size growth
 621         }
 622         else {
 623           // length from edges
 624           double length = 0;
 625           TopTools_MapOfShape tmpMap;
 626           for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() )
 627             if( tmpMap.Add(exp.Current()) )
 628               length += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() ));
 629
 630           if ( ngMesh->GetNSeg() ) {
 631             // we have to multiply length by 2 since for each TopoDS_Edge there
 632             // are double set of NETGEN edges or, in other words, we have to
 633             // divide ngMesh->GetNSeg() on 2.
 634             mparams.maxh = 2*length / ngMesh->GetNSeg();
 635           }
 636           else
 637             mparams.maxh = 1000;
 638           mparams.grading = 0.2; // slow size growth
 639         }
 640         mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
 641         ngMesh->SetGlobalH (mparams.maxh);
 642         netgen::Box<3> bb = occgeo.GetBoundingBox();
 643         bb.Increase (bb.Diam()/20);
 644         ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading);
 645       }
 646       // let netgen compute 2D mesh
 647       startWith = netgen::MESHCONST_MESHSURFACE;
 648       endWith = _optimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
 649       err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
 650       if (err) comment << "Error in netgen::OCCGenerateMesh() at surface mesh generation";
 651     }
 652     // ---------------------
 653     // generate volume mesh
 654     // ---------------------
 655     if (!err && _isVolume)
 656     {
 657       // add ng face descriptors of meshed faces
 658       std::map< int, std::pair<int,int> >::iterator fId_soIds = _faceDescriptors.begin();
 659       for ( ; fId_soIds != _faceDescriptors.end(); ++fId_soIds ) {
 660         int faceID   = fId_soIds->first;
 661         int solidID1 = fId_soIds->second.first;
 662         int solidID2 = fId_soIds->second.second;
 663         ngMesh->AddFaceDescriptor (netgen::FaceDescriptor(faceID, solidID1, solidID2, 0));
 664       }
 665       // pass 3D simple parameters to NETGEN
 666       const NETGENPlugin_SimpleHypothesis_3D* simple3d =
 667         dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
 668       if ( simple3d ) {
 669         if ( double vol = simple3d->GetMaxElementVolume() ) {
 670           // max volume
 671           mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
 672           mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
 673         }
 674         else {
 675           // length from faces
 676           mparams.maxh = ngMesh->AverageH();
 677         }
 678 //      netgen::ARRAY<double> maxhdom;
 679 //      maxhdom.SetSize (occgeo.NrSolids());
 680 //      maxhdom = mparams.maxh;
 681 //      ngMesh->SetMaxHDomain (maxhdom);
 682         ngMesh->SetGlobalH (mparams.maxh);
 683         mparams.grading = 0.4;
 684         ngMesh->CalcLocalH();
 685       }
 686       // let netgen compute 3D mesh
 687       startWith = netgen::MESHCONST_MESHVOLUME;
 688       endWith = _optimize ? netgen::MESHCONST_OPTVOLUME : netgen::MESHCONST_MESHVOLUME;
 689       err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
 690       if (err) comment << "Error in netgen::OCCGenerateMesh()";
 691     }
 692     if (!err && mparams.secondorder > 0)
 693     {
 694       netgen::OCCRefinementSurfaces ref (occgeo);
 695       ref.MakeSecondOrder (*ngMesh);
 696     }
 697   }
 698   catch (netgen::NgException exc)
 699   {
 700     error->myName = err = COMPERR_ALGO_FAILED;
 701     comment << exc.What();
 702   }
 703
 704   int nbNod = ngMesh->GetNP();
 705   int nbSeg = ngMesh->GetNSeg();
 706   int nbFac = ngMesh->GetNSE();
 707   int nbVol = ngMesh->GetNE();
 708
 709   MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
 710           ", nb nodes: " << nbNod <<
 711           ", nb segments: " << nbSeg <<
 712           ", nb faces: " << nbFac <<
 713           ", nb volumes: " << nbVol);
 714
 715   // -----------------------------------------------------------
 716   // Feed back the SMESHDS with the generated Nodes and Elements
 717   // -----------------------------------------------------------
 718
 719   SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
 720   bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
 721   if ( true /*isOK*/ ) // get whatever built
 722   {
 723     // map of nodes assigned to submeshes
 724     NCollection_Map<int> pindMap;
 725     // create and insert nodes into nodeVec
 726     nodeVec.resize( nbNod + 1 );
 727     int i;
 728     for (i = nbInitNod+1; i <= nbNod /*&& isOK*/; ++i )
 729     {
 730       const netgen::MeshPoint& ngPoint = ngMesh->Point(i);
 731       SMDS_MeshNode* node = NULL;
 732       bool newNodeOnVertex = false;
 733       TopoDS_Vertex aVert;
 734       if (i-nbInitNod <= occgeo.vmap.Extent())
 735       {
 736         // point on vertex
 737         aVert = TopoDS::Vertex(occgeo.vmap(i-nbInitNod));
 738         SMESHDS_SubMesh * submesh = meshDS->MeshElements(aVert);
 739         if (submesh)
 740         {
 741           SMDS_NodeIteratorPtr it = submesh->GetNodes();
 742           if (it->more())
 743           {
 744             node = const_cast<SMDS_MeshNode*> (it->next());
 745             pindMap.Add(i);
 746           }
 747         }
 748         if (!node)
 749           newNodeOnVertex = true;
 750       }
 751       if (!node)
 752         node = meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
 753       if (!node)
 754       {
 755         MESSAGE("Cannot create a mesh node");
 756         if ( !comment.size() ) comment << "Cannot create a mesh node";
 757         nbSeg = nbFac = nbVol = isOK = 0;
 758         break;
 759       }
 760       nodeVec.at(i) = node;
 761       if (newNodeOnVertex)
 762       {
 763         // point on vertex
 764         meshDS->SetNodeOnVertex(node, aVert);
 765         pindMap.Add(i);
 766       }
 767     }
 768
 769     // create mesh segments along geometric edges
 770     NCollection_Map<Link> linkMap;
 771     for (i = nbInitSeg+1; i <= nbSeg/* && isOK*/; ++i )
 772     {
 773       const netgen::Segment& seg = ngMesh->LineSegment(i);
 774       Link link(seg.p1, seg.p2);
 775       if (linkMap.Contains(link))
 776         continue;
 777       linkMap.Add(link);
 778       TopoDS_Edge aEdge;
 779       int pinds[3] = { seg.p1, seg.p2, seg.pmid };
 780       int nbp = 0;
 781       double param2 = 0;
 782       for (int j=0; j < 3; ++j)
 783       {
 784         int pind = pinds[j];
 785         if (pind <= 0) continue;
 786         ++nbp;
 787         double param;
 788         if (j < 2)
 789         {
 790           if (aEdge.IsNull())
 791           {
 792             int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
 793             if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
 794               aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
 795           }
 796           param = seg.epgeominfo[j].dist;
 797           param2 += param;
 798         }
 799         else
 800           param = param2 * 0.5;
 801         if (pind <= nbInitNod || pindMap.Contains(pind))
 802           continue;
 803         if (!aEdge.IsNull())
 804         {
 805           meshDS->SetNodeOnEdge(nodeVec.at(pind), aEdge, param);
 806           pindMap.Add(pind);
 807         }
 808       }
 809       SMDS_MeshEdge* edge;
 810       if (nbp < 3) // second order ?
 811         edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]));
 812       else
 813         edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]),
 814                                 nodeVec.at(pinds[2]));
 815       if (!edge)
 816       {
 817         if ( !comment.size() ) comment << "Cannot create a mesh edge";
 818         MESSAGE("Cannot create a mesh edge");
 819         nbSeg = nbFac = nbVol = isOK = 0;
 820         break;
 821       }
 822       if (!aEdge.IsNull())
 823         meshDS->SetMeshElementOnShape(edge, aEdge);
 824     }
 825
 826     // create mesh faces along geometric faces
 827     for (i = nbInitFac+1; i <= nbFac/* && isOK*/; ++i )
 828     {
 829       const netgen::Element2d& elem = ngMesh->SurfaceElement(i);
 830       int aGeomFaceInd = elem.GetIndex();
 831       TopoDS_Face aFace;
 832       if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
 833         aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
 834       vector<SMDS_MeshNode*> nodes;
 835       for (int j=1; j <= elem.GetNP(); ++j)
 836       {
 837         int pind = elem.PNum(j);
 838         SMDS_MeshNode* node = nodeVec.at(pind);
 839         nodes.push_back(node);
 840         if (pind <= nbInitNod || pindMap.Contains(pind))
 841           continue;
 842         if (!aFace.IsNull())
 843         {
 844           const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
 845           meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
 846           pindMap.Add(pind);
 847         }
 848       }
 849       SMDS_MeshFace* face = NULL;
 850       switch (elem.GetType())
 851       {
 852       case netgen::TRIG:
 853         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
 854         break;
 855       case netgen::QUAD:
 856         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
 857         break;
 858       case netgen::TRIG6:
 859         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
 860         break;
 861       case netgen::QUAD8:
 862         face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
 863                                nodes[4],nodes[7],nodes[5],nodes[6]);
 864         break;
 865       default:
 866         MESSAGE("NETGEN created a face of unexpected type, ignoring");
 867         continue;
 868       }
 869       if (!face)
 870       {
 871         if ( !comment.size() ) comment << "Cannot create a mesh face";
 872         MESSAGE("Cannot create a mesh face");
 873         nbSeg = nbFac = nbVol = isOK = 0;
 874         break;
 875       }
 876       if (!aFace.IsNull())
 877         meshDS->SetMeshElementOnShape(face, aFace);
 878     }
 879
 880     // create tetrahedra
 881     for (i = 1; i <= nbVol/* && isOK*/; ++i)
 882     {
 883       const netgen::Element& elem = ngMesh->VolumeElement(i);
 884       int aSolidInd = elem.GetIndex();
 885       TopoDS_Solid aSolid;
 886       if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
 887         aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
 888       vector<SMDS_MeshNode*> nodes;
 889       for (int j=1; j <= elem.GetNP(); ++j)
 890       {
 891         int pind = elem.PNum(j);
 892         SMDS_MeshNode* node = nodeVec.at(pind);
 893         nodes.push_back(node);
 894         if (pind <= nbInitNod || pindMap.Contains(pind))
 895           continue;
 896         if (!aSolid.IsNull())
 897         {
 898           // point in solid
 899           meshDS->SetNodeInVolume(node, aSolid);
 900           pindMap.Add(pind);
 901         }
 902       }
 903       SMDS_MeshVolume* vol = NULL;
 904       switch (elem.GetType())
 905       {
 906       case netgen::TET:
 907         vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
 908         break;
 909       case netgen::TET10:
 910         vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
 911                                 nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
 912         break;
 913       default:
 914         MESSAGE("NETGEN created a volume of unexpected type, ignoring");
 915         continue;
 916       }
 917       if (!vol)
 918       {
 919         if ( !comment.size() ) comment << "Cannot create a mesh volume";
 920         MESSAGE("Cannot create a mesh volume");
 921         nbSeg = nbFac = nbVol = isOK = 0;
 922         break;
 923       }
 924       if (!aSolid.IsNull())
 925         meshDS->SetMeshElementOnShape(vol, aSolid);
 926     }
 927   }
 928
 929   if ( error->IsOK() && ( !isOK || comment.size() > 0 ))
 930     error->myName = COMPERR_ALGO_FAILED;
 931   if ( !comment.empty() )
 932     error->myComment = comment;
 933
 934   // set bad compute error to subshapes of all failed subshapes shapes
 935   if ( !error->IsOK() && err )
 936   {
 937     for (int i = 1; i <= occgeo.fmap.Extent(); i++) {
 938       int status = occgeo.facemeshstatus[i-1];
 939       if (status == 1 ) continue;
 940       if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.fmap( i ))) {
 941         SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 942         if ( !smError || smError->IsOK() ) {
 943           if ( status == -1 )
 944             smError.reset( new SMESH_ComputeError( error->myName, error->myComment ));
 945           else
 946             smError.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, "Ignored" ));
 947         }
 948       }
 949     }
 950   }
 951
 952   nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
 953   nglib::Ng_Exit();
 954
 955   RemoveTmpFiles();
 956
 957   return error->IsOK();
 958 }
 959
 960 //================================================================================
 961 /*!
 962  * \brief Remove "test.out" and "problemfaces" files in current directory
 963  */
 964 //================================================================================
 965
 966 void NETGENPlugin_Mesher::RemoveTmpFiles()
 967 {
 968   removeFile("test.out");
 969   removeFile("problemfaces");
 970 }
 971
 972
 973 //=============================================================================
 974 /*!
 975  * Evaluate
 976  */
 977 //=============================================================================
 978 bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
 979 {
 980 #ifdef WNT
 981   netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
 982 #else
 983   netgen::MeshingParameters& mparams = netgen::mparam;
 984 #endif
 985
 986
 987   // -------------------------
 988   // Prepare OCC geometry
 989   // -------------------------
 990   netgen::OCCGeometry occgeo;
 991   list< SMESH_subMesh* > meshedSM;
 992   PrepareOCCgeometry( occgeo, _shape, *_mesh, &meshedSM );
 993
 994   bool tooManyElems = false;
 995   const int hugeNb = std::numeric_limits<int>::max() / 100;
 996
 997   // ----------------
 998   // evaluate 1D
 999   // ----------------
1000   // pass 1D simple parameters to NETGEN
1001   int nbs = 0;
1002   if ( _simpleHyp ) {
1003     if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
1004       nbs = nbSeg;
1005       // nb of segments
1006       mparams.segmentsperedge = nbSeg + 0.1;
1007       mparams.maxh = occgeo.boundingbox.Diam();
1008       mparams.grading = 0.01;
1009     }
1010     else {
1011       // segment length
1012       mparams.segmentsperedge = 1;
1013       mparams.maxh = _simpleHyp->GetLocalLength();
1014     }
1015   }
1016   TopTools_DataMapOfShapeInteger EdgesMap;
1017   double fullLen = 0.0;
1018   double fullNbSeg = 0;
1019   for (TopExp_Explorer exp(_shape, TopAbs_EDGE); exp.More(); exp.Next()) {
1020     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
1021     if( EdgesMap.IsBound(E) )
1022       continue;
1023     SMESH_subMesh *sm = _mesh->GetSubMesh(E);
1024     std::vector<int> aVec(SMDSEntity_Last, 0);
1025     double aLen = SMESH_Algo::EdgeLength(E);
1026     fullLen += aLen;
1027     int nb1d = nbs;
1028     tooManyElems = ( aLen/hugeNb > mparams.maxh );
1029     if(nb1d==0 && !tooManyElems) {
1030       nb1d = (int)( aLen/mparams.maxh + 1 );
1031     }
1032     if ( tooManyElems ) // avoid FPE
1033     {
1034       aVec[SMDSEntity_Node] = hugeNb;
1035       aVec[ mparams.secondorder > 0 ? SMDSEntity_Quad_Edge : SMDSEntity_Edge] = hugeNb;
1036     }
1037     else
1038     {
1039       fullNbSeg += nb1d;
1040       if( mparams.secondorder > 0 ) {
1041         aVec[SMDSEntity_Node] = 2*nb1d - 1;
1042         aVec[SMDSEntity_Quad_Edge] = nb1d;
1043       }
1044       else {
1045         aVec[SMDSEntity_Node] = nb1d - 1;
1046         aVec[SMDSEntity_Edge] = nb1d;
1047       }
1048     }
1049     aResMap.insert(std::make_pair(sm,aVec));
1050     EdgesMap.Bind(E,nb1d);
1051   }
1052
1053   // ----------------
1054   // evaluate 2D
1055   // ----------------
1056   if ( _simpleHyp ) {
1057     if ( double area = _simpleHyp->GetMaxElementArea() ) {
1058       // face area
1059       mparams.maxh = sqrt(2. * area/sqrt(3.0));
1060       mparams.grading = 0.4; // moderate size growth
1061     }
1062     else {
1063       // length from edges
1064       mparams.maxh = fullLen/fullNbSeg;
1065       mparams.grading = 0.2; // slow size growth
1066     }
1067     mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
1068   }
1069
1070   for (TopExp_Explorer exp(_shape, TopAbs_FACE); exp.More(); exp.Next())
1071   {
1072     TopoDS_Face F = TopoDS::Face( exp.Current() );
1073     SMESH_subMesh *sm = _mesh->GetSubMesh(F);
1074     GProp_GProps G;
1075     BRepGProp::SurfaceProperties(F,G);
1076     double anArea = G.Mass();
1077     tooManyElems = tooManyElems || ( anArea/hugeNb > mparams.maxh*mparams.maxh );
1078     int nb1d = 0;
1079     if ( !tooManyElems )
1080       for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next())
1081         nb1d += EdgesMap.Find(exp1.Current());
1082
1083     int nbFaces = tooManyElems ? hugeNb : int( 4*anArea / mparams.maxh*mparams.maxh*sqrt(3.));
1084     int nbNodes = tooManyElems ? hugeNb : (( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
1085
1086     std::vector<int> aVec(SMDSEntity_Last, 0);
1087     if( mparams.secondorder > 0 ) {
1088       int nb1d_in = (nbFaces*3 - nb1d) / 2;
1089       aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
1090       aVec[SMDSEntity_Quad_Triangle] = nbFaces;
1091     }
1092     else {
1093       aVec[SMDSEntity_Node] = nbNodes;
1094       aVec[SMDSEntity_Triangle] = nbFaces;
1095     }
1096     aResMap.insert(std::make_pair(sm,aVec));
1097   }
1098
1099   // ----------------
1100   // evaluate 3D
1101   // ----------------
1102   if(_isVolume) {
1103     // pass 3D simple parameters to NETGEN
1104     const NETGENPlugin_SimpleHypothesis_3D* simple3d =
1105       dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
1106     if ( simple3d ) {
1107       if ( double vol = simple3d->GetMaxElementVolume() ) {
1108         // max volume
1109         mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
1110         mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
1111       }
1112       else {
1113         // using previous length from faces
1114       }
1115       mparams.grading = 0.4;
1116     }
1117     GProp_GProps G;
1118     BRepGProp::VolumeProperties(_shape,G);
1119     double aVolume = G.Mass();
1120     double tetrVol = 0.1179*mparams.maxh*mparams.maxh*mparams.maxh;
1121     tooManyElems = tooManyElems || ( aVolume/hugeNb > tetrVol );
1122     int nbVols = tooManyElems ? hugeNb : int(aVolume/tetrVol);
1123     int nb1d_in = int(( nbVols*6 - fullNbSeg ) / 6 );
1124     std::vector<int> aVec(SMDSEntity_Last, 0 );
1125     if ( tooManyElems ) // avoid FPE
1126     {
1127       aVec[SMDSEntity_Node] = hugeNb;
1128       aVec[ mparams.secondorder > 0 ? SMDSEntity_Quad_Tetra : SMDSEntity_Tetra] = hugeNb;
1129     }
1130     else
1131     {
1132       if( mparams.secondorder > 0 ) {
1133         aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in;
1134         aVec[SMDSEntity_Quad_Tetra] = nbVols;
1135       }
1136       else {
1137         aVec[SMDSEntity_Node] = nb1d_in/3 + 1;
1138         aVec[SMDSEntity_Tetra] = nbVols;
1139       }
1140     }
1141     SMESH_subMesh *sm = _mesh->GetSubMesh(_shape);
1142     aResMap.insert(std::make_pair(sm,aVec));
1143   }
1144
1145   return true;
1146 }