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