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