0023184: Bad face tessellation result
[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-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21
22 #include <BRepMesh_FastDiscret.ixx>
23
24 #include <BRepMesh_FastDiscretFace.hxx>
25 #include <BRepMesh_FaceAttribute.hxx>
26 #include <BRepMesh_DataStructureOfDelaun.hxx>
27 #include <BRepMesh_ClassifierPtr.hxx>
28 #include <BRepMesh_GeomTool.hxx>
29 #include <BRepMesh_PairOfPolygon.hxx>
30 #include <BRepMesh_DataMapOfShapePairOfPolygon.hxx>
31 #include <BRepMesh_DataMapIteratorOfDataMapOfShapePairOfPolygon.hxx>
32 #include <Geom_Plane.hxx>
33 #include <GeomAbs_IsoType.hxx>
34 #include <GeomAbs_SurfaceType.hxx>
35 #include <TopAbs.hxx>
36 #include <TColStd_HArray1OfReal.hxx>
37 #include <Precision.hxx>
38
39 #include <BRep_Builder.hxx>
40 #include <BRep_Tool.hxx>
41 #include <Poly_Triangulation.hxx>
42 #include <Poly_PolygonOnTriangulation.hxx>
43 #include <Poly_Connect.hxx>
44 #include <TColStd_SequenceOfInteger.hxx>
45 #include <TColStd_Array1OfInteger.hxx>
46 #include <TColStd_HArray1OfInteger.hxx>
47
48 #include <TColgp_Array1OfPnt.hxx>
49 #include <TColgp_Array1OfPnt2d.hxx>
50 #include <Precision.hxx>
51
52 #include <BRepAdaptor_Curve.hxx>
53 #include <BRepAdaptor_Surface.hxx>
54 #include <BRepAdaptor_HSurface.hxx>
55 #include <BRepTools.hxx>
56 #include <BndLib_Add3dCurve.hxx>
57 #include <BRepBndLib.hxx>
58 #include <Bnd_Box.hxx>
59 #include <TopoDS.hxx>
60 #include <TopExp.hxx>
61 #include <TopExp_Explorer.hxx>
62
63 #include <Geom2d_Curve.hxx>
64
65 #include <TColStd_DataMapOfIntegerInteger.hxx>
66 #include <BRepMesh_ShapeTool.hxx>
67 #include <ElSLib.hxx>
68 #include <Geom_Surface.hxx>
69 #include <Adaptor3d_IsoCurve.hxx>
70 #include <BRepMesh_IndexedMapOfVertex.hxx>
71 #include <Extrema_LocateExtPC.hxx>
72
73 #include <BRepMesh_ListOfXY.hxx>
74 #include <BRepMesh_ListIteratorOfListOfXY.hxx>
75
76 #include <TColStd_Array1OfInteger.hxx>
77 #include <BRepMesh_IDMapOfNodeOfDataStructureOfDelaun.hxx>
78 #include <Standard_ErrorHandler.hxx>
79 #include <Standard_Failure.hxx>
80 //#include <TColStd_DataMapOfInteger.hxx>
81 #include <TColGeom2d_SequenceOfCurve.hxx>
82 #include <TopTools_SequenceOfShape.hxx>
83 #include <NCollection_IncAllocator.hxx>
84
85 #include <BRep_ListIteratorOfListOfPointRepresentation.hxx>
86 #include <BRep_PointRepresentation.hxx>
87 #include <BRep_TVertex.hxx>
88 #include <TColStd_MapOfInteger.hxx>
89 #include <SortTools_ShellSortOfReal.hxx>
90 #include <TCollection_CompareOfReal.hxx>
91
92 #include <TopTools_HArray1OfShape.hxx>
93 #include <TopTools_ListIteratorOfListOfShape.hxx>
94
95 #include <vector>
96
97 #ifdef HAVE_TBB
98   // paralleling using Intel TBB
99   #include <tbb/parallel_for_each.h>
100 #endif
101
102 #define UVDEFLECTION 1.e-05
103
104 inline Standard_Real MaxFaceTol (const TopoDS_Face& theFace)
105 {
106   Standard_Real T, TMax = BRep_Tool::Tolerance(theFace);
107   TopExp_Explorer Ex;
108
109   for (Ex.Init(theFace,TopAbs_EDGE); Ex.More(); Ex.Next())
110   {
111     T = BRep_Tool::Tolerance(TopoDS::Edge(Ex.Current()));
112     if (T > TMax) TMax = T;
113   }
114
115   for (Ex.Init(theFace,TopAbs_VERTEX); Ex.More(); Ex.Next())
116   {
117     T = BRep_Tool::Tolerance(TopoDS::Vertex(Ex.Current()));
118     if (T > TMax) TMax = T;
119   }
120
121   return TMax;
122 }
123
124
125 //=======================================================================
126 //function : BRepMesh_FastDiscret
127 //purpose  : 
128 //=======================================================================
129 BRepMesh_FastDiscret::BRepMesh_FastDiscret(const Standard_Real    theDefle,
130                                            const Standard_Real    theAngl,
131                                            const Bnd_Box&         theBox,
132                                            const Standard_Boolean theWithShare,
133                                            const Standard_Boolean theInshape,
134                                            const Standard_Boolean theRelative,
135                                            const Standard_Boolean theShapetrigu) :
136   myAngle (theAngl),
137   myDeflection (theDefle),
138   myWithShare (theWithShare),
139   myInParallel (Standard_False),
140   myNbLocat (0),
141   myRelative (theRelative),
142   myShapetrigu (theShapetrigu), 
143   myInshape (theInshape)
144 {
145   myAllocator = new NCollection_IncAllocator(64000);
146   if(myRelative)
147     BoxMaxDimension(theBox, myDtotale);
148 }
149
150 //=======================================================================
151 //function : BRepMesh_FastDiscret
152 //purpose  : 
153 //=======================================================================
154
155 BRepMesh_FastDiscret::BRepMesh_FastDiscret(const Standard_Real    theDefle,
156                                            const TopoDS_Shape&    theShape,
157                                            const Bnd_Box&         theBox,
158                                            const Standard_Real    theAngl,
159                                            const Standard_Boolean theWithShare,
160                                            const Standard_Boolean theInshape,
161                                            const Standard_Boolean theRelative,
162                                            const Standard_Boolean theShapetrigu): 
163   myAngle (theAngl),
164   myDeflection (theDefle),
165   myWithShare (theWithShare),
166   myInParallel (Standard_False),
167   myNbLocat (0),
168   myRelative (theRelative),
169   myShapetrigu (theShapetrigu),
170   myInshape (theInshape)
171 {
172   myAllocator = new NCollection_IncAllocator(64000);
173   if(myRelative)
174     BoxMaxDimension(theBox, myDtotale);
175   Perform(theShape);
176 }
177
178 //=======================================================================
179 //function : SetParallel
180 //purpose  :
181 //=======================================================================
182 void BRepMesh_FastDiscret::SetParallel (const Standard_Boolean theInParallel)
183 {
184   myInParallel = theInParallel;
185 }
186
187 //=======================================================================
188 //function : IsParallel
189 //purpose  :
190 //=======================================================================
191 Standard_Boolean BRepMesh_FastDiscret::IsParallel() const
192 {
193   return myInParallel;
194 }
195
196 //=======================================================================
197 //function : BoxMaxDimension
198 //purpose  : 
199 //=======================================================================
200
201 void BRepMesh_FastDiscret::BoxMaxDimension(const Bnd_Box& theBox, Standard_Real& theMaxDim)
202 {
203   if(theBox.IsVoid())
204     return;
205   Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
206   theBox.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
207   theMaxDim = TXmax-TXmin;
208   const Standard_Real dy = TYmax-TYmin;
209   const Standard_Real dz = TZmax-TZmin;
210   if (dy > theMaxDim) theMaxDim = dy;
211   if (dz > theMaxDim) theMaxDim = dz;
212 }
213
214 //=======================================================================
215 //function : RelativeEdgeDeflection
216 //purpose  : 
217 //=======================================================================
218
219 Standard_Real BRepMesh_FastDiscret::RelativeEdgeDeflection(const TopoDS_Edge& theEdge,
220                                                                                                    const Standard_Real theDefle,
221                                                                                                      const Standard_Real theDTotale,
222                                                            Standard_Real& theDefCoef)
223 {
224   theDefCoef = 1.;
225   Standard_Real defedge = theDefle;
226   if(theEdge.IsNull())
227     return defedge;
228
229   Bnd_Box B;
230   BRepBndLib::Add(theEdge, B);
231   BoxMaxDimension(B, defedge);
232             
233   // adjusted in relation to the total size:
234   theDefCoef = theDTotale/(2*defedge);
235   if (theDefCoef < 0.5) theDefCoef = 0.5;
236   if (theDefCoef > 2.) theDefCoef = 2.;
237   defedge = theDefCoef * defedge * theDefle;
238
239   return defedge;
240 }
241
242 //=======================================================================
243 //function : Perform(shape)
244 //purpose  : 
245 //=======================================================================
246
247 void BRepMesh_FastDiscret::Perform(const TopoDS_Shape& theShape)
248 {
249   TopTools_IndexedDataMapOfShapeListOfShape anAncestors;
250   TopExp::MapShapesAndAncestors(theShape, TopAbs_EDGE, TopAbs_FACE, anAncestors);
251   std::vector<TopoDS_Face> aFaces;
252   for (TopExp_Explorer ex(theShape, TopAbs_FACE); ex.More(); ex.Next()) {
253     TopoDS_Face aF = TopoDS::Face(ex.Current());
254     Add(aF, anAncestors);
255     aFaces.push_back(aF);
256   }
257
258   if (myInParallel)
259   {
260   #ifdef HAVE_TBB
261     CreateMutexesForSubShapes(theShape, TopAbs_EDGE);
262     // mesh faces in parallel threads using TBB
263     tbb::parallel_for_each (aFaces.begin(), aFaces.end(), *this);
264   #else
265     // alternative parallelization not yet available
266     for (std::vector<TopoDS_Face>::iterator it(aFaces.begin()); it != aFaces.end(); it++)
267       Process (*it);
268   #endif
269     RemoveAllMutexes();
270   }
271   else
272   {
273     for (std::vector<TopoDS_Face>::iterator it(aFaces.begin()); it != aFaces.end(); it++)
274       Process (*it);
275   }
276 }
277
278
279 //=======================================================================
280 //function : Process
281 //purpose  : 
282 //=======================================================================
283
284 void BRepMesh_FastDiscret::Process(const TopoDS_Face& theFace) const
285 {
286   //cout << "START face " << theFace.TShape().operator->() << endl << flush;
287   Handle(BRepMesh_FaceAttribute) fattribute;
288   if ( GetFaceAttribute (theFace, fattribute) ) 
289   {
290     BRepMesh_FastDiscretFace aTool (GetAngle(), WithShare());
291     aTool.Add (theFace, fattribute, GetMapOfDefEdge(), myMutexProvider);
292   }
293   //cout << "END   face " << theFace.TShape().operator->() << endl << flush;
294 }
295
296 //=======================================================================
297 //function : Add(face)
298 //purpose  : 
299 //=======================================================================
300
301 #define MESH_FAILURE(theface)        \
302 {                                    \
303   myFacestate = BRepMesh_Failure;    \
304   myNottriangulated.Append(theface); \
305   return;                            \
306 }
307
308 void BRepMesh_FastDiscret::Add(const TopoDS_Face& theface,
309                                const TopTools_IndexedDataMapOfShapeListOfShape& theAncestors)
310 {
311 #ifndef DEB_MESH
312   try
313   {
314     OCC_CATCH_SIGNALS
315 #endif
316   TopoDS_Face face = theface;
317   BRepTools::Update(face);
318   face.Orientation(TopAbs_FORWARD);
319   myStructure.Nullify();
320   Handle(NCollection_IncAllocator) anAlloc = Handle(NCollection_IncAllocator)::DownCast(myAllocator);
321   anAlloc->Reset(Standard_False);  
322   myStructure=new BRepMesh_DataStructureOfDelaun(anAlloc);
323
324   Standard_Real aUmin, aVmin, aUmax, aVmax;
325   BRepTools::UVBounds (theface, aUmin, aUmax, aVmin, aVmax);
326   Standard_Real aTolU = (aUmax - aUmin) * UVDEFLECTION;
327   Standard_Real aTolV = (aVmax - aVmin) * UVDEFLECTION;
328   myStructure->Data().SetCellSize ( 14 * aTolU, 14 * aTolV );
329   myStructure->Data().SetTolerance( aTolU, aTolV );
330
331   BRepAdaptor_Surface  BS(face, Standard_False);
332   Handle(BRepAdaptor_HSurface) gFace = new BRepAdaptor_HSurface(BS);
333   
334   GeomAbs_SurfaceType thetype;
335   thetype = BS.GetType();
336
337   gp_Pnt2d uvFirst, uvLast;
338
339   Handle(Poly_Triangulation) T;
340   TopLoc_Location loc;
341
342   if (!myWithShare) {          
343     myVertices.Clear();
344     myEdges.Clear();
345   }
346
347   myVemap.Clear();
348   myLocation2d.Clear();
349   myInternaledges.Clear();
350
351   Standard_Integer i;
352   i = 1;
353
354   Standard_Real defedge, defface;
355   Standard_Integer nbEdge = 0;
356   Standard_Real savangle = myAngle;
357   Standard_Real cdef;
358   Standard_Real maxdef = 2.* MaxFaceTol(theface);
359   defface = 0.;
360
361   if (!myRelative) defface = Max(myDeflection, maxdef);
362
363   TColStd_SequenceOfReal aFSeq, aLSeq;
364   TColGeom2d_SequenceOfCurve aCSeq;
365   TopTools_SequenceOfShape aShSeq;
366
367   TopoDS_Iterator exW(face);
368
369   for (; exW.More(); exW.Next()) {
370     const TopoDS_Shape& aWire = exW.Value();
371     if (aWire.ShapeType() != TopAbs_WIRE)
372       continue;
373     TopoDS_Iterator ex(aWire);
374     for(; ex.More(); ex.Next()) {
375       const TopoDS_Edge& edge = TopoDS::Edge(ex.Value());
376       nbEdge++;
377       if (!myMapdefle.IsBound(edge)) {
378         if (myRelative) {
379           if (myEdges.IsBound(edge)) {
380             const BRepMesh_PairOfPolygon& pair = myEdges.Find(edge);
381             const Handle(Poly_PolygonOnTriangulation)& P = pair.First();
382             defedge = P->Deflection();
383           }
384           else {
385             defedge = RelativeEdgeDeflection(edge, myDeflection, myDtotale, cdef);
386             myAngle = savangle * cdef;
387           }
388           defface = defface + defedge;
389           defface = Max(maxdef, defface);
390         }
391         else defedge = myDeflection;
392     
393         defedge = Max(maxdef, defedge);
394         defedge = Max(UVDEFLECTION , defedge);
395         myMapdefle.Bind(edge, defedge);
396       }
397       else{
398         defedge = myMapdefle(edge);
399         if (myRelative) {defface = defface + defedge; defface = Max(maxdef, defface);}
400       }
401       Standard_Real f1,l1;
402       Handle(Geom2d_Curve) C = BRep_Tool::CurveOnSurface(edge, face, f1, l1);
403       if (C.IsNull()) continue;
404
405       aFSeq.Append(f1);
406       aLSeq.Append(l1);
407       aCSeq.Append(C);
408       aShSeq.Append(edge);
409       Add(edge, face, gFace, C, theAncestors, defedge, f1, l1);
410       myAngle = savangle;
411     }
412   }
413
414   if (nbEdge == 0 || myVemap.Extent() < 3)
415   {
416       MESH_FAILURE(theface);
417   }
418
419   if (myRelative ) defface = defface / nbEdge;
420   else             defface = myDeflection;
421
422   if (myWithShare) defface = Max(maxdef, defface);
423   
424   T = BRep_Tool::Triangulation(face, loc);
425
426   if (!myShapetrigu || T.IsNull()) {
427     
428     Standard_Real xCur, yCur;
429     Standard_Real maxX, minX, maxY, minY;
430     minX=minY=1.e100;
431     maxX=maxY=-1.e100;
432     
433     Standard_Integer ipn = 0;
434     Standard_Integer i1 =1;
435     for (i1 = 1; i1 <= myVemap.Extent(); i1++) {
436       const BRepMesh_Vertex& aV = myStructure->GetNode(myVemap.FindKey(i1));
437       ipn++;
438       xCur=aV.Coord().X();
439       yCur=aV.Coord().Y();
440       minX=Min(xCur, minX);
441       maxX=Max(xCur, maxX);
442       minY=Min(yCur, minY);
443       maxY=Max(yCur, maxY);
444     }
445     Standard_Real myumin = minX;
446     Standard_Real myumax = maxX;
447     Standard_Real myvmin = minY;
448     Standard_Real myvmax = maxY;
449     
450     const Standard_Real umin = BS.FirstUParameter();
451     const Standard_Real umax = BS.LastUParameter();
452     const Standard_Real vmin = BS.FirstVParameter();
453     const Standard_Real vmax = BS.LastVParameter();
454
455     if (myumin < umin || myumax > umax)
456     {
457       if (BS.IsUPeriodic())
458       {
459         if ((myumax - myumin) > (umax - umin))
460         {
461           myumax = myumin + (umax - umin);
462         }
463       }
464       else
465       {
466         if (umin > myumin) myumin = umin;
467         if (umax < myumax) myumax = umax;
468       }
469     }
470
471     if (myvmin < vmin || myvmax > vmax)
472     {
473       if (BS.IsVPeriodic())
474       {
475         if ((myvmax - myvmin) > (vmax - vmin))
476         {
477           myvmax = myvmin + (vmax - vmin);
478         }
479       }
480       else
481       {
482         if (vmin > myvmin) myvmin = vmin;
483         if (vmax < myvmax) myvmax = vmax;
484       }
485     }
486
487     // fast verification of the validity of calculated limits. If wrong, 
488     // sure a problem of pcurve.
489     if (thetype == GeomAbs_BezierSurface &&
490         (myumin < -0.5 || myumax > 1.5 || myvmin < -0.5 || myvmax > 1.5))
491     {
492       MESH_FAILURE(theface);
493     }
494  
495     //define parameters for correct parametrics
496     
497     Standard_Real deltaX = 1.0;
498     Standard_Real deltaY = 1.0;
499     Standard_Integer nbVertices = myVemap.Extent();
500     const Standard_Real tolclass = Precision::PConfusion(); //0.03*Max(myumax-myumin, myvmax-myvmin);
501     
502     BRepMesh_ClassifierPtr classifier ( 
503       new BRepMesh_Classifier(face, tolclass,  myInternaledges, myVemap, 
504                               myStructure, myumin, myumax, myvmin, myvmax) );   
505   
506     myFacestate = classifier->State();
507     if (myFacestate == BRepMesh_SelfIntersectingWire)
508     {
509       Standard_Integer nbmaill = 0;
510       Standard_Real eps = Precision::Confusion();
511       while (nbmaill < 5 && myFacestate != BRepMesh_ReMesh)
512       {
513         nbmaill++;
514         
515         //clear the structure of links
516         myStructure.Nullify();
517         myStructure = new BRepMesh_DataStructureOfDelaun(anAlloc);
518         
519         myVemap.Clear();
520         myLocation2d.Clear();
521         myInternaledges.Clear();
522
523         Standard_Integer j1;
524         for(j1 = 1; j1 <= aShSeq.Length(); j1++)
525         {
526           const TopoDS_Edge& edge = TopoDS::Edge(aShSeq.Value(j1));
527           if (myEdges.IsBound(edge))
528           {
529             myEdges.UnBind(edge);
530             myInternaledges.UnBind(edge);
531           }
532         }
533         
534         
535         for( j1 = 1; j1 <= aShSeq.Length(); j1++)
536         {
537           const TopoDS_Edge& edge = TopoDS::Edge(aShSeq.Value(j1));
538           defedge = myMapdefle(edge) / 3.;
539           defedge = Max(defedge, eps);
540           myMapdefle.Bind(edge, defedge);
541           const Handle(Geom2d_Curve)& C = aCSeq.Value(j1);
542           Add(edge, face, gFace, C, theAncestors, defedge, aFSeq.Value(j1), aLSeq.Value(j1));
543         }
544
545         classifier.Nullify();
546
547         classifier = new BRepMesh_Classifier(face, tolclass, myInternaledges, myVemap,
548                                              myStructure, myumin, myumax, myvmin, myvmax);
549
550         if (classifier->State() == BRepMesh_NoError)
551         {
552           myFacestate = BRepMesh_ReMesh;
553         }
554         nbVertices = myVemap.Extent();
555       }
556     }
557     
558     if (myFacestate != BRepMesh_NoError && myFacestate != BRepMesh_ReMesh)
559     {
560       myNottriangulated.Append(face);
561       classifier.Nullify();
562       return;
563     }
564     
565     // try to find the real length:
566     // akm (bug OCC16) : We must calculate these measures in non-singular
567     //     parts of face. Let's try to compute average value of three
568     //     (umin, (umin+umax)/2, umax), and respectively for v.
569     //                 vvvvv
570     Standard_Real longu = 0.0, longv = 0.0; //, last , first;
571     gp_Pnt P11, P12, P21, P22, P31, P32;
572     
573     Standard_Real du = (myumax-myumin)/20;
574     Standard_Real dv = (myvmax-myvmin)/20;
575     Standard_Real dfuave=(myumin+myumax)/2, dfucur;
576     Standard_Real dfvave=(myvmin+myvmax)/2, dfvcur;
577     // U loop
578     BS.D0 (myumin, myvmin, P11);
579     BS.D0 (myumin, dfvave, P21);
580     BS.D0 (myumin, myvmax, P31);
581     for (i1=1, dfucur=myumin+du; i1 <= 20; i1++, dfucur+=du)  {
582       BS.D0 (dfucur, myvmin, P12);
583       BS.D0 (dfucur, dfvave, P22);
584       BS.D0 (dfucur, myvmax, P32);
585       longu += ( P11.Distance(P12) + P21.Distance(P22) + P31.Distance(P32) );
586       P11 = P12;
587       P21 = P22;
588       P31 = P32;
589     }
590     // V loop
591     BS.D0(myumin, myvmin, P11);
592     BS.D0(dfuave, myvmin, P21);
593     BS.D0(myumax, myvmin, P31);
594     for (i1=1, dfvcur=myvmin+dv; i1 <= 20; i1++, dfvcur+=dv) {    
595       BS.D0 (myumin, dfvcur, P12);
596       BS.D0 (dfuave, dfvcur, P22);
597       BS.D0 (myumax, dfvcur, P32);
598       longv += ( P11.Distance(P12) + P21.Distance(P22) + P31.Distance(P32) );
599       P11 = P12;
600       P21 = P22;
601       P31 = P32;
602     }
603     longu /= 3.;
604     longv /= 3.;
605     // akm (bug OCC16) ^^^^^
606     
607     if (longu <= 1.e-16 || longv <= 1.e-16) {
608       //yes, it is seen!!
609 #ifdef DEB_MESH_CHRONO
610       chMaillEdges.Stop();
611       MESH_CHRONO_TSTOP(thetype);
612 #endif
613       MESH_FAILURE(theface);
614     }
615
616
617     if (thetype == GeomAbs_Torus)  {
618       gp_Torus Tor = BS.Torus();
619       Standard_Real r = Tor.MinorRadius(), R = Tor.MajorRadius();
620       Standard_Real Du, Dv;//, pasu, pasv;
621       
622       Dv = Max(1.0e0 - (defface/r),0.0e0) ;
623       Standard_Real oldDv = 2.0 * ACos (Dv);
624       oldDv = Min(oldDv, myAngle);
625       Dv  =  0.9*oldDv; //TWOTHIRD * oldDv;
626       Dv = oldDv;
627       
628       Standard_Integer nbV = Max((Standard_Integer)((myvmax-myvmin)/Dv), 2);
629       Dv = (myvmax-myvmin)/(nbV+1);
630       
631       Standard_Real ru = R + r;
632       if (ru > 1.e-16) {
633         Du = Max(1.0e0 - (defface/ru),0.0e0);
634         Du  = (2.0 * ACos (Du));
635         Du = Min(Du, myAngle);
636         Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
637         if(aa < gp::Resolution())
638           return; 
639
640         Du = Du * Min(oldDv, Du) / aa;
641       }
642       else Du = Dv;     
643   
644       Standard_Integer nbU = Max((Standard_Integer)((myumax-myumin)/Du), 2);
645       nbU = Max(nbU, (Standard_Integer)(nbV*(myumax-myumin)*R/((myvmax-myvmin)*r)/5.));
646       
647       Du = (myumax-myumin)/(nbU+1);
648       //-- DeltaX and DeltaY are chosen so that to avoid "jumping" 
649       //-- of points on the grid
650       deltaX = Du;
651       deltaY = Dv;
652     }
653     else if (thetype == GeomAbs_Cylinder) {
654       /*Standard_Real aMax = Max(myumax,myvmax);
655       deltaX = 0.01;  //1./aMax; //0.01;
656       deltaY = 1.0;*/
657       gp_Cylinder Cyl = BS.Cylinder();
658       Standard_Real R = Cyl.Radius();
659       // Calculate parameters for iteration in U direction
660       Standard_Real Du = 1.0 - (defface/R);
661       if (Du < 0.0) Du = 0.0;
662       Du = 2.0 * ACos (Du);
663       if (Du > GetAngle()) Du = GetAngle();
664       deltaX = Du/longv;
665       deltaY = 1.;
666     }
667     else {
668       deltaX = (myumax-myumin)/longu;
669       deltaY = (myvmax-myvmin)/longv;
670     }    
671
672     if(!myMapattrib.IsBound(theface))
673     {
674       Handle(BRepMesh_FaceAttribute) aFA = new BRepMesh_FaceAttribute();
675       aFA->GetDefFace() = defface;
676       aFA->GetUMax() = myumax;
677       aFA->GetVMax() = myvmax;
678       aFA->GetUMin() = myumin;
679       aFA->GetVMin() = myvmin;
680       aFA->GetDeltaX() = deltaX;
681       aFA->GetDeltaY() = deltaY;
682       aFA->GetMinX() = minX;
683       aFA->GetMinY() = minY;
684       aFA->GetClassifier() = classifier;
685       myMapattrib.Bind(theface, aFA);
686     }
687
688     //Standard_Integer nbVertices = myVemap.Extent();
689     T = new Poly_Triangulation(nbVertices, 1, Standard_True);
690     TColgp_Array1OfPnt&  Nodes = T->ChangeNodes();
691     TColgp_Array1OfPnt2d& Nodes2d = T->ChangeUVNodes();
692     
693     Standard_Integer index;
694     for (i = 1; i <= nbVertices; i++) {
695       index = myVemap.FindKey(i);
696       Nodes(i) = Pnt(index);
697       Nodes2d(i).SetXY(Vertex(index).Coord());
698     }
699
700     // storage of triangulation in the BRep.
701     //TopLoc_Location loc = face.Location();
702     if (!loc.IsIdentity()) {
703       gp_Trsf tr = loc.Transformation();
704       tr.Invert();
705       for (i = Nodes.Lower(); i <= Nodes.Upper(); i++) 
706         Nodes(i).Transform(tr);
707     }
708
709     BRep_Builder B;
710     B.UpdateFace(face, T);
711
712     BRepMesh_DataMapIteratorOfDataMapOfShapePairOfPolygon It(myInternaledges);
713     for (; It.More(); It.Next()) {
714       const BRepMesh_PairOfPolygon& pair = It.Value();
715       const Handle(Poly_PolygonOnTriangulation)& NOD1 = pair.First();
716       const Handle(Poly_PolygonOnTriangulation)& NOD2 = pair.Last();
717       if ( NOD1 == NOD2 )
718         B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, T, loc);
719       else
720         B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, NOD2, T, loc);
721     }
722   }
723
724 #ifndef DEB_MESH
725   }
726   catch(Standard_Failure)
727   {
728     MESH_FAILURE(theface);
729   }
730 #endif // DEB_MESH
731   myStructure.Nullify();
732 }
733
734 //=======================================================================
735 //function : splitSegment
736 //purpose  :
737 //=======================================================================
738 static void splitSegment( BRepMesh_GeomTool&    theGT, 
739                           const Handle(Geom_Surface)& theSurf, 
740                           const Handle(Geom2d_Curve)& theCurve2d,
741                           const BRepAdaptor_Curve&    theBAC,
742                           const Standard_Real         theSquareEDef,
743                           const Standard_Real         theFirst,
744                           const Standard_Real         theLast,
745                           const Standard_Integer      theNbIter)
746 {
747   //limit ineration depth
748   if(theNbIter > 10)
749     return;
750   gp_Pnt2d uvf, uvl, uvm;
751   gp_Pnt   P3dF, P3dL, midP3d, midP3dFromSurf;
752   Standard_Real midpar;
753   
754   if(Abs(theLast - theFirst) < 2*Precision::PConfusion())
755     return;
756
757   theCurve2d->D0(theFirst, uvf);
758   theCurve2d->D0(theLast, uvl);
759
760   P3dF = theSurf->Value(uvf.X(), uvf.Y());
761   P3dL = theSurf->Value(uvl.X(), uvl.Y());
762   
763   if(P3dF.SquareDistance(P3dL) < theSquareEDef)
764     return;
765
766   uvm = gp_Pnt2d((uvf.XY() + uvl.XY())*0.5);
767   midP3dFromSurf = theSurf->Value(uvm.X(), uvm.Y());
768
769   gp_XYZ aVec = P3dL.XYZ()-P3dF.XYZ();
770   aVec.Normalize();
771
772   gp_XYZ Vec1 = midP3dFromSurf.XYZ() - P3dF.XYZ();
773   Standard_Real aModulus = Vec1.Dot(aVec);
774   gp_XYZ aProj = aVec*aModulus;
775   gp_XYZ aDist = Vec1 - aProj;
776     
777   if(aDist.SquareModulus() < theSquareEDef)
778     return;
779
780   midpar = (theFirst + theLast) * 0.5;
781   theBAC.D0(midpar, midP3d);
782   theGT.AddPoint(midP3d, midpar, Standard_False);
783
784   splitSegment(theGT, theSurf, theCurve2d, theBAC, theSquareEDef, theFirst, midpar, theNbIter+1);
785   splitSegment(theGT, theSurf, theCurve2d, theBAC, theSquareEDef, midpar, theLast,  theNbIter+1); 
786 }
787
788 //=======================================================================
789 //function : Add
790 //purpose  : 
791 //=======================================================================
792 void BRepMesh_FastDiscret::Add( const TopoDS_Edge&                  theEdge, 
793                                 const TopoDS_Face&                  theFace, 
794                                 const Handle(BRepAdaptor_HSurface)& theGFace,
795                                 const Handle(Geom2d_Curve)&         theC2d,
796                                 const TopTools_IndexedDataMapOfShapeListOfShape& theAncestors,
797                                 const Standard_Real                 theDefEdge,
798                                 const Standard_Real                 theFirst,
799                                 const Standard_Real                 theLast)
800 {
801   const TopAbs_Orientation orEdge = theEdge.Orientation();
802   if (orEdge == TopAbs_EXTERNAL) return;
803
804   const Standard_Boolean isEdgeBound = myEdges.IsBound(theEdge);
805
806   if (!isEdgeBound)
807   {
808     if (Update(theEdge, theFace, theC2d, theDefEdge, theFirst, theLast))
809     {
810       return;
811     }
812   }
813
814   TopoDS_Vertex pBegin, pEnd;
815   TopExp::Vertices(theEdge, pBegin, pEnd);
816   if (pBegin.IsNull() || pEnd.IsNull())
817   {
818     return;
819   }
820
821   Standard_Real wFirst, wLast;
822   BRep_Tool::Range(theEdge, theFace, wFirst, wLast);
823
824   gp_Pnt2d uvFirst, uvLast;
825   BRep_Tool::UVPoints(theEdge, theFace, uvFirst, uvLast);
826
827   const Standard_Boolean sameUV = uvFirst.IsEqual(uvLast, Precision::PConfusion());
828
829   //Control vertexes tolerances
830   gp_Pnt pFirst = theGFace->Value(uvFirst.X(), uvFirst.Y());
831   gp_Pnt pLast  = theGFace->Value(uvLast.X(), uvLast.Y());
832
833   Standard_Real mindist = 10. * Max(pFirst.Distance(BRep_Tool::Pnt(pBegin)), 
834                                     pLast.Distance(BRep_Tool::Pnt(pEnd)));
835
836   if(mindist < BRep_Tool::Tolerance(pBegin) ||
837      mindist < BRep_Tool::Tolerance(pEnd) ) mindist = theDefEdge;
838
839   // control of degenerated non-coded edges.
840
841   Standard_Boolean degener = BRep_Tool::Degenerated(theEdge);
842
843   if (!degener)
844   {
845     if (pBegin.IsSame(pEnd))
846     {
847       // calculation of the length of the edge in 3D
848       Standard_Real longueur = 0.0;
849       Standard_Real du = (wLast-wFirst)/20;
850       gp_Pnt P1, P2;
851       BRepAdaptor_Curve BC(theEdge);
852       BC.D0(wFirst, P1);
853       Standard_Real tolV = BRep_Tool::Tolerance(pBegin);
854       Standard_Real tolV2 = 1.2*tolV;
855       for (Standard_Integer l = 1; l <= 20; l++) {
856         BC.D0(wFirst + l*du, P2);
857         longueur += P1.Distance(P2);
858         if (longueur > tolV2) break;
859         P1 = P2;
860       }
861
862       if (longueur < tolV2)
863       {
864         degener = Standard_True;
865       }
866     }
867   }
868
869   // Correct UV points
870   if (sameUV)
871   {
872     // 1. is it really sameUV without being degenerated
873     gp_Pnt2d uvF, uvL;
874     theC2d->D0(theFirst, uvF);
875     theC2d->D0(theLast, uvL);
876     if (!uvFirst.IsEqual(uvF, Precision::PConfusion()))
877     {
878       uvFirst = uvF;
879     }
880     if (!uvLast.IsEqual(uvL, Precision::PConfusion()))
881     {
882       uvLast = uvL;
883     }
884   }
885
886   gp_XY theUV;
887
888   // Process first vertex
889   Standard_Integer ipf;
890   if (myVertices.IsBound(pBegin))
891   {
892     ipf = myVertices.Find(pBegin);
893   }
894   else
895   {
896     if (sameUV && myVertices.IsBound(pEnd))
897     {
898       ipf = myVertices.Find(pEnd);
899     }
900     else
901     {
902       myNbLocat++;
903       ipf = myNbLocat;
904       myLocation3d.Bind(ipf, BRep_Tool::Pnt(pBegin));
905     }
906     myVertices.Bind(pBegin, ipf);
907   }
908   theUV = BRepMesh_FastDiscretFace::FindUV(pBegin, uvFirst, ipf, theGFace, mindist, myLocation2d);
909   BRepMesh_Vertex vf(theUV, ipf, BRepMesh_Frontier);
910   Standard_Integer ivf = myStructure->AddNode(vf);
911
912   // Process last vertex
913   Standard_Integer ipl;
914   if (pEnd.IsSame(pBegin))
915   {
916     ipl = ipf;
917   }
918   else
919   {
920     if (myVertices.IsBound(pEnd))
921     {
922       ipl = myVertices.Find(pEnd);
923     }
924     else
925     {
926       if (sameUV)
927       {
928         ipl = ipf;
929       }
930       else
931       {
932         myNbLocat++;
933         ipl = myNbLocat;
934         myLocation3d.Bind(ipl, BRep_Tool::Pnt(pEnd));
935       }
936       myVertices.Bind(pEnd,ipl);
937     }
938   }
939   theUV = BRepMesh_FastDiscretFace::FindUV(pEnd, uvLast, ipl, theGFace, mindist, myLocation2d);
940   BRepMesh_Vertex vl(theUV, ipl, BRepMesh_Frontier);
941   Standard_Integer ivl= myStructure->AddNode(vl);
942
943   Standard_Integer isvf = myVemap.FindIndex(ivf);
944   if (isvf == 0) isvf = myVemap.Add(ivf);
945   Standard_Integer isvl = myVemap.FindIndex(ivl);
946   if (isvl == 0) isvl = myVemap.Add(ivl);
947   
948   Standard_Real otherdefedge = 0.5*theDefEdge;
949   gp_Pnt2d uv;
950   BRepMesh_Vertex v2;
951   gp_Pnt   P3d;
952
953   Handle(Poly_PolygonOnTriangulation) P1;
954
955   if (!isEdgeBound)
956   {
957     Handle(Poly_PolygonOnTriangulation) P2;
958
959     if (degener)
960     {
961       // creation anew:
962       TColStd_Array1OfInteger Nodes(1, 2), NodInStruct(1, 2);
963       TColStd_Array1OfReal Param(1, 2);
964       
965       NodInStruct(1) = ipf;
966       Nodes(1) = isvf;
967       Param(1) = wFirst;
968       
969       NodInStruct(2) = ipl;
970       Nodes(2) = isvl;
971       Param(2) = wLast;
972
973       P1 = new Poly_PolygonOnTriangulation(Nodes, Param);
974       P2 = new Poly_PolygonOnTriangulation(NodInStruct, Param);
975     }
976     else
977     {
978       if (orEdge == TopAbs_INTERNAL) otherdefedge *= 0.5;
979       
980       BRepAdaptor_Curve cons;
981       Standard_Boolean isSameParam = BRep_Tool::SameParameter(theEdge);
982       if (isSameParam)
983       {
984         cons.Initialize(theEdge);
985       }
986       else
987       {
988         cons.Initialize(theEdge, theFace);
989       }
990
991       TopLoc_Location L;
992       Standard_Integer nbpmin = 2;
993       const GeomAbs_CurveType aCurveType = cons.GetType();
994       if ( aCurveType == GeomAbs_Circle )
995         nbpmin = 4; //OCC287
996
997       BRepMesh_GeomTool GT(cons, wFirst, wLast, 0.5 * myAngle, otherdefedge, nbpmin);
998
999       if ( aCurveType == GeomAbs_BSplineCurve )
1000       {
1001         const Standard_Integer aNbInt = cons.NbIntervals( GeomAbs_C1 );
1002         if ( aNbInt > 0 )
1003         {
1004           TColStd_Array1OfReal anIntervals( 1, aNbInt + 1 );
1005           cons.Intervals( anIntervals, GeomAbs_C1 );
1006           for ( Standard_Integer aIntIt = 1; aIntIt <= aNbInt; ++aIntIt )
1007           {
1008             const Standard_Real& aStartInt = anIntervals.Value( aIntIt );
1009             const Standard_Real& anEndInt  = anIntervals.Value( aIntIt + 1 );
1010
1011             BRepMesh_GeomTool aDetalizator( cons, aStartInt, anEndInt,
1012               0.5 * myAngle, otherdefedge, nbpmin );
1013
1014             Standard_Integer aNbAddNodes = aDetalizator.NbPoints();
1015             for ( Standard_Integer aNodeIt = 1; aNodeIt <= aNbAddNodes; ++aNodeIt )
1016             {
1017               Standard_Real aParam;
1018               gp_Pnt        aPoint3d;
1019               gp_Pnt2d      aPoint2d;
1020               aDetalizator.Value( cons, theGFace, aNodeIt, aParam, aPoint3d, aPoint2d );
1021               GT.AddPoint( aPoint3d, aParam, Standard_False );
1022             }
1023           }
1024         }
1025       }
1026
1027       // PTv, chl/922/G9, Take into account internal vertices
1028       // it is necessary for internal edges, which do not split other edges, by their vertex
1029       TopoDS_Iterator exV(theEdge);
1030       for ( ; exV.More(); exV.Next() )
1031       {
1032         TopoDS_Vertex aIntV = TopoDS::Vertex(exV.Value());
1033         if ( aIntV.Orientation() == TopAbs_INTERNAL )
1034          GT.AddPoint(BRep_Tool::Pnt(aIntV),
1035                      BRep_Tool::Parameter(aIntV, theEdge),
1036                      Standard_True);
1037       }
1038
1039       Standard_Integer i; 
1040       Standard_Integer nbnodes = GT.NbPoints();
1041       //Check deflection in 2d space for improvement of edge tesselation.
1042       if( isSameParam && nbnodes > 1)
1043       {
1044         Standard_Real aSquareEdgeDef = otherdefedge * otherdefedge;
1045         const TopTools_ListOfShape& lf = theAncestors.FindFromKey(theEdge);
1046         TopTools_ListIteratorOfListOfShape itl(lf);
1047         for (; itl.More(); itl.Next()) {
1048           const TopoDS_Face& aFace = TopoDS::Face (itl.Value());          
1049
1050           TopLoc_Location aLoc;
1051           Standard_Real aF, aL;
1052           Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace, aLoc);
1053           const Handle(Standard_Type)& aType = aSurf->DynamicType();
1054           if(aType == STANDARD_TYPE(Geom_Plane))
1055             continue;
1056           Handle(Geom2d_Curve) aCurve2d = BRep_Tool::CurveOnSurface(theEdge, aFace, aF, aL);
1057           if(Abs(aF-wFirst)>Precision::PConfusion()||Abs(aL-wLast)>Precision::PConfusion())
1058             continue;
1059           
1060           gp_Pnt2d uvf;
1061           Standard_Real parf;
1062           nbnodes = GT.NbPoints();
1063           TColStd_Array1OfReal aParamArray(1, nbnodes);
1064           for (i = 1; i <= nbnodes; i++)
1065           {          
1066             GT.Value(cons, theGFace, i, parf, P3d, uvf);            
1067             aParamArray.SetValue(i, parf);
1068           }
1069           for (i = 1; i < nbnodes; i++)
1070           {
1071             splitSegment(GT, aSurf, aCurve2d, cons, aSquareEdgeDef, aParamArray(i), aParamArray(i+1), 1);
1072           }
1073         }
1074       }
1075
1076       // Creation of polygons on triangulation:
1077       Standard_Real puv;
1078       nbnodes = GT.NbPoints();
1079
1080       TColStd_Array1OfInteger Nodes(1, nbnodes);
1081       TColStd_Array1OfInteger NodInStruct(1, nbnodes);
1082       TColStd_Array1OfReal Param(1, nbnodes);
1083       
1084       // processing of the 1st point:
1085       Nodes(1) = isvf;
1086       NodInStruct(1) = ipf;
1087       Param(1) = wFirst;
1088       
1089       // last point:
1090       Nodes(nbnodes) = isvl;
1091       NodInStruct(nbnodes) = ipl;
1092       Param(nbnodes) = wLast;
1093
1094       if(nbnodes > 2)
1095       {
1096         Standard_Integer iv2;
1097         for (i = 2; i < GT.NbPoints(); i++)
1098         {
1099           // Record 3d point
1100           GT.Value(cons, theGFace, i, puv, P3d, uv);
1101           myNbLocat++;
1102           myLocation3d.Bind(myNbLocat, P3d);
1103           NodInStruct(i) = myNbLocat;
1104           // Record 2d point
1105           v2.Initialize(uv.Coord(), myNbLocat, BRepMesh_OnCurve);
1106           iv2=myStructure->AddNode(v2);
1107           
1108           Standard_Integer isv = myVemap.FindIndex(iv2);
1109           if (isv == 0) isv = myVemap.Add(iv2);
1110           Nodes(i) = isv;
1111           NodInStruct(i) = myNbLocat;
1112           Param(i) = puv;
1113     
1114           if (orEdge == TopAbs_FORWARD)
1115             myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Frontier));
1116           else if (orEdge == TopAbs_REVERSED)
1117             myStructure->AddLink(BRepMesh_Edge(iv2, ivf, BRepMesh_Frontier));
1118           else if (orEdge == TopAbs_INTERNAL)
1119             myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Fixed));
1120           ivf = iv2;
1121         }
1122       }
1123
1124       P1 = new Poly_PolygonOnTriangulation(Nodes, Param);
1125       P2 = new Poly_PolygonOnTriangulation(NodInStruct, Param);
1126     }
1127
1128     P2->Deflection(otherdefedge);
1129     BRepMesh_PairOfPolygon pair;
1130     pair.Append(P2);
1131     myEdges.Bind(theEdge, pair);
1132
1133     if (ivf != ivl) {
1134       if (orEdge == TopAbs_FORWARD)
1135         myStructure->AddLink(BRepMesh_Edge(ivf, ivl, BRepMesh_Frontier));
1136       else if (orEdge == TopAbs_REVERSED)
1137         myStructure->AddLink(BRepMesh_Edge(ivl, ivf, BRepMesh_Frontier));
1138       else if (orEdge == TopAbs_INTERNAL)
1139         myStructure->AddLink(BRepMesh_Edge(ivf, ivl, BRepMesh_Fixed));
1140     }
1141       
1142
1143   }
1144   // If this Edge has been already checked and it is not degenerated, 
1145   // the points of the polygon calculated at the first check are retrieved :
1146   else
1147   {
1148     if (degener)
1149     {
1150       // Create 2d polygon for degenerated edge
1151       TColStd_Array1OfInteger Nodes(1, 2);
1152       TColStd_Array1OfReal PPar(1, 2);
1153       
1154       Nodes(1) = isvf;
1155       PPar(1) = wFirst;
1156       
1157       Nodes(2) = isvl;
1158       PPar(2) = wLast;
1159       
1160       P1 = new Poly_PolygonOnTriangulation(Nodes, PPar);
1161     }
1162     else
1163     {
1164       // retrieve the polygone:
1165       const BRepMesh_PairOfPolygon& pair = myEdges.Find(theEdge);
1166       const Handle(Poly_PolygonOnTriangulation)& P = pair.First();
1167       const TColStd_Array1OfInteger& NOD = P->Nodes();
1168       Handle(TColStd_HArray1OfReal) Par = P->Parameters();
1169
1170       // creation anew:
1171       const Standard_Integer nbnodes = NOD.Length();
1172       TColStd_Array1OfInteger Nodes(1, nbnodes);
1173       TColStd_Array1OfReal PPar(1, nbnodes);
1174       Standard_Integer iv2;
1175       
1176       Nodes(1) = isvf;
1177       PPar(1) = wFirst;
1178       
1179       Nodes(nbnodes) = isvl;
1180       PPar(nbnodes) = wLast;
1181       
1182       if (nbnodes > 2)
1183       {
1184         Standard_Integer i;
1185         if (BRep_Tool::SameParameter(theEdge))
1186         {
1187           for (i = 2; i < nbnodes; i++)
1188           {
1189             const Standard_Real puv = Par->Value(i);
1190             theC2d->D0(puv, uv);
1191             v2.Initialize(uv.Coord(), NOD(i), BRepMesh_OnCurve);
1192             iv2 = myStructure->AddNode(v2);
1193   
1194             Standard_Integer isv = myVemap.FindIndex(iv2);
1195             if (isv == 0) isv = myVemap.Add(iv2);
1196             Nodes(i) = isv;
1197             PPar(i) = puv;
1198             
1199             if (orEdge==TopAbs_FORWARD)
1200               myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Frontier));
1201             else if (orEdge == TopAbs_REVERSED)
1202               myStructure->AddLink(BRepMesh_Edge(iv2, ivf, BRepMesh_Frontier));
1203             else if (orEdge == TopAbs_INTERNAL)
1204               myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Fixed));
1205     
1206             ivf = iv2;
1207           }
1208         }
1209         else
1210         {
1211           const Standard_Real wFold = Par->Value(Par->Lower());
1212           const Standard_Real wLold = Par->Value(Par->Upper());
1213
1214           Standard_Real wKoef = 1.;
1215           if ((wFold != wFirst || wLold != wLast) && wLold != wFold)
1216           {
1217             wKoef = (wLast - wFirst) / (wLold - wFold);
1218           }
1219           
1220           BRepAdaptor_Curve cons(theEdge, theFace);
1221           Extrema_LocateExtPC pcos;
1222           pcos.Initialize(cons, cons.FirstParameter(), 
1223               cons.LastParameter(), Precision::PConfusion());
1224
1225           Standard_Real wPrev;
1226           Standard_Real wCur      = wFirst;
1227           Standard_Real wCurFound = wFirst;
1228           for (i = 2; i < nbnodes; i++)
1229           {
1230             P3d = myLocation3d(NOD(i));
1231             // Record 2d point
1232             wPrev = wCur;
1233             wCur  = wFirst + wKoef*(Par->Value(i) - wFold);
1234                   wCurFound += (wCur - wPrev);
1235             pcos.Perform(P3d, wCurFound);
1236             if (pcos.IsDone()) wCurFound = pcos.Point().Parameter();
1237             theC2d->D0(wCurFound, uv);
1238             v2.Initialize(uv.Coord(), NOD(i), BRepMesh_OnCurve);
1239             iv2 = myStructure->AddNode(v2);
1240                   
1241             Standard_Integer isv = myVemap.FindIndex(iv2);
1242             if (isv == 0) isv = myVemap.Add(iv2); 
1243                   Nodes(i) = isv;
1244             PPar(i) = wCurFound;
1245                   
1246             if (orEdge==TopAbs_FORWARD)
1247               myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Frontier));
1248             else if (orEdge == TopAbs_REVERSED)
1249               myStructure->AddLink(BRepMesh_Edge(iv2, ivf, BRepMesh_Frontier));
1250             else if (orEdge == TopAbs_INTERNAL)
1251               myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Fixed));
1252     
1253             ivf = iv2;
1254           }
1255         }
1256       }
1257
1258     P1 = new Poly_PolygonOnTriangulation(Nodes, PPar);
1259     
1260     if (ivf != ivl) {
1261       if (orEdge == TopAbs_FORWARD) 
1262         myStructure->AddLink(BRepMesh_Edge(ivf, ivl, BRepMesh_Frontier));
1263       else if (orEdge == TopAbs_REVERSED)
1264         myStructure->AddLink(BRepMesh_Edge(ivl, ivf, BRepMesh_Frontier));
1265       else if (orEdge == TopAbs_INTERNAL)
1266         myStructure->AddLink(BRepMesh_Edge(ivf, ivl, BRepMesh_Fixed));
1267       }
1268     }
1269   }
1270
1271   P1->Deflection(theDefEdge);
1272   if (myInternaledges.IsBound(theEdge))
1273   {
1274     BRepMesh_PairOfPolygon& pair = myInternaledges.ChangeFind(theEdge);
1275     if (orEdge == TopAbs_REVERSED)
1276       pair.Append(P1);
1277     else
1278       pair.Prepend(P1);
1279   }
1280   else
1281   {
1282     BRepMesh_PairOfPolygon pair1;
1283     pair1.Append(P1);
1284     myInternaledges.Bind(theEdge, pair1);
1285   }
1286 }
1287
1288
1289 //=======================================================================
1290 //function : Update(edge)
1291 //purpose  :
1292 //=======================================================================
1293 Standard_Boolean BRepMesh_FastDiscret::Update(const TopoDS_Edge&          theEdge,
1294                                               const TopoDS_Face&          theFace,
1295                                               const Handle(Geom2d_Curve)& theC2d,
1296                                               const Standard_Real         theDefEdge,
1297                                               const Standard_Real         theFirst,
1298                                               const Standard_Real         theLast)
1299 {
1300   TopLoc_Location Loc;
1301   Handle(Poly_Triangulation) T, TNull;
1302   Handle(Poly_PolygonOnTriangulation) Poly, NullPoly;
1303
1304   Standard_Integer i = 1;
1305   Standard_Boolean found = Standard_False;
1306   do
1307   {
1308     BRep_Tool::PolygonOnTriangulation(theEdge,Poly,T,Loc,i);
1309     i++;
1310     if (!found && !T.IsNull() && T->HasUVNodes() && 
1311         !Poly.IsNull() && Poly->HasParameters())
1312     {
1313       if (Poly->Deflection() <= 1.1*theDefEdge)
1314       {
1315         // 2d vertex indices
1316         TopAbs_Orientation orEdge = theEdge.Orientation();
1317         Standard_Integer iv1, iv2, ivl;
1318         Standard_Integer isv1, isv, isvl;
1319         
1320         // Get range on 2d curve
1321         Standard_Real wFirst, wLast;
1322         BRep_Tool::Range(theEdge, theFace, wFirst, wLast);
1323         
1324         // Get end points on 2d curve
1325         gp_Pnt2d uvFirst, uvLast;
1326         BRep_Tool::UVPoints(theEdge, theFace, uvFirst, uvLast);
1327
1328         // Get vertices
1329         TopoDS_Vertex pBegin, pEnd;
1330         TopExp::Vertices(theEdge,pBegin,pEnd);
1331         
1332         const Standard_Boolean sameUV =
1333           uvFirst.IsEqual(uvLast, Precision::PConfusion());
1334         
1335         //Controle vertice tolerances
1336         BRepAdaptor_Surface  BS(theFace, Standard_False);
1337         Handle(BRepAdaptor_HSurface) gFace = new BRepAdaptor_HSurface(BS);
1338         
1339
1340         gp_Pnt pFirst = gFace->Value(uvFirst.X(), uvFirst.Y());
1341         gp_Pnt pLast  = gFace->Value(uvLast.X(), uvLast.Y());
1342         
1343         Standard_Real mindist = 10. * Max(pFirst.Distance(BRep_Tool::Pnt(pBegin)), 
1344                                           pLast.Distance(BRep_Tool::Pnt(pEnd)));
1345
1346         if (mindist < BRep_Tool::Tolerance(pBegin) ||
1347             mindist < BRep_Tool::Tolerance(pEnd) ) mindist = theDefEdge;
1348         
1349         if (sameUV)
1350         {
1351           // 1. is it really sameUV without being degenerated
1352           gp_Pnt2d uvF, uvL;
1353           theC2d->D0(theFirst, uvF);
1354           theC2d->D0(theLast, uvL);
1355           if (!uvFirst.IsEqual(uvF, Precision::PConfusion())) {
1356             uvFirst = uvF;
1357           }
1358           if (!uvLast.IsEqual(uvL, Precision::PConfusion())) {
1359             uvLast = uvL; 
1360           }
1361         }
1362
1363         const TColgp_Array1OfPnt&      Nodes   = T->Nodes();
1364         const TColStd_Array1OfInteger& Indices = Poly->Nodes();
1365         Handle(TColStd_HArray1OfReal)  Param   = Poly->Parameters();
1366         
1367         const Standard_Integer nbnodes = Indices.Length();
1368         TColStd_Array1OfInteger NewNodes(1, nbnodes);
1369         TColStd_Array1OfInteger NewNodInStruct(1, nbnodes);
1370         
1371         gp_Pnt P3d;
1372         gp_XY theUV;
1373         
1374         // Process first vertex
1375         Standard_Integer ipf;
1376         if (myVertices.IsBound(pBegin))
1377         {
1378           ipf = myVertices.Find(pBegin);
1379         }
1380         else
1381         {
1382           if (sameUV && myVertices.IsBound(pEnd))
1383           {
1384               ipf = myVertices.Find(pEnd);
1385           }
1386           else
1387           {
1388             P3d = Nodes(Indices(1));
1389             if (!Loc.IsIdentity())
1390               P3d.Transform(Loc.Transformation());
1391             myNbLocat++;
1392             myLocation3d.Bind(myNbLocat,P3d);
1393             ipf = myNbLocat;
1394           }
1395           myVertices.Bind(pBegin,ipf);
1396         }
1397         NewNodInStruct(1) = ipf;
1398         theUV = BRepMesh_FastDiscretFace::FindUV(pBegin, uvFirst, ipf, gFace, mindist, myLocation2d);
1399         BRepMesh_Vertex vf(theUV,ipf,BRepMesh_Frontier);
1400         iv1 = myStructure->AddNode(vf);
1401         isv1 = myVemap.FindIndex(iv1);
1402         if (isv1 == 0) isv1 = myVemap.Add(iv1);
1403         NewNodes(1) = isv1;
1404         
1405         // Process last vertex
1406         Standard_Integer ipl;
1407         if (pEnd.IsSame(pBegin))
1408         {
1409           ipl = ipf;
1410         }
1411         else
1412         {
1413           if (myVertices.IsBound(pEnd))
1414           {
1415             ipl = myVertices.Find(pEnd);
1416           }
1417           else
1418           {
1419             if (sameUV)
1420             {
1421               ipl = ipf;
1422               ivl = iv1;
1423               isv1 = isv1;
1424             }
1425             else
1426             {
1427               myNbLocat++;
1428               myLocation3d.Bind(myNbLocat,Nodes(Indices(nbnodes)).Transformed(Loc.Transformation()));
1429               ipl = myNbLocat;
1430             }
1431             myVertices.Bind(pEnd,ipl);
1432           }
1433         }
1434         NewNodInStruct(nbnodes) = ipl;
1435         theUV = BRepMesh_FastDiscretFace::FindUV(pEnd, uvLast, ipl, gFace, mindist, myLocation2d);
1436         BRepMesh_Vertex vl(theUV,ipl,BRepMesh_Frontier);
1437         
1438         ivl = myStructure->AddNode(vl);
1439         isvl = myVemap.FindIndex(ivl);
1440         if (isvl == 0) isvl = myVemap.Add(ivl);
1441         
1442         NewNodes(nbnodes) = isvl;
1443   
1444         gp_Pnt2d uv;
1445         BRepMesh_Vertex v;
1446   
1447         if (BRep_Tool::SameParameter(theEdge))
1448         {
1449           for (i = 2; i < Indices.Length(); i++)
1450           {
1451             // Record 3d point
1452             P3d = Nodes(Indices(i));
1453             if (!Loc.IsIdentity())
1454               P3d.Transform(Loc.Transformation());
1455             myNbLocat++;
1456             myLocation3d.Bind(myNbLocat, P3d);
1457             NewNodInStruct(i) = myNbLocat;
1458             // Record 2d point
1459             uv = theC2d->Value(Param->Value(i));
1460             v.Initialize(uv.Coord(), myNbLocat, BRepMesh_Frontier);
1461             iv2 = myStructure->AddNode(v);
1462             isv = myVemap.FindIndex(iv2);
1463             if (isv == 0) isv = myVemap.Add(iv2);
1464             NewNodes(i) = isv;
1465             
1466             //add links
1467             if (orEdge == TopAbs_FORWARD)
1468               myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Frontier));
1469             else if (orEdge == TopAbs_REVERSED)
1470               myStructure->AddLink(BRepMesh_Edge(iv2,iv1,BRepMesh_Frontier));
1471             else if (orEdge == TopAbs_INTERNAL)
1472               myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Fixed));
1473             iv1 = iv2;  
1474           }
1475           
1476           // last point
1477           if (iv1 != ivl) {
1478             if (orEdge == TopAbs_FORWARD)
1479               myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Frontier));
1480             else if (orEdge == TopAbs_REVERSED)
1481               myStructure->AddLink(BRepMesh_Edge(ivl,iv1,BRepMesh_Frontier));
1482             else if (orEdge == TopAbs_INTERNAL)
1483               myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Fixed));
1484           }
1485           
1486           
1487         }
1488         else
1489         {
1490           const Standard_Real wFold = Param->Value(Param->Lower());
1491           const Standard_Real wLold = Param->Value(Param->Upper());
1492           
1493           Standard_Real wKoef = 1.;
1494           if ((wFold != wFirst || wLold != wLast) && wLold != wFold)
1495           {
1496             wKoef = (wLast - wFirst) / (wLold - wFold);
1497           }
1498     
1499           BRepAdaptor_Curve cons(theEdge, theFace);
1500           Extrema_LocateExtPC pcos;
1501           pcos.Initialize(cons, cons.FirstParameter(), 
1502                           cons.LastParameter(), Precision::PConfusion());
1503           
1504           Standard_Real wPrev;
1505           Standard_Real wCur      = wFirst;
1506           Standard_Real wCurFound = wFirst;
1507           for (i = 2; i < Indices.Length(); i++)
1508           {
1509             // Record 3d point
1510             P3d = Nodes(Indices(i));
1511             if (!Loc.IsIdentity())
1512               P3d.Transform(Loc.Transformation());
1513             myNbLocat++;
1514             myLocation3d.Bind(myNbLocat, P3d);
1515             NewNodInStruct(i) = myNbLocat;
1516             // Record 2d point
1517             wPrev = wCur;
1518             wCur  = wFirst + wKoef*(Param->Value(i) - wFold);
1519             wCurFound += (wCur - wPrev);
1520             pcos.Perform(P3d, wCurFound);
1521             if (pcos.IsDone()) wCurFound = pcos.Point().Parameter();
1522             theC2d->D0(wCurFound, uv);
1523             v.Initialize(uv.Coord(), myNbLocat, BRepMesh_Frontier);
1524             iv2 = myStructure->AddNode(v);
1525             isv = myVemap.FindIndex(iv2);
1526             if (isv == 0) isv = myVemap.Add(iv2);
1527             NewNodes(i) = isv;
1528
1529             
1530             //add links
1531             if (orEdge == TopAbs_FORWARD)
1532               myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Frontier));
1533             else if (orEdge == TopAbs_REVERSED)
1534               myStructure->AddLink(BRepMesh_Edge(iv2,iv1,BRepMesh_Frontier));
1535             else if (orEdge == TopAbs_INTERNAL)
1536               myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Fixed));
1537             iv1 = iv2;              
1538           }
1539           
1540           // last point
1541           if (iv1 != ivl) {
1542             if (orEdge == TopAbs_FORWARD)
1543               myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Frontier));
1544             else if (orEdge == TopAbs_REVERSED)
1545               myStructure->AddLink(BRepMesh_Edge(ivl,iv1,BRepMesh_Frontier));
1546             else if (orEdge == TopAbs_INTERNAL)
1547               myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Fixed));
1548           }
1549         }
1550         
1551         Handle(Poly_PolygonOnTriangulation) P1 =
1552           new Poly_PolygonOnTriangulation(NewNodes, Param->Array1());
1553         P1->Deflection(theDefEdge);
1554         if (myInternaledges.IsBound(theEdge))
1555         {
1556           BRepMesh_PairOfPolygon& pair = myInternaledges.ChangeFind(theEdge);
1557           if (theEdge.Orientation() == TopAbs_REVERSED)
1558             pair.Append(P1);
1559           else
1560             pair.Prepend(P1);
1561         }
1562         else
1563         {
1564           BRepMesh_PairOfPolygon pair1;
1565           pair1.Append(P1);
1566           myInternaledges.Bind(theEdge, pair1);
1567         }
1568         
1569         Handle(Poly_PolygonOnTriangulation) P2 =
1570           new Poly_PolygonOnTriangulation(NewNodInStruct, Param->Array1());
1571         P2->Deflection(theDefEdge);
1572         BRepMesh_PairOfPolygon pair;
1573         pair.Append(P2);
1574         myEdges.Bind(theEdge, pair);
1575         
1576         found = Standard_True;
1577       }
1578       else
1579       {
1580         BRep_Builder B;
1581         B.UpdateEdge(theEdge,NullPoly,T,Loc);
1582         B.UpdateFace(theFace,TNull);
1583       }
1584     }
1585     else if (!T.IsNull() && !T->HasUVNodes())
1586     {
1587       BRep_Builder B;
1588       B.UpdateEdge(theEdge,NullPoly,T,Loc);
1589       B.UpdateFace(theFace,TNull);
1590     }
1591   }
1592   while (!Poly.IsNull());
1593
1594   return found;
1595 }
1596
1597 //=======================================================================
1598 //function : output
1599 //purpose  : 
1600 //=======================================================================
1601 Standard_Integer BRepMesh_FastDiscret::NbTriangles() const
1602 {
1603   return myStructure->NbElements();
1604 }
1605
1606 //=======================================================================
1607 //function : Triangle
1608 //purpose  : 
1609 //=======================================================================
1610
1611 const BRepMesh_Triangle& BRepMesh_FastDiscret::Triangle
1612   (const Standard_Integer Index) const
1613 {
1614   return myStructure->GetElement(Index);
1615 }
1616
1617 //=======================================================================
1618 //function : NbEdges
1619 //purpose  : 
1620 //=======================================================================
1621
1622 Standard_Integer BRepMesh_FastDiscret::NbEdges() const
1623 {
1624   return myStructure->NbLinks();
1625 }
1626
1627 //=======================================================================
1628 //function : Edge
1629 //purpose  : 
1630 //=======================================================================
1631
1632 const BRepMesh_Edge& BRepMesh_FastDiscret::Edge(const Standard_Integer Index) const
1633 {
1634   return myStructure->GetLink(Index);
1635 }
1636
1637 //=======================================================================
1638 //function : NbVertices
1639 //purpose  : 
1640 //=======================================================================
1641
1642 Standard_Integer BRepMesh_FastDiscret::NbVertices() const
1643 {
1644   return myStructure->NbNodes();
1645 }
1646
1647 //=======================================================================
1648 //function : Vertex
1649 //purpose  : 
1650 //=======================================================================
1651
1652 const BRepMesh_Vertex& BRepMesh_FastDiscret::Vertex
1653   (const Standard_Integer Index) const
1654 {
1655   return myStructure->GetNode(Index);
1656 }
1657
1658 //=======================================================================
1659 //function : Pnt
1660 //purpose  : 
1661 //=======================================================================
1662
1663 const gp_Pnt& BRepMesh_FastDiscret::Pnt(const Standard_Integer Index) const
1664 {
1665   return myLocation3d(myStructure->GetNode(Index).Location3d());
1666 }
1667
1668 //=======================================================================
1669 //function : VerticesOfDomain
1670 //purpose  : 
1671 //=======================================================================
1672
1673 void BRepMesh_FastDiscret::VerticesOfDomain(BRepMesh_MapOfInteger&  Indices) const 
1674
1675   Indices.Clear();
1676   
1677   // recuperate from the map of edges.
1678   const BRepMesh_MapOfInteger& edmap = myStructure->LinkOfDomain();
1679
1680   // iterator on edges.
1681   BRepMesh_MapOfInteger::Iterator iter(edmap);
1682   
1683   Standard_Integer ind_edge;
1684   for (iter.Reset(); iter.More(); iter.Next()) {
1685     ind_edge = iter.Key();
1686     const BRepMesh_Edge& Ed = Edge(ind_edge);
1687     Indices.Add(Ed.FirstNode());
1688     Indices.Add(Ed.LastNode());
1689   }
1690 }
1691
1692 BRepMesh_Status BRepMesh_FastDiscret::CurrentFaceStatus() const
1693 {
1694   return myFacestate;
1695 }
1696
1697 //=======================================================================
1698 //function : GetFaceAttribute
1699 //purpose  : 
1700 //=======================================================================
1701
1702 Standard_Boolean BRepMesh_FastDiscret::GetFaceAttribute(const TopoDS_Face&              theFace,
1703                                                         Handle(BRepMesh_FaceAttribute)& theFattrib) const
1704 {
1705   if(myMapattrib.IsBound(theFace))
1706   {
1707     theFattrib = myMapattrib(theFace);
1708     return Standard_True;
1709   }
1710   return Standard_False;
1711 }
1712
1713 //=======================================================================
1714 //function : RemoveFaceAttribute
1715 //purpose  : 
1716 //=======================================================================
1717
1718 void BRepMesh_FastDiscret::RemoveFaceAttribute(const TopoDS_Face& theFace)
1719 {
1720   if(myMapattrib.IsBound(theFace))
1721     myMapattrib.UnBind(theFace);
1722 }
1723
1724 //=======================================================================
1725 //function : CreateMutexesForSubShapes
1726 //purpose  : 
1727 //=======================================================================
1728 void BRepMesh_FastDiscret::CreateMutexesForSubShapes(const TopoDS_Shape& theShape,
1729                                                      const TopAbs_ShapeEnum theType)
1730 {
1731   myMutexProvider.CreateMutexesForSubShapes(theShape, theType);
1732 }
1733
1734 //=======================================================================
1735 //function : RemoveAllMutexes
1736 //purpose  : 
1737 //=======================================================================
1738 void BRepMesh_FastDiscret::RemoveAllMutexes()
1739 {
1740   myMutexProvider.RemoveAllMutexes();
1741 }