src/NETGENPlugin/NETGENPlugin_NETGEN_2D_ONLY.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 // File      : NETGENPlugin_NETGEN_2D_ONLY.cxx
  21 // Author    : Edward AGAPOV (OCC)
  22 // Project   : SALOME
  23 //
  24 #include "NETGENPlugin_NETGEN_2D_ONLY.hxx"
  25 #include "NETGENPlugin_Hypothesis_2D.hxx"
  26
  27 #include <SMDS_MeshElement.hxx>
  28 #include <SMDS_MeshNode.hxx>
  29 #include <SMESHDS_Mesh.hxx>
  30 #include <SMESH_Comment.hxx>
  31 #include <SMESH_Gen.hxx>
  32 #include <SMESH_Mesh.hxx>
  33 #include <SMESH_MesherHelper.hxx>
  34 #include <SMESH_subMesh.hxx>
  35 #include <StdMeshers_FaceSide.hxx>
  36 #include <StdMeshers_LengthFromEdges.hxx>
  37 #include <StdMeshers_MaxElementArea.hxx>
  38 #include <StdMeshers_QuadranglePreference.hxx>
  39 #include <StdMeshers_ViscousLayers2D.hxx>
  40
  41 #include <Precision.hxx>
  42 #include <Standard_ErrorHandler.hxx>
  43 #include <Standard_Failure.hxx>
  44
  45 #include <utilities.h>
  46
  47 #include <GEOMUtils.hxx>
  48
  49 #include <list>
  50 #include <vector>
  51 #include <limits>
  52
  53 /*
  54   Netgen include files
  55 */
  56 namespace nglib {
  57 #include <nglib.h>
  58 }
  59 #ifndef OCCGEOMETRY
  60 #define OCCGEOMETRY
  61 #endif
  62 #include <occgeom.hpp>
  63 #include <meshing.hpp>
  64
  65 #include "SMESH_Octree.hxx"
  66
  67 //// Used for node projection in curve
  68 #include <BRepAdaptor_Curve.hxx>
  69 #include <BRep_Builder.hxx>
  70 #include <BRep_Tool.hxx>
  71 #include <BndLib_Add3dCurve.hxx>
  72 #include <GCPnts_TangentialDeflection.hxx>
  73 #include <ShapeAnalysis_Curve.hxx>
  74 #include <TopExp.hxx>
  75 #include <TopExp_Explorer.hxx>
  76 #include <TopoDS.hxx>
  77 #include <TopoDS_Compound.hxx>
  78 #include <TopoDS_Edge.hxx>
  79 #include <TopoDS_Vertex.hxx>
  80
  81 //#include <meshtype.hpp>
  82 namespace netgen {
  83   NETGENPLUGIN_DLL_HEADER
  84   extern MeshingParameters mparam;
  85 #ifdef NETGEN_V5
  86   extern void OCCSetLocalMeshSize(OCCGeometry & geom, Mesh & mesh);
  87 #endif
  88 }
  89
  90 using namespace std;
  91 using namespace netgen;
  92 using namespace nglib;
  93
  94 namespace // copied class from StdMesher_Importe_1D
  95 {
  96   /*!
  97   * \brief Compute point position on a curve. Use octree to fast reject far points
  98   */
  99   class CurveProjector : public SMESH_Octree
 100   {
 101   public:
 102     CurveProjector( const TopoDS_Edge& edge, double enlarge );
 103
 104     bool IsOnCurve( const gp_XYZ& point, double & distance2, double & u );
 105
 106     bool IsOut( const gp_XYZ& point ) const { return getBox()->IsOut( point ); }
 107
 108   protected:
 109     CurveProjector() {}
 110     SMESH_Octree* newChild() const { return new CurveProjector; }
 111     void          buildChildrenData();
 112     Bnd_B3d*      buildRootBox();
 113
 114   private:
 115     struct CurveSegment : public Bnd_B3d
 116     {
 117       double _chord, _chord2, _length2;
 118       gp_Pnt _pFirst, _pLast;
 119       gp_Lin _line;
 120       Handle(Geom_Curve) _curve;
 121
 122       CurveSegment() {}
 123       void Init( const gp_Pnt& pf, const gp_Pnt& pl,
 124                  double uf, double ul, double tol, Handle(Geom_Curve)& curve );
 125       bool IsOn( const gp_XYZ& point, double & distance2, double & u );
 126       bool IsInContact( const Bnd_B3d& bb );
 127     };
 128     std::vector< CurveSegment > _segments;
 129   };
 130
 131     //===============================================================================
 132   /*!
 133    * \brief Create an octree of curve segments
 134    */
 135   //================================================================================
 136
 137   CurveProjector::CurveProjector( const TopoDS_Edge& edge, double enlarge )
 138     :SMESH_Octree( 0 )
 139   {
 140     double f,l;
 141     Handle(Geom_Curve) curve = BRep_Tool::Curve( edge, f, l );
 142     double curDeflect = 0.3; // Curvature deflection
 143     double angDeflect = 1e+100; // Angular deflection - don't control chordal error
 144     GCPnts_TangentialDeflection div( BRepAdaptor_Curve( edge ), angDeflect, curDeflect );
 145     _segments.resize( div.NbPoints() - 1 );
 146     for ( int i = 1; i < div.NbPoints(); ++i )
 147       try {
 148         _segments[ i - 1 ].Init( div.Value( i ),     div.Value( i+1 ),
 149                                  div.Parameter( i ), div.Parameter( i+1 ),
 150                                  enlarge, curve );
 151       }
 152       catch ( Standard_Failure& ) {
 153         _segments.resize( _segments.size() - 1 );
 154         --i;
 155       }
 156     if ( _segments.size() < 3 )
 157       myIsLeaf = true;
 158
 159     compute();
 160
 161     if ( _segments.size() == 1 )
 162       myBox->Enlarge( enlarge );
 163   }
 164
 165   //================================================================================
 166   /*!
 167    * \brief Return the maximal box
 168    */
 169   //================================================================================
 170
 171   Bnd_B3d* CurveProjector::buildRootBox()
 172   {
 173     Bnd_B3d* box = new Bnd_B3d;
 174     for ( size_t i = 0; i < _segments.size(); ++i )
 175       box->Add( _segments[i] );
 176     return box;
 177   }
 178
 179   //================================================================================
 180   /*!
 181    * \brief Redistribute segments among children
 182    */
 183   //================================================================================
 184
 185   void CurveProjector::buildChildrenData()
 186   {
 187     bool allIn = true;
 188     for ( size_t i = 0; i < _segments.size(); ++i )
 189     {
 190       for (int j = 0; j < 8; j++)
 191       {
 192         if ( _segments[i].IsInContact( *myChildren[j]->getBox() ))
 193           ((CurveProjector*)myChildren[j])->_segments.push_back( _segments[i]);
 194         else
 195           allIn = false;
 196       }
 197     }
 198     if ( allIn && _segments.size() < 3 )
 199     {
 200       myIsLeaf = true;
 201       for (int j = 0; j < 8; j++)
 202         static_cast<CurveProjector*>( myChildren[j])->myIsLeaf = true;
 203     }
 204     else
 205     {
 206       SMESHUtils::FreeVector( _segments ); // = _segments.clear() + free memory
 207
 208       for (int j = 0; j < 8; j++)
 209       {
 210         CurveProjector* child = static_cast<CurveProjector*>( myChildren[j]);
 211         if ( child->_segments.size() < 3 )
 212           child->myIsLeaf = true;
 213       }
 214     }
 215   }
 216
 217   //================================================================================
 218   /*!
 219    * \brief Return true if a point is close to the curve
 220    *  \param [in] point - the point
 221    *  \param [out] distance2 - distance to the curve
 222    *  \param [out] u - parameter on the curve
 223    *  \return bool - is the point is close to the curve
 224    */
 225   //================================================================================
 226
 227   bool CurveProjector::IsOnCurve( const gp_XYZ& point, double & distance2, double & u )
 228   {
 229     if ( getBox()->IsOut( point ))
 230       return false;
 231
 232     bool ok = false;
 233     double dist2, param;
 234     distance2 = Precision::Infinite();
 235
 236     if ( isLeaf() )
 237     {
 238       for ( size_t i = 0; i < _segments.size(); ++i )
 239         if ( !_segments[i].IsOut( point ) &&
 240              _segments[i].IsOn( point, dist2, param ) &&
 241              dist2 < distance2 )
 242         {
 243           distance2 = dist2;
 244           u         = param;
 245           ok        = true;
 246         }
 247       return ok;
 248     }
 249     else
 250     {
 251       for (int i = 0; i < 8; i++)
 252         if (((CurveProjector*) myChildren[i])->IsOnCurve( point, dist2, param ) &&
 253             dist2 < distance2 )
 254         {
 255           distance2 = dist2;
 256           u         = param;
 257           ok        = true;
 258         }
 259     }
 260     return ok;
 261   }
 262
 263   //================================================================================
 264   /*!
 265    * \brief Initialize
 266    */
 267   //================================================================================
 268
 269   void CurveProjector::CurveSegment::Init(const gp_Pnt&       pf,
 270                                           const gp_Pnt&       pl,
 271                                           const double        uf,
 272                                           const double        ul,
 273                                           const double        tol,
 274                                           Handle(Geom_Curve)& curve )
 275   {
 276     _pFirst  = pf;
 277     _pLast   = pl;
 278     _curve   = curve;
 279     _length2 = pf.SquareDistance( pl );
 280     _line.SetLocation( pf );
 281     _line.SetDirection( gp_Vec( pf, pl ));
 282     _chord2  = Max( _line.     SquareDistance( curve->Value( uf + 0.25 * ( ul - uf ))),
 283                     Max( _line.SquareDistance( curve->Value( uf + 0.5  * ( ul - uf ))),
 284                          _line.SquareDistance( curve->Value( uf + 0.75 * ( ul - uf )))));
 285     _chord2 *= ( 1.05 * 1.05 ); // +5%
 286     _chord2  = Max( tol, _chord2 );
 287     _chord   = Sqrt( _chord2 );
 288
 289     Bnd_Box bb;
 290     BndLib_Add3dCurve::Add( GeomAdaptor_Curve( curve, uf, ul ), tol, bb );
 291     Add( bb.CornerMin() );
 292     Add( bb.CornerMax() );
 293   }
 294
 295   //================================================================================
 296   /*!
 297    * \brief Return true if a point is close to the curve segment
 298    *  \param [in] point - the point
 299    *  \param [out] distance2 - distance to the curve
 300    *  \param [out] u - parameter on the curve
 301    *  \return bool - is the point is close to the curve segment
 302    */
 303   //================================================================================
 304
 305   bool CurveProjector::CurveSegment::IsOn( const gp_XYZ& point, double & distance2, double & u )
 306   {
 307     distance2 = _line.SquareDistance( point );
 308     if ( distance2 > _chord2 )
 309       return false;
 310
 311     // check if the point projection falls into the segment range
 312     {
 313       gp_Vec edge( _pFirst, _pLast );
 314       gp_Vec n1p ( _pFirst, point  );
 315       u = ( edge * n1p ) / _length2; // param [0,1] on the edge
 316       if ( u < 0. )
 317       {
 318         if ( _pFirst.SquareDistance( point ) > _chord2 )
 319           return false;
 320       }
 321       else if ( u > 1. )
 322       {
 323         if ( _pLast.SquareDistance( point ) > _chord2 )
 324           return false;
 325       }
 326     }
 327     gp_Pnt proj;
 328     distance2 = ShapeAnalysis_Curve().Project( _curve, point, Precision::Confusion(),
 329                                                proj, u, false );
 330     distance2 *= distance2;
 331     return true;
 332   }
 333
 334   //================================================================================
 335   /*!
 336    * \brief Check if the segment is in contact with a box
 337    */
 338   //================================================================================
 339
 340   bool CurveProjector::CurveSegment::IsInContact( const Bnd_B3d& bb )
 341   {
 342     if ( bb.IsOut( _line.Position(), /*isRay=*/true, _chord ))
 343       return false;
 344
 345     gp_Ax1 axRev = _line.Position().Reversed();
 346     axRev.SetLocation( _pLast );
 347     return !bb.IsOut( axRev, /*isRay=*/true, _chord );
 348   }
 349
 350 }
 351
 352 //=============================================================================
 353 /*!
 354  *
 355  */
 356 //=============================================================================
 357
 358 NETGENPlugin_NETGEN_2D_ONLY::NETGENPlugin_NETGEN_2D_ONLY(int        hypId,
 359                                                          SMESH_Gen* gen)
 360   : SMESH_2D_Algo(hypId, gen)
 361 {
 362   _name = "NETGEN_2D_ONLY";
 363
 364   _shapeType = (1 << TopAbs_FACE);// 1 bit /shape type
 365   _onlyUnaryInput = false; // treat all FACEs at once
 366
 367   _compatibleHypothesis.push_back("MaxElementArea");
 368   _compatibleHypothesis.push_back("LengthFromEdges");
 369   _compatibleHypothesis.push_back("QuadranglePreference");
 370   _compatibleHypothesis.push_back("NETGEN_Parameters_2D");
 371   _compatibleHypothesis.push_back("ViscousLayers2D");
 372
 373   _hypMaxElementArea       = 0;
 374   _hypLengthFromEdges      = 0;
 375   _hypQuadranglePreference = 0;
 376   _hypParameters           = 0;
 377 }
 378
 379 //=============================================================================
 380 /*!
 381  *
 382  */
 383 //=============================================================================
 384
 385 NETGENPlugin_NETGEN_2D_ONLY::~NETGENPlugin_NETGEN_2D_ONLY()
 386 {
 387   //MESSAGE("NETGENPlugin_NETGEN_2D_ONLY::~NETGENPlugin_NETGEN_2D_ONLY");
 388 }
 389
 390 //=============================================================================
 391 /*!
 392  *
 393  */
 394 //=============================================================================
 395
 396 bool NETGENPlugin_NETGEN_2D_ONLY::CheckHypothesis (SMESH_Mesh&         aMesh,
 397                                                    const TopoDS_Shape& aShape,
 398                                                    Hypothesis_Status&  aStatus)
 399 {
 400   _hypMaxElementArea = 0;
 401   _hypLengthFromEdges = 0;
 402   _hypQuadranglePreference = 0;
 403   _hypParameters = 0;
 404   _progressByTic = -1;
 405
 406   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape, false);
 407
 408   if (hyps.empty())
 409   {
 410     aStatus = HYP_OK; //SMESH_Hypothesis::HYP_MISSING;
 411     return true;  // (PAL13464) can work with no hypothesis, LengthFromEdges is default one
 412   }
 413
 414   aStatus = HYP_MISSING;
 415
 416   bool hasVL = false;
 417   list<const SMESHDS_Hypothesis*>::const_iterator ith;
 418   for (ith = hyps.begin(); ith != hyps.end(); ++ith )
 419   {
 420     const SMESHDS_Hypothesis* hyp = (*ith);
 421
 422     string hypName = hyp->GetName();
 423
 424     if      ( hypName == "MaxElementArea")
 425       _hypMaxElementArea = static_cast<const StdMeshers_MaxElementArea*> (hyp);
 426     else if ( hypName == "LengthFromEdges" )
 427       _hypLengthFromEdges = static_cast<const StdMeshers_LengthFromEdges*> (hyp);
 428     else if ( hypName == "QuadranglePreference" )
 429       _hypQuadranglePreference = static_cast<const StdMeshers_QuadranglePreference*>(hyp);
 430     else if ( hypName == "NETGEN_Parameters_2D" )
 431       _hypParameters = static_cast<const NETGENPlugin_Hypothesis_2D*>(hyp);
 432     else if ( hypName == StdMeshers_ViscousLayers2D::GetHypType() )
 433       hasVL = true;
 434     else {
 435       aStatus = HYP_INCOMPATIBLE;
 436       return false;
 437     }
 438   }
 439
 440   int nbHyps = bool(_hypMaxElementArea) + bool(_hypLengthFromEdges) + bool(_hypParameters );
 441   if ( nbHyps > 1 )
 442     aStatus = HYP_CONCURRENT;
 443   else if ( hasVL )
 444     error( StdMeshers_ViscousLayers2D::CheckHypothesis( aMesh, aShape, aStatus ));
 445   else
 446     aStatus = HYP_OK;
 447
 448   if ( aStatus == HYP_OK && _hypParameters && _hypQuadranglePreference )
 449   {
 450     aStatus = HYP_INCOMPAT_HYPS;
 451     return error(SMESH_Comment("\"") << _hypQuadranglePreference->GetName()
 452                  << "\" and \"" << _hypParameters->GetName()
 453                  << "\" are incompatible hypotheses");
 454   }
 455
 456   return ( aStatus == HYP_OK );
 457 }
 458
 459 // namespace
 460 // {
 461 //   void limitSize( netgen::Mesh* ngMesh,
 462 //                   const double  maxh )
 463 //   {
 464 //     // get bnd box
 465 //     netgen::Point3d pmin, pmax;
 466 //     ngMesh->GetBox( pmin, pmax, 0 );
 467 //     const double dx = pmax.X() - pmin.X();
 468 //     const double dy = pmax.Y() - pmin.Y();
 469 //     const double dz = pmax.Z() - pmin.Z();
 470
 471 //     const int nbX = Max( 2, int( dx / maxh * 3 ));
 472 //     const int nbY = Max( 2, int( dy / maxh * 3 ));
 473 //     const int nbZ = Max( 2, int( dz / maxh * 3 ));
 474
 475 //     if ( ! & ngMesh->LocalHFunction() )
 476 //       ngMesh->SetLocalH( pmin, pmax, 0.1 );
 477
 478 //     netgen::Point3d p;
 479 //     for ( int i = 0; i <= nbX; ++i )
 480 //     {
 481 //       p.X() = pmin.X() +  i * dx / nbX;
 482 //       for ( int j = 0; j <= nbY; ++j )
 483 //       {
 484 //         p.Y() = pmin.Y() +  j * dy / nbY;
 485 //         for ( int k = 0; k <= nbZ; ++k )
 486 //         {
 487 //           p.Z() = pmin.Z() +  k * dz / nbZ;
 488 //           ngMesh->RestrictLocalH( p, maxh );
 489 //         }
 490 //       }
 491 //     }
 492 //   }
 493 // }
 494
 495
 496
 497 /**
 498  * @brief MapSegmentsToEdges.
 499  * @remark To feed 1D segments not associated to any geometry we need:
 500  *          1) For each face, check all segments that are in the face (use ShapeAnalysis_Surface class)
 501  *          2) Check to which edge the segment below to, use the copied [from StdMesher_Importe_1D] CurveProjector class
 502  *          3) Create new netgen segments with the (u,v) parameters obtained from the ShapeAnalysis_Surface projector
 503  *          4) also define the 'param' value of the nodes relative to the edges obtained from CurveProjector
 504  *          5) Add the new netgen segments IN ORDER into the netgen::mesh data structure to form a closed chain
 505  *          6) Beware with the occ::edge orientation
 506  *
 507  * @param aMesh Mesh file (containing 1D elements)
 508  * @param aShape Shape file (BREP or STEP format)
 509  * @param ngLib netgenlib library wrapper
 510  * @param nodeVec vector of nodes used internally to feed smesh aMesh after computation
 511  * @param premeshedNodes map of prmeshed nodes and the smesh nodeID associate to it
 512  * @param newNetgenCoordinates map of 3D coordinate of new points created by netgen
 513  * @param newNetgenElements map of triangular or quadrangular elements ID and the nodes defining the 2D element
 514  * @return false
 515  */
 516 bool NETGENPlugin_NETGEN_2D_ONLY::MapSegmentsToEdges(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape, NETGENPlugin_NetgenLibWrapper &ngLib,
 517                                                       vector< const SMDS_MeshNode* >& nodeVec, std::map<int,const SMDS_MeshNode*>& premeshedNodes,
 518                                                       std::map<int,std::vector<double>>& newNetgenCoordinates, std::map<int,std::vector<smIdType>>& newNetgenElements )
 519 {
 520   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 521   SMESH_MesherHelper helper(aMesh);
 522   helper.SetElementsOnShape( true );
 523   const int numberOfPremeshedNodes = aMesh.NbNodes();
 524   TopTools_IndexedMapOfShape faces;
 525   TopExp::MapShapes( aShape, TopAbs_FACE, faces );
 526   int err = 0;
 527   for ( int i = 1; i <= faces.Size(); ++i )
 528   {
 529     int numOfEdges = 0;
 530     int totalEdgeLenght = 0;
 531
 532     TopoDS_Face face = TopoDS::Face(faces( i ) );
 533     Bnd_Box FaceBox;
 534     GEOMUtils::PreciseBoundingBox( face, FaceBox );
 535     if ( face.Orientation() != TopAbs_FORWARD && face.Orientation() != TopAbs_REVERSED )
 536       face.Orientation( TopAbs_FORWARD );
 537
 538     // Set the occgeom to be meshed and some ngMesh parameteres
 539     netgen::OCCGeometry occgeom;
 540     netgen::Mesh * ngMesh =  (netgen::Mesh*) ngLib._ngMesh;
 541     ngMesh->DeleteMesh();
 542
 543     occgeom.shape = face;
 544     occgeom.fmap.Add( face );
 545     occgeom.CalcBoundingBox();
 546     occgeom.facemeshstatus.SetSize(1);
 547     occgeom.facemeshstatus = 0;
 548     occgeom.face_maxh_modified.SetSize(1);
 549     occgeom.face_maxh_modified = 0;
 550     occgeom.face_maxh.SetSize(1);
 551     occgeom.face_maxh = netgen::mparam.maxh;
 552
 553     // Set the face descriptor
 554     const int solidID = 0, faceID = 1; /*always 1 because faces are meshed one by one*/
 555     if ( ngMesh->GetNFD() < 1 )
 556       ngMesh->AddFaceDescriptor( netgen::FaceDescriptor( faceID, solidID, solidID, 0 ));
 557
 558     ngMesh->SetGlobalH ( netgen::mparam.maxh );
 559     ngMesh->SetMinimalH( netgen::mparam.minh );
 560     Box<3> bb = occgeom.GetBoundingBox();
 561     bb.Increase (bb.Diam()/10);
 562     ngMesh->SetLocalH (bb.PMin(), bb.PMax(), netgen::mparam.grading);
 563     // end set the occgeom to be meshed and some ngMesh parameteres
 564
 565     Handle(ShapeAnalysis_Surface) sprojector = new ShapeAnalysis_Surface( BRep_Tool::Surface( face ));
 566     double tol = BRep_Tool::MaxTolerance( face, TopAbs_FACE );
 567     gp_Pnt surfPnt(0,0,0);
 568
 569     map<const SMDS_MeshNode*, int > node2ngID;
 570     map<int, const SMDS_MeshNode* > ng2smesh;
 571     // Figure out which edge is this onde in!
 572     TopTools_IndexedMapOfShape edges;
 573     TopExp::MapShapes( face, TopAbs_EDGE, edges );
 574     TopoDS_Edge meshingEdge;
 575     // Check wich nodes are in this face!
 576     SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge);
 577     std::unique_ptr<CurveProjector> myCurveProjector;
 578     while ( iteratorElem->more() ) // loop on elements on a geom face
 579     {
 580       // check mesh face
 581       const SMDS_MeshElement* elem = iteratorElem->next();
 582       const SMDS_MeshNode* node0 = elem->GetNode( 0 );
 583       const SMDS_MeshNode* node1 = elem->GetNode( 1 );
 584       SMESH_NodeXYZ nXYZ0( node0 );
 585       SMESH_NodeXYZ nXYZ1( node1 );
 586       double segmentLength = ( nXYZ0 - nXYZ1 ).Modulus();
 587       nXYZ0 += nXYZ1;
 588       nXYZ0 /= 2.0;
 589       if( FaceBox.IsOut( nXYZ0 ) )
 590         continue;
 591
 592       gp_XY uv = sprojector->ValueOfUV( nXYZ0, tol ).XY();
 593       surfPnt = sprojector->Value( uv );
 594       double dist = surfPnt.Distance( nXYZ0 );
 595
 596       // used for the curve projector of edges
 597       double geomTol = Precision::Confusion();
 598       if ( dist < tol /*element is on face*/ )
 599       {
 600         numOfEdges++;
 601         totalEdgeLenght += segmentLength;
 602         int occEdgeIdFound = -1;
 603         for ( int edgeId = 1; edgeId <= edges.Size(); ++edgeId ) /*find in which edge the node is placed*/
 604         {
 605           meshingEdge = TopoDS::Edge(edges( edgeId ));
 606           myCurveProjector = std::unique_ptr<CurveProjector>( new CurveProjector(meshingEdge, geomTol) );
 607           if ( myCurveProjector->IsOut( nXYZ0 ) /*keep searching*/)
 608             continue;
 609           else
 610           {
 611             occEdgeIdFound = edgeId;
 612             break;
 613           }
 614         }
 615
 616         netgen::Segment seg;
 617         for ( size_t nodeId = 0; nodeId < 2; nodeId++)
 618         {
 619           const SMDS_MeshNode* node = elem->GetNode( nodeId );
 620           int ngId = ngMesh->GetNP() + 1;
 621           ngId = node2ngID.insert( make_pair( node, ngId )).first->second;
 622           if ( ngId > ngMesh->GetNP() /* mean it is a new node to be add to the mesh*/)
 623           {
 624             // Restric size of mesh based on the edge dimension
 625             {
 626               SMESH_NodeXYZ nXYZ( node );
 627               netgen::Point3d pi(nXYZ.X(), nXYZ.Y(), nXYZ.Z());
 628               ngMesh->RestrictLocalH( pi, segmentLength );
 629             }
 630
 631             netgen::MeshPoint mp( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
 632             ngMesh->AddPoint ( mp, 1, netgen::EDGEPOINT );
 633             ng2smesh.insert( make_pair(ngId, node) );
 634             premeshedNodes.insert( make_pair( (int)node->GetID(), node ) );
 635           }
 636           seg[nodeId] = ngId;                // ng node id
 637           SMESH_NodeXYZ nXYZ( node );
 638           gp_XY uv    = sprojector->ValueOfUV( nXYZ, tol ).XY();
 639
 640           // Compute the param (relative distance) of the node in relation the edge
 641           // fundamental for netgen working properly
 642           double dist2, param;
 643           if ( myCurveProjector->IsOnCurve( nXYZ, dist2, param ) )
 644           {
 645             seg.epgeominfo[ nodeId ].dist = param;
 646             seg.epgeominfo[ nodeId ].u    = uv.X();
 647             seg.epgeominfo[ nodeId ].v    = uv.Y();
 648             seg.epgeominfo[ nodeId ].edgenr = occEdgeIdFound;
 649           }
 650         }
 651
 652         if ( meshingEdge.Orientation() != TopAbs_FORWARD )
 653         {
 654           swap(seg[0], seg[1]);
 655           swap(seg.epgeominfo[0], seg.epgeominfo[1] );
 656         }
 657
 658         seg.edgenr = ngMesh->GetNSeg() + 1; // netgen segment id
 659         seg.si     = faceID;
 660         ngMesh->AddSegment( seg );
 661       } // end if edge is on the face
 662     } // end iteration on elements
 663
 664     // set parameters from _hypLengthFromEdges if needed
 665     if ( !_hypParameters && _hypLengthFromEdges && numOfEdges )
 666     {
 667       netgen::mparam.maxh = totalEdgeLenght / numOfEdges;
 668       if ( netgen::mparam.maxh < DBL_MIN )
 669         netgen::mparam.maxh = occgeom.GetBoundingBox().Diam();
 670
 671       occgeom.face_maxh = netgen::mparam.maxh;
 672       ngMesh->SetGlobalH ( netgen::mparam.maxh );
 673     }
 674     // end set parameters
 675
 676     ngMesh->CalcSurfacesOfNode();
 677     const int startWith = MESHCONST_MESHSURFACE;
 678     const int endWith   = MESHCONST_OPTSURFACE;
 679     err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh);
 680
 681     // Ng_Mesh * ngMeshptr = (Ng_Mesh*) ngLib._ngMesh;
 682     // int NetgenNodes = Ng_GetNP(ngMeshptr);
 683     // int NetgenSeg_2D = Ng_GetNSeg_2D(ngMeshptr);
 684     // int NetgenFaces = Ng_GetNSE(ngMeshptr);
 685     // std::cout << "\n";
 686     // std::cout << "Number of nodes, seg, faces, vols: " << NetgenNodes << ", " << NetgenSeg_2D << ", " << NetgenFaces << "\n";
 687     // std::cout << "err from netgen computation: " << err << "\n";
 688
 689     // ----------------------------------------------------
 690     // Fill the SMESHDS with the generated nodes and faces
 691     // ----------------------------------------------------
 692     nodeVec.clear();
 693     FillNodesAndElements( aMesh, helper, ngMesh, nodeVec, ng2smesh, newNetgenCoordinates, newNetgenElements, numberOfPremeshedNodes );
 694   } // Face iteration
 695
 696   return (bool) err;
 697 }
 698
 699 std::tuple<bool,bool> NETGENPlugin_NETGEN_2D_ONLY::SetParameteres( SMESH_Mesh& aMesh, const TopoDS_Shape& aShape,
 700                                                                   NETGENPlugin_Mesher& aMesher, netgen::Mesh * ngMeshes,
 701                                                                   netgen::OCCGeometry& occgeoComm, bool isSubMeshSupported )
 702 {
 703   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 704
 705   aMesher.SetParameters( _hypParameters ); // _hypParameters -> netgen::mparam
 706   if ( _hypMaxElementArea )
 707   {
 708     netgen::mparam.maxh = sqrt( 2. * _hypMaxElementArea->GetMaxArea() / sqrt(3.0) );
 709   }
 710   if ( _hypQuadranglePreference )
 711     netgen::mparam.quad = true;
 712
 713   // local size is common for all FACEs in aShape?
 714   const bool isCommonLocalSize = ( !_hypLengthFromEdges && !_hypMaxElementArea && netgen::mparam.uselocalh );
 715   const bool isDefaultHyp = ( !_hypLengthFromEdges && !_hypMaxElementArea && !_hypParameters );
 716
 717   if ( isCommonLocalSize ) // compute common local size in ngMeshes[0]
 718   {
 719     aMesher.PrepareOCCgeometry( occgeoComm, aShape, aMesh );//, meshedSM );
 720
 721     // local size set at MESHCONST_ANALYSE step depends on
 722     // minh, face_maxh, grading and curvaturesafety; find minh if not set by the user
 723     if ( !_hypParameters || netgen::mparam.minh < DBL_MIN )
 724     {
 725       if ( !_hypParameters )
 726         netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam() / 3.;
 727       netgen::mparam.minh = aMesher.GetDefaultMinSize( aShape, netgen::mparam.maxh );
 728     }
 729     // set local size depending on curvature and NOT closeness of EDGEs
 730 #ifdef NETGEN_V6
 731     const double factor = 2; //netgen::occparam.resthcloseedgefac;
 732 #else
 733     const double factor = netgen::occparam.resthcloseedgefac;
 734     netgen::occparam.resthcloseedgeenable = false;
 735     netgen::occparam.resthcloseedgefac = 1.0 + netgen::mparam.grading;
 736 #endif
 737     occgeoComm.face_maxh = netgen::mparam.maxh;
 738 #ifdef NETGEN_V6
 739     netgen::OCCParameters occparam;
 740     netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes, netgen::mparam, occparam );
 741 #else
 742     netgen::OCCSetLocalMeshSize( occgeoComm, *ngMeshes );
 743 #endif
 744     occgeoComm.emap.Clear();
 745     occgeoComm.vmap.Clear();
 746
 747     if ( isSubMeshSupported )
 748     {
 749       // set local size according to size of existing segments
 750       TopTools_IndexedMapOfShape edgeMap;
 751       TopExp::MapShapes( aMesh.GetShapeToMesh(), TopAbs_EDGE, edgeMap );
 752       for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
 753       {
 754         const TopoDS_Shape& edge = edgeMap( iE );
 755         if ( SMESH_Algo::isDegenerated( TopoDS::Edge( edge )))
 756           continue;
 757         SMESHDS_SubMesh* smDS = meshDS->MeshElements( edge );
 758         if ( !smDS ) continue;
 759         SMDS_ElemIteratorPtr segIt = smDS->GetElements();
 760         while ( segIt->more() )
 761         {
 762           const SMDS_MeshElement* seg = segIt->next();
 763           SMESH_TNodeXYZ n1 = seg->GetNode(0);
 764           SMESH_TNodeXYZ n2 = seg->GetNode(1);
 765           gp_XYZ p = 0.5 * ( n1 + n2 );
 766           netgen::Point3d pi(p.X(), p.Y(), p.Z());
 767           ngMeshes->RestrictLocalH( pi, factor * ( n1 - n2 ).Modulus() );
 768         }
 769       }
 770     }
 771     else
 772     {
 773       SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Edge);
 774       while ( iteratorElem->more() ) // loop on elements on a geom face
 775       {
 776         const SMDS_MeshElement* elem = iteratorElem->next();
 777         const SMDS_MeshNode* node0 = elem->GetNode( 0 );
 778         const SMDS_MeshNode* node1 = elem->GetNode( 1 );
 779         SMESH_NodeXYZ nXYZ0( node0 );
 780         SMESH_NodeXYZ nXYZ1( node1 );
 781         double segmentLength = ( nXYZ0 - nXYZ1 ).Modulus();
 782         gp_XYZ p = 0.5 * ( nXYZ0 + nXYZ1 );
 783         netgen::Point3d pi(p.X(), p.Y(), p.Z());
 784         ngMeshes->RestrictLocalH( pi, factor * segmentLength );
 785       }
 786     }
 787
 788     // set local size defined on shapes
 789     aMesher.SetLocalSize( occgeoComm, *ngMeshes );
 790     aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMeshes );
 791     try {
 792       ngMeshes->LoadLocalMeshSize( mparam.meshsizefilename );
 793     } catch (NgException & ex) {
 794       throw error( COMPERR_BAD_PARMETERS, ex.What() );
 795     }
 796   }
 797
 798   return std::make_tuple( isCommonLocalSize, isDefaultHyp );
 799 }
 800
 801 bool NETGENPlugin_NETGEN_2D_ONLY::ComputeMaxhOfFace( TopoDS_Face& Face, NETGENPlugin_Mesher& aMesher, TSideVector& wires,
 802                                                       netgen::OCCGeometry& occgeoComm, bool isDefaultHyp, bool isCommonLocalSize )
 803 {
 804   size_t nbWires = wires.size();
 805   if ( !_hypParameters )
 806   {
 807     double edgeLength = 0;
 808     if (_hypLengthFromEdges )
 809     {
 810       // compute edgeLength as an average segment length
 811       smIdType nbSegments = 0;
 812       for ( size_t iW = 0; iW < nbWires; ++iW )
 813       {
 814         edgeLength += wires[ iW ]->Length();
 815         nbSegments += wires[ iW ]->NbSegments();
 816       }
 817       if ( nbSegments )
 818         edgeLength /= double( nbSegments );
 819       netgen::mparam.maxh = edgeLength;
 820     }
 821     else if ( isDefaultHyp )
 822     {
 823       // set edgeLength by a longest segment
 824       double maxSeg2 = 0;
 825       for ( size_t iW = 0; iW < nbWires; ++iW )
 826       {
 827         const UVPtStructVec& points = wires[ iW ]->GetUVPtStruct();
 828         if ( points.empty() )
 829           return error( COMPERR_BAD_INPUT_MESH );
 830         gp_Pnt pPrev = SMESH_TNodeXYZ( points[0].node );
 831         for ( size_t i = 1; i < points.size(); ++i )
 832         {
 833           gp_Pnt p = SMESH_TNodeXYZ( points[i].node );
 834           maxSeg2 = Max( maxSeg2, p.SquareDistance( pPrev ));
 835           pPrev = p;
 836         }
 837       }
 838       edgeLength = sqrt( maxSeg2 ) * 1.05;
 839       netgen::mparam.maxh = edgeLength;
 840     }
 841     if ( netgen::mparam.maxh < DBL_MIN )
 842       netgen::mparam.maxh = occgeoComm.GetBoundingBox().Diam();
 843
 844     if ( !isCommonLocalSize )
 845     {
 846       netgen::mparam.minh = aMesher.GetDefaultMinSize( Face, netgen::mparam.maxh );
 847     }
 848   }
 849   return true;
 850 }
 851
 852 void NETGENPlugin_NETGEN_2D_ONLY::FillNodesAndElements( SMESH_Mesh& aMesh, SMESH_MesherHelper& helper, netgen::Mesh * ngMesh,
 853                                                         vector< const SMDS_MeshNode* >& nodeVec, map<int, const SMDS_MeshNode* >& ng2smesh,
 854                                                         std::map<int,std::vector<double>>& newNetgenCoordinates,
 855                                                         std::map<int,std::vector<smIdType>>& newNetgenElements, const int numberOfPremeshedNodes )
 856 {
 857   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 858
 859   int nbNodes = ngMesh->GetNP();
 860   int nbFaces = ngMesh->GetNSE();
 861   nodeVec.resize( nbNodes + 1, 0 );
 862   // map to index local node numeration to global numeration used by Remote mesher to write the final global resulting mesh
 863   std::map<smIdType,smIdType> local2globalMap;
 864
 865   smIdType myNewCoordinateCounter = newNetgenCoordinates.size() > 0 ? newNetgenCoordinates.rbegin()->first + 1: numberOfPremeshedNodes+1;
 866   int myNewFaceCounter = newNetgenElements.size() > 0 ? newNetgenElements.rbegin()->first + 1 : 1;
 867   // add nodes
 868   for ( int ngID = 1; ngID <= nbNodes; ++ngID )
 869   {
 870     const MeshPoint& ngPoint = ngMesh->Point( ngID );
 871     // Check if ngPoint is not already present because was in the premeshed mesh boundary
 872     if ( ng2smesh.count( ngID ) == 0 )
 873     {
 874       std::vector<double> netgenCoordinates = {ngPoint(0), ngPoint(1), ngPoint(2)};
 875       SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
 876       nodeVec[ ngID ] = node;
 877       newNetgenCoordinates.insert( make_pair( myNewCoordinateCounter, std::move(netgenCoordinates)) );
 878       local2globalMap.insert( std::make_pair( node->GetID(), myNewCoordinateCounter ) );
 879       myNewCoordinateCounter++;
 880     }
 881     else
 882     {
 883       nodeVec[ ngID ] = ng2smesh[ ngID ];
 884       local2globalMap.insert( std::make_pair( nodeVec[ ngID ]->GetID(), nodeVec[ ngID ]->GetID() ) );
 885     }
 886   }
 887   // create faces
 888   int i,j;
 889   vector<const SMDS_MeshNode*> nodes;
 890   for ( i = 1; i <= nbFaces ; ++i )
 891   {
 892     const Element2d& elem = ngMesh->SurfaceElement(i);
 893     nodes.resize( elem.GetNP() );
 894     for (j=1; j <= elem.GetNP(); ++j)
 895     {
 896       int pind = elem.PNum(j);
 897       if ( pind < 1 )
 898         break;
 899       nodes[ j-1 ] = nodeVec[ pind ];
 900     }
 901     if ( j > elem.GetNP() )
 902     {
 903       std::vector<smIdType> netgenCoordinates = { local2globalMap[nodes[0]->GetID()], local2globalMap[nodes[1]->GetID()], local2globalMap[nodes[2]->GetID()] };
 904       newNetgenElements.insert( std::make_pair( myNewFaceCounter, std::move( netgenCoordinates ) ) );
 905       helper.AddFace(nodes[0],nodes[1],nodes[2]);
 906       myNewFaceCounter++;
 907     }
 908   }
 909 }
 910
 911 void NETGENPlugin_NETGEN_2D_ONLY::FillNodesAndElements( SMESH_Mesh& aMesh, SMESH_MesherHelper& helper, netgen::Mesh * ngMesh, vector< const SMDS_MeshNode* >& nodeVec, int faceId )
 912 {
 913   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 914
 915   int nbNodes = ngMesh->GetNP();
 916   int nbFaces = ngMesh->GetNSE();
 917
 918   int nbInputNodes = (int) nodeVec.size()-1;
 919   nodeVec.resize( nbNodes+1, 0 );
 920
 921   // add nodes
 922   for ( int ngID = nbInputNodes + 1; ngID <= nbNodes; ++ngID )
 923   {
 924     const MeshPoint& ngPoint = ngMesh->Point( ngID );
 925     SMDS_MeshNode * node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
 926     nodeVec[ ngID ] = node;
 927   }
 928
 929   // create faces
 930   int i,j;
 931   vector<const SMDS_MeshNode*> nodes;
 932   for ( i = 1; i <= nbFaces ; ++i )
 933   {
 934     const Element2d& elem = ngMesh->SurfaceElement(i);
 935     nodes.resize( elem.GetNP() );
 936     for (j=1; j <= elem.GetNP(); ++j)
 937     {
 938       int pind = elem.PNum(j);
 939       if ( pind < 1 )
 940         break;
 941       nodes[ j-1 ] = nodeVec[ pind ];
 942       if ( nodes[ j-1 ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_3DSPACE )
 943       {
 944         const PointGeomInfo& pgi = elem.GeomInfoPi(j);
 945         meshDS->SetNodeOnFace( nodes[ j-1 ], faceId, pgi.u, pgi.v);
 946       }
 947     }
 948     if ( j > elem.GetNP() )
 949     {
 950       if ( elem.GetType() == TRIG )
 951         helper.AddFace(nodes[0],nodes[1],nodes[2]);
 952       else
 953         helper.AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
 954     }
 955   }
 956 }
 957
 958
 959 //=============================================================================
 960 /*!
 961  *Here we are going to use the NETGEN mesher
 962  */
 963 //=============================================================================
 964
 965 bool NETGENPlugin_NETGEN_2D_ONLY::Compute(SMESH_Mesh&         aMesh,
 966                                           const TopoDS_Shape& aShape)
 967 {
 968   netgen::multithread.terminate = 0;
 969   //netgen::multithread.task = "Surface meshing";
 970
 971   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 972   SMESH_MesherHelper helper(aMesh);
 973   helper.SetElementsOnShape( true );
 974
 975   NETGENPlugin_NetgenLibWrapper ngLib;
 976   ngLib._isComputeOk = false;
 977
 978   netgen::Mesh   ngMeshNoLocSize;
 979   netgen::Mesh * ngMeshes[2] = { (netgen::Mesh*) ngLib._ngMesh,  & ngMeshNoLocSize };
 980
 981   netgen::OCCGeometry occgeoComm;
 982
 983   // min / max sizes are set as follows:
 984   // if ( _hypParameters )
 985   //    min and max are defined by the user
 986   // else if ( _hypLengthFromEdges )
 987   //    min = aMesher.GetDefaultMinSize()
 988   //    max = average segment len of a FACE
 989   // else if ( _hypMaxElementArea )
 990   //    min = aMesher.GetDefaultMinSize()
 991   //    max = f( _hypMaxElementArea )
 992   // else
 993   //    min = aMesher.GetDefaultMinSize()
 994   //    max = max segment len of a FACE
 995
 996   NETGENPlugin_Mesher aMesher( &aMesh, aShape, /*isVolume=*/false);
 997   auto options = SetParameteres( aMesh, aShape, aMesher, ngMeshes[0], occgeoComm );
 998   const bool isCommonLocalSize  = std::get<0>( options );
 999   const bool isDefaultHyp       = std::get<1>( options );
1000   const bool toOptimize = _hypParameters ? _hypParameters->GetOptimize() : true;
1001
1002   netgen::mparam.uselocalh = toOptimize; // restore as it is used at surface optimization
1003
1004   // ==================
1005   // Loop on all FACEs
1006   // ==================
1007
1008   vector< const SMDS_MeshNode* > nodeVec;
1009
1010   TopExp_Explorer fExp( aShape, TopAbs_FACE );
1011   for ( int iF = 0; fExp.More(); fExp.Next(), ++iF )
1012   {
1013     TopoDS_Face F = TopoDS::Face( fExp.Current() /*.Oriented( TopAbs_FORWARD )*/);
1014     int    faceID = meshDS->ShapeToIndex( F );
1015     SMESH_ComputeErrorPtr& faceErr = aMesh.GetSubMesh( F )->GetComputeError();
1016
1017     _quadraticMesh = helper.IsQuadraticSubMesh( F );
1018     const bool ignoreMediumNodes = _quadraticMesh;
1019
1020     // build viscous layers if required
1021     if ( F.Orientation() != TopAbs_FORWARD &&
1022          F.Orientation() != TopAbs_REVERSED )
1023       F.Orientation( TopAbs_FORWARD ); // avoid pb with TopAbs_INTERNAL
1024     SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
1025     if ( !proxyMesh )
1026       continue;
1027
1028     // ------------------------
1029     // get all EDGEs of a FACE
1030     // ------------------------
1031     TSideVector wires =
1032       StdMeshers_FaceSide::GetFaceWires( F, aMesh, ignoreMediumNodes, faceErr, &helper, proxyMesh );
1033     if ( faceErr && !faceErr->IsOK() )
1034       continue;
1035     size_t nbWires = wires.size();
1036     if ( nbWires == 0 )
1037     {
1038       faceErr.reset
1039         ( new SMESH_ComputeError
1040           ( COMPERR_ALGO_FAILED, "Problem in StdMeshers_FaceSide::GetFaceWires()" ));
1041       continue;
1042     }
1043     if ( wires[0]->NbSegments() < 3 ) // ex: a circle with 2 segments
1044     {
1045       faceErr.reset
1046         ( new SMESH_ComputeError
1047           ( COMPERR_BAD_INPUT_MESH, SMESH_Comment("Too few segments: ")<<wires[0]->NbSegments()) );
1048       continue;
1049     }
1050
1051     // ----------------------
1052     // compute maxh of a FACE
1053     // ----------------------
1054
1055     bool setMaxh = ComputeMaxhOfFace( F, aMesher, wires, occgeoComm, isDefaultHyp, isCommonLocalSize );
1056     if (!setMaxh)
1057       return setMaxh;
1058
1059     // prepare occgeom
1060     netgen::OCCGeometry occgeom;
1061     occgeom.shape = F;
1062     occgeom.fmap.Add( F );
1063     occgeom.CalcBoundingBox();
1064     occgeom.facemeshstatus.SetSize(1);
1065     occgeom.facemeshstatus = 0;
1066     occgeom.face_maxh_modified.SetSize(1);
1067     occgeom.face_maxh_modified = 0;
1068     occgeom.face_maxh.SetSize(1);
1069     occgeom.face_maxh = netgen::mparam.maxh;
1070
1071     // -------------------------
1072     // Fill netgen mesh
1073     // -------------------------
1074
1075     // MESHCONST_ANALYSE step may lead to a failure, so we make an attempt
1076     // w/o MESHCONST_ANALYSE at the second loop
1077     int err = 0;
1078     enum { LOC_SIZE, NO_LOC_SIZE };
1079     int iLoop = isCommonLocalSize ? 0 : 1;
1080     for ( ; iLoop < 2; iLoop++ )
1081     {
1082       //bool isMESHCONST_ANALYSE = false;
1083       InitComputeError();
1084
1085       netgen::Mesh * ngMesh = ngMeshes[ iLoop ];
1086       ngMesh->DeleteMesh();
1087
1088       if ( iLoop == NO_LOC_SIZE )
1089       {
1090         ngMesh->SetGlobalH ( mparam.maxh );
1091         ngMesh->SetMinimalH( mparam.minh );
1092         Box<3> bb = occgeom.GetBoundingBox();
1093         bb.Increase (bb.Diam()/10);
1094         ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparam.grading);
1095         aMesher.SetLocalSize( occgeom, *ngMesh );
1096         aMesher.SetLocalSizeForChordalError( occgeoComm, *ngMesh );
1097         try {
1098           ngMesh->LoadLocalMeshSize( mparam.meshsizefilename );
1099         } catch (NgException & ex) {
1100           return error( COMPERR_BAD_PARMETERS, ex.What() );
1101         }
1102       }
1103
1104       nodeVec.clear();
1105       faceErr = aMesher.AddSegmentsToMesh( *ngMesh, occgeom, wires, helper, nodeVec,
1106                                            /*overrideMinH=*/!_hypParameters);
1107       if ( faceErr && !faceErr->IsOK() )
1108         break;
1109
1110       //if ( !isCommonLocalSize )
1111       //limitSize( ngMesh, mparam.maxh * 0.8);
1112
1113       // -------------------------
1114       // Generate surface mesh
1115       // -------------------------
1116
1117       const int startWith = MESHCONST_MESHSURFACE;
1118       const int endWith   = toOptimize ? MESHCONST_OPTSURFACE : MESHCONST_MESHSURFACE;
1119
1120       SMESH_Comment str;
1121       try {
1122         OCC_CATCH_SIGNALS;
1123
1124         err = ngLib.GenerateMesh(occgeom, startWith, endWith, ngMesh);
1125
1126         if ( netgen::multithread.terminate )
1127           return false;
1128         if ( err )
1129           str << "Error in netgen::OCCGenerateMesh() at " << netgen::multithread.task;
1130       }
1131       catch (Standard_Failure& ex)
1132       {
1133         err = 1;
1134         str << "Exception in  netgen::OCCGenerateMesh()"
1135             << " at " << netgen::multithread.task
1136             << ": " << ex.DynamicType()->Name();
1137         if ( ex.GetMessageString() && strlen( ex.GetMessageString() ))
1138           str << ": " << ex.GetMessageString();
1139       }
1140       catch (...) {
1141         err = 1;
1142         str << "Exception in  netgen::OCCGenerateMesh()"
1143             << " at " << netgen::multithread.task;
1144       }
1145       if ( err )
1146       {
1147         if ( aMesher.FixFaceMesh( occgeom, *ngMesh, 1 ))
1148           break;
1149         if ( iLoop == LOC_SIZE )
1150         {
1151           netgen::mparam.minh = netgen::mparam.maxh;
1152           netgen::mparam.maxh = 0;
1153           for ( size_t iW = 0; iW < wires.size(); ++iW )
1154           {
1155             StdMeshers_FaceSidePtr wire = wires[ iW ];
1156             const vector<UVPtStruct>& uvPtVec = wire->GetUVPtStruct();
1157             for ( size_t iP = 1; iP < uvPtVec.size(); ++iP )
1158             {
1159               SMESH_TNodeXYZ   p( uvPtVec[ iP ].node );
1160               netgen::Point3d np( p.X(),p.Y(),p.Z());
1161               double segLen = p.Distance( uvPtVec[ iP-1 ].node );
1162               double   size = ngMesh->GetH( np );
1163               netgen::mparam.minh = Min( netgen::mparam.minh, size );
1164               netgen::mparam.maxh = Max( netgen::mparam.maxh, segLen );
1165             }
1166           }
1167           //cerr << "min " << netgen::mparam.minh << " max " << netgen::mparam.maxh << endl;
1168           netgen::mparam.minh *= 0.9;
1169           netgen::mparam.maxh *= 1.1;
1170           continue;
1171         }
1172         else
1173         {
1174           faceErr.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, str ));
1175         }
1176       }
1177
1178       // ----------------------------------------------------
1179       // Fill the SMESHDS with the generated nodes and faces
1180       // ----------------------------------------------------
1181       FillNodesAndElements( aMesh, helper, ngMesh, nodeVec, faceID );
1182
1183       break;
1184     } // two attempts
1185   } // loop on FACEs
1186
1187   return true;
1188 }
1189
1190 void NETGENPlugin_NETGEN_2D_ONLY::CancelCompute()
1191 {
1192   SMESH_Algo::CancelCompute();
1193   netgen::multithread.terminate = 1;
1194 }
1195
1196 //================================================================================
1197 /*!
1198  * \brief Return progress of Compute() [0.,1]
1199  */
1200 //================================================================================
1201
1202 double NETGENPlugin_NETGEN_2D_ONLY::GetProgress() const
1203 {
1204   return -1;
1205   // const char* task1 = "Surface meshing";
1206   // //const char* task2 = "Optimizing surface";
1207   // double& progress = const_cast<NETGENPlugin_NETGEN_2D_ONLY*>( this )->_progress;
1208   // if ( _progressByTic < 0. &&
1209   //      strncmp( netgen::multithread.task, task1, 3 ) == 0 )
1210   // {
1211   //   progress = Min( 0.25, SMESH_Algo::GetProgressByTic() ); // [0, 0.25]
1212   // }
1213   // else //if ( strncmp( netgen::multithread.task, task2, 3 ) == 0)
1214   // {
1215   //   if ( _progressByTic < 0 )
1216   //   {
1217   //     NETGENPlugin_NETGEN_2D_ONLY* me = (NETGENPlugin_NETGEN_2D_ONLY*) this;
1218   //     me->_progressByTic = 0.25 / (_progressTic+1);
1219   //   }
1220   //   const_cast<NETGENPlugin_NETGEN_2D_ONLY*>( this )->_progressTic++;
1221   //   progress = Max( progress, _progressByTic * _progressTic );
1222   // }
1223   // //cout << netgen::multithread.task << " " << _progressTic << endl;
1224   // return Min( progress, 0.99 );
1225 }
1226
1227 //=============================================================================
1228 /*!
1229  *
1230  */
1231 //=============================================================================
1232
1233 bool NETGENPlugin_NETGEN_2D_ONLY::Evaluate(SMESH_Mesh& aMesh,
1234                                            const TopoDS_Shape& aShape,
1235                                            MapShapeNbElems& aResMap)
1236 {
1237   TopoDS_Face F = TopoDS::Face(aShape);
1238   if(F.IsNull())
1239     return false;
1240
1241   // collect info from edges
1242   smIdType nb0d = 0, nb1d = 0;
1243   bool IsQuadratic = false;
1244   bool IsFirst = true;
1245   double fullLen = 0.0;
1246   TopTools_MapOfShape tmpMap;
1247   for (TopExp_Explorer exp(F, TopAbs_EDGE); exp.More(); exp.Next()) {
1248     TopoDS_Edge E = TopoDS::Edge(exp.Current());
1249     if( tmpMap.Contains(E) )
1250       continue;
1251     tmpMap.Add(E);
1252     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
1253     MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
1254     if( anIt==aResMap.end() ) {
1255       SMESH_subMesh *sm = aMesh.GetSubMesh(F);
1256       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1257       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
1258       return false;
1259     }
1260     std::vector<smIdType> aVec = (*anIt).second;
1261     nb0d += aVec[SMDSEntity_Node];
1262     nb1d += std::max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
1263     double aLen = SMESH_Algo::EdgeLength(E);
1264     fullLen += aLen;
1265     if(IsFirst) {
1266       IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
1267       IsFirst = false;
1268     }
1269   }
1270   tmpMap.Clear();
1271
1272   // compute edge length
1273   double ELen = 0;
1274   if (( _hypLengthFromEdges ) || ( !_hypLengthFromEdges && !_hypMaxElementArea )) {
1275     if ( nb1d > 0 )
1276       ELen = fullLen / double( nb1d );
1277   }
1278   if ( _hypMaxElementArea ) {
1279     double maxArea = _hypMaxElementArea->GetMaxArea();
1280     ELen = sqrt(2. * maxArea/sqrt(3.0));
1281   }
1282   GProp_GProps G;
1283   BRepGProp::SurfaceProperties(F,G);
1284   double anArea = G.Mass();
1285
1286   const int hugeNb = numeric_limits<int>::max()/10;
1287   if ( anArea / hugeNb > ELen*ELen )
1288   {
1289     SMESH_subMesh *sm = aMesh.GetSubMesh(F);
1290     SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1291     smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated.\nToo small element length",this));
1292     return false;
1293   }
1294   smIdType nbFaces = (smIdType) ( anArea / ( ELen*ELen*sqrt(3.) / 4 ) );
1295   smIdType nbNodes = (smIdType) ( ( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
1296   std::vector<smIdType> aVec(SMDSEntity_Last);
1297   for(smIdType i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i]=0;
1298   if( IsQuadratic ) {
1299     aVec[SMDSEntity_Node] = nbNodes;
1300     aVec[SMDSEntity_Quad_Triangle] = nbFaces;
1301   }
1302   else {
1303     aVec[SMDSEntity_Node] = nbNodes;
1304     aVec[SMDSEntity_Triangle] = nbFaces;
1305   }
1306   SMESH_subMesh *sm = aMesh.GetSubMesh(F);
1307   aResMap.insert(std::make_pair(sm,aVec));
1308
1309   return true;
1310 }