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