src/BRepMesh/BRepMesh_FastDiscretFace.cxx

   1 // Created by: Ekaterina SMIRNOVA
   2 // Copyright (c) 2008-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 #include <BRepMesh_FastDiscretFace.hxx>
  16
  17 #include <BRepMesh_PairOfPolygon.hxx>
  18 #include <BRepMesh_ShapeTool.hxx>
  19 #include <Poly_PolygonOnTriangulation.hxx>
  20 #include <Poly_Triangulation.hxx>
  21
  22 #include <BRepAdaptor_Surface.hxx>
  23 #include <BRepAdaptor_HSurface.hxx>
  24 #include <BRepAdaptor_Curve.hxx>
  25 #include <Adaptor3d_IsoCurve.hxx>
  26
  27 #include <BRep_ListIteratorOfListOfPointRepresentation.hxx>
  28 #include <BRep_PointRepresentation.hxx>
  29 #include <BRep_TVertex.hxx>
  30 #include <BRep_Tool.hxx>
  31
  32 #include <GeomLib.hxx>
  33 #include <Geom_Surface.hxx>
  34 #include <Geom_BSplineSurface.hxx>
  35 #include <GCPnts_TangentialDeflection.hxx>
  36 #include <GCPnts_AbscissaPoint.hxx>
  37
  38 #include <Standard_ErrorHandler.hxx>
  39 #include <Standard_Failure.hxx>
  40 #include <SortTools_ShellSortOfReal.hxx>
  41 #include <TCollection_CompareOfReal.hxx>
  42 #include <TColStd_Array1OfReal.hxx>
  43 #include <TColStd_ListOfInteger.hxx>
  44 #include <TColStd_SequenceOfReal.hxx>
  45 #include <TColStd_Array1OfInteger.hxx>
  46 #include <TColStd_HArray1OfReal.hxx>
  47 #include <TColgp_Array1OfPnt2d.hxx>
  48 #include <TopTools_DataMapOfShapeReal.hxx>
  49
  50 #include <TopExp_Explorer.hxx>
  51 #include <TopoDS.hxx>
  52 #include <TopoDS_Vertex.hxx>
  53 #include <TopExp.hxx>
  54
  55 #include <NCollection_Map.hxx>
  56 #include <Bnd_Box2d.hxx>
  57
  58 #include <algorithm>
  59
  60 IMPLEMENT_STANDARD_HANDLE (BRepMesh_FastDiscretFace, Standard_Transient)
  61 IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_FastDiscretFace, Standard_Transient)
  62
  63 static Standard_Real FUN_CalcAverageDUV(TColStd_Array1OfReal& P, const Standard_Integer PLen)
  64 {
  65   Standard_Integer i, j, n = 0;
  66   Standard_Real p, result = 0.;
  67
  68   for(i = 1; i <= PLen; i++)
  69   {
  70     // Sort
  71     for(j = i + 1; j <= PLen; j++)
  72     {
  73       if(P(i) > P(j))
  74       {
  75         p = P(i);
  76         P(i) = P(j);
  77         P(j) = p;
  78       }
  79     }
  80     // Accumulate
  81     if (i != 1)
  82     {
  83       p = Abs(P(i) - P(i-1));
  84       if(p > 1.e-7)
  85       {
  86         result += p;
  87         n++;
  88       }
  89     }
  90   }
  91   return (n? (result / (Standard_Real) n) : -1.);
  92 }
  93
  94 namespace
  95 {
  96
  97   //! Auxiliary class used to extract geometrical parameters of fixed TopoDS_Vertex.
  98   class FixedVExplorer
  99   {
 100   public:
 101
 102     DEFINE_STANDARD_ALLOC
 103
 104     FixedVExplorer(const TopoDS_Vertex& theVertex)
 105       : myVertex(theVertex)
 106     {
 107     }
 108
 109     const TopoDS_Vertex& Vertex() const
 110     {
 111       return myVertex;
 112     }
 113
 114     Standard_Boolean IsSameUV() const
 115     {
 116       return Standard_False;
 117     }
 118
 119     TopoDS_Vertex SameVertex() const
 120     {
 121       return TopoDS_Vertex();
 122     }
 123
 124     gp_Pnt Point() const
 125     {
 126       return BRep_Tool::Pnt(myVertex);
 127     }
 128
 129   private:
 130
 131     void operator =(const FixedVExplorer& /*theOther*/)
 132     {
 133     }
 134
 135   private:
 136     const TopoDS_Vertex& myVertex;
 137   };
 138 }
 139
 140
 141 //=======================================================================
 142 //function : BRepMesh_FastDiscretFace
 143 //purpose  :
 144 //=======================================================================
 145 BRepMesh_FastDiscretFace::BRepMesh_FastDiscretFace(
 146   const Standard_Real    theAngle,
 147   const Standard_Real    theMinSize,
 148   const Standard_Boolean isInternalVerticesMode,
 149   const Standard_Boolean isControlSurfaceDeflection)
 150 : myAngle(theAngle),
 151   myInternalVerticesMode(isInternalVerticesMode),
 152   myMinSize(theMinSize),
 153   myIsControlSurfaceDeflection(isControlSurfaceDeflection)
 154 {
 155 }
 156
 157 //=======================================================================
 158 //function : Perform
 159 //purpose  :
 160 //=======================================================================
 161 void BRepMesh_FastDiscretFace::Perform(const Handle(BRepMesh_FaceAttribute)& theAttribute)
 162 {
 163   add(theAttribute);
 164   commitSurfaceTriangulation();
 165 }
 166
 167 //=======================================================================
 168 //function : initDataStructure
 169 //purpose  :
 170 //=======================================================================
 171 void BRepMesh_FastDiscretFace::initDataStructure()
 172 {
 173   const Standard_Real aTolU = myAttribute->ToleranceU();
 174   const Standard_Real aTolV = myAttribute->ToleranceV();
 175   const Standard_Real uCellSize = 14.0 * aTolU;
 176   const Standard_Real vCellSize = 14.0 * aTolV;
 177
 178   const Standard_Real deltaX = myAttribute->GetDeltaX();
 179   const Standard_Real deltaY = myAttribute->GetDeltaY();
 180
 181   Handle(NCollection_IncAllocator) aAllocator =
 182     new NCollection_IncAllocator(BRepMesh::MEMORY_BLOCK_SIZE_HUGE);
 183   myStructure = new BRepMesh_DataStructureOfDelaun(aAllocator);
 184   myStructure->Data()->SetCellSize ( uCellSize / deltaX, vCellSize / deltaY);
 185   myStructure->Data()->SetTolerance( aTolU     / deltaX, aTolV     / deltaY);
 186
 187   myAttribute->ChangeStructure() = myStructure;
 188   myAttribute->ChangeSurfacePoints() = new BRepMesh::DMapOfIntegerPnt;
 189   myAttribute->ChangeSurfaceVertices()= new BRepMesh::DMapOfVertexInteger;
 190
 191   // Check the necessity to fill the map of parameters
 192   const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
 193   GeomAbs_SurfaceType thetype = gFace->GetType();
 194   const Standard_Boolean isBSpline = (thetype == GeomAbs_BezierSurface ||
 195                                       thetype == GeomAbs_BSplineSurface);
 196   const Standard_Boolean useUVParam = (thetype == GeomAbs_Torus ||isBSpline);
 197
 198   myUParam.Clear();
 199   myVParam.Clear();
 200
 201   // essai de determination de la longueur vraie:
 202   // akm (bug OCC16) : We must calculate these measures in non-singular
 203   //     parts of face. Let`s try to compute average value of three
 204   //     (umin, (umin+umax)/2, umax), and respectively for v.
 205   //                 vvvvv
 206   //Standard_Real longu = 0.0, longv = 0.0; //, last , first;
 207   //gp_Pnt P11, P12, P21, P22, P31, P32;
 208   BRepMesh::HVectorOfVertex& aBoundaryNodes = myAttribute->ChangeMeshNodes();
 209   BRepMesh::VectorOfVertex::Iterator aNodesIt(*aBoundaryNodes);
 210   for (; aNodesIt.More(); aNodesIt.Next())
 211   {
 212     BRepMesh_Vertex& aNode = aNodesIt.ChangeValue();
 213     gp_XY aPnt2d = aNode.Coord();
 214
 215     if (useUVParam)
 216     {
 217       myUParam.Add(aPnt2d.X());
 218       myVParam.Add(aPnt2d.Y());
 219     }
 220
 221     aNode.ChangeCoord() = myAttribute->Scale(aPnt2d, Standard_True);
 222     myStructure->AddNode(aNode, Standard_True);
 223   }
 224   aBoundaryNodes.Nullify();
 225
 226   if (isBSpline)
 227   {
 228     const Standard_Real aRange[2][2] = {
 229       {myAttribute->GetUMin(), myAttribute->GetUMax()},
 230       {myAttribute->GetVMin(), myAttribute->GetVMax()}
 231     };
 232
 233     const GeomAbs_Shape aContinuity = GeomAbs_CN;
 234     for (Standard_Integer i = 0; i < 2; ++i)
 235     {
 236       const Standard_Boolean isU = (i == 0);
 237       const Standard_Integer aIntervalsNb = isU ?
 238         gFace->NbUIntervals(aContinuity) :
 239         gFace->NbVIntervals(aContinuity);
 240
 241       BRepMesh::IMapOfReal& aParams = isU ? myUParam : myVParam;
 242       if (aIntervalsNb < aParams.Size())
 243         continue;
 244
 245       TColStd_Array1OfReal aIntervals(1, aIntervalsNb + 1);
 246       if (isU)
 247         gFace->UIntervals(aIntervals, aContinuity);
 248       else
 249         gFace->VIntervals(aIntervals, aContinuity);
 250
 251       for (Standard_Integer j = 1; j <= aIntervals.Upper(); ++j)
 252       {
 253         const Standard_Real aParam = aIntervals(j);
 254         if (aParam > aRange[i][0] && aParam < aRange[i][1])
 255           aParams.Add(aParam);
 256       }
 257     }
 258   }
 259
 260   ////////////////////////////////////////////////////////////
 261   //add internal vertices after self-intersection check
 262   if ( myInternalVerticesMode )
 263   {
 264     TopExp_Explorer anExplorer(myAttribute->Face(), TopAbs_VERTEX, TopAbs_EDGE);
 265     for ( ; anExplorer.More(); anExplorer.Next() )
 266       add(TopoDS::Vertex(anExplorer.Current()));
 267   }
 268
 269   const BRepMesh::HDMapOfShapePairOfPolygon& aEdges = myAttribute->ChangeInternalEdges();
 270   TopExp_Explorer aWireIt(myAttribute->Face(), TopAbs_WIRE);
 271   for (; aWireIt.More(); aWireIt.Next())
 272   {
 273     TopExp_Explorer aEdgeIt(aWireIt.Current(), TopAbs_EDGE);
 274     for (; aEdgeIt.More(); aEdgeIt.Next())
 275     {
 276       const TopoDS_Edge& aEdge = TopoDS::Edge(aEdgeIt.Current());
 277       BRepMesh_PairOfPolygon aPair;
 278       if (!aEdges->Find(aEdge, aPair))
 279         continue;
 280
 281       TopAbs_Orientation aOri = aEdge.Orientation();
 282       const Handle(Poly_PolygonOnTriangulation)& aPolygon =
 283         aOri == TopAbs_REVERSED ? aPair.Last() : aPair.First();
 284
 285       const TColStd_Array1OfInteger& aIndices = aPolygon->Nodes();
 286       const Standard_Integer aNodesNb = aPolygon->NbNodes();
 287
 288       Standard_Integer aPrevId = aIndices(1);
 289       for (Standard_Integer i = 2; i <= aNodesNb; ++i)
 290       {
 291         const Standard_Integer aCurId = aIndices(i);
 292         addLinkToMesh(aPrevId, aCurId, aOri);
 293         aPrevId = aCurId;
 294       }
 295     }
 296   }
 297 }
 298
 299 //=======================================================================
 300 //function : addLinkToMesh
 301 //purpose  :
 302 //=======================================================================
 303 void BRepMesh_FastDiscretFace::addLinkToMesh(
 304   const Standard_Integer   theFirstNodeId,
 305   const Standard_Integer   theLastNodeId,
 306   const TopAbs_Orientation theOrientation)
 307 {
 308   if (theOrientation == TopAbs_FORWARD)
 309     myStructure->AddLink(BRepMesh_Edge(theFirstNodeId, theLastNodeId, BRepMesh_Frontier));
 310   else if (theOrientation == TopAbs_REVERSED)
 311     myStructure->AddLink(BRepMesh_Edge(theLastNodeId, theFirstNodeId, BRepMesh_Frontier));
 312   else if (theOrientation == TopAbs_INTERNAL)
 313     myStructure->AddLink(BRepMesh_Edge(theFirstNodeId, theLastNodeId, BRepMesh_Fixed));
 314 }
 315
 316 //=======================================================================
 317 //function : Add
 318 //purpose  :
 319 //=======================================================================
 320 void BRepMesh_FastDiscretFace::add(const Handle(BRepMesh_FaceAttribute)& theAttribute)
 321 {
 322   if (!theAttribute->IsValid() || theAttribute->ChangeMeshNodes()->IsEmpty())
 323     return;
 324
 325   myAttribute = theAttribute;
 326   initDataStructure();
 327
 328   BRepMesh::HIMapOfInteger& aVertexEdgeMap = myAttribute->ChangeVertexEdgeMap();
 329   Standard_Integer nbVertices = aVertexEdgeMap->Extent();
 330   BRepMesh::Array1OfInteger tabvert_corr(1, nbVertices);
 331   for ( Standard_Integer i = 1; i <= nbVertices; ++i )
 332     tabvert_corr(i) = i;
 333
 334   BRepMesh_Delaun trigu(myStructure, tabvert_corr);
 335
 336   //removed all free edges from triangulation
 337   const Standard_Integer nbLinks = myStructure->NbLinks();
 338   for( Standard_Integer i = 1; i <= nbLinks; i++ )
 339   {
 340     if( myStructure->ElementsConnectedTo(i).Extent() < 1 )
 341     {
 342       BRepMesh_Edge& anEdge = (BRepMesh_Edge&)trigu.GetEdge(i);
 343       if ( anEdge.Movability() == BRepMesh_Deleted )
 344         continue;
 345
 346       anEdge.SetMovability(BRepMesh_Free);
 347       myStructure->RemoveLink(i);
 348     }
 349   }
 350
 351   const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
 352   GeomAbs_SurfaceType thetype = gFace->GetType();
 353
 354   Standard_Boolean rajout =
 355     (thetype == GeomAbs_Sphere || thetype == GeomAbs_Torus);
 356
 357   // Check the necessity to fill the map of parameters
 358   const Standard_Boolean useUVParam = (thetype == GeomAbs_Torus         ||
 359                                        thetype == GeomAbs_BezierSurface ||
 360                                        thetype == GeomAbs_BSplineSurface);
 361
 362   const Standard_Real umax = myAttribute->GetUMax();
 363   const Standard_Real umin = myAttribute->GetUMin();
 364   const Standard_Real vmax = myAttribute->GetVMax();
 365   const Standard_Real vmin = myAttribute->GetVMin();
 366
 367   Standard_Boolean isaline =
 368     ((umax - umin) < Precision::PConfusion() ||
 369      (vmax - vmin) < Precision::PConfusion());
 370
 371   Standard_Real aDef = -1;
 372   if ( !isaline && myStructure->ElementsOfDomain().Extent() > 0 )
 373   {
 374     BRepMesh::ListOfVertex aNewVertices;
 375     if (!rajout)
 376     {
 377       aDef = control(aNewVertices, trigu, Standard_True);
 378       rajout = (aDef > myAttribute->GetDefFace() || aDef < 0.);
 379     }
 380
 381     if (!rajout && useUVParam)
 382     {
 383       rajout = (myVParam.Extent() > 2 &&
 384         (gFace->IsUClosed() || gFace->IsVClosed()));
 385     }
 386
 387     if (rajout)
 388     {
 389       insertInternalVertices(aNewVertices, trigu);
 390
 391       //control internal points
 392       if (myIsControlSurfaceDeflection)
 393         aDef = control(aNewVertices, trigu, Standard_False);
 394     }
 395   }
 396
 397   //modify myStructure back
 398   BRepMesh::HVectorOfVertex& aMeshNodes = myStructure->Data()->ChangeVertices();
 399   for ( Standard_Integer i = 1; i <= myStructure->NbNodes(); ++i )
 400   {
 401     BRepMesh_Vertex& aNode = aMeshNodes->ChangeValue(i - 1);
 402     aNode.ChangeCoord() = myAttribute->Scale(aNode.Coord(), Standard_False);
 403
 404     const BRepMesh::ListOfInteger& alist = myStructure->LinksConnectedTo(i);
 405     // Register internal nodes used in triangulation
 406     if (!alist.IsEmpty() && aVertexEdgeMap->FindIndex(i) == 0)
 407       aVertexEdgeMap->Add(i);
 408   }
 409
 410   if (!(aDef < 0.))
 411     myAttribute->SetDefFace(aDef);
 412 }
 413
 414 //=======================================================================
 415 //function : addVerticesToMesh
 416 //purpose  :
 417 //=======================================================================
 418 Standard_Boolean BRepMesh_FastDiscretFace::addVerticesToMesh(
 419   const BRepMesh::ListOfVertex& theVertices,
 420   BRepMesh_Delaun&              theMeshBuilder)
 421 {
 422   if (theVertices.IsEmpty())
 423     return Standard_False;
 424
 425   BRepMesh::Array1OfVertexOfDelaun aArrayOfNewVertices(1, theVertices.Extent());
 426   BRepMesh::ListOfVertex::Iterator aVertexIt(theVertices);
 427   for (Standard_Integer aVertexId = 0; aVertexIt.More(); aVertexIt.Next())
 428     aArrayOfNewVertices(++aVertexId) = aVertexIt.Value();
 429
 430   theMeshBuilder.AddVertices(aArrayOfNewVertices);
 431   return Standard_True;
 432 }
 433
 434 //=======================================================================
 435 //function : insertInternalVertices
 436 //purpose  :
 437 //=======================================================================
 438 static void filterParameters(const BRepMesh::IMapOfReal& theParams,
 439                              const Standard_Real         theMinDist,
 440                              const Standard_Real         theFilterDist,
 441                              BRepMesh::SequenceOfReal&   theResult)
 442 {
 443   // Sort sequence of parameters
 444   BRepMesh::SequenceOfReal aParamTmp;
 445   Standard_Integer aParamLength = 1;
 446   const Standard_Integer anInitLen = theParams.Extent();
 447
 448   TColStd_Array1OfReal aParamArray(1, anInitLen);
 449   Standard_Integer j;
 450   for (j = 1; j <= anInitLen; j++)
 451     aParamArray(j) = theParams(j);
 452
 453   TCollection_CompareOfReal aCompare;
 454   SortTools_ShellSortOfReal::Sort(aParamArray, aCompare);
 455
 456   // mandadory pre-filtering using the first (minimal) filter value
 457   Standard_Real aP1, aP2;
 458   aP1 = aParamArray(1);
 459   aParamTmp.Append(aP1);
 460   for (j = 2; j <= anInitLen; j++)
 461   {
 462     aP2 = aParamArray(j);
 463     if ((aP2-aP1) > theMinDist)
 464     {
 465       aParamTmp.Append(aP2);
 466       aP1 = aP2;
 467       aParamLength++;
 468     }
 469   }
 470
 471   //add last point if required
 472   if(aParamArray(anInitLen)-theParams(aParamLength) > theMinDist)
 473   {
 474     aParamTmp.Append(aParamArray(anInitLen));
 475     aParamLength++;
 476   }
 477
 478   //perform filtering on series
 479   Standard_Real aLastAdded, aLastCandidate;
 480   Standard_Boolean isCandidateDefined = Standard_False;
 481   aLastAdded = aParamTmp.First();
 482   aLastCandidate = aLastAdded;
 483   theResult.Append(aParamTmp.First());
 484
 485   for(j=2;j<aParamTmp.Length();j++)
 486   {
 487     Standard_Real aVal = aParamTmp.Value(j);
 488     if(aVal-aLastAdded > theFilterDist)
 489     {
 490       //adds the parameter
 491       if(isCandidateDefined) {
 492         aLastAdded = aLastCandidate;
 493         isCandidateDefined = Standard_False;
 494         j--;
 495       }
 496       else
 497       {
 498         aLastAdded = aVal;
 499       }
 500       theResult.Append(aLastAdded);
 501       continue;
 502     }
 503
 504     aLastCandidate = aVal;
 505     isCandidateDefined = Standard_True;
 506   }
 507   theResult.Append(aParamTmp.Last());
 508
 509   if( theResult.Length() == 2 )
 510   {
 511     Standard_Real    dist  = theResult.Last() - theResult.First();
 512     Standard_Integer nbint = (Standard_Integer)((dist / theFilterDist) + 0.5);
 513
 514     if( nbint > 1 )
 515     {
 516       //Five points more is maximum
 517       if( nbint > 5 )
 518       {
 519         nbint = 5;
 520       }
 521
 522       Standard_Integer i;
 523       Standard_Real dU = dist / nbint;
 524       for( i = 1; i < nbint; i++ )
 525       {
 526         theResult.InsertAfter(i, theResult.First()+i*dU);
 527       }
 528     }
 529   }
 530 }
 531
 532 void BRepMesh_FastDiscretFace::insertInternalVertices(
 533   BRepMesh::ListOfVertex&  theNewVertices,
 534   BRepMesh_Delaun&         theMeshBuilder)
 535 {
 536   const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
 537   switch (gFace->GetType())
 538   {
 539   case GeomAbs_Sphere:
 540       insertInternalVerticesSphere(theNewVertices);
 541       break;
 542
 543   case GeomAbs_Cylinder:
 544       insertInternalVerticesCylinder(theNewVertices);
 545       break;
 546
 547   case GeomAbs_Cone:
 548     insertInternalVerticesCone(theNewVertices);
 549     break;
 550
 551   case GeomAbs_Torus:
 552     insertInternalVerticesTorus(theNewVertices);
 553     break;
 554
 555   case GeomAbs_BezierSurface:
 556   case GeomAbs_BSplineSurface:
 557     insertInternalVerticesBSpline(theNewVertices);
 558     break;
 559
 560   default:
 561     insertInternalVerticesOther(theNewVertices);
 562     break;
 563   }
 564
 565   addVerticesToMesh(theNewVertices, theMeshBuilder);
 566 }
 567
 568 //=======================================================================
 569 //function : insertInternalVerticesSphere
 570 //purpose  :
 571 //=======================================================================
 572 void BRepMesh_FastDiscretFace::insertInternalVerticesSphere(
 573   BRepMesh::ListOfVertex& theNewVertices)
 574 {
 575   Standard_Real aRange[] = {
 576     myAttribute->GetVMin(), myAttribute->GetVMax(),
 577     myAttribute->GetUMin(), myAttribute->GetUMax()
 578   };
 579
 580   gp_Sphere aSphere = myAttribute->Surface()->Sphere();
 581
 582   // Calculate parameters for iteration in V direction
 583   Standard_Real aStep = 0.7 * GCPnts_TangentialDeflection::ArcAngularStep(
 584     aSphere.Radius(), myAttribute->GetDefFace(), myAngle, myMinSize);
 585
 586   Standard_Real aDd[2] = {aStep, aStep};
 587   Standard_Real aPasMax[2] = {0., 0.};
 588   for (Standard_Integer i = 0; i < 2; ++i)
 589   {
 590     const Standard_Real aMax = aRange[2 * i + 1];
 591     const Standard_Real aDiff = aMax - aRange[2 * i + 0];
 592     aDd[i] = aDiff / ((Standard_Integer)(aDiff / aDd[i]) + 1);
 593     aPasMax[i] = aMax - Precision::PConfusion();
 594   }
 595
 596   const Standard_Real aHalfDu = aDd[1] * 0.5;
 597   Standard_Boolean Shift = Standard_False;
 598   Standard_Real aPasV = aRange[0] + aDd[0];
 599   for (; aPasV < aPasMax[0]; aPasV += aDd[0])
 600   {
 601     Shift = !Shift;
 602     const Standard_Real d = (Shift) ? aHalfDu : 0.;
 603     Standard_Real aPasU = aRange[2] + d;
 604     for (; aPasU < aPasMax[1]; aPasU += aDd[1])
 605     {
 606       tryToInsertAnalyticVertex(gp_Pnt2d(aPasU, aPasV), aSphere, theNewVertices);
 607     }
 608   }
 609 }
 610
 611 //=======================================================================
 612 //function : insertInternalVerticesCylinder
 613 //purpose  :
 614 //=======================================================================
 615 void BRepMesh_FastDiscretFace::insertInternalVerticesCylinder(
 616   BRepMesh::ListOfVertex& theNewVertices)
 617 {
 618   const Standard_Real umax = myAttribute->GetUMax();
 619   const Standard_Real umin = myAttribute->GetUMin();
 620   const Standard_Real vmax = myAttribute->GetVMax();
 621   const Standard_Real vmin = myAttribute->GetVMin();
 622
 623   gp_Cylinder aCylinder = myAttribute->Surface()->Cylinder();
 624   const Standard_Real aRadius = aCylinder.Radius();
 625
 626   // Calculate parameters for iteration in U direction
 627   Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep(
 628     aRadius, myAttribute->GetDefFace(), myAngle, myMinSize);
 629
 630   const Standard_Real su = umax - umin;
 631   const Standard_Integer nbU = (Standard_Integer)(su / Du);
 632   Du = su / (nbU + 1);
 633
 634   // Calculate parameters for iteration in V direction
 635   const Standard_Real sv = vmax - vmin;
 636   Standard_Integer nbV = (Standard_Integer)(nbU*sv / (su*aRadius));
 637   nbV = Min(nbV, 100 * nbU);
 638   Standard_Real Dv = sv / (nbV + 1);
 639
 640   Standard_Real pasu, pasv, pasvmax = vmax - Dv*0.5, pasumax = umax - Du*0.5;
 641   for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv)
 642   {
 643     for (pasu = umin + Du; pasu < pasumax; pasu += Du)
 644     {
 645       tryToInsertAnalyticVertex(gp_Pnt2d(pasu, pasv), aCylinder, theNewVertices);
 646     }
 647   }
 648 }
 649
 650 //=======================================================================
 651 //function : insertInternalVerticesCone
 652 //purpose  :
 653 //=======================================================================
 654 void BRepMesh_FastDiscretFace::insertInternalVerticesCone(
 655   BRepMesh::ListOfVertex& theNewVertices)
 656 {
 657   const Standard_Real umax = myAttribute->GetUMax();
 658   const Standard_Real umin = myAttribute->GetUMin();
 659   const Standard_Real vmax = myAttribute->GetVMax();
 660   const Standard_Real vmin = myAttribute->GetVMin();
 661
 662   gp_Cone aCone = myAttribute->Surface()->Cone();
 663   Standard_Real RefR = aCone.RefRadius();
 664   Standard_Real SAng = aCone.SemiAngle();
 665   Standard_Real aRadius = Max(Abs(RefR + vmin*Sin(SAng)), Abs(RefR + vmax*Sin(SAng)));
 666
 667   Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep(
 668     aRadius, myAttribute->GetDefFace(), myAngle, myMinSize);
 669
 670   Standard_Real Dv, pasu, pasv;
 671   Standard_Integer nbU = (Standard_Integer)((umax - umin) / Du);
 672   Standard_Integer nbV = (Standard_Integer)(nbU*(vmax - vmin) / ((umax - umin)*aRadius));
 673   Du = (umax - umin) / (nbU + 1);
 674   Dv = (vmax - vmin) / (nbV + 1);
 675
 676   Standard_Real pasvmax = vmax - Dv*0.5, pasumax = umax - Du*0.5;
 677   for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv)
 678   {
 679     for (pasu = umin + Du; pasu < pasumax; pasu += Du)
 680     {
 681       tryToInsertAnalyticVertex(gp_Pnt2d(pasu, pasv), aCone, theNewVertices);
 682     }
 683   }
 684 }
 685
 686 //=======================================================================
 687 //function : insertInternalVerticesTorus
 688 //purpose  :
 689 //=======================================================================
 690 void BRepMesh_FastDiscretFace::insertInternalVerticesTorus(
 691   BRepMesh::ListOfVertex& theNewVertices)
 692 {
 693   const Standard_Real umax     = myAttribute->GetUMax();
 694   const Standard_Real umin     = myAttribute->GetUMin();
 695   const Standard_Real vmax     = myAttribute->GetVMax();
 696   const Standard_Real vmin     = myAttribute->GetVMin();
 697   const Standard_Real deltaX   = myAttribute->GetDeltaX();
 698   const Standard_Real deltaY   = myAttribute->GetDeltaY();
 699   const Standard_Real aDefFace = myAttribute->GetDefFace();
 700
 701   gp_Torus T = myAttribute->Surface()->Torus();
 702
 703   Standard_Boolean insert;
 704   Standard_Integer i, j, ParamULength, ParamVLength;
 705   Standard_Real pp, pasu, pasv;
 706   Standard_Real r = T.MinorRadius(), R = T.MajorRadius();
 707
 708   BRepMesh::SequenceOfReal ParamU, ParamV;
 709
 710   Standard_Real oldDv = GCPnts_TangentialDeflection::ArcAngularStep(
 711     r, aDefFace, myAngle, myMinSize);
 712
 713   Standard_Real Dv = 0.9*oldDv; //TWOTHIRD * oldDv;
 714   Dv = oldDv;
 715
 716   Standard_Real Du;
 717   Standard_Integer nbV = Max((Standard_Integer)((vmax - vmin) / Dv), 2);
 718   Dv = (vmax - vmin) / (nbV + 1);
 719   Standard_Real ru = R + r;
 720   if (ru > 1.e-16)
 721   {
 722     Du = GCPnts_TangentialDeflection::ArcAngularStep(
 723       ru, aDefFace, myAngle, myMinSize);
 724
 725     Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
 726     if (aa < gp::Resolution())
 727       return;
 728     Du *= Min(oldDv, Du) / aa;
 729   }
 730   else Du = Dv;
 731
 732   Standard_Integer nbU = Max((Standard_Integer)((umax - umin) / Du), 2);
 733   nbU = Max(nbU, (int)(nbV*(umax - umin)*R / ((vmax - vmin)*r) / 5.));
 734   Du = (umax - umin) / (nbU + 1);
 735
 736   if (R < r)
 737   {
 738     // As the points of edges are returned.
 739     // in this case, the points are not representative.
 740
 741     //-- Choose DeltaX and DeltaY so that to avoid skipping points on the grid
 742     for (i = 0; i <= nbU; i++) ParamU.Append(umin + i* Du);
 743   }//R<r
 744   else //U if R > r
 745   {
 746     //--ofv: U
 747     // Number of mapped U parameters
 748     const Standard_Integer LenU = myUParam.Extent();
 749     // Fill array of U parameters
 750     TColStd_Array1OfReal Up(1, LenU);
 751     for (j = 1; j <= LenU; j++) Up(j) = myUParam(j);
 752
 753     // Calculate DU, leave array of parameters
 754     Standard_Real aDU = FUN_CalcAverageDUV(Up, LenU);
 755     aDU = Max(aDU, Abs(umax - umin) / (Standard_Real)nbU / 2.);
 756     Standard_Real dUstd = Abs(umax - umin) / (Standard_Real)LenU;
 757     if (aDU > dUstd) dUstd = aDU;
 758     // Add U parameters
 759     for (j = 1; j <= LenU; j++)
 760     {
 761       pp = Up(j);
 762       insert = Standard_True;
 763       ParamULength = ParamU.Length();
 764       for (i = 1; i <= ParamULength && insert; i++)
 765       {
 766         insert = (Abs(ParamU.Value(i) - pp) > (0.5*dUstd));
 767       }
 768       if (insert) ParamU.Append(pp);
 769     }
 770   }
 771
 772   //--ofv: V
 773   // Number of mapped V parameters
 774   const Standard_Integer LenV = myVParam.Extent();
 775   // Fill array of V parameters
 776   TColStd_Array1OfReal Vp(1, LenV);
 777   for (j = 1; j <= LenV; j++) Vp(j) = myVParam(j);
 778   // Calculate DV, sort array of parameters
 779   Standard_Real aDV = FUN_CalcAverageDUV(Vp, LenV);
 780   aDV = Max(aDV, Abs(vmax - vmin) / (Standard_Real)nbV / 2.);
 781
 782   Standard_Real dVstd = Abs(vmax - vmin) / (Standard_Real)LenV;
 783   if (aDV > dVstd) dVstd = aDV;
 784   // Add V parameters
 785   for (j = 1; j <= LenV; j++)
 786   {
 787     pp = Vp(j);
 788
 789     insert = Standard_True;
 790     ParamVLength = ParamV.Length();
 791     for (i = 1; i <= ParamVLength && insert; i++)
 792     {
 793       insert = (Abs(ParamV.Value(i) - pp) > (dVstd*2. / 3.));
 794     }
 795     if (insert) ParamV.Append(pp);
 796   }
 797
 798   Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
 799   Standard_Real uminnew = umin + deltaY*0.1;
 800   Standard_Real vminnew = vmin + deltaX*0.1;
 801   Standard_Real umaxnew = umax - deltaY*0.1;
 802   Standard_Real vmaxnew = vmax - deltaX*0.1;
 803
 804   for (i = 1; i <= Lu; i++)
 805   {
 806     pasu = ParamU.Value(i);
 807     if (pasu >= uminnew && pasu < umaxnew)
 808     {
 809       for (j = 1; j <= Lv; j++)
 810       {
 811         pasv = ParamV.Value(j);
 812         if (pasv >= vminnew && pasv < vmaxnew)
 813         {
 814           tryToInsertAnalyticVertex(gp_Pnt2d(pasu, pasv), T, theNewVertices);
 815         }
 816       }
 817     }
 818   }
 819 }
 820
 821 //=======================================================================
 822 //function : insertInternalVerticesBSpline
 823 //purpose  :
 824 //=======================================================================
 825 void BRepMesh_FastDiscretFace::insertInternalVerticesBSpline(
 826   BRepMesh::ListOfVertex& theNewVertices)
 827 {
 828   const Standard_Real aRange[2][2] = {
 829       { myAttribute->GetUMax(), myAttribute->GetUMin() },
 830       { myAttribute->GetVMax(), myAttribute->GetVMin() }
 831   };
 832
 833   const Standard_Real aDelta[2] = {
 834     myAttribute->GetDeltaX(),
 835     myAttribute->GetDeltaY()
 836   };
 837
 838   const Standard_Real                 aDefFace = myAttribute->GetDefFace();
 839   const Handle(BRepAdaptor_HSurface)& gFace    = myAttribute->Surface();
 840
 841   BRepMesh::SequenceOfReal aParams[2];
 842   for (Standard_Integer i = 0; i < 2; ++i)
 843   {
 844     Standard_Boolean isU = (i == 0);
 845     Standard_Real aRes = isU ?
 846       gFace->UResolution(aDefFace) :
 847       gFace->VResolution(aDefFace);
 848
 849     // Sort and filter sequence of parameters
 850     Standard_Real aMinDiff = Precision::PConfusion();
 851     if (aDelta[i] < 1.)
 852       aMinDiff /= aDelta[i];
 853
 854     aMinDiff = Max(myMinSize, aMinDiff);
 855
 856     Standard_Real aRangeDiff = aRange[i][0] - aRange[i][1];
 857     Standard_Real aDiffMaxLim = 0.1 * aRangeDiff;
 858     Standard_Real aDiffMinLim = Max(0.005 * aRangeDiff, 2. * aRes);
 859     Standard_Real aDiff = Max(myMinSize, Min(aDiffMaxLim, aDiffMinLim));
 860     filterParameters(isU ? myUParam : myVParam, aMinDiff, aDiff, aParams[i]);
 861   }
 862
 863   // check intermediate isolines
 864   Handle(Geom_Surface) aBSpline;
 865   GeomAbs_SurfaceType  aSurfType = gFace->GetType();
 866   if (aSurfType == GeomAbs_BezierSurface)
 867     aBSpline = gFace->Bezier();
 868   else if (aSurfType == GeomAbs_BSplineSurface)
 869     aBSpline = gFace->BSpline();
 870
 871   const BRepMesh::HClassifier& aClassifier = myAttribute->ChangeClassifier();
 872
 873   // precision for compare square distances
 874   const Standard_Real aPrecision   = Precision::Confusion();
 875   const Standard_Real aSqPrecision = aPrecision * aPrecision;
 876   for (Standard_Integer k = 0; k < 2; ++k)
 877   {
 878     BRepMesh::SequenceOfReal& aParams1 = aParams[k];
 879     BRepMesh::SequenceOfReal& aParams2 = aParams[(k + 1) % 2];
 880     const Standard_Boolean isU = (k == 0);
 881     Standard_Integer aStartIndex, aEndIndex;
 882     if (isU)
 883     {
 884       aStartIndex = 1;
 885       aEndIndex   = aParams1.Length();
 886     }
 887     else
 888     {
 889       aStartIndex = 2;
 890       aEndIndex   = aParams1.Length() - 1;
 891     }
 892
 893     for (Standard_Integer i = aStartIndex; i <= aEndIndex; ++i)
 894     {
 895       const Standard_Real aParam1 = aParams1(i);
 896       GeomAdaptor_Curve aIso(isU ?
 897         aBSpline->UIso(aParam1) : aBSpline->VIso(aParam1));
 898
 899       Standard_Real aPrevParam2 = aParams2(1);
 900       gp_Pnt aPrevPnt2 = aIso.Value(aPrevParam2);
 901       for (Standard_Integer j = 2; j <= aParams2.Length();)
 902       {
 903         Standard_Real aParam2 = aParams2(j);
 904         gp_Pnt aPnt2 = aIso.Value(aParam2);
 905         Standard_Real aMidParam = 0.5 * (aPrevParam2 + aParam2);
 906         gp_Pnt aMidPnt = aIso.Value(aMidParam);
 907
 908         // 23.03.2010 skl for OCC21645 - change precision for comparison
 909         Standard_Real aDist;
 910         if (aPrevPnt2.SquareDistance(aPnt2) > aSqPrecision)
 911         {
 912           gp_Lin aLin(aPrevPnt2, gp_Dir(gp_Vec(aPrevPnt2, aPnt2)));
 913           aDist = aLin.Distance(aMidPnt);
 914         }
 915         else
 916           aDist = aPrevPnt2.Distance(aMidPnt);
 917
 918         if (aDist > aDefFace && aDist > myMinSize)
 919         {
 920           // insertion
 921           aParams2.InsertBefore(j, aMidParam);
 922         }
 923         else
 924         {
 925           //put regular grig for normals
 926           gp_Pnt2d aStPnt1, aStPnt2;
 927           if (isU)
 928           {
 929             aStPnt1 = gp_Pnt2d(aParam1, aPrevParam2);
 930             aStPnt2 = gp_Pnt2d(aParam1, aMidParam);
 931           }
 932           else
 933           {
 934             aStPnt1 = gp_Pnt2d(aPrevParam2, aParam1);
 935             aStPnt2 = gp_Pnt2d(aMidParam,   aParam1);
 936           }
 937
 938           gp_Dir N1(0, 0, 1), N2(0, 0, 1);
 939           Standard_Boolean aSt1 = GeomLib::NormEstim(aBSpline, aStPnt1, aPrecision, N1);
 940           Standard_Boolean aSt2 = GeomLib::NormEstim(aBSpline, aStPnt2, aPrecision, N2);
 941
 942           Standard_Real aAngle = N2.Angle(N1);
 943           if (aSt1 < 1 && aSt2 < 1 && aAngle > myAngle)
 944           {
 945             Standard_Real aLen = GCPnts_AbscissaPoint::Length(aIso,
 946               aPrevParam2, aMidParam, aDefFace);
 947
 948             if (aLen > myMinSize)
 949             {
 950               // insertion
 951               aParams2.InsertBefore(j, aMidParam);
 952               continue;
 953             }
 954           }
 955
 956           aPrevParam2 = aParam2;
 957           aPrevPnt2   = aPnt2;
 958
 959           ++j;
 960         }
 961       }
 962     }
 963   }
 964
 965   // insert nodes of the regular grid
 966   for (Standard_Integer i = 1; i <= aParams[0].Length(); ++i)
 967   {
 968     const Standard_Real aParam1 = aParams[0].Value(i);
 969     for (Standard_Integer j = 1; j <= aParams[1].Length(); ++j)
 970     {
 971       gp_Pnt2d aPnt2d(aParam1, aParams[1].Value(j));
 972
 973       // Classify intersection point
 974       if (aClassifier->Perform(aPnt2d) == TopAbs_IN)
 975       {
 976         gp_Pnt aPnt;
 977         gFace->D0(aPnt2d.X(), aPnt2d.Y(), aPnt);
 978         insertVertex(aPnt, aPnt2d.Coord(), theNewVertices);
 979       }
 980     }
 981   }
 982 }
 983
 984 //=======================================================================
 985 //function : insertInternalVerticesOther
 986 //purpose  :
 987 //=======================================================================
 988 void BRepMesh_FastDiscretFace::insertInternalVerticesOther(
 989   BRepMesh::ListOfVertex& theNewVertices)
 990 {
 991   const Standard_Real aAngle = myAngle;//0.35;
 992   const Standard_Real aRange[2][2] = {
 993       { myAttribute->GetUMax(), myAttribute->GetUMin() },
 994       { myAttribute->GetVMax(), myAttribute->GetVMin() }
 995   };
 996
 997   const Standard_Real                 aDefFace = myAttribute->GetDefFace();
 998   const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
 999
1000   BRepMesh::SequenceOfReal aParams[2];
1001   const Standard_Integer aIsoPointsNb = 11;
1002   for (Standard_Integer k = 0; k < 2; ++k)
1003   {
1004     const Standard_Boolean isU = (k == 0);
1005     const GeomAbs_IsoType  aIsoType = isU ? GeomAbs_IsoU : GeomAbs_IsoV;
1006     const Standard_Integer aOtherParamsIndex = (k + 1) % 2;
1007     const Standard_Real (&aRange1)[2] = aRange[k];
1008     const Standard_Real (&aRange2)[2] = aRange[aOtherParamsIndex];
1009
1010     GCPnts_TangentialDeflection aDiscretIso[aIsoPointsNb];
1011     const Standard_Real aStepWidth = (aRange1[0] - aRange1[1]) / aIsoPointsNb;
1012
1013     // Find the most curved Iso.
1014     Standard_Integer aMaxIndex = 1, aMaxPointsNb = 0;
1015     for (Standard_Integer aIsoIt = 0; aIsoIt < aIsoPointsNb; ++aIsoIt)
1016     {
1017       Standard_Real aParam = aRange1[1] + aIsoIt * aStepWidth;
1018       Adaptor3d_IsoCurve aIso(gFace, aIsoType, aParam);
1019
1020       Standard_Real aFirstParam = Max(aRange2[1], aIso.FirstParameter());
1021       Standard_Real aLastParam  = Min(aRange2[0], aIso.LastParameter());
1022
1023       aDiscretIso[aIsoIt].Initialize(aIso, aFirstParam, aLastParam,
1024         aAngle, 0.7 * aDefFace, 2, Precision::PConfusion(), myMinSize);
1025
1026       const Standard_Integer aPointsNb = aDiscretIso[aIsoIt].NbPoints();
1027       if (aPointsNb > aMaxPointsNb)
1028       {
1029         aMaxPointsNb = aPointsNb;
1030         aMaxIndex    = aIsoIt;
1031       }
1032     }
1033
1034     BRepMesh::SequenceOfReal& aParams2 = aParams[aOtherParamsIndex];
1035     GCPnts_TangentialDeflection& aDIso = aDiscretIso[aMaxIndex];
1036     for (Standard_Integer i = 1; i <= aMaxPointsNb; ++i)
1037       aParams2.Append(aDIso.Parameter(i));
1038
1039     if (aParams2.Length() == 2)
1040       aParams2.InsertAfter(1, 0.5 * (aRange2[1] + aRange2[0]));
1041   }
1042
1043   Adaptor3d_IsoCurve aIsoV;
1044   aIsoV.Load(gFace);
1045
1046   const BRepMesh::HClassifier& aClassifier = myAttribute->ChangeClassifier();
1047   const Standard_Integer aUPointsNb = aParams[0].Length();
1048   const Standard_Integer aVPointsNb = aParams[1].Length();
1049   for (Standard_Integer i = 2; i < aVPointsNb; ++i)
1050   {
1051     const Standard_Real aV = aParams[1].Value(i);
1052     aIsoV.Load(GeomAbs_IsoV, aV);
1053
1054     for (Standard_Integer j = 2; j < aUPointsNb; ++j)
1055     {
1056       const Standard_Real aU = aParams[0].Value(j);
1057
1058       const gp_Pnt2d aNewPoint(aU, aV);
1059       if (aClassifier->Perform(aNewPoint) == TopAbs_IN)
1060         insertVertex(aIsoV.Value(aU), aNewPoint.Coord(), theNewVertices);
1061     }
1062   }
1063 }
1064
1065 //=======================================================================
1066 //function : checkDeflectionAndInsert
1067 //purpose  :
1068 //=======================================================================
1069 Standard_Boolean BRepMesh_FastDiscretFace::checkDeflectionAndInsert(
1070   const gp_Pnt&              thePnt3d,
1071   const gp_XY&               theUV,
1072   const Standard_Boolean     isDeflectionCheckOnly,
1073   const Standard_Real        theTriangleDeflection,
1074   const Standard_Real        theFaceDeflection,
1075   const BRepMesh_CircleTool& theCircleTool,
1076   BRepMesh::ListOfVertex&    theVertices,
1077   Standard_Real&             theMaxTriangleDeflection)
1078 {
1079   if (theTriangleDeflection > theMaxTriangleDeflection)
1080     theMaxTriangleDeflection = theTriangleDeflection;
1081
1082   if (theTriangleDeflection < theFaceDeflection)
1083     return Standard_True;
1084
1085   if (myMinSize > Precision::Confusion())
1086   {
1087     // Iterator in the list of indexes of circles containing the node
1088     BRepMesh::ListOfInteger& aCirclesList =
1089       const_cast<BRepMesh_CircleTool&>(theCircleTool).Select(
1090       myAttribute->Scale(theUV, Standard_True));
1091
1092     Handle(NCollection_IncAllocator) aAllocator =
1093       new NCollection_IncAllocator(BRepMesh::MEMORY_BLOCK_SIZE_HUGE);
1094     BRepMesh::MapOfInteger aUsedNodes(10, aAllocator);
1095     BRepMesh::ListOfInteger::Iterator aCircleIt(aCirclesList);
1096     for (; aCircleIt.More(); aCircleIt.Next())
1097     {
1098       const BRepMesh_Triangle& aTriangle =
1099         myStructure->GetElement(aCircleIt.Value());
1100
1101       Standard_Integer aNodes[3];
1102       myStructure->ElementNodes(aTriangle, aNodes);
1103
1104       for (Standard_Integer i = 0; i < 3; ++i)
1105       {
1106         const Standard_Integer aNodeId = aNodes[i];
1107         if (aUsedNodes.Contains(aNodeId))
1108           continue;
1109
1110         aUsedNodes.Add(aNodeId);
1111         const BRepMesh_Vertex& aNode = myStructure->GetNode(aNodeId);
1112         const gp_Pnt& aPoint = myAttribute->GetPoint(aNode);
1113
1114         if (thePnt3d.SquareDistance(aPoint) < myMinSize * myMinSize)
1115           return Standard_True;
1116       }
1117     }
1118   }
1119
1120   if (isDeflectionCheckOnly)
1121     return Standard_False;
1122
1123   insertVertex(thePnt3d, theUV, theVertices);
1124   return Standard_True;
1125 }
1126
1127 //=======================================================================
1128 //function : control
1129 //purpose  :
1130 //=======================================================================
1131 Standard_Real BRepMesh_FastDiscretFace::control(
1132   BRepMesh::ListOfVertex&  theNewVertices,
1133   BRepMesh_Delaun&         theTrigu,
1134   const Standard_Boolean   theIsFirst)
1135 {
1136   Standard_Integer aTrianglesNb = myStructure->ElementsOfDomain().Extent();
1137   if (aTrianglesNb < 1)
1138     return -1.0;
1139
1140   //IMPORTANT: Constants used in calculations
1141   const Standard_Real MinimalArea2d     = 1.e-9;
1142   const Standard_Real MinimalSqLength3d = 1.e-12;
1143   const Standard_Real aSqDefFace = myAttribute->GetDefFace() * myAttribute->GetDefFace();
1144
1145   const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
1146
1147   Handle(Geom_Surface) aBSpline;
1148   GeomAbs_SurfaceType  aSurfType = gFace->GetType();
1149   if (aSurfType == GeomAbs_BezierSurface)
1150     aBSpline = gFace->Bezier();
1151   else if (aSurfType == GeomAbs_BSplineSurface)
1152     aBSpline = gFace->BSpline();
1153
1154   NCollection_DataMap<Standard_Integer, gp_Dir> aNorMap;
1155   BRepMesh::MapOfIntegerInteger                 aStatMap;
1156   NCollection_Map<BRepMesh_OrientedEdge>        aCouples(3 * aTrianglesNb);
1157   const BRepMesh_CircleTool& aCircles = theTrigu.Circles();
1158
1159   // Perform refinement passes
1160   // Define the number of iterations
1161   Standard_Integer       aIterationsNb = 11;
1162   const Standard_Integer aPassesNb = (theIsFirst ? 1 : aIterationsNb);
1163   // Initialize stop condition
1164   Standard_Real aMaxSqDef = -1.;
1165   Standard_Integer aPass = 1, aInsertedNb = 1;
1166   Standard_Boolean isAllDegenerated = Standard_False;
1167   for (; aPass <= aPassesNb && aInsertedNb && !isAllDegenerated; ++aPass)
1168   {
1169     theNewVertices.Clear();
1170
1171     // Reset stop condition
1172     aInsertedNb      = 0;
1173     aMaxSqDef        = -1.;
1174     isAllDegenerated = Standard_True;
1175
1176     aTrianglesNb = myStructure->ElementsOfDomain().Extent();
1177     if (aTrianglesNb < 1)
1178       break;
1179
1180     // Iterate on current triangles
1181     const BRepMesh::MapOfInteger& aTriangles = myStructure->ElementsOfDomain();
1182     BRepMesh::MapOfInteger::Iterator aTriangleIt(aTriangles);
1183     for (; aTriangleIt.More(); aTriangleIt.Next())
1184     {
1185       const Standard_Integer aTriangleId = aTriangleIt.Key();
1186       const BRepMesh_Triangle& aCurrentTriangle = myStructure->GetElement(aTriangleId);
1187
1188       if (aCurrentTriangle.Movability() == BRepMesh_Deleted)
1189         continue;
1190
1191       Standard_Integer v[3];
1192       myStructure->ElementNodes(aCurrentTriangle, v);
1193
1194       Standard_Integer e[3];
1195       Standard_Boolean o[3];
1196       aCurrentTriangle.Edges(e, o);
1197
1198       gp_XY xy[3];
1199       gp_XYZ p[3];
1200       Standard_Boolean m[3];
1201       for (Standard_Integer i = 0; i < 3; ++i)
1202       {
1203         m[i] = (myStructure->GetLink(e[i]).Movability() == BRepMesh_Frontier);
1204
1205         const BRepMesh_Vertex& aVertex = myStructure->GetNode(v[i]);
1206         xy[i] = myAttribute->Scale(aVertex.Coord(), Standard_False);
1207         p [i] = myAttribute->GetPoint(aVertex).Coord();
1208       }
1209
1210       gp_XYZ aLinkVec[3];
1211       Standard_Boolean isDegeneratedTri = Standard_False;
1212       for (Standard_Integer i = 0; i < 3 && !isDegeneratedTri; ++i)
1213       {
1214         aLinkVec[i] = p[(i + 1) % 3] - p[i];
1215         isDegeneratedTri = (aLinkVec[i].SquareModulus() < MinimalSqLength3d);
1216       }
1217
1218       if (isDegeneratedTri)
1219         continue;
1220
1221       isAllDegenerated = Standard_False;
1222
1223       // Check triangle area in 2d
1224       if (Abs((xy[1]-xy[0])^(xy[2]-xy[1])) < MinimalArea2d)
1225         continue;
1226
1227       // Check triangle normal
1228       gp_Pnt pDef;
1229       Standard_Real aSqDef = -1.;
1230       Standard_Boolean isSkipped = Standard_False;
1231       gp_XYZ normal(aLinkVec[0] ^ aLinkVec[1]);
1232       try
1233       {
1234         OCC_CATCH_SIGNALS
1235
1236         normal.Normalize();
1237
1238         // Check deflection on triangle
1239         gp_XY mi2d = (xy[0] + xy[1] + xy[2]) / 3.0;
1240         gFace->D0(mi2d.X(), mi2d.Y(), pDef);
1241         aSqDef = Abs(normal * (pDef.XYZ() - p[0]));
1242         aSqDef *= aSqDef;
1243
1244         isSkipped = !checkDeflectionAndInsert(pDef, mi2d, theIsFirst,
1245           aSqDef, aSqDefFace, aCircles, theNewVertices, aMaxSqDef);
1246
1247         if (isSkipped)
1248           break;
1249       }
1250       catch (Standard_Failure)
1251       {
1252         continue;
1253       }
1254
1255       for (Standard_Integer i = 0; i < 3 && !isSkipped; ++i)
1256       {
1257         if (m[i]) // is a boundary
1258           continue;
1259
1260         Standard_Integer j = (i + 1) % 3;
1261         // Check if this link was already processed
1262         Standard_Integer aFirstVertex, aLastVertex;
1263         if (v[i] < v[j])
1264         {
1265           aFirstVertex = v[i];
1266           aLastVertex = v[j];
1267         }
1268         else
1269         {
1270           aFirstVertex = v[j];
1271           aLastVertex = v[i];
1272         }
1273
1274         if (aCouples.Add(BRepMesh_OrientedEdge(aFirstVertex, aLastVertex)))
1275         {
1276           // Check deflection on edge 1
1277           gp_XY mi2d = (xy[i] + xy[j]) * 0.5;
1278           gFace->D0(mi2d.X(), mi2d.Y(), pDef);
1279           gp_Lin aLin(p[i], gp_Vec(p[i], p[j]));
1280           aSqDef = aLin.SquareDistance(pDef);
1281
1282           isSkipped = !checkDeflectionAndInsert(pDef, mi2d, theIsFirst,
1283             aSqDef, aSqDefFace, aCircles, theNewVertices, aMaxSqDef);
1284         }
1285       }
1286
1287       if (isSkipped)
1288         break;
1289
1290       //check normal on bsplines
1291       if (theIsFirst && !aBSpline.IsNull())
1292       {
1293         gp_Dir N[3] = { gp::DZ(), gp::DZ(), gp::DZ() };
1294         Standard_Integer aSt[3];
1295
1296         for (Standard_Integer i = 0; i < 3; ++i)
1297         {
1298           if (aNorMap.IsBound(v[i]))
1299           {
1300             aSt[i] = aStatMap.Find(v[i]);
1301             N[i] = aNorMap.Find(v[i]);
1302           }
1303           else
1304           {
1305             aSt[i] = GeomLib::NormEstim(aBSpline, gp_Pnt2d(xy[i]), Precision::Confusion(), N[i]);
1306             aStatMap.Bind(v[i], aSt[i]);
1307             aNorMap.Bind(v[i], N[i]);
1308           }
1309         }
1310
1311         Standard_Real aAngle[3];
1312         for (Standard_Integer i = 0; i < 3; ++i)
1313           aAngle[i] = N[(i + 1) % 3].Angle(N[i]);
1314
1315         if (aSt[0] < 1 && aSt[1] < 1 && aSt[2] < 1)
1316         {
1317           if (aAngle[0] > myAngle || aAngle[1] > myAngle || aAngle[2] > myAngle)
1318           {
1319             aMaxSqDef = -1.;
1320             break;
1321           }
1322         }
1323       }
1324     }
1325
1326     if (theIsFirst)
1327       continue;
1328
1329     if (addVerticesToMesh(theNewVertices, theTrigu))
1330       ++aInsertedNb;
1331   }
1332
1333   return (aMaxSqDef < 0) ? aMaxSqDef : Sqrt(aMaxSqDef);
1334 }
1335
1336 //=======================================================================
1337 //function : add
1338 //purpose  :
1339 //=======================================================================
1340 void BRepMesh_FastDiscretFace::add(const TopoDS_Vertex& theVertex)
1341 {
1342   if (theVertex.Orientation() != TopAbs_INTERNAL)
1343     return;
1344
1345   try
1346   {
1347     OCC_CATCH_SIGNALS
1348
1349     gp_Pnt2d aPnt2d = BRep_Tool::Parameters(theVertex, myAttribute->Face());
1350     // check UV values for internal vertices
1351     if (myAttribute->ChangeClassifier()->Perform(aPnt2d) != TopAbs_IN)
1352       return;
1353
1354     NCollection_Handle<FixedVExplorer> aFixedVExplorer = new FixedVExplorer(theVertex);
1355     Standard_Integer aIndex = myAttribute->GetVertexIndex(aFixedVExplorer);
1356     gp_XY anUV = BRepMesh_ShapeTool::FindUV(aIndex, aPnt2d,
1357       theVertex, BRep_Tool::Tolerance(theVertex), myAttribute);
1358
1359     Standard_Integer aTmpId1, aTmpId2;
1360     anUV = myAttribute->Scale(anUV, Standard_True);
1361     myAttribute->AddNode(aIndex, anUV, BRepMesh_Fixed, aTmpId1, aTmpId2);
1362   }
1363   catch (Standard_Failure)
1364   {
1365   }
1366 }
1367
1368 //=======================================================================
1369 //function : insertVertex
1370 //purpose  :
1371 //=======================================================================
1372 void BRepMesh_FastDiscretFace::insertVertex(
1373   const gp_Pnt&           thePnt3d,
1374   const gp_XY&            theUV,
1375   BRepMesh::ListOfVertex& theVertices)
1376 {
1377   Standard_Integer aNbLocat = myAttribute->LastPointId();
1378   myAttribute->ChangeSurfacePoints()->Bind(++aNbLocat, thePnt3d);
1379
1380   gp_XY aPnt2d  = myAttribute->Scale(theUV, Standard_True);
1381   BRepMesh_Vertex aVertex(aPnt2d, aNbLocat, BRepMesh_Free);
1382   theVertices.Append(aVertex);
1383 }
1384
1385 //=======================================================================
1386 //function : commitSurfaceTriangulation
1387 //purpose  :
1388 //=======================================================================
1389 void BRepMesh_FastDiscretFace::commitSurfaceTriangulation()
1390 {
1391   if (myAttribute.IsNull() || !myAttribute->IsValid())
1392     return;
1393
1394   const TopoDS_Face& aFace = myAttribute->Face();
1395   BRepMesh_ShapeTool::NullifyFace(aFace);
1396
1397   Handle(BRepMesh_DataStructureOfDelaun)& aStructure = myAttribute->ChangeStructure();
1398   const BRepMesh::MapOfInteger&           aTriangles = aStructure->ElementsOfDomain();
1399
1400   if (aTriangles.IsEmpty())
1401   {
1402     myAttribute->SetStatus(BRepMesh_Failure);
1403     return;
1404   }
1405
1406   BRepMesh::HIMapOfInteger& aVetrexEdgeMap = myAttribute->ChangeVertexEdgeMap();
1407
1408   // Store triangles
1409   Standard_Integer aVerticesNb  = aVetrexEdgeMap->Extent();
1410   Standard_Integer aTrianglesNb = aTriangles.Extent();
1411   Handle(Poly_Triangulation) aNewTriangulation =
1412     new Poly_Triangulation(aVerticesNb, aTrianglesNb, Standard_True);
1413
1414   Poly_Array1OfTriangle& aPolyTrianges = aNewTriangulation->ChangeTriangles();
1415
1416   Standard_Integer aTriangeId = 1;
1417   BRepMesh::MapOfInteger::Iterator aTriIt(aTriangles);
1418   for (; aTriIt.More(); aTriIt.Next())
1419   {
1420     const BRepMesh_Triangle& aCurElem = aStructure->GetElement(aTriIt.Key());
1421
1422     Standard_Integer aNode[3];
1423     aStructure->ElementNodes(aCurElem, aNode);
1424
1425     Standard_Integer aNodeId[3];
1426     for (Standard_Integer i = 0; i < 3; ++i)
1427       aNodeId[i] = aVetrexEdgeMap->FindIndex(aNode[i]);
1428
1429     aPolyTrianges(aTriangeId++).Set(aNodeId[0], aNodeId[1], aNodeId[2]);
1430   }
1431
1432   // Store mesh nodes
1433   TColgp_Array1OfPnt&   aNodes   = aNewTriangulation->ChangeNodes();
1434   TColgp_Array1OfPnt2d& aNodes2d = aNewTriangulation->ChangeUVNodes();
1435
1436   for (Standard_Integer i = 1; i <= aVerticesNb; ++i)
1437   {
1438     Standard_Integer       aVertexId = aVetrexEdgeMap->FindKey(i);
1439     const BRepMesh_Vertex& aVertex   = aStructure->GetNode(aVertexId);
1440     const gp_Pnt&          aPoint    = myAttribute->GetPoint(aVertex);
1441
1442     aNodes(i)   = aPoint;
1443     aNodes2d(i) = aVertex.Coord();
1444   }
1445
1446   aNewTriangulation->Deflection(myAttribute->GetDefFace());
1447   BRepMesh_ShapeTool::AddInFace(aFace, aNewTriangulation);
1448
1449   // Delete unused data
1450   myAttribute->ChangeStructure().Nullify();
1451   myAttribute->ChangeSurfacePoints().Nullify();
1452   myAttribute->ChangeSurfaceVertices().Nullify();
1453
1454   myAttribute->ChangeClassifier().Nullify();
1455   myAttribute->ChangeLocation2D().Nullify();
1456   myAttribute->ChangeVertexEdgeMap().Nullify();
1457 }