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