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