src/NETGENPlugin/NETGENPlugin_Remesher_2D.cxx

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