src/NETGENPlugin/NETGENPlugin_Remesher_2D.cxx

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