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