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