src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx

   1 // Copyright (C) 2007-2022  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, or (at your option) any later version.
  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 // TODO: remove use of netgen_param
  36 #include "NETGENPlugin_DriverParam.hxx"
  37
  38 #include <SMDS_MeshElement.hxx>
  39 #include <SMDS_MeshNode.hxx>
  40 #include <SMESHDS_Mesh.hxx>
  41 #include <SMESH_Comment.hxx>
  42 #include <SMESH_ControlsDef.hxx>
  43 #include <SMESH_Gen.hxx>
  44 #include <SMESH_Mesh.hxx>
  45 #include <SMESH_MeshEditor.hxx>
  46 #include <SMESH_MesherHelper.hxx>
  47 #include <SMESH_subMesh.hxx>
  48 #include <StdMeshers_MaxElementVolume.hxx>
  49 #include <StdMeshers_QuadToTriaAdaptor.hxx>
  50 #include <StdMeshers_ViscousLayers.hxx>
  51 #include <SMESH_subMesh.hxx>
  52
  53
  54 #include <BRepGProp.hxx>
  55 #include <BRep_Tool.hxx>
  56 #include <GProp_GProps.hxx>
  57 #include <TopExp.hxx>
  58 #include <TopExp_Explorer.hxx>
  59 #include <TopTools_ListIteratorOfListOfShape.hxx>
  60 #include <TopoDS.hxx>
  61
  62 #include <Standard_Failure.hxx>
  63 #include <Standard_ErrorHandler.hxx>
  64
  65 #include <utilities.h>
  66
  67 #include <list>
  68 #include <vector>
  69 #include <map>
  70
  71 #include <cstdlib>
  72
  73 /*
  74   Netgen include files
  75 */
  76
  77 #ifndef OCCGEOMETRY
  78 #define OCCGEOMETRY
  79 #endif
  80 #include <occgeom.hpp>
  81
  82 #ifdef NETGEN_V5
  83 #include <ngexception.hpp>
  84 #endif
  85 #ifdef NETGEN_V6
  86 #include <core/exception.hpp>
  87 #endif
  88
  89 namespace nglib {
  90 #include <nglib.h>
  91 }
  92 namespace netgen {
  93
  94   NETGENPLUGIN_DLL_HEADER
  95   extern MeshingParameters mparam;
  96
  97   NETGENPLUGIN_DLL_HEADER
  98   extern volatile multithreadt multithread;
  99 }
 100 using namespace nglib;
 101 using namespace std;
 102
 103 //=============================================================================
 104 /*!
 105  *
 106  */
 107 //=============================================================================
 108
 109 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, SMESH_Gen* gen)
 110   : SMESH_3D_Algo(hypId, gen)
 111 {
 112   _name = "NETGEN_3D";
 113   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
 114   _compatibleHypothesis.push_back("MaxElementVolume");
 115   _compatibleHypothesis.push_back("NETGEN_Parameters");
 116   _compatibleHypothesis.push_back("ViscousLayers");
 117
 118   _maxElementVolume = 0.;
 119
 120   _hypMaxElementVolume = NULL;
 121   _hypParameters = NULL;
 122   _viscousLayersHyp = NULL;
 123
 124   _requireShape = false; // can work without shape
 125 }
 126
 127 //=============================================================================
 128 /*!
 129  *
 130  */
 131 //=============================================================================
 132
 133 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
 134 {
 135 }
 136
 137 //=============================================================================
 138 /*!
 139  *
 140  */
 141 //=============================================================================
 142
 143 bool NETGENPlugin_NETGEN_3D::CheckHypothesis (SMESH_Mesh&         aMesh,
 144                                               const TopoDS_Shape& aShape,
 145                                               Hypothesis_Status&  aStatus)
 146 {
 147   _hypMaxElementVolume = NULL;
 148   _hypParameters = NULL;
 149   _viscousLayersHyp = NULL;
 150   _maxElementVolume = DBL_MAX;
 151
 152   // for correct work of GetProgress():
 153   //netgen::multithread.percent = 0.;
 154   //netgen::multithread.task = "Volume meshing";
 155   _progressByTic = -1.;
 156
 157   list<const SMESHDS_Hypothesis*>::const_iterator itl;
 158   //const SMESHDS_Hypothesis* theHyp;
 159
 160   const list<const SMESHDS_Hypothesis*>& hyps =
 161     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
 162   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
 163   if ( h == hyps.end())
 164   {
 165     aStatus = SMESH_Hypothesis::HYP_OK;
 166     return true;  // can work with no hypothesis
 167   }
 168
 169   aStatus = HYP_OK;
 170   for ( ; h != hyps.end(); ++h )
 171   {
 172     if ( !_hypMaxElementVolume )
 173       _hypMaxElementVolume = dynamic_cast< const StdMeshers_MaxElementVolume*> ( *h );
 174     if ( !_viscousLayersHyp ) // several _viscousLayersHyp's allowed
 175       _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
 176     if ( ! _hypParameters )
 177       _hypParameters = dynamic_cast< const NETGENPlugin_Hypothesis*> ( *h );
 178
 179     if ( *h != _hypMaxElementVolume &&
 180          *h != _viscousLayersHyp &&
 181          *h != _hypParameters &&
 182          !dynamic_cast< const StdMeshers_ViscousLayers*>(*h)) // several VL hyps allowed
 183       aStatus = HYP_INCOMPATIBLE;
 184   }
 185   if ( _hypMaxElementVolume && _hypParameters )
 186     aStatus = HYP_INCOMPATIBLE;
 187   else if ( aStatus == HYP_OK && _viscousLayersHyp )
 188     error( _viscousLayersHyp->CheckHypothesis( aMesh, aShape, aStatus ));
 189
 190   if ( _hypMaxElementVolume )
 191     _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
 192
 193   return aStatus == HYP_OK;
 194 }
 195
 196
 197
 198 //=============================================================================
 199 /*!
 200  *Here we are going to use the NETGEN mesher
 201  */
 202 //=============================================================================
 203
 204 /**
 205  * @brief Get an iterator on the Surface element with their orientation
 206  *
 207  */
 208
 209 bool NETGENPlugin_NETGEN_3D::getSurfaceElements(
 210     SMESH_Mesh&         aMesh,
 211     const TopoDS_Shape& aShape,
 212     SMESH_ProxyMesh::Ptr proxyMesh,
 213     NETGENPlugin_Internals &internals,
 214     SMESH_MesherHelper &helper,
 215     std::map<const SMDS_MeshElement*, tuple<bool, bool>>& listElements
 216 )
 217 {
 218   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 219   TopAbs_ShapeEnum mainType = aMesh.GetShapeToMesh().ShapeType();
 220   bool checkReverse = ( mainType == TopAbs_COMPOUND || mainType == TopAbs_COMPSOLID );
 221
 222   for ( TopExp_Explorer exFa( aShape, TopAbs_FACE ); exFa.More(); exFa.Next())
 223   {
 224     const TopoDS_Shape& aShapeFace = exFa.Current();
 225     int faceID = meshDS->ShapeToIndex( aShapeFace );
 226     bool isInternalFace = internals.isInternalShape( faceID );
 227     bool isRev = false;
 228     if ( checkReverse && !isInternalFace &&
 229           helper.NbAncestors(aShapeFace, aMesh, aShape.ShapeType()) > 1 )
 230       // IsReversedSubMesh() can work wrong on strongly curved faces,
 231       // so we use it as less as possible
 232       isRev = helper.IsReversedSubMesh( TopoDS::Face( aShapeFace ));
 233
 234     const SMESHDS_SubMesh * aSubMeshDSFace = proxyMesh->GetSubMesh( aShapeFace );
 235     if ( !aSubMeshDSFace ) continue;
 236
 237     SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
 238     if ( _quadraticMesh &&
 239           dynamic_cast< const SMESH_ProxyMesh::SubMesh*>( aSubMeshDSFace ))
 240     {
 241       // add medium nodes of proxy triangles to helper (#16843)
 242       while ( iteratorElem->more() )
 243         helper.AddTLinks( static_cast< const SMDS_MeshFace* >( iteratorElem->next() ));
 244
 245       iteratorElem = aSubMeshDSFace->GetElements();
 246     }
 247     while(iteratorElem->more()){
 248       const SMDS_MeshElement* elem = iteratorElem->next();
 249       // check mesh face
 250       if ( !elem ){
 251         return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 252       }
 253       if ( elem->NbCornerNodes() != 3 ){
 254         return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
 255       }
 256       listElements[elem] = tuple(isRev, isInternalFace);
 257     }
 258   }
 259
 260   return false;
 261 }
 262
 263
 264 bool NETGENPlugin_NETGEN_3D::computeFillNgMesh(
 265   SMESH_Mesh&         aMesh,
 266   const TopoDS_Shape& aShape,
 267   vector< const SMDS_MeshNode* > &nodeVec,
 268   NETGENPlugin_NetgenLibWrapper &ngLib,
 269   SMESH_MesherHelper &helper,
 270   int &Netgen_NbOfNodes)
 271 {
 272   netgen::multithread.terminate = 0;
 273   netgen::multithread.task = "Volume meshing";
 274   _progressByTic = -1.;
 275
 276   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 277
 278   _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
 279   helper.SetElementsOnShape( true );
 280
 281   Netgen_NbOfNodes = 0;
 282   double Netgen_point[3];
 283   int Netgen_triangle[3];
 284
 285   Ng_Mesh * Netgen_mesh = (Ng_Mesh*)ngLib._ngMesh;
 286
 287   {
 288     const int invalid_ID = -1;
 289
 290     SMESH::Controls::Area areaControl;
 291     SMESH::Controls::TSequenceOfXYZ nodesCoords;
 292
 293     // maps nodes to ng ID
 294     typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 295     typedef TNodeToIDMap::value_type                     TN2ID;
 296     TNodeToIDMap nodeToNetgenID;
 297
 298     // find internal shapes
 299     NETGENPlugin_Internals internals( aMesh, aShape, /*is3D=*/true );
 300
 301     // ---------------------------------
 302     // Feed the Netgen with surface mesh
 303     // ---------------------------------
 304     bool isRev=false;
 305     bool isInternalFace=false;
 306
 307     SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
 308     if ( _viscousLayersHyp )
 309     {
 310       netgen::multithread.percent = 3;
 311       proxyMesh = _viscousLayersHyp->Compute( aMesh, aShape );
 312       if ( !proxyMesh )
 313         return false;
 314     }
 315     if ( aMesh.NbQuadrangles() > 0 )
 316     {
 317       netgen::multithread.percent = 6;
 318       StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
 319       Adaptor->Compute(aMesh,aShape,proxyMesh.get());
 320       proxyMesh.reset( Adaptor );
 321     }
 322
 323     std::map<const SMDS_MeshElement*, tuple<bool, bool>> listElements;
 324     bool ret = getSurfaceElements(aMesh, aShape, proxyMesh, internals, helper, listElements);
 325     if(ret)
 326       return ret;
 327
 328     for ( auto const& [elem, info] : listElements ) // loop on elements on a geom face
 329     {
 330       isRev = get<0>(info);
 331       isInternalFace = get<1>(info);
 332       // Add nodes of triangles and triangles them-selves to netgen mesh
 333
 334       // add three nodes of triangle
 335       bool hasDegen = false;
 336       for ( int iN = 0; iN < 3; ++iN )
 337       {
 338         const SMDS_MeshNode* node = elem->GetNode( iN );
 339         const int shapeID = node->getshapeId();
 340         if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
 341               helper.IsDegenShape( shapeID ))
 342         {
 343           // ignore all nodes on degeneraged edge and use node on its vertex instead
 344           TopoDS_Shape vertex = TopoDS_Iterator( meshDS->IndexToShape( shapeID )).Value();
 345           node = SMESH_Algo::VertexNode( TopoDS::Vertex( vertex ), meshDS );
 346           hasDegen = true;
 347         }
 348         int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
 349         if ( ngID == invalid_ID )
 350         {
 351           ngID = ++Netgen_NbOfNodes;
 352           Netgen_point [ 0 ] = node->X();
 353           Netgen_point [ 1 ] = node->Y();
 354           Netgen_point [ 2 ] = node->Z();
 355           Ng_AddPoint(Netgen_mesh, Netgen_point);
 356         }
 357         Netgen_triangle[ isRev ? 2-iN : iN ] = ngID;
 358       }
 359       // add triangle
 360       if ( hasDegen && (Netgen_triangle[0] == Netgen_triangle[1] ||
 361                         Netgen_triangle[0] == Netgen_triangle[2] ||
 362                         Netgen_triangle[2] == Netgen_triangle[1] ))
 363         continue;
 364
 365       Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 366
 367       if ( isInternalFace && !proxyMesh->IsTemporary( elem ))
 368       {
 369         swap( Netgen_triangle[1], Netgen_triangle[2] );
 370         Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 371       }
 372     } // loop on elements on a face
 373
 374     // insert old nodes into nodeVec
 375     nodeVec.resize( nodeToNetgenID.size() + 1, 0 );
 376     TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 377     for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 378       nodeVec[ n_id->second ] = n_id->first;
 379     nodeToNetgenID.clear();
 380
 381     if ( internals.hasInternalVertexInSolid() )
 382     {
 383       netgen::OCCGeometry occgeo;
 384       NETGENPlugin_Mesher::AddIntVerticesInSolids( occgeo,
 385                                                    (netgen::Mesh&) *Netgen_mesh,
 386                                                    nodeVec,
 387                                                    internals);
 388     }
 389   }
 390   Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
 391   return false;
 392 }
 393
 394 bool NETGENPlugin_NETGEN_3D::computePrepareParam(
 395   SMESH_Mesh&         aMesh,
 396   NETGENPlugin_NetgenLibWrapper &ngLib,
 397   netgen::OCCGeometry &occgeo,
 398   SMESH_MesherHelper &helper,
 399   int &endWith)
 400
 401 {
 402   netgen::multithread.terminate = 0;
 403
 404   netgen::Mesh* ngMesh = ngLib._ngMesh;
 405
 406   NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
 407
 408
 409   if ( _hypParameters )
 410   {
 411     aMesher.SetParameters( _hypParameters );
 412
 413     if ( !_hypParameters->GetLocalSizesAndEntries().empty() ||
 414          !_hypParameters->GetMeshSizeFile().empty() )
 415     {
 416       if ( ! &ngMesh->LocalHFunction() )
 417       {
 418         netgen::Point3d pmin, pmax;
 419         ngMesh->GetBox( pmin, pmax, 0 );
 420         ngMesh->SetLocalH( pmin, pmax, _hypParameters->GetGrowthRate() );
 421       }
 422       aMesher.SetLocalSize( occgeo, *ngMesh );
 423
 424       try {
 425         ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
 426       } catch (netgen::NgException & ex) {
 427         return error( COMPERR_BAD_PARMETERS, ex.What() );
 428       }
 429     }
 430     if ( !_hypParameters->GetOptimize() )
 431       endWith = netgen::MESHCONST_MESHVOLUME;
 432   }
 433   else if ( _hypMaxElementVolume )
 434   {
 435     netgen::mparam.maxh = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
 436     // limitVolumeSize( ngMesh, mparam.maxh ); // result is unpredictable
 437   }
 438   else if ( aMesh.HasShapeToMesh() )
 439   {
 440     aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh );
 441     netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2;
 442   }
 443   else
 444   {
 445     netgen::Point3d pmin, pmax;
 446     ngMesh->GetBox (pmin, pmax);
 447     netgen::mparam.maxh = Dist(pmin, pmax)/2;
 448   }
 449
 450   if ( !_hypParameters && aMesh.HasShapeToMesh() )
 451   {
 452     netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh );
 453   }
 454   return false;
 455 }
 456
 457 bool NETGENPlugin_NETGEN_3D::computeRunMesher(
 458   netgen::OCCGeometry &occgeo,
 459   vector< const SMDS_MeshNode* > &nodeVec,
 460   netgen::Mesh* ngMesh,
 461   NETGENPlugin_NetgenLibWrapper &ngLib,
 462   int &startWith, int &endWith)
 463 {
 464   int err = 1;
 465
 466   try
 467   {
 468     OCC_CATCH_SIGNALS;
 469
 470     ngLib.CalcLocalH(ngMesh);
 471     err = ngLib.GenerateMesh(occgeo, startWith, endWith);
 472
 473     if(netgen::multithread.terminate)
 474       return false;
 475     if ( err ){
 476       error(SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task);
 477     }
 478   }
 479   catch (Standard_Failure& ex)
 480   {
 481     SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
 482     str << " at " << netgen::multithread.task
 483         << ": " << ex.DynamicType()->Name();
 484     if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
 485       str << ": " << ex.GetMessageString();
 486     error(str);
 487   }
 488   catch (netgen::NgException& exc)
 489   {
 490     SMESH_Comment str("NgException");
 491     if ( strlen( netgen::multithread.task ) > 0 )
 492       str << " at " << netgen::multithread.task;
 493     str << ": " << exc.What();
 494     error(str);
 495   }
 496   catch (...)
 497   {
 498     SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
 499     if ( strlen( netgen::multithread.task ) > 0 )
 500       str << " at " << netgen::multithread.task;
 501     error(str);
 502   }
 503
 504   if ( err )
 505   {
 506     SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
 507     if ( ce && ce->HasBadElems() ){
 508       error( ce );
 509     }
 510   }
 511
 512   return false;
 513 }
 514
 515 bool NETGENPlugin_NETGEN_3D::computeFillMesh(
 516   vector< const SMDS_MeshNode* > &nodeVec,
 517   NETGENPlugin_NetgenLibWrapper &ngLib,
 518   SMESH_MesherHelper &helper,
 519   int &Netgen_NbOfNodes
 520   )
 521 {
 522   Ng_Mesh* Netgen_mesh = ngLib.ngMesh();
 523
 524   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 525   int Netgen_NbOfTetra    = Ng_GetNE(Netgen_mesh);
 526
 527   bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
 528   if ( isOK )
 529   {
 530     double Netgen_point[3];
 531     int    Netgen_tetrahedron[4];
 532
 533     // create and insert new nodes into nodeVec
 534     nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 );
 535     int nodeIndex = Netgen_NbOfNodes + 1;
 536     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 537     {
 538       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 539       nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], Netgen_point[1], Netgen_point[2]);
 540     }
 541
 542     // create tetrahedrons
 543     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 544     {
 545       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 546       try
 547       {
 548         helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
 549                           nodeVec.at( Netgen_tetrahedron[1] ),
 550                           nodeVec.at( Netgen_tetrahedron[2] ),
 551                           nodeVec.at( Netgen_tetrahedron[3] ));
 552       }
 553       catch (...)
 554       {
 555       }
 556     }
 557   }
 558   return false;
 559 }
 560
 561 bool NETGENPlugin_NETGEN_3D::Compute(
 562   SMESH_Mesh&         aMesh,
 563   const TopoDS_Shape& aShape)
 564 {
 565   // vector of nodes in which node index == netgen ID
 566   vector< const SMDS_MeshNode* > nodeVec;
 567   NETGENPlugin_NetgenLibWrapper ngLib;
 568   SMESH_MesherHelper helper(aMesh);
 569   int startWith = netgen::MESHCONST_MESHVOLUME;
 570   int endWith   = netgen::MESHCONST_OPTVOLUME;
 571   int Netgen_NbOfNodes;
 572
 573   computeFillNgMesh(aMesh, aShape, nodeVec, ngLib, helper, Netgen_NbOfNodes);
 574
 575   netgen::OCCGeometry occgeo;
 576   computePrepareParam(aMesh, ngLib, occgeo, helper, endWith);
 577   computeRunMesher(occgeo, nodeVec, ngLib._ngMesh, ngLib, startWith, endWith);
 578
 579   computeFillMesh(nodeVec, ngLib, helper, Netgen_NbOfNodes);
 580
 581   return false;
 582
 583 }
 584
 585 //================================================================================
 586 /*!
 587  * \brief set parameters and generate the volume mesh
 588  */
 589 //================================================================================
 590
 591 bool NETGENPlugin_NETGEN_3D::compute(SMESH_Mesh&                     aMesh,
 592                                      SMESH_MesherHelper&             helper,
 593                                      vector< const SMDS_MeshNode* >& nodeVec,
 594                                      NETGENPlugin_NetgenLibWrapper&  ngLib)
 595 {
 596   netgen::multithread.terminate = 0;
 597
 598   netgen::Mesh* ngMesh = ngLib._ngMesh;
 599   Ng_Mesh* Netgen_mesh = ngLib.ngMesh();
 600   int Netgen_NbOfNodes = Ng_GetNP( Netgen_mesh );
 601
 602   int startWith = netgen::MESHCONST_MESHVOLUME;
 603   int endWith   = netgen::MESHCONST_OPTVOLUME;
 604   int err = 1;
 605
 606   NETGENPlugin_Mesher aMesher( &aMesh, helper.GetSubShape(), /*isVolume=*/true );
 607   netgen::OCCGeometry occgeo;
 608
 609   if ( _hypParameters )
 610   {
 611     aMesher.SetParameters( _hypParameters );
 612
 613     if ( !_hypParameters->GetLocalSizesAndEntries().empty() ||
 614          !_hypParameters->GetMeshSizeFile().empty() )
 615     {
 616       if ( ! &ngMesh->LocalHFunction() )
 617       {
 618         netgen::Point3d pmin, pmax;
 619         ngMesh->GetBox( pmin, pmax, 0 );
 620         ngMesh->SetLocalH( pmin, pmax, _hypParameters->GetGrowthRate() );
 621       }
 622       aMesher.SetLocalSize( occgeo, *ngMesh );
 623
 624       try {
 625         ngMesh->LoadLocalMeshSize( netgen::mparam.meshsizefilename );
 626       } catch (netgen::NgException & ex) {
 627         return error( COMPERR_BAD_PARMETERS, ex.What() );
 628       }
 629     }
 630     if ( !_hypParameters->GetOptimize() )
 631       endWith = netgen::MESHCONST_MESHVOLUME;
 632   }
 633   else if ( _hypMaxElementVolume )
 634   {
 635     netgen::mparam.maxh = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
 636     // limitVolumeSize( ngMesh, mparam.maxh ); // result is unpredictable
 637   }
 638   else if ( aMesh.HasShapeToMesh() )
 639   {
 640     aMesher.PrepareOCCgeometry( occgeo, helper.GetSubShape(), aMesh );
 641     netgen::mparam.maxh = occgeo.GetBoundingBox().Diam()/2;
 642   }
 643   else
 644   {
 645     netgen::Point3d pmin, pmax;
 646     ngMesh->GetBox (pmin, pmax);
 647     netgen::mparam.maxh = Dist(pmin, pmax)/2;
 648   }
 649
 650   if ( !_hypParameters && aMesh.HasShapeToMesh() )
 651   {
 652     netgen::mparam.minh = aMesher.GetDefaultMinSize( helper.GetSubShape(), netgen::mparam.maxh );
 653   }
 654
 655   try
 656   {
 657     OCC_CATCH_SIGNALS;
 658
 659     ngLib.CalcLocalH(ngMesh);
 660     err = ngLib.GenerateMesh(occgeo, startWith, endWith);
 661
 662     if(netgen::multithread.terminate)
 663       return false;
 664     if ( err )
 665       error(SMESH_Comment("Error in netgen::OCCGenerateMesh() at ") << netgen::multithread.task);
 666   }
 667   catch (Standard_Failure& ex)
 668   {
 669     SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
 670     str << " at " << netgen::multithread.task
 671         << ": " << ex.DynamicType()->Name();
 672     if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
 673       str << ": " << ex.GetMessageString();
 674     error(str);
 675   }
 676   catch (netgen::NgException& exc)
 677   {
 678     SMESH_Comment str("NgException");
 679     if ( strlen( netgen::multithread.task ) > 0 )
 680       str << " at " << netgen::multithread.task;
 681     str << ": " << exc.What();
 682     error(str);
 683   }
 684   catch (...)
 685   {
 686     SMESH_Comment str("Exception in  netgen::OCCGenerateMesh()");
 687     if ( strlen( netgen::multithread.task ) > 0 )
 688       str << " at " << netgen::multithread.task;
 689     error(str);
 690   }
 691
 692   int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 693   int Netgen_NbOfTetra    = Ng_GetNE(Netgen_mesh);
 694
 695   // -------------------------------------------------------------------
 696   // Feed back the SMESHDS with the generated Nodes and Volume Elements
 697   // -------------------------------------------------------------------
 698
 699   if ( err )
 700   {
 701     SMESH_ComputeErrorPtr ce = NETGENPlugin_Mesher::ReadErrors(nodeVec);
 702     if ( ce && ce->HasBadElems() )
 703       error( ce );
 704   }
 705
 706   bool isOK = ( /*status == NG_OK &&*/ Netgen_NbOfTetra > 0 );// get whatever built
 707   if ( isOK )
 708   {
 709     double Netgen_point[3];
 710     int    Netgen_tetrahedron[4];
 711
 712     // create and insert new nodes into nodeVec
 713     nodeVec.resize( Netgen_NbOfNodesNew + 1, 0 );
 714     int nodeIndex = Netgen_NbOfNodes + 1;
 715     for ( ; nodeIndex <= Netgen_NbOfNodesNew; ++nodeIndex )
 716     {
 717       Ng_GetPoint( Netgen_mesh, nodeIndex, Netgen_point );
 718       nodeVec.at(nodeIndex) = helper.AddNode(Netgen_point[0], Netgen_point[1], Netgen_point[2]);
 719     }
 720
 721     // create tetrahedrons
 722     for ( int elemIndex = 1; elemIndex <= Netgen_NbOfTetra; ++elemIndex )
 723     {
 724       Ng_GetVolumeElement(Netgen_mesh, elemIndex, Netgen_tetrahedron);
 725       try
 726       {
 727         helper.AddVolume (nodeVec.at( Netgen_tetrahedron[0] ),
 728                           nodeVec.at( Netgen_tetrahedron[1] ),
 729                           nodeVec.at( Netgen_tetrahedron[2] ),
 730                           nodeVec.at( Netgen_tetrahedron[3] ));
 731       }
 732       catch (...)
 733       {
 734       }
 735     }
 736   }
 737
 738   return !err;
 739 }
 740
 741 //================================================================================
 742 /*!
 743  * \brief Compute tetrahedral mesh from 2D mesh without geometry
 744  */
 745 //================================================================================
 746
 747 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh&         aMesh,
 748                                      SMESH_MesherHelper* aHelper)
 749 {
 750   const int invalid_ID = -1;
 751
 752   netgen::multithread.terminate = 0;
 753   _progressByTic = -1.;
 754
 755   SMESH_MesherHelper::MType MeshType = aHelper->IsQuadraticMesh();
 756   if ( MeshType == SMESH_MesherHelper::COMP )
 757     return error( COMPERR_BAD_INPUT_MESH,
 758                   SMESH_Comment("Mesh with linear and quadratic elements given"));
 759
 760   aHelper->SetIsQuadratic( MeshType == SMESH_MesherHelper::QUADRATIC );
 761
 762   // ---------------------------------
 763   // Feed the Netgen with surface mesh
 764   // ---------------------------------
 765
 766   int Netgen_NbOfNodes = 0;
 767   double Netgen_point[3];
 768   int Netgen_triangle[3];
 769
 770   NETGENPlugin_NetgenLibWrapper ngLib;
 771   Ng_Mesh * Netgen_mesh = ngLib.ngMesh();
 772
 773   SMESH_ProxyMesh::Ptr proxyMesh( new SMESH_ProxyMesh( aMesh ));
 774   if ( aMesh.NbQuadrangles() > 0 )
 775   {
 776     StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
 777     Adaptor->Compute(aMesh);
 778     proxyMesh.reset( Adaptor );
 779
 780     if ( aHelper->IsQuadraticMesh() )
 781     {
 782       SMDS_ElemIteratorPtr fIt = proxyMesh->GetFaces();
 783       while( fIt->more())
 784         aHelper->AddTLinks( static_cast< const SMDS_MeshFace* >( fIt->next() ));
 785     }
 786   }
 787
 788   // maps nodes to ng ID
 789   typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
 790   typedef TNodeToIDMap::value_type                     TN2ID;
 791   TNodeToIDMap nodeToNetgenID;
 792
 793   SMDS_ElemIteratorPtr fIt = proxyMesh->GetFaces();
 794   while( fIt->more())
 795   {
 796     // check element
 797     const SMDS_MeshElement* elem = fIt->next();
 798     if ( !elem )
 799       return error( COMPERR_BAD_INPUT_MESH, "Null element encounters");
 800     if ( elem->NbCornerNodes() != 3 )
 801       return error( COMPERR_BAD_INPUT_MESH, "Not triangle element encounters");
 802
 803     // add three nodes of triangle
 804     for ( int iN = 0; iN < 3; ++iN )
 805     {
 806       const SMDS_MeshNode* node = elem->GetNode( iN );
 807       int& ngID = nodeToNetgenID.insert(TN2ID( node, invalid_ID )).first->second;
 808       if ( ngID == invalid_ID )
 809       {
 810         ngID = ++Netgen_NbOfNodes;
 811         Netgen_point [ 0 ] = node->X();
 812         Netgen_point [ 1 ] = node->Y();
 813         Netgen_point [ 2 ] = node->Z();
 814         Ng_AddPoint(Netgen_mesh, Netgen_point);
 815       }
 816       Netgen_triangle[ iN ] = ngID;
 817     }
 818     Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 819   }
 820   proxyMesh.reset(); // delete tmp faces
 821
 822   // vector of nodes in which node index == netgen ID
 823   vector< const SMDS_MeshNode* > nodeVec ( nodeToNetgenID.size() + 1 );
 824   // insert old nodes into nodeVec
 825   TNodeToIDMap::iterator n_id = nodeToNetgenID.begin();
 826   for ( ; n_id != nodeToNetgenID.end(); ++n_id )
 827     nodeVec.at( n_id->second ) = n_id->first;
 828   nodeToNetgenID.clear();
 829
 830   // -------------------------
 831   // Generate the volume mesh
 832   // -------------------------
 833
 834   return ( ngLib._isComputeOk = compute( aMesh, *aHelper, nodeVec, ngLib ));
 835 }
 836
 837 void NETGENPlugin_NETGEN_3D::CancelCompute()
 838 {
 839   SMESH_Algo::CancelCompute();
 840   netgen::multithread.terminate = 1;
 841 }
 842
 843 //================================================================================
 844 /*!
 845  * \brief Return Compute progress
 846  */
 847 //================================================================================
 848
 849 double NETGENPlugin_NETGEN_3D::GetProgress() const
 850 {
 851   double res;
 852   const char* volMeshing = "Volume meshing";
 853   const char* dlnMeshing = "Delaunay meshing";
 854   const double meshingRatio = 0.15;
 855   const_cast<NETGENPlugin_NETGEN_3D*>( this )->_progressTic++;
 856
 857   if ( _progressByTic < 0. &&
 858        ( strncmp( netgen::multithread.task, dlnMeshing, 3 ) == 0 ||
 859          strncmp( netgen::multithread.task, volMeshing, 3 ) == 0 ))
 860   {
 861     res = 0.001 + meshingRatio * netgen::multithread.percent / 100.;
 862   }
 863   else // different otimizations
 864   {
 865     if ( _progressByTic < 0. )
 866       ((NETGENPlugin_NETGEN_3D*)this)->_progressByTic = meshingRatio / _progressTic;
 867     res = _progressByTic * _progressTic;
 868   }
 869   return Min ( res, 0.98 );
 870 }
 871
 872 //=============================================================================
 873 /*!
 874  *
 875  */
 876 //=============================================================================
 877
 878 bool NETGENPlugin_NETGEN_3D::Evaluate(SMESH_Mesh& aMesh,
 879                                       const TopoDS_Shape& aShape,
 880                                       MapShapeNbElems& aResMap)
 881 {
 882   smIdType nbtri = 0, nbqua = 0;
 883   double fullArea = 0.0;
 884   for (TopExp_Explorer expF(aShape, TopAbs_FACE); expF.More(); expF.Next()) {
 885     TopoDS_Face F = TopoDS::Face( expF.Current() );
 886     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
 887     MapShapeNbElemsItr anIt = aResMap.find(sm);
 888     if( anIt==aResMap.end() ) {
 889       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 890       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 891       return false;
 892     }
 893     std::vector<smIdType> aVec = (*anIt).second;
 894     nbtri += std::max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
 895     nbqua += std::max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 896     GProp_GProps G;
 897     BRepGProp::SurfaceProperties(F,G);
 898     double anArea = G.Mass();
 899     fullArea += anArea;
 900   }
 901
 902   // collect info from edges
 903   smIdType nb0d_e = 0, nb1d_e = 0;
 904   bool IsQuadratic = false;
 905   bool IsFirst = true;
 906   TopTools_MapOfShape tmpMap;
 907   for (TopExp_Explorer expF(aShape, TopAbs_EDGE); expF.More(); expF.Next()) {
 908     TopoDS_Edge E = TopoDS::Edge(expF.Current());
 909     if( tmpMap.Contains(E) )
 910       continue;
 911     tmpMap.Add(E);
 912     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(expF.Current());
 913     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
 914     if( anIt==aResMap.end() ) {
 915       SMESH_ComputeErrorPtr& smError = aSubMesh->GetComputeError();
 916       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
 917                                             "Submesh can not be evaluated",this));
 918       return false;
 919     }
 920     std::vector<smIdType> aVec = (*anIt).second;
 921     nb0d_e += aVec[SMDSEntity_Node];
 922     nb1d_e += std::max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
 923     if(IsFirst) {
 924       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 925       IsFirst = false;
 926     }
 927   }
 928   tmpMap.Clear();
 929
 930   double ELen_face = sqrt(2.* ( fullArea/double(nbtri+nbqua*2) ) / sqrt(3.0) );
 931   double ELen_vol = pow( 72, 1/6. ) * pow( _maxElementVolume, 1/3. );
 932   double ELen = Min(ELen_vol,ELen_face*2);
 933
 934   GProp_GProps G;
 935   BRepGProp::VolumeProperties(aShape,G);
 936   double aVolume = G.Mass();
 937   double tetrVol = 0.1179*ELen*ELen*ELen;
 938   double CoeffQuality = 0.9;
 939   smIdType nbVols = (smIdType)( aVolume/tetrVol/CoeffQuality );
 940   smIdType nb1d_f = (nbtri*3 + nbqua*4 - nb1d_e) / 2;
 941   smIdType nb1d_in = (nbVols*6 - nb1d_e - nb1d_f ) / 5;
 942   std::vector<smIdType> aVec(SMDSEntity_Last);
 943   for(smIdType i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
 944   if( IsQuadratic ) {
 945     aVec[SMDSEntity_Node] = nb1d_in/6 + 1 + nb1d_in;
 946     aVec[SMDSEntity_Quad_Tetra] = nbVols - nbqua*2;
 947     aVec[SMDSEntity_Quad_Pyramid] = nbqua;
 948   }
 949   else {
 950     aVec[SMDSEntity_Node] = nb1d_in/6 + 1;
 951     aVec[SMDSEntity_Tetra] = nbVols - nbqua*2;
 952     aVec[SMDSEntity_Pyramid] = nbqua;
 953   }
 954   SMESH_subMesh *sm = aMesh.GetSubMesh(aShape);
 955   aResMap.insert(std::make_pair(sm,aVec));
 956
 957   return true;
 958 }
 959
 960