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