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