src/NETGENPlugin/NETGENPlugin_Remesher_2D.cxx

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