0026106: BRepMesh - revision of data model
[occt.git] / src / BRepMesh / BRepMesh_FastDiscret.cxx
1 // Created on: 1996-02-27
2 // Created by: Ekaterina SMIRNOVA
3 // Copyright (c) 1996-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #include <BRepMesh_FastDiscret.hxx>
18
19 #include <BRepMesh_WireChecker.hxx>
20 #include <BRepMesh_FastDiscretFace.hxx>
21 #include <BRepMesh_FaceAttribute.hxx>
22 #include <BRepMesh_DataStructureOfDelaun.hxx>
23 #include <BRepMesh_GeomTool.hxx>
24 #include <BRepMesh_PairOfPolygon.hxx>
25 #include <BRepMesh_Classifier.hxx>
26 #include <BRepMesh_EdgeParameterProvider.hxx>
27 #include <BRepMesh_IEdgeTool.hxx>
28 #include <BRepMesh_EdgeTessellator.hxx>
29 #include <BRepMesh_EdgeTessellationExtractor.hxx>
30
31 #include <BRepAdaptor_Curve.hxx>
32 #include <BRepAdaptor_Surface.hxx>
33 #include <BRepAdaptor_HSurface.hxx>
34
35 #include <Bnd_Box.hxx>
36 #include <BRepTools.hxx>
37 #include <BRepBndLib.hxx>
38 #include <BndLib_Add3dCurve.hxx>
39 #include <Poly_Triangulation.hxx>
40 #include <Poly_PolygonOnTriangulation.hxx>
41
42 #include <Precision.hxx>
43 #include <Geom2d_Curve.hxx>
44 #include <Geom2dAdaptor_HCurve.hxx>
45 #include <Geom_Surface.hxx>
46 #include <Geom_Plane.hxx>
47 #include <GeomAbs_SurfaceType.hxx>
48 #include <Extrema_LocateExtPC.hxx>
49
50 #include <TColStd_Array1OfInteger.hxx>
51 #include <TColStd_Array1OfCharacter.hxx>
52 #include <TColStd_HArray1OfReal.hxx>
53 #include <TColgp_Array1OfPnt2d.hxx>
54 #include <TColGeom2d_SequenceOfCurve.hxx>
55
56 #include <TopTools_SequenceOfShape.hxx>
57 #include <TopTools_ListIteratorOfListOfShape.hxx>
58
59 #include <TopAbs.hxx>
60 #include <TopoDS.hxx>
61 #include <TopoDS_Vertex.hxx>
62 #include <TopExp.hxx>
63 #include <TopExp_Explorer.hxx>
64
65 #include <OSD_Parallel.hxx>
66
67 #include <Standard_ErrorHandler.hxx>
68 #include <Standard_Failure.hxx>
69 #include <NCollection_IncAllocator.hxx>
70
71 #include <BRep_ListIteratorOfListOfPointRepresentation.hxx>
72 #include <BRep_PointRepresentation.hxx>
73
74 #include <vector>
75
76 IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_FastDiscret,Standard_Transient)
77
78 #define UVDEFLECTION 1.e-05
79
80 //=======================================================================
81 //function : BRepMesh_FastDiscret
82 //purpose  : 
83 //=======================================================================
84 BRepMesh_FastDiscret::BRepMesh_FastDiscret( const Bnd_Box&         theBox,
85                                             const BRepMesh_FastDiscret::Parameters& theParams)
86    :
87   myMapdefle(1000, new NCollection_IncAllocator()),
88   myBoundaryVertices(new BRepMesh::DMapOfVertexInteger),
89   myBoundaryPoints(new BRepMesh::DMapOfIntegerPnt),
90   myParameters(theParams),
91   myDtotale(0.)
92
93   if ( myParameters.Relative )
94     BRepMesh_ShapeTool::BoxMaxDimension(theBox, myDtotale);
95 }
96
97 //=======================================================================
98 //function : InitSharedFaces
99 //purpose  : 
100 //=======================================================================
101 void BRepMesh_FastDiscret::InitSharedFaces(const TopoDS_Shape& theShape)
102 {
103   TopExp::MapShapesAndAncestors(theShape, TopAbs_EDGE, TopAbs_FACE, mySharedFaces);
104 }
105
106 //=======================================================================
107 //function : Perform(shape)
108 //purpose  : 
109 //=======================================================================
110 void BRepMesh_FastDiscret::Perform(const TopoDS_Shape& theShape)
111 {
112   InitSharedFaces(theShape);
113
114   std::vector<TopoDS_Face> aFaces;
115   TopExp_Explorer anExplorer(theShape, TopAbs_FACE);
116   for (; anExplorer.More(); anExplorer.Next())
117   {
118     const TopoDS_Face& aFace = TopoDS::Face(anExplorer.Current());
119     Add(aFace);
120     aFaces.push_back(aFace);
121   }
122
123   OSD_Parallel::ForEach(aFaces.begin(), aFaces.end(), *this, !myParameters.InParallel, (Standard_Integer )aFaces.size());
124 }
125
126
127 //=======================================================================
128 //function : Process
129 //purpose  : 
130 //=======================================================================
131 void BRepMesh_FastDiscret::Process(const TopoDS_Face& theFace) const
132 {
133   Handle(BRepMesh_FaceAttribute) anAttribute;
134   if (GetFaceAttribute(theFace, anAttribute))
135   {
136     try
137     {
138       OCC_CATCH_SIGNALS
139
140       BRepMesh_FastDiscretFace aTool(myParameters.Angle, myParameters.MinSize, 
141         myParameters.InternalVerticesMode, myParameters.ControlSurfaceDeflection);
142       aTool.Perform(anAttribute);
143     }
144     catch (Standard_Failure)
145     {
146       anAttribute->SetStatus(BRepMesh_Failure);
147     }
148   }
149 }
150
151 //=======================================================================
152 //function : resetDataStructure
153 //purpose  : 
154 //=======================================================================
155 void BRepMesh_FastDiscret::resetDataStructure()
156 {
157   Handle(NCollection_IncAllocator) aAllocator;
158   if (myAttribute->ChangeStructure().IsNull())
159     aAllocator = new NCollection_IncAllocator(BRepMesh::MEMORY_BLOCK_SIZE_HUGE);
160   else
161     aAllocator = myAttribute->ChangeStructure()->Allocator();
162
163   myAttribute->Clear();
164   aAllocator->Reset(Standard_False);
165   Handle(BRepMesh_DataStructureOfDelaun) aStructure = 
166     new BRepMesh_DataStructureOfDelaun(aAllocator);
167
168   const Standard_Real aTolU = myAttribute->ToleranceU();
169   const Standard_Real aTolV = myAttribute->ToleranceV();
170   const Standard_Real uCellSize = 14.0 * aTolU;
171   const Standard_Real vCellSize = 14.0 * aTolV;
172
173   aStructure->Data()->SetCellSize ( uCellSize, vCellSize);
174   aStructure->Data()->SetTolerance( aTolU    , aTolV    );
175
176   myAttribute->ChangeStructure() = aStructure;
177 }
178
179 //=======================================================================
180 //function : Add(face)
181 //purpose  : 
182 //=======================================================================
183 Standard_Integer BRepMesh_FastDiscret::Add(const TopoDS_Face& theFace)
184 {
185   myAttribute.Nullify();
186   GetFaceAttribute(theFace, myAttribute, Standard_True);
187
188   try
189   {
190     OCC_CATCH_SIGNALS
191
192     // Initialize face attributes
193     if (!myAttribute->IsInitialized ())
194       myAttribute->SetFace (theFace, myParameters.AdaptiveMin);
195     
196     BRepMesh::HIMapOfInteger&            aVertexEdgeMap = myAttribute->ChangeVertexEdgeMap();
197     BRepMesh::HDMapOfShapePairOfPolygon& aInternalEdges = myAttribute->ChangeInternalEdges();
198
199     resetDataStructure();
200
201     Standard_Real defedge = myParameters.Deflection;
202     Standard_Integer nbEdge = 0;
203     Standard_Real savangle = myParameters.Angle;
204     Standard_Real cdef;
205     Standard_Real maxdef = 2.* BRepMesh_ShapeTool::MaxFaceTolerance(theFace);
206
207     Standard_Real defface = 0.;
208     if (!myParameters.Relative)
209     {
210       defedge = Max(UVDEFLECTION, defedge);
211       defface = Max(myParameters.Deflection, maxdef);
212     }
213
214     const TopoDS_Face&                  aFace = myAttribute->Face();
215     for (TopoDS_Iterator aWireIt(aFace); aWireIt.More(); aWireIt.Next())
216     {
217       for (TopoDS_Iterator aEdgeIt(aWireIt.Value()); aEdgeIt.More(); aEdgeIt.Next(), ++nbEdge)
218       {
219         const TopoDS_Edge& aEdge = TopoDS::Edge(aEdgeIt.Value());
220         if (aEdge.IsNull())
221           continue;
222         if (myParameters.Relative)
223         {
224           if (!myMapdefle.IsBound(aEdge))
225           {
226             if (myEdges.IsBound(aEdge))
227             {
228               const BRepMesh_PairOfPolygon& aPair = myEdges.Find(aEdge);
229               const Handle(Poly_PolygonOnTriangulation)& aPolygon = aPair.First();
230               defedge = aPolygon->Deflection();
231             }
232             else
233             {
234               defedge = BRepMesh_ShapeTool::RelativeEdgeDeflection(
235                 aEdge, myParameters.Deflection, myDtotale, cdef);
236
237               myParameters.Angle = savangle * cdef;
238             }
239           }
240           else
241           {
242             defedge = myMapdefle(aEdge);
243           }
244
245           defface += defedge;
246           defface = Max(maxdef, defface);
247
248           if (!myMapdefle.IsBound(aEdge))
249           {
250             defedge = Max(UVDEFLECTION, defedge);
251             myMapdefle.Bind(aEdge, defedge);
252           }
253         }
254         else
255         {
256           if (!myMapdefle.IsBound(aEdge))
257           {
258             myMapdefle.Bind(aEdge, defedge);
259           }
260         }
261
262         Standard_Real aFirstParam, aLastParam;
263         Handle(Geom2d_Curve) aCurve2d = 
264           BRep_Tool::CurveOnSurface(aEdge, aFace, aFirstParam, aLastParam);
265
266         if (aCurve2d.IsNull())
267           continue;
268         Handle(Geom2dAdaptor_HCurve) aPCurve =
269           new Geom2dAdaptor_HCurve(aCurve2d, aFirstParam, aLastParam);
270
271         add(aEdge, aPCurve, defedge);
272         myParameters.Angle = savangle;
273       }
274     }
275
276     if ( nbEdge == 0 || aVertexEdgeMap->Extent() < 3 )
277     {
278       myAttribute->ChangeStructure().Nullify();
279       myAttribute->SetStatus(BRepMesh_Failure);
280       return myAttribute->GetStatus();
281     }
282
283     if ( myParameters.Relative )
284     {
285       defface = defface / nbEdge;
286     }
287     else
288     {
289       defface = myParameters.Deflection;
290     }
291
292     defface = Max(maxdef, defface);
293
294     TopLoc_Location aLoc;
295     Handle(Poly_Triangulation) aTriangulation = BRep_Tool::Triangulation(aFace, aLoc);
296     const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
297
298     if ( aTriangulation.IsNull() )
299     {
300       Standard_Real xCur, yCur;
301       Standard_Real maxX, minX, maxY, minY;
302
303       minX = minY = 1.e100;
304       maxX = maxY =-1.e100;
305
306       Standard_Integer ipn = 0;
307       Standard_Integer i1 = 1;
308       for ( i1 = 1; i1 <= aVertexEdgeMap->Extent(); ++i1 )
309       {
310         const BRepMesh_Vertex& aVertex = 
311           myAttribute->ChangeStructure()->GetNode(aVertexEdgeMap->FindKey(i1));
312
313         ++ipn;
314
315         xCur = aVertex.Coord().X();
316         yCur = aVertex.Coord().Y();
317
318         minX = Min(xCur, minX);
319         maxX = Max(xCur, maxX);
320         minY = Min(yCur, minY);
321         maxY = Max(yCur, maxY);
322       }
323
324       Standard_Real myumin = minX;
325       Standard_Real myumax = maxX;
326       Standard_Real myvmin = minY;
327       Standard_Real myvmax = maxY;
328
329       const Standard_Real umin = gFace->FirstUParameter();
330       const Standard_Real umax = gFace->LastUParameter();
331       const Standard_Real vmin = gFace->FirstVParameter();
332       const Standard_Real vmax = gFace->LastVParameter();
333
334       if (myumin < umin || myumax > umax)
335       {
336         if (gFace->IsUPeriodic())
337         {
338           if ((myumax - myumin) > (umax - umin))
339             myumax = myumin + (umax - umin);
340         }
341         else
342         {
343           if (umin > myumin)
344             myumin = umin;
345
346           if (umax < myumax)
347             myumax = umax;
348         }
349       }
350
351       if (myvmin < vmin || myvmax > vmax)
352       {
353         if (gFace->IsVPeriodic())
354         {
355           if ((myvmax - myvmin) > (vmax - vmin))
356             myvmax = myvmin + (vmax - vmin);
357         }
358         else
359         {
360           if ( vmin > myvmin )
361             myvmin = vmin;
362
363           if (vmax < myvmax)
364             myvmax = vmax;
365         }
366       }
367
368       GeomAbs_SurfaceType aSurfType = gFace->GetType();
369       // Fast verification of the validity of calculated limits.
370       // If wrong, sure a problem of pcurve.
371       if (aSurfType == GeomAbs_BezierSurface &&
372          (myumin < -0.5 || myumax > 1.5 || myvmin < -0.5 || myvmax > 1.5) )
373       {
374         myAttribute->ChangeStructure().Nullify();
375         myAttribute->SetStatus(BRepMesh_Failure);
376         return myAttribute->GetStatus();
377       }
378
379       //define parameters for correct parametrics
380       Standard_Real deltaX = 1.0;
381       Standard_Real deltaY = 1.0;
382
383       {
384         Standard_Real aTolU, aTolV;
385         myAttribute->ChangeStructure()->Data()->GetTolerance(aTolU, aTolV);
386         const Standard_Real aTol = Sqrt(aTolU * aTolU + aTolV * aTolV);
387
388         BRepMesh::HClassifier& aClassifier = myAttribute->ChangeClassifier();
389         BRepMesh_WireChecker aDFaceChecker(aFace, aTol, aInternalEdges, 
390           aVertexEdgeMap, myAttribute->ChangeStructure(),
391           myumin, myumax, myvmin, myvmax, myParameters.InParallel );
392
393         aDFaceChecker.ReCompute(aClassifier);
394         BRepMesh_Status aCheckStatus = aDFaceChecker.Status();
395
396         if (aCheckStatus == BRepMesh_SelfIntersectingWire)
397         {
398           Standard_Integer nbmaill = 0;
399           Standard_Real eps = Precision::Confusion();
400           while (nbmaill < 5 && aCheckStatus != BRepMesh_ReMesh)
401           {
402             ++nbmaill;
403
404             resetDataStructure();
405
406             for (TopoDS_Iterator aWireIt(aFace); aWireIt.More(); aWireIt.Next())
407             {
408               for (TopoDS_Iterator aEdgeIt(aWireIt.Value()); aEdgeIt.More(); aEdgeIt.Next(), ++nbEdge)
409               {
410                 const TopoDS_Edge& anEdge = TopoDS::Edge(aEdgeIt.Value());
411                 if (anEdge.IsNull())
412                   continue;
413                 if (myEdges.IsBound(anEdge))
414                   myEdges.UnBind(anEdge);
415
416                 defedge = Max(myMapdefle(anEdge) / 3.0, eps);
417                 myMapdefle.Bind(anEdge, defedge);
418
419                 Standard_Real aFirstParam, aLastParam;
420                 Handle(Geom2d_Curve) aCurve2d =
421                   BRep_Tool::CurveOnSurface(anEdge, aFace, aFirstParam, aLastParam);
422                 if (aCurve2d.IsNull())
423                   continue;
424
425                 Handle(Geom2dAdaptor_HCurve) aPCurve =
426                   new Geom2dAdaptor_HCurve(aCurve2d, aFirstParam, aLastParam);
427                 add(anEdge, aPCurve, defedge);
428               }
429             }
430
431             aDFaceChecker.ReCompute(aClassifier);
432             if (aDFaceChecker.Status() == BRepMesh_NoError)
433               aCheckStatus = BRepMesh_ReMesh;
434           }
435         }
436
437         myAttribute->SetStatus(aCheckStatus);
438         if (!myAttribute->IsValid())
439         {
440           myAttribute->ChangeStructure().Nullify();
441           return myAttribute->GetStatus();
442         }
443       }
444
445       // try to find the real length:
446       // akm (bug OCC16) : We must calculate these measures in non-singular
447       //     parts of face. Let's try to compute average value of three
448       //     (umin, (umin+umax)/2, umax), and respectively for v.
449       //                 vvvvv
450       Standard_Real longu = 0.0, longv = 0.0; //, last , first;
451       gp_Pnt P11, P12, P21, P22, P31, P32;
452
453       Standard_Real du = 0.05 * ( myumax - myumin );
454       Standard_Real dv = 0.05 * ( myvmax - myvmin );
455       Standard_Real dfuave = 0.5 * ( myumin + myumax );
456       Standard_Real dfvave = 0.5 * ( myvmin + myvmax );
457       Standard_Real dfucur, dfvcur;
458
459       // U loop
460       gFace->D0(myumin, myvmin, P11);
461       gFace->D0(myumin, dfvave, P21);
462       gFace->D0(myumin, myvmax, P31);
463       for (i1=1, dfucur=myumin+du; i1 <= 20; i1++, dfucur+=du)
464       {
465         gFace->D0(dfucur, myvmin, P12);
466         gFace->D0(dfucur, dfvave, P22);
467         gFace->D0(dfucur, myvmax, P32);
468         longu += ( P11.Distance(P12) + P21.Distance(P22) + P31.Distance(P32) );
469         P11 = P12;
470         P21 = P22;
471         P31 = P32;
472       }
473
474       // V loop
475       gFace->D0(myumin, myvmin, P11);
476       gFace->D0(dfuave, myvmin, P21);
477       gFace->D0(myumax, myvmin, P31);
478       for (i1=1, dfvcur=myvmin+dv; i1 <= 20; i1++, dfvcur+=dv)
479       {
480         gFace->D0(myumin, dfvcur, P12);
481         gFace->D0(dfuave, dfvcur, P22);
482         gFace->D0(myumax, dfvcur, P32);
483         longv += ( P11.Distance(P12) + P21.Distance(P22) + P31.Distance(P32) );
484         P11 = P12;
485         P21 = P22;
486         P31 = P32;
487       }
488
489       longu /= 3.;
490       longv /= 3.;
491       // akm (bug OCC16) ^^^^^
492
493       if (longu <= 1.e-16 || longv <= 1.e-16)
494       {
495         //yes, it is seen!!
496         myAttribute->ChangeStructure().Nullify();
497         myAttribute->SetStatus(BRepMesh_Failure);
498         return myAttribute->GetStatus();
499       }
500
501
502       if (aSurfType == GeomAbs_Torus)
503       {
504         gp_Torus Tor = gFace->Torus();
505         Standard_Real r = Tor.MinorRadius(), R = Tor.MajorRadius();
506         Standard_Real Du, Dv;//, pasu, pasv;
507
508         Dv = Max(1.0e0 - (defface/r),0.0e0) ;
509         Standard_Real oldDv = 2.0 * ACos (Dv);
510         oldDv = Min(oldDv, myParameters.Angle);
511         Dv  =  0.9*oldDv; //TWOTHIRD * oldDv;
512         Dv = oldDv;
513
514         Standard_Integer nbV = Max((Standard_Integer)((myvmax-myvmin)/Dv), 2);
515         Dv = (myvmax-myvmin)/(nbV+1);
516
517         Standard_Real ru = R + r;
518         if ( ru > 1.e-16 )
519         {
520           Du = Max(1.0e0 - (defface/ru),0.0e0);
521           Du  = (2.0 * ACos (Du));
522           Du = Min(Du, myParameters.Angle);
523           Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
524
525           if (aa < gp::Resolution())
526           {
527             myAttribute->ChangeStructure().Nullify();
528             return myAttribute->GetStatus();
529           }
530
531           Du = Du * Min(oldDv, Du) / aa;
532         }
533         else
534         {
535           Du = Dv;
536         }
537
538         Standard_Integer nbU = Max((Standard_Integer)((myumax-myumin)/Du), 2);
539         nbU = Max(nbU, (Standard_Integer)(nbV*(myumax-myumin)*R/((myvmax-myvmin)*r)/5.));
540       
541         Du = (myumax-myumin)/(nbU+1);
542         //-- DeltaX and DeltaY are chosen so that to avoid "jumping" 
543         //-- of points on the grid
544         deltaX = Du;
545         deltaY = Dv;
546       }
547       else if (aSurfType == GeomAbs_Cylinder)
548       {
549         gp_Cylinder Cyl = gFace->Cylinder();
550         Standard_Real R = Cyl.Radius();
551
552         // Calculate parameters for iteration in U direction
553         Standard_Real Du = 1.0 - (defface/R);
554         if (Du < 0.0)
555           Du = 0.0;
556
557         Du = 2.0 * ACos (Du);
558         if (Du > myParameters.Angle)
559           Du = myParameters.Angle;
560
561         deltaX = Du / longv;
562         deltaY = 1.;
563       }
564       else
565       {
566         deltaX = (myumax-myumin)/longu;
567         deltaY = (myvmax-myvmin)/longv;
568       }
569
570       // Restore face attribute
571       myAttribute->SetDefFace(defface);
572       myAttribute->SetUMax(myumax);
573       myAttribute->SetVMax(myvmax);
574       myAttribute->SetUMin(myumin);
575       myAttribute->SetVMin(myvmin);
576       myAttribute->SetDeltaX(deltaX);
577       myAttribute->SetDeltaY(deltaY);
578     }
579
580     myAttribute->ChangeMeshNodes() = 
581       myAttribute->ChangeStructure()->Data()->Vertices();
582   }
583   catch(Standard_Failure)
584   {
585     myAttribute->SetStatus(BRepMesh_Failure);
586   }
587
588   myAttribute->ChangeStructure().Nullify();
589   return myAttribute->GetStatus();
590 }
591
592 //=======================================================================
593 //function : getEdgeAttributes
594 //purpose  :
595 //=======================================================================
596 Standard_Boolean BRepMesh_FastDiscret::getEdgeAttributes(
597   const TopoDS_Edge&                      theEdge,
598   const Handle(Geom2dAdaptor_HCurve)&     thePCurve,
599   const Standard_Real                     theDefEdge,
600   BRepMesh_FastDiscret::EdgeAttributes&   theAttributes) const
601 {
602   EdgeAttributes& aEAttr = theAttributes;
603
604   // Get vertices
605   TopExp::Vertices(theEdge, aEAttr.FirstVertex, aEAttr.LastVertex);
606   if (aEAttr.FirstVertex.IsNull() || aEAttr.LastVertex.IsNull())
607     return Standard_False;
608
609   // Get range on 2d curve
610   const TopoDS_Face& aFace = myAttribute->Face();
611   BRep_Tool::Range(theEdge, aFace, aEAttr.FirstParam, aEAttr.LastParam);
612
613   // Get end points on 2d curve
614   BRep_Tool::UVPoints(theEdge, aFace, aEAttr.FirstUV, aEAttr.LastUV);
615
616   aEAttr.IsSameUV =
617     aEAttr.FirstUV.IsEqual(aEAttr.LastUV, Precision::PConfusion());
618   if (aEAttr.IsSameUV)
619   {
620     // 1. is it really sameUV without being degenerated
621     gp_Pnt2d uvF, uvL;
622     thePCurve->D0(thePCurve->FirstParameter(), uvF);
623     thePCurve->D0(thePCurve->LastParameter(),  uvL);
624
625     if (!aEAttr.FirstUV.IsEqual(uvF, Precision::PConfusion()))
626       aEAttr.FirstUV = uvF;
627
628     if (!aEAttr.LastUV.IsEqual(uvL, Precision::PConfusion()))
629       aEAttr.LastUV = uvL;
630
631     aEAttr.IsSameUV =
632       aEAttr.FirstUV.IsEqual(aEAttr.LastUV, Precision::PConfusion());
633   }
634
635   //Control tolerance of vertices
636   const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
637   gp_Pnt pFirst = gFace->Value(aEAttr.FirstUV.X(), aEAttr.FirstUV.Y());
638   gp_Pnt pLast  = gFace->Value(aEAttr.LastUV.X(),  aEAttr.LastUV.Y());
639
640   Standard_Real aSqDist = pFirst.SquareDistance(BRep_Tool::Pnt(aEAttr.FirstVertex));
641   aSqDist = Max(aSqDist, pLast.SquareDistance(BRep_Tool::Pnt(aEAttr.LastVertex)));
642
643   aEAttr.Deflection = Max(theDefEdge, BRep_Tool::Tolerance(theEdge));
644   aEAttr.Deflection = Max(aEAttr.Deflection, sqrt(aSqDist));
645
646   return Standard_True;
647 }
648
649 //=======================================================================
650 //function : registerEdgeVertices
651 //purpose  : 
652 //=======================================================================
653 void BRepMesh_FastDiscret::registerEdgeVertices(
654   BRepMesh_FastDiscret::EdgeAttributes& theAttributes,
655   Standard_Integer&                     ipf,
656   Standard_Integer&                     ivf,
657   Standard_Integer&                     isvf,
658   Standard_Integer&                     ipl,
659   Standard_Integer&                     ivl,
660   Standard_Integer&                     isvl)
661 {
662   EdgeAttributes& aEAttr = theAttributes;
663   if (aEAttr.FirstVExtractor.IsNull())
664   {
665     // Use edge geometry to produce tesselation.
666     aEAttr.FirstVExtractor = 
667       new TopoDSVExplorer(aEAttr.FirstVertex, aEAttr.IsSameUV, aEAttr.LastVertex);
668   }
669
670   if (aEAttr.LastVExtractor.IsNull())
671   {
672     // Use edge geometry to produce tesselation.
673     aEAttr.LastVExtractor = 
674       new TopoDSVExplorer(aEAttr.LastVertex, aEAttr.IsSameUV, aEAttr.FirstVertex);
675   }
676
677   // Process first vertex
678   ipf = myAttribute->GetVertexIndex(aEAttr.FirstVExtractor, Standard_True);
679   Standard_Real aMinDist = 2. * BRep_Tool::Tolerance(aEAttr.FirstVertex);
680   gp_XY aTmpUV1 = BRepMesh_ShapeTool::FindUV(ipf, aEAttr.FirstUV, aMinDist, myAttribute);
681
682   myAttribute->AddNode(ipf, aTmpUV1, BRepMesh_Frontier, ivf, isvf);
683
684   // Process last vertex
685   ipl = aEAttr.LastVertex.IsSame(aEAttr.FirstVertex) ? ipf :
686     myAttribute->GetVertexIndex(aEAttr.LastVExtractor, Standard_True);
687   aMinDist = 2. * BRep_Tool::Tolerance(aEAttr.LastVertex);
688   gp_XY aTmpUV2 = BRepMesh_ShapeTool::FindUV(ipl, aEAttr.LastUV, aMinDist, myAttribute);
689
690   myAttribute->AddNode(ipl, aTmpUV2, BRepMesh_Frontier, ivl, isvl);
691
692   // Update edge deflection
693   const Handle(BRepAdaptor_HSurface)& gFace = myAttribute->Surface();
694   gp_Pnt aPntE1 = gFace->Value(aEAttr.FirstUV.X(), aEAttr.FirstUV.Y());
695   gp_Pnt aPntFound1 = gFace->Value(aTmpUV1.X(), aTmpUV1.Y());
696   Standard_Real aSqDist = aPntE1.SquareDistance(aPntFound1);
697   gp_Pnt aPntE2 = gFace->Value(aEAttr.LastUV.X(), aEAttr.LastUV.Y());
698   gp_Pnt aPntFound2 = gFace->Value(aTmpUV2.X(), aTmpUV2.Y());
699   aSqDist = Max(aSqDist, aPntE2.SquareDistance(aPntFound2));
700   aEAttr.Deflection = Max(aEAttr.Deflection, sqrt(aSqDist));
701 }
702
703 //=======================================================================
704 //function : add
705 //purpose  : 
706 //=======================================================================
707 void BRepMesh_FastDiscret::add(
708   const TopoDS_Edge&                      theEdge,
709   const Handle(Geom2dAdaptor_HCurve)&     thePCurve,
710   const Standard_Real                     theDefEdge)
711 {
712   const TopAbs_Orientation orEdge = theEdge.Orientation();
713   if (orEdge == TopAbs_EXTERNAL)
714     return;
715
716   EdgeAttributes aEAttr;
717   if (!getEdgeAttributes(theEdge, thePCurve, theDefEdge, aEAttr))
718     return;
719
720   if (!myEdges.IsBound(theEdge))
721   {
722     update(theEdge, thePCurve, theDefEdge, aEAttr);
723     return;
724   }
725
726   Standard_Integer ipf, ivf, isvf, ipl, ivl, isvl;
727   registerEdgeVertices(aEAttr, ipf, ivf, isvf, ipl, ivl, isvl);
728
729   // If this Edge has been already checked and it is not degenerated, 
730   // the points of the polygon calculated at the first check are retrieved :
731
732   // retrieve the polygone:
733   const BRepMesh_PairOfPolygon&              aPair    = myEdges.Find(theEdge);
734   const Handle(Poly_PolygonOnTriangulation)& aPolygon = aPair.First();
735   const TColStd_Array1OfInteger&             aNodes   = aPolygon->Nodes();
736   Handle(TColStd_HArray1OfReal)              aParams  = aPolygon->Parameters();
737
738   // creation anew:
739   const Standard_Integer  aNodesNb = aNodes.Length();
740   TColStd_Array1OfInteger aNewNodes (1, aNodesNb);
741   TColStd_Array1OfReal    aNewParams(1, aNodesNb);
742
743   aNewNodes (1) = isvf;
744   aNewParams(1) = aEAttr.FirstParam;
745
746   aNewNodes (aNodesNb) = isvl;
747   aNewParams(aNodesNb) = aEAttr.LastParam;
748
749   const TopoDS_Face& aFace = myAttribute->Face();
750   if (!BRepMesh_ShapeTool::IsDegenerated(theEdge, aFace))
751   {
752     BRepMesh_EdgeParameterProvider aProvider(theEdge, aFace, aParams);
753     for (Standard_Integer i = 2; i < aNodesNb; ++i)
754     {
755       const Standard_Integer aPointId = aNodes(i);
756       const gp_Pnt& aPnt = myBoundaryPoints->Find(aPointId);
757
758       const Standard_Real aParam = aProvider.Parameter(i, aPnt);
759       gp_Pnt2d aUV = thePCurve->Value(aParam);
760
761       Standard_Integer iv2, isv;
762       myAttribute->AddNode(aPointId, aUV.Coord(), BRepMesh_OnCurve, iv2, isv);
763
764       aNewNodes (i) = isv;
765       aNewParams(i) = aParam;
766     }
767   }
768
769   Handle(Poly_PolygonOnTriangulation) P1 = 
770     new Poly_PolygonOnTriangulation(aNewNodes, aNewParams);
771
772   storePolygon(theEdge, P1, aEAttr.Deflection);
773 }
774
775 //=======================================================================
776 //function : update(edge)
777 //purpose  :
778 //=======================================================================
779 void BRepMesh_FastDiscret::update(
780   const TopoDS_Edge&                                theEdge,
781   const Handle(Geom2dAdaptor_HCurve)&               theC2d,
782   const Standard_Real                               theDefEdge,
783   BRepMesh_FastDiscret::EdgeAttributes&             theAttributes)
784 {
785   EdgeAttributes& aEAttr = theAttributes;
786
787   const TopoDS_Face& aFace = myAttribute->Face();
788   Handle(BRepMesh_IEdgeTool) aEdgeTool;
789   // Try to find existing tessellation.
790   for (Standard_Integer i = 1; aEdgeTool.IsNull(); ++i)
791   {
792     TopLoc_Location aLoc;
793     Handle(Poly_Triangulation) aTriangulation;
794     Handle(Poly_PolygonOnTriangulation) aPolygon;
795     BRep_Tool::PolygonOnTriangulation(theEdge, aPolygon, aTriangulation, aLoc, i);
796
797     if (aPolygon.IsNull())
798       break;
799
800     if (aTriangulation.IsNull() || !aPolygon->HasParameters())
801       continue;
802
803     if (aPolygon->Deflection() > 1.1 * theDefEdge)
804       continue;
805
806     const TColgp_Array1OfPnt&      aNodes   = aTriangulation->Nodes();
807     const TColStd_Array1OfInteger& aIndices = aPolygon->Nodes();
808     Handle(TColStd_HArray1OfReal)  aParams  = aPolygon->Parameters();
809
810     aEAttr.FirstVExtractor = new PolyVExplorer(aEAttr.FirstVertex, 
811       aEAttr.IsSameUV, aEAttr.LastVertex, aIndices(1), aNodes, aLoc);
812
813     aEAttr.LastVExtractor = new PolyVExplorer(aEAttr.LastVertex, 
814       aEAttr.IsSameUV, aEAttr.FirstVertex, aIndices(aIndices.Length()), aNodes, aLoc);
815
816     aEdgeTool = new BRepMesh_EdgeTessellationExtractor(theEdge, theC2d, 
817       aFace, aTriangulation, aPolygon, aLoc);
818   }
819
820   if (aEdgeTool.IsNull())
821   {
822     aEdgeTool = new BRepMesh_EdgeTessellator(theEdge, myAttribute, 
823       mySharedFaces, theDefEdge, myParameters.Angle, myParameters.MinSize);
824   }
825
826   Standard_Integer ipf, ivf, isvf, ipl, ivl, isvl;
827   registerEdgeVertices(aEAttr, ipf, ivf, isvf, ipl, ivl, isvl);
828
829   Handle(Poly_PolygonOnTriangulation) P1, P2;
830   if (BRepMesh_ShapeTool::IsDegenerated(theEdge, aFace))
831   {
832     // two nodes
833     Standard_Integer aNewNodesArr[] = {isvf, isvl};
834     Standard_Integer aNewNodInStructArr[] = {ipf, ipl};
835     Standard_Real aNewParamsArr[] = {aEAttr.FirstParam, aEAttr.LastParam};
836     TColStd_Array1OfInteger aNewNodes      (aNewNodesArr[0], 1, 2);
837     TColStd_Array1OfInteger aNewNodInStruct(aNewNodInStructArr[0], 1, 2);
838     TColStd_Array1OfReal    aNewParams     (aNewParamsArr[0], 1, 2);
839
840     P1 = new Poly_PolygonOnTriangulation(aNewNodes,       aNewParams);
841     P2 = new Poly_PolygonOnTriangulation(aNewNodInStruct, aNewParams);
842   }
843   else
844   {
845     const Standard_Integer  aNodesNb = aEdgeTool->NbPoints();
846     // Allocate the memory for arrays aNewNodesVec, aNewNodesInStructVec, aNewParamsVec
847     // only once using the buffer aBuf.
848     TColStd_Array1OfCharacter aBuf(1, aNodesNb * (2*sizeof(Standard_Integer) + sizeof(Standard_Real)));
849     TColStd_Array1OfInteger aNewNodesVec(*reinterpret_cast<const Standard_Integer*>
850       (&aBuf(1)), 1, aNodesNb);
851     TColStd_Array1OfInteger aNewNodesInStructVec(*reinterpret_cast<const Standard_Integer*>
852       (&aBuf(1 + aNodesNb*sizeof(Standard_Integer))), 1, aNodesNb);
853     TColStd_Array1OfReal    aNewParamsVec(*reinterpret_cast<const Standard_Real*>
854       (&aBuf(1 + aNodesNb*2*sizeof(Standard_Integer))), 1, aNodesNb);
855
856     Standard_Integer aNodesCount = 1;
857     aNewNodesInStructVec(aNodesCount) = ipf;
858     aNewNodesVec        (aNodesCount) = isvf;
859     aNewParamsVec       (aNodesCount) = aEAttr.FirstParam;
860
861     ++aNodesCount;
862     Standard_Integer aPrevNodeId  = ivf;
863     Standard_Integer aLastPointId = myAttribute->LastPointId();
864     for (Standard_Integer i = 2; i < aNodesNb; ++i)
865     {
866       gp_Pnt        aPnt;
867       gp_Pnt2d      aUV;
868       Standard_Real aParam;
869       if (!aEdgeTool->Value(i, aParam, aPnt, aUV))
870         continue;
871
872       // Imitate index of 3d point in order to not to add points to map without necessity.
873       Standard_Integer iv2, isv;
874       myAttribute->AddNode(aLastPointId + 1, aUV.Coord(), BRepMesh_Frontier, iv2, isv);
875       if (aPrevNodeId == iv2)
876         continue;
877
878       // Ok, now we can add point to the map.
879       myBoundaryPoints->Bind (++aLastPointId, aPnt);
880
881       aNewNodesInStructVec(aNodesCount) = aLastPointId;
882       aNewNodesVec        (aNodesCount) = isv;
883       aNewParamsVec       (aNodesCount) = aParam;
884
885       ++aNodesCount;
886       aPrevNodeId = iv2;
887     }
888
889     aNewNodesInStructVec(aNodesCount) = ipl;
890     aNewNodesVec        (aNodesCount) = isvl;
891     aNewParamsVec       (aNodesCount) = aEAttr.LastParam;
892
893     TColStd_Array1OfInteger aNewNodes      (aNewNodesVec.First (),        1, aNodesCount);
894     TColStd_Array1OfInteger aNewNodInStruct(aNewNodesInStructVec.First(), 1, aNodesCount);
895     TColStd_Array1OfReal    aNewParams     (aNewParamsVec.First(),        1, aNodesCount);
896
897     P1 = new Poly_PolygonOnTriangulation(aNewNodes,       aNewParams);
898     P2 = new Poly_PolygonOnTriangulation(aNewNodInStruct, aNewParams);
899   }
900
901   storePolygon(theEdge, P1, aEAttr.Deflection);
902   storePolygonSharedData(theEdge, P2, aEAttr.Deflection);
903 }
904
905 //=======================================================================
906 //function : storeInternalPolygon
907 //purpose  : 
908 //=======================================================================
909 void BRepMesh_FastDiscret::storePolygon(
910   const TopoDS_Edge&                         theEdge,
911   Handle(Poly_PolygonOnTriangulation)&       thePolygon,
912   const Standard_Real                        theDeflection)
913 {
914   thePolygon->Deflection(theDeflection);
915   BRepMesh::HDMapOfShapePairOfPolygon& aInternalEdges = myAttribute->ChangeInternalEdges();
916   if (aInternalEdges->IsBound(theEdge))
917   {
918     BRepMesh_PairOfPolygon& aPair = aInternalEdges->ChangeFind(theEdge);
919     if (theEdge.Orientation() == TopAbs_REVERSED)
920       aPair.Append(thePolygon);
921     else
922       aPair.Prepend(thePolygon);
923
924     return;
925   }
926
927   BRepMesh_PairOfPolygon aPair;
928   aPair.Append(thePolygon);
929   aInternalEdges->Bind(theEdge, aPair);
930 }
931
932 //=======================================================================
933 //function : storePolygonSharedData
934 //purpose  : 
935 //=======================================================================
936 void BRepMesh_FastDiscret::storePolygonSharedData(
937   const TopoDS_Edge&                         theEdge,
938   Handle(Poly_PolygonOnTriangulation)&       thePolygon,
939   const Standard_Real                        theDeflection)
940 {
941   thePolygon->Deflection(theDeflection);
942   BRepMesh_PairOfPolygon aPair;
943   aPair.Append(thePolygon);
944   myEdges.Bind(theEdge, aPair);
945 }
946
947 //=======================================================================
948 //function : GetFaceAttribute
949 //purpose  : 
950 //=======================================================================
951 Standard_Boolean BRepMesh_FastDiscret::GetFaceAttribute(
952   const TopoDS_Face&              theFace,
953   Handle(BRepMesh_FaceAttribute)& theAttribute,
954   const Standard_Boolean          isForceCreate) const
955 {
956   if (myAttributes.IsBound(theFace))
957   {
958     theAttribute = myAttributes(theFace);
959     return Standard_True;
960   }
961   else if (isForceCreate)
962   {
963     theAttribute = new BRepMesh_FaceAttribute(myBoundaryVertices, myBoundaryPoints);
964     myAttributes.Bind(theFace, theAttribute);
965   }
966
967   return Standard_False;
968 }
969
970 //=======================================================================
971 //function : RemoveFaceAttribute
972 //purpose  : 
973 //=======================================================================
974 void BRepMesh_FastDiscret::RemoveFaceAttribute(const TopoDS_Face& theFace)
975 {
976   if (myAttributes.IsBound(theFace))
977     myAttributes.UnBind(theFace);
978 }