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