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