src/NETGENPlugin/NETGENPlugin_Remesher_2D.cxx

   1 // Copyright (C) 2007-2016  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 // File      : NETGENPlugin_Remesher_2D.cxx
  23 // Created   : Thu Sep 21 16:48:46 2017
  24 // Author    : Edward AGAPOV (eap)
  25 //
  26
  27 #include "NETGENPlugin_Remesher_2D.hxx"
  28
  29 #include "NETGENPlugin_Mesher.hxx"
  30 #include "NETGENPlugin_Hypothesis_2D.hxx"
  31
  32 #include <SMDS_SetIterator.hxx>
  33 #include <SMESHDS_Mesh.hxx>
  34 #include <SMESH_ControlsDef.hxx>
  35 #include <SMESH_Gen.hxx>
  36 #include <SMESH_MeshAlgos.hxx>
  37 #include <SMESH_MesherHelper.hxx>
  38 #include <SMESH_subMesh.hxx>
  39
  40 #include <Bnd_B3d.hxx>
  41 #include <Precision.hxx>
  42
  43 #include <occgeom.hpp>
  44 #include <meshing.hpp>
  45 #include <stlgeom.hpp>
  46 //#include <stltool.hpp>
  47
  48 #include <boost/container/flat_set.hpp>
  49
  50 using namespace nglib;
  51
  52 namespace netgen {
  53
  54 #if defined(NETGEN_V5) && defined(WIN32)
  55    DLL_HEADER
  56  #endif
  57    extern STLParameters stlparam;
  58 }
  59 namespace nglib
  60 {
  61 #if defined(NETGEN_V5) && defined(WIN32)
  62         DLL_HEADER
  63 #endif
  64   extern netgen::Array<netgen::Point<3> > readedges;
  65 }
  66
  67 namespace
  68 {
  69   //=============================================================================
  70   /*!
  71    * \brief Fill holes in the mesh, since netgen can remesh only a closed shell mesh.
  72    *        At destruction, remove triangles filling the holes
  73    */
  74   class HoleFiller
  75   {
  76   public:
  77     HoleFiller( SMESH_Mesh& meshDS );
  78     ~HoleFiller();
  79     void AddHoleBorders( Ng_STL_Geometry * ngStlGeo );
  80     void KeepHole() { myHole.clear(); myCapElems.clear(); }
  81     void ClearCapElements() { myCapElems.clear(); }
  82
  83   private:
  84     SMESHDS_Mesh*                          myMeshDS;
  85     std::vector< std::vector< gp_XYZ > >   myHole;      // initial border nodes
  86     std::vector< gp_XYZ >                  myInHolePos; // position inside each hole
  87     std::vector< const SMDS_MeshElement* > myCapElems;  // elements closing holes
  88   };
  89
  90   //================================================================================
  91   /*!
  92    * \brief Fill holes in the mesh
  93    */
  94   //================================================================================
  95
  96   HoleFiller::HoleFiller( SMESH_Mesh& theMesh ):
  97     myMeshDS( theMesh.GetMeshDS() )
  98   {
  99     SMESH_MeshEditor editor( &theMesh );
 100     {
 101       // // merge nodes
 102       // const double tol = Max( 0.1 * netgen::mparam.minh, Precision::Confusion() );
 103       // TIDSortedNodeSet allNodes;
 104       // SMESH_MeshEditor::TListOfListOfNodes equalNodes;
 105       // editor.FindCoincidentNodes( allNodes, tol, equalNodes, true );
 106       // editor.MergeNodes( equalNodes, /*noHoles=*/false );
 107     }
 108
 109     // find holes
 110     SMESH_MeshAlgos::TFreeBorderVec holes;
 111     bool isManifold = true, isGoodOri = true;
 112     SMESH_MeshAlgos::FindFreeBorders( *myMeshDS, holes, /*closedOnly=*/true,
 113                                       &isManifold, &isGoodOri );
 114
 115     if ( !isManifold )
 116     {
 117       // set bad faces into a compute error
 118       const char* text = "Non-manifold mesh. Only manifold mesh can be re-meshed";
 119       SMESH_BadInputElements* error =
 120         new SMESH_BadInputElements( myMeshDS, COMPERR_BAD_INPUT_MESH, text );
 121       SMESH::Controls::MultiConnection2D fun;
 122       fun.SetMesh( myMeshDS );
 123       SMDS_ElemIteratorPtr fIt = myMeshDS->elementsIterator( SMDSAbs_Face );
 124       while ( fIt->more() )
 125       {
 126         const SMDS_MeshElement* f = fIt->next();
 127         if ( fun.GetValue( f->GetID() ) > 2 )
 128           error->myBadElements.push_back( f );
 129       }
 130       theMesh.GetSubMesh( theMesh.GetShapeToMesh() )->GetComputeError().reset( error );
 131       throw SALOME_Exception( text );
 132     }
 133
 134     // don't want to sew coincident borders
 135     if ( !holes.empty() )
 136     {
 137       // define tolerance
 138       double tol, len, sumLen = 0, minLen = 1e100;
 139       int     nbSeg = 0;
 140       for ( size_t i = 0; i < holes.size(); ++i )
 141       {
 142         nbSeg += holes[i].size();
 143         SMESH_NodeXYZ p1 = holes[i][0];
 144         for ( size_t iP = 1; iP < holes[i].size(); ++iP )
 145         {
 146           SMESH_NodeXYZ p2 = holes[i][iP];
 147           len = ( p1 - p2 ).Modulus();
 148           sumLen += len;
 149           minLen = Min( minLen, len );
 150           p1 = p2;
 151         }
 152       }
 153       double avgLen = sumLen / nbSeg;
 154       if ( minLen > 1e-5 * avgLen )
 155         tol = 0.1 * minLen; // minLen is not degenerate
 156       else
 157         tol = 0.1 * avgLen;
 158
 159       SMESH_MeshAlgos::CoincidentFreeBorders freeBords;
 160       SMESH_MeshAlgos::FindCoincidentFreeBorders( *myMeshDS, tol, freeBords );
 161       if ( !freeBords._coincidentGroups.empty() )
 162       {
 163         const char* text = "Can't re-meshed a mesh with coincident free edges";
 164         SMESH_BadInputElements* error =
 165           new SMESH_BadInputElements( myMeshDS, COMPERR_BAD_INPUT_MESH, text );
 166         for ( size_t i = 0; i < freeBords._borders.size(); ++i )
 167           error->myBadElements.insert( error->myBadElements.end(),
 168                                        freeBords._borders[i].begin(),
 169                                        freeBords._borders[i].end() );
 170         theMesh.GetSubMesh( theMesh.GetShapeToMesh() )->GetComputeError().reset( error );
 171         throw SALOME_Exception( text );
 172       }
 173     }
 174
 175     // fill holes
 176     myHole.resize( holes.size() );
 177     myInHolePos.resize( holes.size() );
 178     std::vector<const SMDS_MeshElement*> newFaces;
 179     for ( size_t i = 0; i < holes.size(); ++i )
 180     {
 181       newFaces.clear();
 182       SMESH_MeshAlgos::FillHole( holes[i], *myMeshDS, newFaces );
 183
 184       // keep data to be able to remove hole filling faces after remeshing
 185       if ( !newFaces.empty() )
 186       {
 187         myHole[i].resize( holes[i].size() );
 188         for ( size_t iP = 0; iP < holes[i].size(); ++iP )
 189           myHole[i][iP] = SMESH_NodeXYZ( holes[i][iP] );
 190
 191         myInHolePos[i] = ( SMESH_NodeXYZ( newFaces[0]->GetNode(0)) +
 192                            SMESH_NodeXYZ( newFaces[0]->GetNode(1)) +
 193                            SMESH_NodeXYZ( newFaces[0]->GetNode(2)) ) / 3.;
 194         myCapElems.insert( myCapElems.end(), newFaces.begin(), newFaces.end() );
 195         // unmark to be able to remove them if meshing is canceled
 196         // for ( size_t iF = 0; iF < newFaces.size(); ++iF )
 197         //   newFaces[iF]->setIsMarked( false );
 198       }
 199     }
 200     // fix orientation
 201     if ( !isGoodOri )
 202     {
 203       SMDS_ElemIteratorPtr fIt = myMeshDS->elementsIterator( SMDSAbs_Face );
 204       while ( fIt->more() )
 205       {
 206         const SMDS_MeshElement* f = fIt->next();
 207         gp_XYZ normal;
 208         if ( SMESH_MeshAlgos::FaceNormal( f, normal ))
 209         {
 210           TIDSortedElemSet allFaces;
 211           editor.Reorient2D( allFaces, normal, f );
 212           break;
 213         }
 214       }
 215     }
 216   }
 217   //================================================================================
 218   /*!
 219    * \brief Add hole borders to be kept in a new mesh
 220    */
 221   //================================================================================
 222
 223   void HoleFiller::AddHoleBorders( Ng_STL_Geometry * ngStlGeo )
 224   {
 225     nglib::readedges.SetSize(0);
 226
 227     for ( size_t i = 0; i < myHole.size(); ++i )
 228       for ( size_t iP = 1; iP < myHole[i].size(); ++iP )
 229       {
 230         Ng_STL_AddEdge( ngStlGeo,
 231                         myHole[i][iP-1].ChangeData(),
 232                         myHole[i][iP-0].ChangeData() );
 233       }
 234   }
 235   //================================================================================
 236   /*!
 237    * \brief Remove triangles filling the holes
 238    */
 239   //================================================================================
 240
 241   HoleFiller::~HoleFiller()
 242   {
 243     if ( myMeshDS->NbNodes() < 3 )
 244       return;
 245
 246     if ( !myCapElems.empty() ) // old mesh not removed; simply remove myCapElems
 247     {
 248       for ( size_t i = 0; i < myCapElems.size(); ++i )
 249         myMeshDS->RemoveFreeElement( myCapElems[i], /*sm=*/0 );
 250       return;
 251     }
 252
 253     bool hasOrphanNodes = true;
 254
 255     const double tol = Max( 1e-3 * netgen::mparam.minh, Precision::Confusion() );
 256
 257     for ( size_t i = 0; i < myHole.size(); ++i )
 258     {
 259       std::vector< gp_XYZ >& borderPnt = myHole[i];
 260       const gp_XYZ&          inHolePos = myInHolePos[i];
 261       if ( borderPnt.empty() ) continue;
 262       borderPnt.pop_back(); // first point repeated at end
 263
 264       // mark all nodes located on the hole border
 265
 266       // new nodeSearcher for each hole, otherwise it contains removed nodes for i > 0
 267       SMESHUtils::Deleter< SMESH_NodeSearcher > nodeSearcher;
 268       if ( hasOrphanNodes )
 269       {
 270         std::vector< const SMDS_MeshNode* > sharedNodes;
 271         sharedNodes.reserve( myMeshDS->NbNodes() );
 272         SMDS_NodeIteratorPtr nIt = myMeshDS->nodesIterator();
 273         while ( nIt->more() )
 274         {
 275           const SMDS_MeshNode* n = nIt->next();
 276           if ( n->NbInverseElements() )
 277             sharedNodes.push_back( n );
 278         }
 279         hasOrphanNodes = ((int) sharedNodes.size() < myMeshDS->NbNodes() );
 280         SMDS_ElemIteratorPtr elemIt( new SMDS_NodeVectorElemIterator( sharedNodes.begin(),
 281                                                                       sharedNodes.end() ));
 282         nodeSearcher._obj = SMESH_MeshAlgos::GetNodeSearcher( elemIt );
 283       }
 284       else
 285       {
 286         nodeSearcher._obj = SMESH_MeshAlgos::GetNodeSearcher( *myMeshDS );
 287       }
 288
 289       std::vector< const SMDS_MeshElement* > edgesToRemove;
 290       edgesToRemove.reserve( borderPnt.size() );
 291
 292       // look for a border point coincident with a node
 293       size_t iP = 0;
 294       SMESH_NodeXYZ bordNode1;
 295       for ( ; iP < borderPnt.size(); ++iP )
 296       {
 297         bordNode1 = nodeSearcher->FindClosestTo( borderPnt[iP] );
 298         if (( bordNode1 - borderPnt[iP] ).SquareModulus() < tol * tol )
 299           break;
 300       }
 301       ++iP;
 302       bordNode1._node->setIsMarked( true );
 303
 304       // find the rest nodes located on the hole border
 305       boost::container::flat_set< const SMDS_MeshNode* > checkedNodes;
 306       gp_XYZ p1 = bordNode1;
 307       for ( size_t j = 0; j < borderPnt.size()+1; ++j,  iP = ( iP+1 ) % borderPnt.size() )
 308       {
 309         // among nodes surrounding bordNode1 find one most close to vec12
 310         gp_XYZ vec12 = borderPnt[iP] - p1;
 311         bool pntReached = false; // last found node is at iP
 312         while ( !pntReached )
 313         {
 314           const SMDS_MeshNode* bordNode = bordNode1._node;
 315           SMDS_ElemIteratorPtr fIt = bordNode->GetInverseElementIterator( SMDSAbs_Face );
 316           double minArea = 1e100;
 317           checkedNodes.clear();
 318           checkedNodes.insert( bordNode );
 319           while ( fIt->more() )
 320           {
 321             const SMDS_MeshElement* f = fIt->next();
 322             for ( int iN = 0, nbN = f->NbNodes(); iN < nbN; ++iN )
 323             {
 324               const SMDS_MeshNode* n = f->GetNode( iN );
 325               if ( !checkedNodes.insert( n ).second )
 326                 continue;
 327               SMESH_NodeXYZ pn = n;
 328               gp_XYZ vecPN = pn - bordNode1;
 329               if ( vecPN * vec12 <= 0 )
 330                 continue;
 331               gp_XYZ vec1N = pn - p1;
 332               double     a = vec12.CrossSquareMagnitude( vec1N );
 333               if ( a < minArea )
 334               {
 335                 bordNode = n;
 336                 minArea = a;
 337               }
 338             }
 339             if ( minArea < std::numeric_limits<double>::min() )
 340               break;
 341           }
 342           if ( bordNode == bordNode1._node )
 343             return; // bug in the loop above
 344
 345           SMESH_NodeXYZ bordNode2 = bordNode;
 346           gp_XYZ            vec1N = bordNode2 - p1;
 347           double u = ( vec12 * vec1N ) / vec12.SquareModulus(); // param [0,1] of bordNode on vec12
 348           if ( u < 1 + tol )
 349           {
 350             bordNode->setIsMarked( true );
 351             //cout << bordNode->GetID() << " ";
 352
 353             if ( const SMDS_MeshElement* edge = myMeshDS->FindEdge( bordNode1._node, bordNode ))
 354               edgesToRemove.push_back( edge );
 355             else
 356               edgesToRemove.push_back( myMeshDS->AddEdge( bordNode1._node, bordNode ));
 357             edgesToRemove.back()->setIsMarked( true );
 358
 359             if ( minArea > std::numeric_limits<double>::min() &&
 360                  minArea / vec12.SquareModulus() > tol * tol )
 361             {
 362               // node is far from the border, move it
 363               gp_XYZ p = p1 + u * vec12;
 364               myMeshDS->MoveNode( bordNode, p.X(), p.Y(), p.Z() );
 365             }
 366             bordNode1 = bordNode2;
 367           }
 368           //else -- there must be another border point between bordNode1 and bordNode
 369           pntReached = ( u > 1 - tol );
 370         }
 371         p1 = borderPnt[iP];
 372
 373       }
 374       //cout << endl << endl;
 375
 376       // remove all faces starting from inHolePos
 377
 378       // get a starting face
 379       std::vector< const SMDS_MeshNode* >     nodesToRemove;
 380       std::vector< const SMDS_MeshElement* >  facesToRemove;
 381       const SMDS_MeshNode* inHoleNode = nodeSearcher->FindClosestTo( inHolePos );
 382       if ( inHoleNode && ! inHoleNode->isMarked() )
 383       {
 384         SMDS_ElemIteratorPtr fIt = inHoleNode->GetInverseElementIterator( SMDSAbs_Face );
 385         while ( fIt->more() )
 386           facesToRemove.push_back( fIt->next() );
 387       }
 388       else
 389       {
 390         SMESHUtils::Deleter< SMESH_ElementSearcher > faceSearcher
 391           ( SMESH_MeshAlgos::GetElementSearcher( *myMeshDS ));
 392         if ( const SMDS_MeshElement* f = faceSearcher->FindClosestTo( inHolePos, SMDSAbs_Face ))
 393           facesToRemove.push_back( f );
 394         else
 395           continue;
 396       }
 397       for ( size_t iF = 0; iF < facesToRemove.size(); ++iF )
 398         facesToRemove[iF]->setIsMarked( true );
 399
 400       // remove faces and nodes
 401       TIDSortedElemSet elemSet, avoidSet;
 402       const SMDS_MeshElement* e;
 403       while ( !facesToRemove.empty() )
 404       {
 405         const SMDS_MeshElement* inHoleFace = facesToRemove.back();
 406         facesToRemove.pop_back();
 407
 408         // add adjacent faces into facesToRemove
 409         for ( int iN = 0, nbN = inHoleFace->NbNodes(); iN < nbN; ++iN )
 410         {
 411           const SMDS_MeshNode* n1 = inHoleFace->GetNode( iN );
 412           if ( !n1->isMarked() )
 413           {
 414             SMDS_ElemIteratorPtr eIt = n1->GetInverseElementIterator();
 415             while ( eIt->more() )
 416             {
 417               e = eIt->next();
 418               if ( e->GetType() == SMDSAbs_Face )
 419               {
 420                 if ( !e->isMarked() )
 421                   facesToRemove.push_back( e );
 422                 e->setIsMarked( true );
 423               }
 424               else if ( e->GetType() == SMDSAbs_Edge )
 425               {
 426                 myMeshDS->RemoveFreeElement( e, 0, /*fromGroups=*/false );
 427               }
 428             }
 429             if ( n1->NbInverseElements() == 1 )
 430               nodesToRemove.push_back( n1 );
 431           }
 432           else
 433           {
 434             const SMDS_MeshNode* n2 = inHoleFace->GetNodeWrap( iN+1 );
 435             if (( n2->isMarked() ) &&
 436                 ( !(e = myMeshDS->FindEdge( n1, n2 )) || !e->isMarked() )) // n1-n2 not hole border
 437             {
 438               if ( e ) // remove edge
 439                 myMeshDS->RemoveFreeElement( e, 0, /*fromGroups=*/false );
 440               avoidSet.clear();
 441               avoidSet.insert( inHoleFace );
 442               if (( e = SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet )))
 443               {
 444                 if ( !e->isMarked() )
 445                   facesToRemove.push_back( e );
 446                 e->setIsMarked( true );
 447               }
 448             }
 449           }
 450         }
 451         myMeshDS->RemoveFreeElement( inHoleFace, 0, /*fromGroups=*/false );
 452
 453         for ( size_t iN = 0; iN < nodesToRemove.size(); ++iN )
 454           myMeshDS->RemoveFreeNode( nodesToRemove[iN], 0, /*fromGroups=*/false );
 455         nodesToRemove.clear();
 456       }
 457
 458       // remove edges from the hole border
 459       // for ( size_t iE = 0; iE < edgesToRemove.size(); ++iE )
 460       //   myMeshDS->RemoveFreeElement( edgesToRemove[iE], 0, /*fromGroups=*/false );
 461
 462     } // loop on holes
 463
 464     return;
 465   } // ~HoleFiller()
 466
 467 } // namespace
 468
 469 //=============================================================================
 470 /*!
 471  * Constructor
 472  */
 473 //=============================================================================
 474
 475 NETGENPlugin_Remesher_2D::NETGENPlugin_Remesher_2D(int hypId, SMESH_Gen* gen)
 476   : SMESH_2D_Algo(hypId, gen)
 477 {
 478   _name = "NETGEN_Remesher_2D";
 479   _shapeType = (1 << TopAbs_FACE); // 1 bit /shape type
 480   _compatibleHypothesis.push_back("NETGEN_RemesherParameters_2D");
 481   _requireShape = false;
 482
 483   _hypothesis = 0;
 484 }
 485
 486 //=============================================================================
 487 /*!
 488  * Check assigned hypotheses
 489  */
 490 //=============================================================================
 491
 492 bool NETGENPlugin_Remesher_2D::CheckHypothesis (SMESH_Mesh&         theMesh,
 493                                                 const TopoDS_Shape& theShape,
 494                                                 Hypothesis_Status&  theStatus)
 495 {
 496   _hypothesis = 0;
 497
 498   // can work with no hypothesis
 499   theStatus = SMESH_Hypothesis::HYP_OK;
 500
 501   const list<const SMESHDS_Hypothesis*>& hyps =
 502     GetUsedHypothesis( theMesh, theShape, /*skipAux=*/true );
 503
 504   switch ( hyps.size() ) {
 505   case 0:
 506     break;
 507   case 1:
 508     _hypothesis = hyps.front();
 509     break;
 510   default:
 511     theStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 512   }
 513
 514   return theStatus == SMESH_Hypothesis::HYP_OK;
 515 }
 516
 517 //=============================================================================
 518 /*!
 519  * Compute mesh on an input mesh
 520  */
 521 //=============================================================================
 522
 523 bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh&         theMesh,
 524                                        SMESH_MesherHelper* theHelper)
 525 {
 526   if ( theMesh.NbFaces() == 0 )
 527     return !error( COMPERR_WARNING, "No faces in input mesh");
 528
 529   NETGENPlugin_Mesher mesher( &theMesh, theMesh.GetShapeToMesh(), /*isVol=*/false);
 530   NETGENPlugin_NetgenLibWrapper ngLib;
 531   netgen::Mesh *        ngMesh = (netgen::Mesh*) ngLib._ngMesh;
 532   Ng_STL_Geometry *   ngStlGeo = Ng_STL_NewGeometry();
 533   netgen::STLTopology* stlTopo = (netgen::STLTopology*) ngStlGeo;
 534   netgen::multithread.terminate = 0;
 535
 536   const NETGENPlugin_RemesherHypothesis_2D* hyp =
 537     dynamic_cast<const NETGENPlugin_RemesherHypothesis_2D*>( _hypothesis );
 538   mesher.SetParameters( hyp );// for holeFiller
 539
 540   SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
 541   HoleFiller holeFiller( theMesh );
 542   //theHelper->SetIsQuadratic( theMesh.NbFaces( ORDER_QUADRATIC ));
 543
 544   // fill ngStlGeo with triangles
 545   SMDS_ElemIteratorPtr fIt = meshDS->elementsIterator( SMDSAbs_Face );
 546   while ( fIt->more() )
 547   {
 548     const SMDS_MeshElement* f = fIt->next();
 549     SMESH_NodeXYZ n1 = f->GetNode( 0 );
 550     SMESH_NodeXYZ n2 = f->GetNode( 1 );
 551     SMESH_NodeXYZ n3 = f->GetNode( 2 );
 552     Ng_STL_AddTriangle( ngStlGeo,
 553                         n1.ChangeData(),
 554                         n2.ChangeData(),
 555                         n3.ChangeData() );
 556     if ( f->NbNodes() > 3 )
 557     {
 558       n2.Set( f->GetNode( 3 ));
 559       Ng_STL_AddTriangle( ngStlGeo,
 560                           n1.ChangeData(),
 561                           n3.ChangeData(),
 562                           n2.ChangeData());
 563     }
 564   }
 565   // add edges
 566   holeFiller.AddHoleBorders( ngStlGeo );
 567
 568   // init stl DS
 569   Ng_Result ng_res = Ng_STL_InitSTLGeometry( ngStlGeo );
 570   if ( ng_res != NG_OK )
 571   {
 572 #ifdef _DEBUG_
 573     holeFiller.KeepHole();
 574 #endif
 575     std::string txt = "Error Initialising the STL Geometry";
 576     if ( !stlTopo->GetStatusText().empty() )
 577       txt += ". " + stlTopo->GetStatusText();
 578     return error( COMPERR_BAD_INPUT_MESH, txt );
 579   }
 580
 581   Ng_Meshing_Parameters ngParams;
 582   ng_res = Ng_STL_MakeEdges( ngStlGeo, ngLib._ngMesh, &ngParams );
 583   if ( ng_res != NG_OK )
 584     return error( "Error in Edge Meshing" );
 585
 586   // set parameters
 587   if ( hyp )
 588   {
 589     ngParams.maxh              = hyp->GetMaxSize();
 590     ngParams.minh              = hyp->GetMinSize();
 591     ngParams.meshsize_filename = (char*) hyp->GetMeshSizeFile().c_str();
 592     ngParams.quad_dominated    = hyp->GetQuadAllowed();
 593     netgen::stlparam.yangle    = hyp->GetRidgeAngle();
 594     mesher.SetParameters( hyp );
 595   }
 596   else
 597   {
 598     double diagSize = Dist( stlTopo->GetBoundingBox().PMin(), stlTopo->GetBoundingBox().PMax());
 599     netgen::mparam.maxh = diagSize / GetGen()->GetBoundaryBoxSegmentation();
 600     netgen::mparam.minh = netgen::mparam.maxh;
 601   }
 602
 603   // TODO: expose stlparam.resth* to the user
 604   // netgen::stlparam.resthcloseedgeenable = 0; // Restrict H due to close edges
 605   // netgen::stlparam.resthlinelengthenable = 0; // Restrict H due to line-length
 606   // netgen::stlparam.resthatlasenable = 0;
 607   // //netgen::stlparam.resthchartdistenable = 0;
 608
 609   double h = netgen::mparam.maxh;
 610   ngMesh->SetGlobalH( h );
 611   ngMesh->SetMinimalH( netgen::mparam.minh );
 612   ngMesh->SetLocalH( stlTopo->GetBoundingBox().PMin() - netgen::Vec3d(h, h, h),
 613                      stlTopo->GetBoundingBox().PMax() + netgen::Vec3d(h, h, h),
 614                      netgen::mparam.grading );
 615   ngMesh->LoadLocalMeshSize( ngParams.meshsize_filename );
 616
 617   netgen::OCCGeometry occgeo;
 618   mesher.SetLocalSize( occgeo, *ngMesh );
 619
 620   // meshing
 621   try
 622   {
 623     ng_res = Ng_STL_GenerateSurfaceMesh( ngStlGeo, ngLib._ngMesh, &ngParams );
 624   }
 625   catch (netgen::NgException & ex)
 626   {
 627     if ( netgen::multithread.terminate )
 628       return false;
 629   }
 630   if ( ng_res != NG_OK )
 631     return error( "Error in Surface Meshing" );
 632
 633   int nbN = ngMesh->GetNP();
 634   int nbE = ngMesh->GetNSeg();
 635   int nbF = ngMesh->GetNSE();
 636   if ( nbF == 0 )
 637     return error( "Error in Surface Meshing" );
 638
 639   // remove existing mesh
 640   holeFiller.ClearCapElements();
 641   SMDS_ElemIteratorPtr eIt = meshDS->elementsIterator();
 642   while ( eIt->more() )
 643     meshDS->RemoveFreeElement( eIt->next(), /*sm=*/0 );
 644   SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator();
 645   while ( nIt->more() )
 646     meshDS->RemoveFreeNode( nIt->next(), /*sm=*/0 );
 647
 648   // retrieve new mesh
 649
 650   // add nodes
 651   std::vector< const SMDS_MeshNode* > newNodes( nbN+1 );
 652   for ( int i = 1; i <= nbN; ++i )
 653   {
 654     const netgen::MeshPoint& p = ngMesh->Point(i);
 655     newNodes[i] = meshDS->AddNode( p(0),p(1),p(2) );
 656   }
 657
 658   // add edges
 659   std::vector<const SMDS_MeshNode*> nodes(4);
 660   for ( int i = 1; i <= nbE; ++i )
 661   {
 662     const netgen::Segment& seg = ngMesh->LineSegment(i);
 663     nodes.clear();
 664     for ( int j = 0; j < 2; ++j )
 665     {
 666       size_t pind = seg.pnums[j];
 667       if ( pind > 0 && pind < newNodes.size() )
 668         nodes.push_back( newNodes[ pind ]);
 669       else
 670         break;
 671     }
 672     if ( nodes.size() == 2 && !meshDS->FindEdge( nodes[0], nodes[1] ))
 673       meshDS->AddEdge( nodes[0], nodes[1] );
 674   }
 675
 676   // add faces
 677   for ( int i = 1; i <= nbF; ++i )
 678   {
 679     const netgen::Element2d& elem = ngMesh->SurfaceElement(i);
 680     nodes.clear();
 681     for ( int j = 1; j <= elem.GetNP(); ++j )
 682     {
 683       size_t pind = elem.PNum(j);
 684       if ( pind > 0 && pind < newNodes.size() )
 685         nodes.push_back( newNodes[ pind ]);
 686       else
 687         break;
 688     }
 689     switch( nodes.size() )
 690     {
 691     case 3: meshDS->AddFace( nodes[0], nodes[1], nodes[2] ); break;
 692     case 4: meshDS->AddFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
 693     }
 694   }
 695
 696   // as we don't assign the new triangles to a shape (the pseudo-shape),
 697   // to avoid their removal at hypothesis modification,
 698   // we mark the shape as always computed to avoid the error messages
 699   // that no elements assigned to the shape
 700   theMesh.GetSubMesh( theHelper->GetSubShape() )->SetIsAlwaysComputed( true );
 701
 702   return true;
 703 }
 704
 705 //=============================================================================
 706 /*!
 707  * Do not compute mesh on geometry
 708  */
 709 //=============================================================================
 710
 711 bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh&         theMesh,
 712                                        const TopoDS_Shape& theShape)
 713 {
 714   return false;
 715 }
 716
 717 //=============================================================================
 718 /*!
 719  * Terminate Compute()
 720  */
 721 //=============================================================================
 722
 723 void NETGENPlugin_Remesher_2D::CancelCompute()
 724 {
 725   SMESH_Algo::CancelCompute();
 726   netgen::multithread.terminate = 1;
 727 }
 728
 729 //================================================================================
 730 /*!
 731  * \brief Return progress of Compute() [0.,1]
 732  */
 733 //================================================================================
 734
 735 double NETGENPlugin_Remesher_2D::GetProgress() const
 736 {
 737   return netgen::multithread.percent / 100.;
 738 }
 739
 740 //=============================================================================
 741 /*!
 742  *
 743  */
 744 //=============================================================================
 745
 746 bool NETGENPlugin_Remesher_2D::Evaluate(SMESH_Mesh&         aMesh,
 747                                         const TopoDS_Shape& aShape,
 748                                         MapShapeNbElems& aResMap)
 749 {
 750   return false;
 751 }