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