0031031: Incorrect result is returned from BRepPrimAPI_MakePrism::Generated()
[occt.git] / src / BRepLib / BRepLib.cxx
1 // Created on: 1993-12-15
2 // Created by: Remi LEQUETTE
3 // Copyright (c) 1993-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 //pmn 26/09/97 Add parameters of approximation in BuildCurve3d
18 //  Modified by skv - Thu Jun  3 12:39:19 2004 OCC5898
19
20 #include <Adaptor3d_CurveOnSurface.hxx>
21 #include <AdvApprox_ApproxAFunction.hxx>
22 #include <AppParCurves_MultiBSpCurve.hxx>
23 #include <AppParCurves_MultiCurve.hxx>
24 #include <Approx_CurvilinearParameter.hxx>
25 #include <Approx_SameParameter.hxx>
26 #include <Bnd_Box.hxx>
27 #include <BRep_Builder.hxx>
28 #include <BRep_CurveRepresentation.hxx>
29 #include <BRep_GCurve.hxx>
30 #include <BRep_ListIteratorOfListOfCurveRepresentation.hxx>
31 #include <BRep_ListOfCurveRepresentation.hxx>
32 #include <BRepCheck.hxx>
33 #include <BRep_TEdge.hxx>
34 #include <BRep_TFace.hxx>
35 #include <BRep_Tool.hxx>
36 #include <BRep_TVertex.hxx>
37 #include <BRepAdaptor_HCurve.hxx>
38 #include <BRepAdaptor_HCurve2d.hxx>
39 #include <BRepAdaptor_HSurface.hxx>
40 #include <BRepAdaptor_Surface.hxx>
41 #include <BRepBndLib.hxx>
42 #include <BRepClass3d_SolidClassifier.hxx>
43 #include <BRepLib.hxx>
44 #include <BRepLib_MakeFace.hxx>
45 #include <BSplCLib.hxx>
46 #include <ElSLib.hxx>
47 #include <Extrema_LocateExtPC.hxx>
48 #include <GCPnts_QuasiUniformDeflection.hxx>
49 #include <Geom2d_BSplineCurve.hxx>
50 #include <Geom2d_Curve.hxx>
51 #include <Geom2d_TrimmedCurve.hxx>
52 #include <Geom2dAdaptor.hxx>
53 #include <Geom2dAdaptor_Curve.hxx>
54 #include <Geom2dAdaptor_HCurve.hxx>
55 #include <Geom2dConvert.hxx>
56 #include <Geom_BSplineCurve.hxx>
57 #include <Geom_BSplineSurface.hxx>
58 #include <Geom_Curve.hxx>
59 #include <Geom_Plane.hxx>
60 #include <Geom_RectangularTrimmedSurface.hxx>
61 #include <Geom_Surface.hxx>
62 #include <Geom_TrimmedCurve.hxx>
63 #include <GeomAdaptor_Curve.hxx>
64 #include <GeomAdaptor_HCurve.hxx>
65 #include <GeomAdaptor_HSurface.hxx>
66 #include <GeomAdaptor_Surface.hxx>
67 #include <GeomLib.hxx>
68 #include <GeomLProp_SLProps.hxx>
69 #include <gp.hxx>
70 #include <gp_Ax2.hxx>
71 #include <gp_Pln.hxx>
72 #include <Poly_PolygonOnTriangulation.hxx>
73 #include <Poly_Triangulation.hxx>
74 #include <Precision.hxx>
75 #include <ProjLib_ProjectedCurve.hxx>
76 #include <Standard_ErrorHandler.hxx>
77 #include <Standard_Real.hxx>
78 #include <TColgp_Array1OfPnt.hxx>
79 #include <TColgp_Array1OfPnt2d.hxx>
80 #include <TColStd_Array1OfReal.hxx>
81 #include <TColStd_MapOfTransient.hxx>
82 #include <TopExp.hxx>
83 #include <TopExp_Explorer.hxx>
84 #include <TopoDS.hxx>
85 #include <TopoDS_Edge.hxx>
86 #include <TopoDS_Face.hxx>
87 #include <TopoDS_Shape.hxx>
88 #include <TopoDS_Solid.hxx>
89 #include <TopoDS_Vertex.hxx>
90 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
91 #include <TopTools_ListIteratorOfListOfShape.hxx>
92 #include <TopTools_MapOfShape.hxx>
93 #include <TShort_HArray1OfShortReal.hxx>
94 #include <TColgp_Array1OfXY.hxx>
95 #include <BRepTools_ReShape.hxx>
96 #include <TopTools_DataMapOfShapeReal.hxx>
97 #include <TopoDS_LockedShape.hxx>
98
99 #include <algorithm>
100
101 // TODO - not thread-safe static variables
102 static Standard_Real thePrecision = Precision::Confusion();     
103 static Handle(Geom_Plane) thePlane;
104
105 static void InternalUpdateTolerances(const TopoDS_Shape& theOldShape,
106   const Standard_Boolean IsVerifyTolerance, const Standard_Boolean IsMutableInput, BRepTools_ReShape& theReshaper);
107
108 //=======================================================================
109 // function: BRepLib_ComparePoints
110 // purpose:  implementation of IsLess() function for two points
111 //=======================================================================
112 struct BRepLib_ComparePoints {
113   bool operator()(const gp_Pnt& theP1, const gp_Pnt& theP2)
114   {
115     for (Standard_Integer i = 1; i <= 3; ++i) {
116       if (theP1.Coord(i) < theP2.Coord(i)) {
117         return Standard_True;
118       }
119       else if (theP1.Coord(i) > theP2.Coord(i)) {
120         return Standard_False;
121       }
122     }
123     return Standard_False;
124   }
125 };
126
127
128 //=======================================================================
129 //function : Precision
130 //purpose  : 
131 //=======================================================================
132
133 void BRepLib::Precision(const Standard_Real P)
134 {
135   thePrecision = P;
136 }
137
138 //=======================================================================
139 //function : Precision
140 //purpose  : 
141 //=======================================================================
142
143 Standard_Real  BRepLib::Precision()
144 {
145   return thePrecision;
146 }
147
148 //=======================================================================
149 //function : Plane
150 //purpose  : 
151 //=======================================================================
152
153 void  BRepLib::Plane(const Handle(Geom_Plane)& P)
154 {
155   thePlane = P;
156 }
157
158
159 //=======================================================================
160 //function : Plane
161 //purpose  : 
162 //=======================================================================
163
164 const Handle(Geom_Plane)&  BRepLib::Plane()
165 {
166   if (thePlane.IsNull()) thePlane = new Geom_Plane(gp::XOY());
167   return thePlane;
168 }
169 //=======================================================================
170 //function : CheckSameRange
171 //purpose  : 
172 //=======================================================================
173
174 Standard_Boolean  BRepLib::CheckSameRange(const TopoDS_Edge& AnEdge,
175   const Standard_Real Tolerance) 
176 {
177   Standard_Boolean  IsSameRange = Standard_True,
178     first_time_in = Standard_True ;
179
180   BRep_ListIteratorOfListOfCurveRepresentation an_Iterator
181     ((*((Handle(BRep_TEdge)*)&AnEdge.TShape()))->ChangeCurves());
182
183   Standard_Real first, last;
184   Standard_Real current_first =0., current_last =0. ;
185   Handle(BRep_GCurve) geometric_representation_ptr ;
186
187   while (IsSameRange && an_Iterator.More()) {
188     geometric_representation_ptr =
189       Handle(BRep_GCurve)::DownCast(an_Iterator.Value());
190     if (!geometric_representation_ptr.IsNull()) {
191
192       first = geometric_representation_ptr->First();
193       last =  geometric_representation_ptr->Last();
194       if (first_time_in ) {
195         current_first = first ;
196         current_last = last   ;
197         first_time_in = Standard_False ;
198       }
199       else {
200         IsSameRange = (Abs(current_first - first) <= Tolerance) 
201           && (Abs(current_last -last) <= Tolerance ) ;
202       }
203     }
204     an_Iterator.Next() ;
205   }
206   return IsSameRange ;
207 }
208
209 //=======================================================================
210 //function : SameRange
211 //purpose  : 
212 //=======================================================================
213
214 void BRepLib::SameRange(const TopoDS_Edge& AnEdge,
215   const Standard_Real Tolerance) 
216 {
217   BRep_ListIteratorOfListOfCurveRepresentation an_Iterator
218     ((*((Handle(BRep_TEdge)*)&AnEdge.TShape()))->ChangeCurves());
219
220   Handle(Geom2d_Curve) Curve2dPtr, Curve2dPtr2, NewCurve2dPtr, NewCurve2dPtr2;
221   TopLoc_Location LocalLoc ;
222
223   Standard_Boolean first_time_in = Standard_True,
224     has_curve,
225     has_closed_curve ;
226   Handle(BRep_GCurve) geometric_representation_ptr ;
227   Standard_Real first,
228     current_first,
229     last,
230     current_last ;
231
232   const Handle(Geom_Curve) C = BRep_Tool::Curve(AnEdge,
233     LocalLoc,
234     current_first,
235     current_last);
236   if (!C.IsNull()) {
237     first_time_in = Standard_False ;
238   }
239
240   while (an_Iterator.More()) {
241     geometric_representation_ptr =
242       Handle(BRep_GCurve)::DownCast(an_Iterator.Value());
243     if (! geometric_representation_ptr.IsNull()) {
244       has_closed_curve =
245         has_curve = Standard_False ;
246       first = geometric_representation_ptr->First();
247       last =  geometric_representation_ptr->Last();
248       if (geometric_representation_ptr->IsCurveOnSurface()) {
249         Curve2dPtr = geometric_representation_ptr->PCurve() ; 
250         has_curve = Standard_True ;
251       }
252       if (geometric_representation_ptr->IsCurveOnClosedSurface()) {
253         Curve2dPtr2 = geometric_representation_ptr->PCurve2() ;
254         has_closed_curve = Standard_True ;
255       }
256       if (has_curve || has_closed_curve) {
257         if (first_time_in) {
258           current_first = first ;
259           current_last = last ;
260           first_time_in = Standard_False ;
261         }
262
263         if (Abs(first - current_first) > Precision::Confusion() ||
264           Abs(last - current_last) > Precision::Confusion() )
265         {
266           if (has_curve)
267           {
268             GeomLib::SameRange(Tolerance,
269               Curve2dPtr,
270               geometric_representation_ptr->First(),
271               geometric_representation_ptr->Last(),
272               current_first,
273               current_last,
274               NewCurve2dPtr);
275             geometric_representation_ptr->PCurve(NewCurve2dPtr) ;
276           }
277           if (has_closed_curve)
278           {
279             GeomLib::SameRange(Tolerance,
280               Curve2dPtr2,
281               geometric_representation_ptr->First(),
282               geometric_representation_ptr->Last(),
283               current_first,
284               current_last,
285               NewCurve2dPtr2);
286             geometric_representation_ptr->PCurve2(NewCurve2dPtr2) ;
287           }
288         }
289       }
290     }
291     an_Iterator.Next() ;
292   }
293   BRep_Builder B;
294   B.Range(TopoDS::Edge(AnEdge),
295     current_first,
296     current_last) ;
297
298   B.SameRange(AnEdge,
299     Standard_True) ;
300 }
301
302 //=======================================================================
303 //function : EvaluateMaxSegment
304 //purpose  : return MaxSegment to pass in approximation, if MaxSegment==0 provided
305 //=======================================================================
306
307 static Standard_Integer evaluateMaxSegment(const Standard_Integer aMaxSegment,
308   const Adaptor3d_CurveOnSurface& aCurveOnSurface)
309 {
310   if (aMaxSegment != 0) return aMaxSegment;
311
312   Handle(Adaptor3d_HSurface) aSurf   = aCurveOnSurface.GetSurface();
313   Handle(Adaptor2d_HCurve2d) aCurv2d = aCurveOnSurface.GetCurve();
314
315   Standard_Real aNbSKnots = 0, aNbC2dKnots = 0;
316
317   if (aSurf->GetType() == GeomAbs_BSplineSurface) {
318     Handle(Geom_BSplineSurface) aBSpline = aSurf->BSpline();
319     aNbSKnots = Max(aBSpline->NbUKnots(), aBSpline->NbVKnots());
320   }
321   if (aCurv2d->GetType() == GeomAbs_BSplineCurve) {
322     aNbC2dKnots = aCurv2d->NbKnots();
323   }
324   Standard_Integer aReturn = (Standard_Integer) (  30 + Max(aNbSKnots, aNbC2dKnots) ) ;
325   return aReturn;
326 }
327
328 //=======================================================================
329 //function : BuildCurve3d
330 //purpose  : 
331 //=======================================================================
332
333 Standard_Boolean  BRepLib::BuildCurve3d(const TopoDS_Edge& AnEdge,
334   const Standard_Real Tolerance,
335   const GeomAbs_Shape Continuity,
336   const Standard_Integer MaxDegree,
337   const Standard_Integer MaxSegment)
338 {
339   Standard_Integer //ErrorCode,
340     //                   ReturnCode = 0,
341     ii,
342     //                   num_knots,
343     jj;
344
345   TopLoc_Location LocalLoc,L[2],LC;
346   Standard_Real f,l,fc,lc, first[2], last[2],
347     tolerance,
348     max_deviation,
349     average_deviation ;
350   Handle(Geom2d_Curve) Curve2dPtr, Curve2dArray[2]  ;
351   Handle(Geom_Surface) SurfacePtr, SurfaceArray[2]  ;
352
353   Standard_Integer not_done ;
354   // if the edge has a 3d curve returns true
355
356
357   const Handle(Geom_Curve) C = BRep_Tool::Curve(AnEdge,LocalLoc,f,l);
358   if (!C.IsNull()) 
359     return Standard_True;
360   //
361   // this should not exists but UpdateEdge makes funny things 
362   // if the edge is not same range 
363   //
364   if (! CheckSameRange(AnEdge,
365     Precision::Confusion())) {
366       SameRange(AnEdge,
367         Tolerance) ;
368   }
369
370
371
372   // search a curve on a plane
373   Handle(Geom_Surface) S;
374   Handle(Geom2d_Curve) PC;
375   Standard_Integer i = 0;
376   Handle(Geom_Plane) P;
377   not_done = 1 ;
378
379   while (not_done) {
380     i++;
381     BRep_Tool::CurveOnSurface(AnEdge,PC,S,LocalLoc,f,l,i);
382     Handle(Geom_RectangularTrimmedSurface) RT = 
383       Handle(Geom_RectangularTrimmedSurface)::DownCast(S);
384     if ( RT.IsNull()) {
385       P = Handle(Geom_Plane)::DownCast(S);
386     }
387     else {
388       P = Handle(Geom_Plane)::DownCast(RT->BasisSurface());
389     }
390     not_done = ! (S.IsNull() || !P.IsNull()) ;
391   }
392   if (! P.IsNull()) {
393     // compute the 3d curve
394     gp_Ax2 axes = P->Position().Ax2();
395     Handle(Geom_Curve) C3d = GeomLib::To3d(axes,PC);
396     if (C3d.IsNull())
397       return Standard_False;
398     // update the edge
399     Standard_Real First, Last;
400
401     BRep_Builder B;
402     B.UpdateEdge(AnEdge,C3d,LocalLoc,0.0e0);
403     BRep_Tool::Range(AnEdge, S, LC, First, Last);
404     B.Range(AnEdge, First, Last); //Do not forget 3D range.(PRO6412)
405   }
406   else {
407     //
408     // compute the 3d curve using existing surface
409     //
410     fc = f ;
411     lc = l ;
412     if (!BRep_Tool::Degenerated(AnEdge)) {
413       jj = 0 ;
414       for (ii = 0 ; ii < 3 ; ii++ ) {
415         BRep_Tool::CurveOnSurface(TopoDS::Edge(AnEdge),
416           Curve2dPtr,
417           SurfacePtr,
418           LocalLoc,
419           fc,
420           lc,
421           ii) ;
422
423         if (!Curve2dPtr.IsNull() && jj < 2){
424           Curve2dArray[jj] = Curve2dPtr ;
425           SurfaceArray[jj] = SurfacePtr ;
426           L[jj] = LocalLoc ;
427           first[jj] = fc ;
428           last[jj] = lc ;
429           jj += 1 ;
430         }
431       }
432       f = first[0] ;
433       l = last[0] ;
434       Curve2dPtr = Curve2dArray[0] ;
435       SurfacePtr = SurfaceArray[0] ;
436
437       Geom2dAdaptor_Curve     AnAdaptor3dCurve2d (Curve2dPtr, f, l) ;
438       GeomAdaptor_Surface     AnAdaptor3dSurface (SurfacePtr) ;
439       Handle(Geom2dAdaptor_HCurve) AnAdaptor3dCurve2dPtr =
440         new Geom2dAdaptor_HCurve(AnAdaptor3dCurve2d) ;
441       Handle(GeomAdaptor_HSurface) AnAdaptor3dSurfacePtr =
442         new GeomAdaptor_HSurface (AnAdaptor3dSurface) ;
443       Adaptor3d_CurveOnSurface  CurveOnSurface( AnAdaptor3dCurve2dPtr,
444         AnAdaptor3dSurfacePtr) ;
445
446       Handle(Geom_Curve) NewCurvePtr ;
447
448       GeomLib::BuildCurve3d(Tolerance,
449         CurveOnSurface,
450         f,
451         l,
452         NewCurvePtr,
453         max_deviation,
454         average_deviation,
455         Continuity,
456         MaxDegree,
457         evaluateMaxSegment(MaxSegment,CurveOnSurface)) ;
458       BRep_Builder B;   
459       tolerance = BRep_Tool::Tolerance(AnEdge) ;
460       //Patch
461       //max_deviation = Max(tolerance, max_deviation) ;
462       max_deviation = Max( tolerance, Tolerance );
463       if (NewCurvePtr.IsNull())
464         return Standard_False;
465       B.UpdateEdge(TopoDS::Edge(AnEdge),
466         NewCurvePtr,
467         L[0],
468         max_deviation) ;
469       if (jj == 1 ) {
470         //
471         // if there is only one curve on surface attached to the edge
472         // than it can be qualified sameparameter
473         //
474         B.SameParameter(TopoDS::Edge(AnEdge),
475           Standard_True) ;
476       }
477     }
478     else {
479       return Standard_False ;
480     }
481
482   }         
483   return Standard_True;
484 }
485 //=======================================================================
486 //function : BuildCurves3d
487 //purpose  : 
488 //=======================================================================
489
490 Standard_Boolean  BRepLib::BuildCurves3d(const TopoDS_Shape& S) 
491
492 {
493   return BRepLib::BuildCurves3d(S,
494     1.0e-5) ;
495 }
496
497 //=======================================================================
498 //function : BuildCurves3d
499 //purpose  : 
500 //=======================================================================
501
502 Standard_Boolean  BRepLib::BuildCurves3d(const TopoDS_Shape& S,
503   const Standard_Real Tolerance,
504   const GeomAbs_Shape Continuity,
505   const Standard_Integer MaxDegree,
506   const Standard_Integer MaxSegment)
507 {
508   Standard_Boolean boolean_value,
509     ok = Standard_True;
510   TopTools_MapOfShape a_counter ;
511   TopExp_Explorer ex(S,TopAbs_EDGE);
512
513   while (ex.More()) {
514     if (a_counter.Add(ex.Current())) {
515       boolean_value = 
516         BuildCurve3d(TopoDS::Edge(ex.Current()),
517         Tolerance, Continuity,
518         MaxDegree, MaxSegment);
519       ok = ok && boolean_value ;
520     }
521     ex.Next();
522   }
523   return ok;
524 }
525 //=======================================================================
526 //function : UpdateEdgeTolerance
527 //purpose  : 
528 //=======================================================================
529
530 Standard_Boolean  BRepLib::UpdateEdgeTol(const TopoDS_Edge& AnEdge,
531   const Standard_Real MinToleranceRequested,
532   const Standard_Real MaxToleranceToCheck)
533 {
534
535   Standard_Integer curve_on_surface_index,
536     curve_index,
537     not_done,
538     has_closed_curve,
539     has_curve,
540     jj,
541     ii,
542     geom_reference_curve_flag = 0,
543     max_sampling_points = 90,
544     min_sampling_points = 30 ;
545
546   Standard_Real factor = 100.0e0,
547     //     sampling_array[2],
548     safe_factor = 1.4e0,
549     current_last,
550     current_first,
551     max_distance,
552     coded_edge_tolerance,
553     edge_tolerance = 0.0e0 ;
554   Handle(TColStd_HArray1OfReal) parameters_ptr ;
555   Handle(BRep_GCurve) geometric_representation_ptr ;
556
557   if (BRep_Tool::Degenerated(AnEdge)) return Standard_False ;
558   coded_edge_tolerance = BRep_Tool::Tolerance(AnEdge) ;
559   if (coded_edge_tolerance > MaxToleranceToCheck) return Standard_False ;
560
561   const Handle(BRep_TEdge)& TE = *((Handle(BRep_TEdge)*)&AnEdge.TShape());
562   BRep_ListOfCurveRepresentation& list_curve_rep = TE->ChangeCurves() ;
563   BRep_ListIteratorOfListOfCurveRepresentation an_iterator(list_curve_rep),
564     second_iterator(list_curve_rep) ;
565   Handle(Geom2d_Curve) curve2d_ptr, new_curve2d_ptr;
566   Handle(Geom_Surface) surface_ptr ;
567   TopLoc_Location local_location ;
568   GCPnts_QuasiUniformDeflection  a_sampler ;
569   GeomAdaptor_Curve  geom_reference_curve ;
570   Adaptor3d_CurveOnSurface  curve_on_surface_reference ; 
571   Handle(Geom_Curve) C = BRep_Tool::Curve(AnEdge,
572     local_location,
573     current_first,
574     current_last);
575   curve_on_surface_index = -1 ;
576   if (!C.IsNull()) {
577     if (! local_location.IsIdentity()) {
578       C = Handle(Geom_Curve)::
579         DownCast(C-> Transformed(local_location.Transformation()) ) ;
580     }
581     geom_reference_curve.Load(C) ;
582     geom_reference_curve_flag = 1 ;
583     a_sampler.Initialize(geom_reference_curve,
584       MinToleranceRequested * factor,
585       current_first,
586       current_last) ;
587   }
588   else {
589     not_done = 1 ;
590     curve_on_surface_index = 0 ;  
591
592     while (not_done && an_iterator.More()) {
593       geometric_representation_ptr =
594         Handle(BRep_GCurve)::DownCast(second_iterator.Value());
595       if (!geometric_representation_ptr.IsNull() 
596         && geometric_representation_ptr->IsCurveOnSurface()) {
597           curve2d_ptr = geometric_representation_ptr->PCurve() ;
598           local_location = geometric_representation_ptr->Location() ;
599           current_first = geometric_representation_ptr->First();
600           //first = geometric_representation_ptr->First();
601           current_last =  geometric_representation_ptr->Last();
602           // must be inverted 
603           //
604           if (! local_location.IsIdentity() ) {
605             surface_ptr = Handle(Geom_Surface)::
606               DownCast( geometric_representation_ptr->Surface()->
607               Transformed(local_location.Transformation()) ) ;
608           }
609           else {
610             surface_ptr = 
611               geometric_representation_ptr->Surface() ;
612           }
613           not_done = 0 ;
614       }
615       curve_on_surface_index += 1 ;
616     }
617     Geom2dAdaptor_Curve     AnAdaptor3dCurve2d (curve2d_ptr) ;
618     GeomAdaptor_Surface     AnAdaptor3dSurface (surface_ptr) ;
619     Handle(Geom2dAdaptor_HCurve) AnAdaptor3dCurve2dPtr =
620       new Geom2dAdaptor_HCurve(AnAdaptor3dCurve2d) ;
621     Handle(GeomAdaptor_HSurface) AnAdaptor3dSurfacePtr =
622       new GeomAdaptor_HSurface (AnAdaptor3dSurface) ;
623     curve_on_surface_reference.Load (AnAdaptor3dCurve2dPtr, AnAdaptor3dSurfacePtr);
624     a_sampler.Initialize(curve_on_surface_reference,
625       MinToleranceRequested * factor,
626       current_first,
627       current_last) ;
628   }
629   TColStd_Array1OfReal   sampling_parameters(1,a_sampler.NbPoints()) ;
630   for (ii = 1 ; ii <= a_sampler.NbPoints() ; ii++) {
631     sampling_parameters(ii) = a_sampler.Parameter(ii) ;
632   }
633   if (a_sampler.NbPoints() < min_sampling_points) {
634     GeomLib::DensifyArray1OfReal(min_sampling_points,
635       sampling_parameters,
636       parameters_ptr) ;
637   }
638   else if (a_sampler.NbPoints() > max_sampling_points) {
639     GeomLib::RemovePointsFromArray(max_sampling_points,
640       sampling_parameters,
641       parameters_ptr) ; 
642   }
643   else {
644     jj = 1 ;
645     parameters_ptr =
646       new TColStd_HArray1OfReal(1,sampling_parameters.Length()) ;
647     for (ii = sampling_parameters.Lower() ; ii <= sampling_parameters.Upper() ; ii++) {
648       parameters_ptr->ChangeArray1()(jj) =
649         sampling_parameters(ii) ;
650       jj +=1 ;
651     }
652   }
653
654   curve_index = 0 ;
655
656   while (second_iterator.More()) {
657     geometric_representation_ptr =
658       Handle(BRep_GCurve)::DownCast(second_iterator.Value());
659     if (! geometric_representation_ptr.IsNull() && 
660       curve_index != curve_on_surface_index) {
661         has_closed_curve =
662           has_curve = Standard_False ;
663         //      first = geometric_representation_ptr->First();
664         //      last =  geometric_representation_ptr->Last();
665         local_location = geometric_representation_ptr->Location() ;
666         if (geometric_representation_ptr->IsCurveOnSurface()) {
667           curve2d_ptr = geometric_representation_ptr->PCurve() ; 
668           has_curve = Standard_True ;
669         }
670         if (geometric_representation_ptr->IsCurveOnClosedSurface()) {
671           curve2d_ptr = geometric_representation_ptr->PCurve2() ;
672           has_closed_curve = Standard_True ;
673         }
674
675         if (has_curve ||
676           has_closed_curve) {
677             if (! local_location.IsIdentity() ) {
678               surface_ptr = Handle(Geom_Surface)::
679                 DownCast( geometric_representation_ptr->Surface()->
680                 Transformed(local_location.Transformation()) ) ;
681             }
682             else {
683               surface_ptr = 
684                 geometric_representation_ptr->Surface() ;
685             }
686             Geom2dAdaptor_Curve     an_adaptor_curve2d (curve2d_ptr) ;
687             GeomAdaptor_Surface     an_adaptor_surface(surface_ptr) ;
688             Handle(Geom2dAdaptor_HCurve) an_adaptor_curve2d_ptr =
689               new Geom2dAdaptor_HCurve(an_adaptor_curve2d) ;
690             Handle(GeomAdaptor_HSurface) an_adaptor_surface_ptr =
691               new GeomAdaptor_HSurface (an_adaptor_surface) ;
692             Adaptor3d_CurveOnSurface a_curve_on_surface(an_adaptor_curve2d_ptr,
693               an_adaptor_surface_ptr) ;
694
695             if (BRep_Tool::SameParameter(AnEdge)) {
696
697               GeomLib::EvalMaxParametricDistance(a_curve_on_surface,
698                 geom_reference_curve,
699                 MinToleranceRequested,
700                 parameters_ptr->Array1(),
701                 max_distance) ;
702             }
703             else if (geom_reference_curve_flag) {
704               GeomLib::EvalMaxDistanceAlongParameter(a_curve_on_surface,
705                 geom_reference_curve,
706                 MinToleranceRequested,
707                 parameters_ptr->Array1(),
708                 max_distance) ;
709             }
710             else {
711
712               GeomLib::EvalMaxDistanceAlongParameter(a_curve_on_surface,
713                 curve_on_surface_reference,
714                 MinToleranceRequested,
715                 parameters_ptr->Array1(),
716                 max_distance) ;
717             }
718             max_distance *= safe_factor ;
719             edge_tolerance = Max(max_distance, edge_tolerance) ;
720         }
721
722
723     }
724     curve_index += 1 ;
725     second_iterator.Next() ; 
726   }
727
728   TE->Tolerance(edge_tolerance);
729   return Standard_True ;
730
731 }
732 //=======================================================================
733 //function : UpdateEdgeTolerance
734 //purpose  : 
735 //=======================================================================
736 Standard_Boolean BRepLib::UpdateEdgeTolerance(const TopoDS_Shape& S,
737   const Standard_Real MinToleranceRequested,
738   const Standard_Real MaxToleranceToCheck) 
739 {
740   TopExp_Explorer ex(S,TopAbs_EDGE);
741   TopTools_MapOfShape  a_counter ;
742
743   Standard_Boolean     return_status = Standard_False,
744     local_flag ;
745
746   while (ex.More()) {
747     if (a_counter.Add(ex.Current())) {
748       local_flag =
749         BRepLib::UpdateEdgeTol(TopoDS::Edge(ex.Current()),
750         MinToleranceRequested,
751         MaxToleranceToCheck) ;
752       if (local_flag && ! return_status) {
753         return_status = Standard_True ;
754       }
755     }
756     ex.Next();
757   }
758   return return_status ;
759 }
760
761 //=======================================================================
762 //function : GetEdgeTol
763 //purpose  : 
764 //=======================================================================
765 static void GetEdgeTol(const TopoDS_Edge& theEdge,
766   const TopoDS_Face& theFace, Standard_Real& theEdTol)
767 {
768   TopLoc_Location L;
769   const Handle(Geom_Surface)& S = BRep_Tool::Surface(theFace,L);
770   TopLoc_Location l = L.Predivided(theEdge.Location());
771
772   const Handle(BRep_TEdge)& TE = *((Handle(BRep_TEdge)*)&theEdge.TShape());
773   BRep_ListIteratorOfListOfCurveRepresentation itcr(TE->Curves());
774
775   while (itcr.More()) {
776     const Handle(BRep_CurveRepresentation)& cr = itcr.Value();
777     if(cr->IsCurveOnSurface(S,l)) return;
778     itcr.Next();
779   }
780
781   Handle(Geom_Plane) GP;
782   Handle(Geom_RectangularTrimmedSurface) GRTS;
783   GRTS = Handle(Geom_RectangularTrimmedSurface)::DownCast(S);
784   if(!GRTS.IsNull())
785     GP = Handle(Geom_Plane)::DownCast(GRTS->BasisSurface());
786   else
787     GP = Handle(Geom_Plane)::DownCast(S);
788
789   Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve();
790   Handle(GeomAdaptor_HSurface) HS = new GeomAdaptor_HSurface();
791
792   TopLoc_Location LC;
793   Standard_Real First, Last;
794   GeomAdaptor_Curve& GAC = HC->ChangeCurve();
795   GAC.Load(BRep_Tool::Curve(theEdge,LC,First,Last));
796   LC = L.Predivided(LC);
797
798   if (!LC.IsIdentity()) {
799     GP = Handle(Geom_Plane)::DownCast(
800       GP->Transformed(LC.Transformation()));
801   }
802   GeomAdaptor_Surface& GAS = HS->ChangeSurface();
803   GAS.Load(GP);
804
805   ProjLib_ProjectedCurve Proj(HS,HC);
806   Handle(Geom2d_Curve) pc = Geom2dAdaptor::MakeCurve(Proj);
807
808   gp_Pln pln = GAS.Plane();
809   Standard_Real d2 = 0.;
810   Standard_Integer nn = 23;
811   Standard_Real unsurnn = 1./nn;
812   for(Standard_Integer i = 0; i <= nn; i++){
813     Standard_Real t = unsurnn*i;
814     Standard_Real u = First*(1.-t) + Last*t;
815     gp_Pnt Pc3d = HC->Value(u);
816     gp_Pnt2d p2d = pc->Value(u);
817     gp_Pnt Pcons = ElSLib::Value(p2d.X(),p2d.Y(),pln);
818     Standard_Real eps = Max(Pc3d.XYZ().SquareModulus(), Pcons.XYZ().SquareModulus());
819     eps = Epsilon(eps);
820     Standard_Real temp = Pc3d.SquareDistance(Pcons);
821     if(temp <= eps)
822     {
823       temp = 0.;
824     }
825     if(temp > d2) d2 = temp;
826   }
827   d2 = 1.5*sqrt(d2);
828   theEdTol = d2;
829 }
830
831 //=======================================================================
832 //function : UpdTolMap
833 //purpose  : Update map ShToTol (shape to tolerance)
834 //=======================================================================
835 static void UpdTolMap(const TopoDS_Shape& theSh, Standard_Real theNewTol, 
836   TopTools_DataMapOfShapeReal& theShToTol)
837 {
838   TopAbs_ShapeEnum aSt = theSh.ShapeType();
839   Standard_Real aShTol;
840   if (aSt == TopAbs_VERTEX)
841     aShTol = BRep_Tool::Tolerance(TopoDS::Vertex(theSh));
842   else if (aSt == TopAbs_EDGE)
843     aShTol = BRep_Tool::Tolerance(TopoDS::Edge(theSh));
844   else
845     return;
846   //
847   if (theNewTol > aShTol)
848   {
849     const Standard_Real* anOldtol = theShToTol.Seek(theSh);
850     if (!anOldtol)
851       theShToTol.Bind(theSh, theNewTol);
852     else
853       theShToTol(theSh) = Max(*anOldtol, theNewTol);
854   }
855 }
856
857 //=======================================================================
858 //function : UpdShTol
859 //purpose  : Update vertices/edges/faces according to ShToTol map (create copies of necessary)
860 //=======================================================================
861 static void UpdShTol(const TopTools_DataMapOfShapeReal& theShToTol,
862   const Standard_Boolean IsMutableInput, BRepTools_ReShape& theReshaper,
863   Standard_Boolean theVForceUpdate)
864 {
865   BRep_Builder aB;
866   TopTools_DataMapIteratorOfDataMapOfShapeReal SHToTolit(theShToTol);
867   for (;SHToTolit.More();SHToTolit.Next())
868   {
869     const TopoDS_Shape& aSh = SHToTolit.Key();
870     Standard_Real aTol = SHToTolit.Value();
871     //
872     TopoDS_Shape aNsh;
873     const TopoDS_Shape& aVsh = theReshaper.Value(aSh);
874     Standard_Boolean UseOldSh = IsMutableInput || theReshaper.IsNewShape(aSh) || !aVsh.IsSame(aSh);
875     if (UseOldSh)
876       aNsh = aVsh;
877     else
878     {
879       aNsh = aSh.EmptyCopied();
880       //add subshapes from the original shape
881       TopoDS_Iterator sit(aSh);
882       for (;sit.More();sit.Next())
883         aB.Add(aNsh, sit.Value());
884       //
885       aNsh.Free(aSh.Free());
886       aNsh.Checked(aSh.Checked());
887       aNsh.Orientable(aSh.Orientable());
888       aNsh.Closed(aSh.Closed());
889       aNsh.Infinite(aSh.Infinite());
890       aNsh.Convex(aSh.Convex());
891       //
892     }
893     //
894     switch (aSh.ShapeType())
895     {
896     case TopAbs_FACE: 
897       {
898         aB.UpdateFace(TopoDS::Face(aNsh), aTol); 
899         break;
900       }
901     case TopAbs_EDGE: 
902       {
903         aB.UpdateEdge(TopoDS::Edge(aNsh), aTol);   
904         break;
905       }
906     case TopAbs_VERTEX: 
907       {
908         const Handle(BRep_TVertex)& aTV = *((Handle(BRep_TVertex)*)&aNsh.TShape());
909         //
910         if(aTV->Locked())
911           throw TopoDS_LockedShape("BRep_Builder::UpdateVertex");
912         //
913         if (theVForceUpdate)
914           aTV->Tolerance(aTol);
915         else
916           aTV->UpdateTolerance(aTol);
917         aTV->Modified(Standard_True);
918         break;
919       }
920     default: 
921       break;
922     }
923     //
924     if (!UseOldSh)
925       theReshaper.Replace(aSh, aNsh);
926   }
927 }
928
929 //=======================================================================
930 //function : InternalSameParameter
931 //purpose  : 
932 //=======================================================================
933 static void InternalSameParameter(const TopoDS_Shape& theSh, BRepTools_ReShape& theReshaper,
934   const Standard_Real theTol, const Standard_Boolean IsForced, const Standard_Boolean IsMutableInput ) 
935 {
936   TopExp_Explorer ex(theSh,TopAbs_EDGE);
937   TopTools_MapOfShape  Done;
938   BRep_Builder aB;
939   TopTools_DataMapOfShapeReal aShToTol;
940  
941   while (ex.More()) 
942   {
943     const TopoDS_Edge& aCE = TopoDS::Edge(ex.Current());
944     if (Done.Add(aCE))
945     {
946       TopoDS_Edge aNE = TopoDS::Edge(theReshaper.Value(aCE));
947       Standard_Boolean UseOldEdge = IsMutableInput || theReshaper.IsNewShape(aCE) || !aNE.IsSame(aCE);
948       if (IsForced && (BRep_Tool::SameRange(aCE) || BRep_Tool::SameParameter(aCE)))
949       {
950         if (!UseOldEdge)
951         {
952           aNE = TopoDS::Edge(aCE.EmptyCopied());
953           TopoDS_Iterator sit(aCE);
954           for (;sit.More();sit.Next())
955             aB.Add(aNE, sit.Value());
956           theReshaper.Replace(aCE, aNE);
957           UseOldEdge = Standard_True;
958         }
959         aB.SameRange(aNE, Standard_False);
960         aB.SameParameter(aNE, Standard_False);
961       }
962       Standard_Real aNewTol = -1;
963       TopoDS_Edge aResEdge = BRepLib::SameParameter(aNE, theTol, aNewTol, UseOldEdge);
964       if (!UseOldEdge && !aResEdge.IsNull())
965         //NE have been empty-copied
966         theReshaper.Replace(aNE, aResEdge);
967       if (aNewTol > 0)
968       {
969         TopoDS_Vertex aV1, aV2;
970         TopExp::Vertices(aCE,aV1,aV2);
971         if (!aV1.IsNull())
972           UpdTolMap(aV1, aNewTol, aShToTol);
973         if (!aV2.IsNull()) 
974           UpdTolMap(aV2, aNewTol, aShToTol);
975       }
976     }
977     ex.Next();
978   }
979  
980   Done.Clear();
981   BRepAdaptor_Surface BS;
982   for(ex.Init(theSh,TopAbs_FACE); ex.More(); ex.Next()){
983     const TopoDS_Face& curface = TopoDS::Face(ex.Current());
984     if(!Done.Add(curface)) continue;
985     BS.Initialize(curface);
986     if(BS.GetType() != GeomAbs_Plane) continue;
987     TopExp_Explorer ex2;
988     for(ex2.Init(curface,TopAbs_EDGE); ex2.More(); ex2.Next()){
989       const TopoDS_Edge& E = TopoDS::Edge(ex2.Current());
990       TopoDS_Shape aNe = theReshaper.Value(E);
991       Standard_Real aNewEtol = -1;
992       GetEdgeTol(TopoDS::Edge(aNe), curface, aNewEtol);
993       if (aNewEtol >= 0) //not equal to -1
994         UpdTolMap(E, aNewEtol, aShToTol);
995     }
996   }
997  
998   //
999   UpdShTol(aShToTol, IsMutableInput, theReshaper, Standard_False );
1000
1001   InternalUpdateTolerances(theSh, Standard_False, IsMutableInput, theReshaper );
1002 }
1003
1004 //================================================================
1005 //function : SameParameter
1006 //WARNING  : New spec DUB LBO 9/9/97.
1007 //  Recode in the edge the best tolerance found, 
1008 //  for vertex extremities it is required to find something else
1009 //================================================================
1010 void  BRepLib::SameParameter(const TopoDS_Shape& S,
1011   const Standard_Real Tolerance,
1012   const Standard_Boolean forced)
1013 {
1014   BRepTools_ReShape reshaper;
1015   InternalSameParameter( S, reshaper, Tolerance, forced, Standard_True);
1016 }
1017
1018 //=======================================================================
1019 //function : SameParameter
1020 //purpose  : 
1021 //=======================================================================
1022 void BRepLib::SameParameter(const TopoDS_Shape& S, BRepTools_ReShape& theReshaper,
1023   const Standard_Real Tolerance, const Standard_Boolean forced ) 
1024 {
1025   InternalSameParameter( S, theReshaper, Tolerance, forced, Standard_False);
1026 }
1027
1028 //=======================================================================
1029 //function : EvalTol
1030 //purpose  : 
1031 //=======================================================================
1032 static Standard_Boolean EvalTol(const Handle(Geom2d_Curve)& pc,
1033   const Handle(Geom_Surface)& s,
1034   const GeomAdaptor_Curve&    gac,
1035   const Standard_Real         tol,
1036   Standard_Real&              tolbail)
1037 {
1038   Standard_Integer ok = 0;
1039   Standard_Real f = gac.FirstParameter();
1040   Standard_Real l = gac.LastParameter();
1041   Extrema_LocateExtPC Projector;
1042   Projector.Initialize(gac,f,l,tol);
1043   Standard_Real u,v;
1044   gp_Pnt p;
1045   tolbail = tol;
1046   for(Standard_Integer i = 1; i <= 5; i++){
1047     Standard_Real t = i/6.;
1048     t = (1.-t) * f + t * l;
1049     pc->Value(t).Coord(u,v);
1050     p = s->Value(u,v);
1051     Projector.Perform(p,t);
1052     if (Projector.IsDone()) {
1053       Standard_Real dist2 = Projector.SquareDistance();
1054       if(dist2 > tolbail * tolbail) tolbail = sqrt (dist2);
1055       ok++;
1056     }
1057   }
1058   return (ok > 2);
1059 }
1060
1061 //=======================================================================
1062 //function : ComputeTol
1063 //purpose  : 
1064 //=======================================================================
1065 static Standard_Real ComputeTol(const Handle(Adaptor3d_HCurve)& c3d,
1066   const Handle(Adaptor2d_HCurve2d)& c2d,
1067   const Handle(Adaptor3d_HSurface)& surf,
1068   const Standard_Integer        nbp)
1069
1070 {
1071
1072   TColStd_Array1OfReal dist(1,nbp+10);
1073   dist.Init(-1.);
1074
1075   //Adaptor3d_CurveOnSurface  cons(c2d,surf);
1076   Standard_Real uf = surf->FirstUParameter(), ul = surf->LastUParameter(),
1077                 vf = surf->FirstVParameter(), vl = surf->LastVParameter();
1078   Standard_Real du = 0.01 * (ul - uf), dv = 0.01 * (vl - vf);
1079   Standard_Boolean isUPeriodic = surf->IsUPeriodic(), isVPeriodic = surf->IsVPeriodic();
1080   Standard_Real DSdu = 1./surf->UResolution(1.), DSdv = 1./surf->VResolution(1.);
1081   Standard_Real d2 = 0.;
1082   Standard_Real first = c3d->FirstParameter();
1083   Standard_Real last  = c3d->LastParameter();
1084   Standard_Real dapp = -1.;
1085   for (Standard_Integer i = 0; i <= nbp; ++i)
1086   {
1087     const Standard_Real t = IntToReal(i)/IntToReal(nbp);
1088     const Standard_Real u = first*(1.-t) + last*t;
1089     gp_Pnt Pc3d = c3d->Value(u);
1090     gp_Pnt2d Puv = c2d->Value(u);
1091     if(!isUPeriodic)
1092     {
1093       if(Puv.X() < uf - du)
1094       {
1095         dapp = Max(dapp, DSdu * (uf - Puv.X()));
1096         continue;
1097       }
1098       else if(Puv.X() > ul + du)
1099       {
1100         dapp = Max(dapp, DSdu * (Puv.X() - ul));
1101         continue;
1102       }
1103     }
1104     if(!isVPeriodic)
1105     {
1106       if(Puv.Y() < vf - dv)
1107       {
1108         dapp = Max(dapp, DSdv * (vf - Puv.Y()));
1109         continue;
1110       }
1111       else if(Puv.Y() > vl + dv)
1112       {
1113         dapp = Max(dapp, DSdv * (Puv.Y() - vl));
1114         continue;
1115       }
1116     }
1117     gp_Pnt Pcons = surf->Value(Puv.X(), Puv.Y());
1118     if (Precision::IsInfinite(Pcons.X()) ||
1119         Precision::IsInfinite(Pcons.Y()) ||
1120         Precision::IsInfinite(Pcons.Z()))
1121     {
1122       d2 = Precision::Infinite();
1123       break;
1124     }
1125     Standard_Real temp = Pc3d.SquareDistance(Pcons);
1126
1127     dist(i+1) = temp;
1128
1129     d2 = Max (d2, temp);
1130   }
1131
1132   if(Precision::IsInfinite(d2))
1133   {
1134     return d2;
1135   }
1136
1137   d2 = Sqrt(d2);
1138   if(dapp > d2)
1139   {
1140     return dapp;
1141   }
1142
1143   Standard_Boolean ana = Standard_False;
1144   Standard_Real D2 = 0;
1145   Standard_Integer N1 = 0;
1146   Standard_Integer N2 = 0;
1147   Standard_Integer N3 = 0;
1148
1149   for (Standard_Integer i = 1; i<= nbp+10; ++i)
1150   {
1151     if (dist(i) > 0)
1152     {
1153       if (dist(i) < 1.0)
1154       {
1155         ++N1;
1156       }
1157       else
1158       {
1159         ++N2;
1160       }
1161     }
1162   }
1163
1164   if (N1 > N2 && N2 != 0)
1165   {
1166     N3 = 100*N2/(N1+N2);
1167   }
1168   if (N3 < 10 && N3 != 0)
1169   {
1170     ana = Standard_True;
1171     for (Standard_Integer i = 1; i <= nbp+10; ++i)
1172     {
1173       if (dist(i) > 0 && dist(i) < 1.0)
1174       {
1175         D2 = Max (D2, dist(i));
1176       }
1177     }
1178   }
1179
1180   //d2 = 1.5*sqrt(d2);
1181   d2 = (!ana) ? 1.5 * d2 : 1.5*sqrt(D2);
1182   d2 = Max (d2, 1.e-7);
1183   return d2;
1184 }
1185
1186 //=======================================================================
1187 //function : GetCurve3d
1188 //purpose  : 
1189 //=======================================================================
1190 static void GetCurve3d(const TopoDS_Edge& theEdge, Handle(Geom_Curve)& theC3d, Standard_Real& theF3d, 
1191   Standard_Real& theL3d, TopLoc_Location& theLoc3d, BRep_ListOfCurveRepresentation& theCList)
1192 {
1193   const Handle(BRep_TEdge)& aTE = *((Handle(BRep_TEdge)*) &theEdge.TShape());
1194   theCList = aTE->ChangeCurves(); // current function (i.e. GetCurve3d()) will not change any of this curves
1195   BRep_ListIteratorOfListOfCurveRepresentation anIt(theCList);
1196   Standard_Boolean NotDone = Standard_True;
1197   while (NotDone && anIt.More()) {
1198     Handle(BRep_GCurve) GCurve = Handle(BRep_GCurve)::DownCast(anIt.Value());
1199     if (!GCurve.IsNull() && GCurve->IsCurve3D()) {
1200       theC3d = GCurve->Curve3D() ;
1201       theF3d = GCurve->First();
1202       theL3d = GCurve->Last();
1203       theLoc3d = GCurve->Location() ;
1204       NotDone = Standard_False;
1205     } 
1206     anIt.Next() ;
1207   }
1208 }
1209
1210 //=======================================================================
1211 //function : UpdateVTol
1212 //purpose  : 
1213 //=======================================================================
1214 void UpdateVTol(const TopoDS_Vertex theV1, const TopoDS_Vertex& theV2, Standard_Real theTol)
1215 {
1216   BRep_Builder aB;
1217   if (!theV1.IsNull())
1218     aB.UpdateVertex(theV1,theTol);
1219   if (!theV2.IsNull())
1220     aB.UpdateVertex(theV2,theTol);
1221 }
1222
1223 //=======================================================================
1224 //function : SameParameter
1225 //purpose  : 
1226 //=======================================================================
1227 void BRepLib::SameParameter(const TopoDS_Edge& theEdge,
1228   const Standard_Real theTolerance)
1229 {
1230   Standard_Real aNewTol = -1;
1231   SameParameter(theEdge, theTolerance, aNewTol, Standard_True);
1232   if (aNewTol > 0)
1233   {
1234     TopoDS_Vertex aV1, aV2;
1235     TopExp::Vertices(theEdge,aV1,aV2);
1236     UpdateVTol(aV1, aV2, aNewTol);
1237   }
1238 }
1239
1240 //=======================================================================
1241 //function : SameParameter
1242 //purpose  : 
1243 //=======================================================================
1244 TopoDS_Edge BRepLib::SameParameter(const TopoDS_Edge& theEdge,
1245   const Standard_Real theTolerance, Standard_Real& theNewTol, Standard_Boolean IsUseOldEdge)
1246 {
1247   if (BRep_Tool::SameParameter(theEdge)) 
1248     return TopoDS_Edge();
1249   Standard_Real f3d =0.,l3d =0.;
1250   TopLoc_Location L3d;
1251   Handle(Geom_Curve) C3d;
1252   BRep_ListOfCurveRepresentation CList;
1253   GetCurve3d(theEdge, C3d, f3d, l3d, L3d, CList);
1254   if(C3d.IsNull()) 
1255     return TopoDS_Edge();
1256
1257   BRep_Builder B;
1258   TopoDS_Edge aNE;
1259   Handle(BRep_TEdge) aNTE;
1260   if (IsUseOldEdge)
1261   {
1262     aNE = theEdge;
1263     aNTE = *((Handle(BRep_TEdge)*) &theEdge.TShape());
1264   }
1265   else
1266   {
1267     aNE = TopoDS::Edge(theEdge.EmptyCopied()); //will be modified a little bit later, so copy anyway  
1268     GetCurve3d(aNE, C3d, f3d, l3d, L3d, CList); //C3d pointer and CList will be differ after copying
1269     aNTE = *((Handle(BRep_TEdge)*) &aNE.TShape());
1270     TopoDS_Iterator sit(theEdge);
1271     for (;sit.More();sit.Next()) //add vertices from old edge to the new ones
1272       B.Add(aNE, sit.Value());
1273   }
1274
1275   BRep_ListIteratorOfListOfCurveRepresentation It(CList);
1276
1277   const Standard_Integer NCONTROL = 22;
1278
1279   Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve();
1280   Handle(Geom2dAdaptor_HCurve) HC2d = new Geom2dAdaptor_HCurve();
1281   Handle(GeomAdaptor_HSurface) HS = new GeomAdaptor_HSurface();
1282   GeomAdaptor_Curve& GAC = HC->ChangeCurve();
1283   Geom2dAdaptor_Curve& GAC2d = HC2d->ChangeCurve2d();
1284   GeomAdaptor_Surface& GAS = HS->ChangeSurface(); 
1285
1286   // modified by NIZHNY-OCC486  Tue Aug 27 17:15:13 2002 :
1287   Standard_Boolean m_TrimmedPeriodical = Standard_False;
1288   Handle(Standard_Type) TheType = C3d->DynamicType();
1289   if( TheType == STANDARD_TYPE(Geom_TrimmedCurve))
1290   {
1291     Handle(Geom_Curve) gtC (Handle(Geom_TrimmedCurve)::DownCast (C3d)->BasisCurve());
1292     m_TrimmedPeriodical = gtC->IsPeriodic();
1293   }
1294   // modified by NIZHNY-OCC486  Tue Aug 27 17:15:17 2002 .
1295
1296
1297   if(!C3d->IsPeriodic()) {
1298     Standard_Real Udeb = C3d->FirstParameter();
1299     Standard_Real Ufin = C3d->LastParameter();
1300     // modified by NIZHNY-OCC486  Tue Aug 27 17:17:14 2002 :
1301     //if (Udeb > f3d) f3d = Udeb;
1302     //if (l3d > Ufin) l3d = Ufin;
1303     if(!m_TrimmedPeriodical)
1304     {
1305       if (Udeb > f3d) f3d = Udeb;
1306       if (l3d > Ufin) l3d = Ufin;
1307     }
1308     // modified by NIZHNY-OCC486  Tue Aug 27 17:17:55 2002 .
1309   }
1310   if(!L3d.IsIdentity()){
1311     C3d = Handle(Geom_Curve)::DownCast(C3d->Transformed(L3d.Transformation()));
1312   }
1313   GAC.Load(C3d,f3d,l3d);
1314
1315   Standard_Real Prec_C3d = BRepCheck::PrecCurve(GAC);
1316
1317   Standard_Boolean IsSameP = 1;
1318   Standard_Real maxdist = 0.;
1319
1320   //  Modified by skv - Thu Jun  3 12:39:19 2004 OCC5898 Begin
1321   Standard_Real anEdgeTol = BRep_Tool::Tolerance(aNE);
1322   //  Modified by skv - Thu Jun  3 12:39:20 2004 OCC5898 End
1323   Standard_Boolean SameRange = BRep_Tool::SameRange(aNE);
1324   Standard_Boolean YaPCu = Standard_False;
1325   const Standard_Real BigError = 1.e10;
1326   It.Initialize(CList);
1327
1328   while (It.More()) {
1329
1330     Standard_Boolean isANA = Standard_False;
1331     Standard_Boolean isBSP = Standard_False;
1332     Handle(BRep_GCurve) GCurve = Handle(BRep_GCurve)::DownCast(It.Value());
1333     Handle(Geom2d_Curve) PC[2];
1334     Handle(Geom_Surface) S;
1335     if (!GCurve.IsNull() && GCurve->IsCurveOnSurface()) {
1336       YaPCu = Standard_True;
1337       PC[0] = GCurve->PCurve();
1338       TopLoc_Location PCLoc = GCurve->Location();
1339       S = GCurve->Surface();
1340       if (!PCLoc.IsIdentity() ) {
1341         S = Handle(Geom_Surface)::DownCast(S->Transformed(PCLoc.Transformation()));
1342       }
1343
1344       GAS.Load(S);
1345       if (GCurve->IsCurveOnClosedSurface()) {
1346         PC[1] = GCurve->PCurve2();
1347       }
1348
1349       // Eval tol2d to compute SameRange
1350       Standard_Real UResol = Max(GAS.UResolution(theTolerance), Precision::PConfusion());
1351       Standard_Real VResol = Max(GAS.VResolution(theTolerance), Precision::PConfusion());
1352       Standard_Real Tol2d  = Min(UResol, VResol);
1353       for(Standard_Integer i = 0; i < 2; i++){
1354         Handle(Geom2d_Curve) curPC = PC[i];
1355         Standard_Boolean updatepc = 0;
1356         if(curPC.IsNull()) break;
1357         if(!SameRange){
1358           GeomLib::SameRange(Tol2d,
1359             PC[i],GCurve->First(),GCurve->Last(),
1360             f3d,l3d,curPC);
1361
1362           updatepc = (curPC != PC[i]);
1363
1364         }
1365         Standard_Boolean goodpc = 1;
1366         GAC2d.Load(curPC,f3d,l3d);
1367
1368         Standard_Real error = ComputeTol(HC, HC2d, HS, NCONTROL);
1369
1370         if(error > BigError)
1371         {
1372           maxdist = error;
1373           break;
1374         }
1375
1376         if(GAC2d.GetType() == GeomAbs_BSplineCurve && 
1377           GAC2d.Continuity() == GeomAbs_C0) {
1378             Handle(Geom2d_BSplineCurve) bs2d = GAC2d.BSpline();
1379             Handle(Geom2d_BSplineCurve) bs2dsov = bs2d;
1380             Standard_Real fC0 = bs2d->FirstParameter(), lC0 = bs2d->LastParameter();
1381             Standard_Boolean repar = Standard_True;
1382             gp_Pnt2d OriginPoint;
1383             bs2d->D0(fC0, OriginPoint);
1384             Geom2dConvert::C0BSplineToC1BSplineCurve(bs2d, Tol2d);
1385             isBSP = Standard_True; 
1386
1387             if(bs2d->IsPeriodic()) { // -------- IFV, Jan 2000
1388               gp_Pnt2d NewOriginPoint;
1389               bs2d->D0(bs2d->FirstParameter(), NewOriginPoint);
1390               if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1391                 Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) {
1392
1393                   TColStd_Array1OfReal Knotbs2d (1, bs2d->NbKnots());
1394                   bs2d->Knots(Knotbs2d);
1395
1396                   for(Standard_Integer Index = 1; Index <= bs2d->NbKnots(); Index++) {
1397                     bs2d->D0(Knotbs2d(Index), NewOriginPoint);
1398                     if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1399                       Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) continue;
1400
1401                     bs2d->SetOrigin(Index);
1402                     break;
1403                   }
1404               }
1405             }
1406
1407             if(bs2d->Continuity() == GeomAbs_C0) {
1408               Standard_Real tolbail;
1409               if(EvalTol(curPC,S,GAC,theTolerance,tolbail)){
1410                 bs2d = bs2dsov;
1411                 Standard_Real UResbail = GAS.UResolution(tolbail);
1412                 Standard_Real VResbail = GAS.VResolution(tolbail);
1413                 Standard_Real Tol2dbail  = Min(UResbail,VResbail);
1414                 bs2d->D0(bs2d->FirstParameter(), OriginPoint); 
1415
1416                 Standard_Integer nbp = bs2d->NbPoles();
1417                 TColgp_Array1OfPnt2d poles(1,nbp);
1418                 bs2d->Poles(poles);
1419                 gp_Pnt2d p = poles(1), p1;
1420                 Standard_Real d = Precision::Infinite();
1421                 for(Standard_Integer ip = 2; ip <= nbp; ip++) {
1422                   p1 = poles(ip);
1423                   d = Min(d,p.SquareDistance(p1));
1424                   p = p1;
1425                 }
1426                 d = sqrt(d)*.1;
1427
1428                 Tol2dbail = Max(Min(Tol2dbail,d),Tol2d);
1429
1430                 Geom2dConvert::C0BSplineToC1BSplineCurve(bs2d,Tol2dbail);
1431
1432                 if(bs2d->IsPeriodic()) { // -------- IFV, Jan 2000
1433                   gp_Pnt2d NewOriginPoint;
1434                   bs2d->D0(bs2d->FirstParameter(), NewOriginPoint);
1435                   if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1436                     Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) {
1437
1438                       TColStd_Array1OfReal Knotbs2d (1, bs2d->NbKnots());
1439                       bs2d->Knots(Knotbs2d);
1440
1441                       for(Standard_Integer Index = 1; Index <= bs2d->NbKnots(); Index++) {
1442                         bs2d->D0(Knotbs2d(Index), NewOriginPoint);
1443                         if(Abs(OriginPoint.X() - NewOriginPoint.X()) > Precision::PConfusion() ||
1444                           Abs(OriginPoint.Y() - NewOriginPoint.Y()) > Precision::PConfusion()    ) continue;
1445
1446                         bs2d->SetOrigin(Index);
1447                         break;
1448                       }
1449                   }
1450                 }
1451
1452
1453                 if(bs2d->Continuity() == GeomAbs_C0) {
1454                   goodpc = 1;
1455                   bs2d = bs2dsov;
1456                   repar = Standard_False;
1457                 }
1458               }
1459               else goodpc = 0;
1460             }
1461
1462             if(goodpc){
1463               if(repar) {
1464                 Standard_Integer NbKnots = bs2d->NbKnots();
1465                 TColStd_Array1OfReal Knots(1,NbKnots);
1466                 bs2d->Knots(Knots);
1467                 //          BSplCLib::Reparametrize(f3d,l3d,Knots);
1468                 BSplCLib::Reparametrize(fC0,lC0,Knots);
1469                 bs2d->SetKnots(Knots);
1470                 GAC2d.Load(bs2d,f3d,l3d);
1471                 curPC = bs2d;
1472                 Standard_Boolean updatepcsov = updatepc;
1473                 updatepc = Standard_True;
1474
1475                 Standard_Real error1 = ComputeTol(HC, HC2d, HS, NCONTROL);
1476                 if(error1 > error) {
1477                   bs2d = bs2dsov;
1478                   GAC2d.Load(bs2d,f3d,l3d);
1479                   curPC = bs2d;
1480                   updatepc = updatepcsov;
1481                   isANA = Standard_True;
1482                 }
1483                 else {
1484                   error = error1;
1485                 }
1486               }
1487
1488               //check, if new BSpline "good" or not --------- IFV, Jan of 2000
1489               GeomAbs_Shape cont = bs2d->Continuity();
1490               Standard_Boolean IsBad = Standard_False;
1491
1492               if(cont > GeomAbs_C0 && error > Max(1.e-3,theTolerance)) {
1493                 Standard_Integer NbKnots = bs2d->NbKnots();
1494                 TColStd_Array1OfReal Knots(1,NbKnots);
1495                 bs2d->Knots(Knots);
1496                 Standard_Real critratio = 10.; 
1497                 Standard_Real dtprev = Knots(2) - Knots(1), dtratio = 1.;
1498                 Standard_Real dtmin = dtprev;
1499                 Standard_Real dtcur;
1500                 for(Standard_Integer j = 2; j < NbKnots; j++) {
1501                   dtcur = Knots(j+1) - Knots(j);
1502                   dtmin = Min(dtmin, dtcur);
1503
1504                   if(IsBad) continue;
1505
1506                   if(dtcur > dtprev) dtratio = dtcur/dtprev;
1507                   else dtratio = dtprev/dtcur;
1508                   if(dtratio > critratio) {IsBad = Standard_True;}
1509                   dtprev = dtcur;
1510
1511                 }
1512                 if(IsBad) {
1513                   // To avoid failures in Approx_CurvilinearParameter 
1514                   bs2d->Resolution(Max(1.e-3,theTolerance), dtcur);
1515                   if(dtmin < dtcur) IsBad = Standard_False;
1516                 }
1517               }
1518
1519
1520               if(IsBad ) { //if BSpline "bad", try to reparametrize it
1521                 // by its curve length
1522
1523                 //            GeomAbs_Shape cont = bs2d->Continuity();
1524                 if(cont > GeomAbs_C2) cont = GeomAbs_C2;
1525                 Standard_Integer maxdeg = bs2d->Degree();
1526                 if(maxdeg == 1) maxdeg = 14;
1527                 Approx_CurvilinearParameter AppCurPar(HC2d, HS, Max(1.e-3,theTolerance),
1528                   cont, maxdeg, 10);
1529                 if(AppCurPar.IsDone() || AppCurPar.HasResult()) {
1530                   bs2d = AppCurPar.Curve2d1();
1531                   GAC2d.Load(bs2d,f3d,l3d);
1532                   curPC = bs2d;
1533
1534                   if(Abs(bs2d->FirstParameter() - fC0) > Tol2d ||
1535                     Abs(bs2d->LastParameter() - lC0) > Tol2d    ) {
1536                       Standard_Integer NbKnots = bs2d->NbKnots();
1537                       TColStd_Array1OfReal Knots(1,NbKnots);
1538                       bs2d->Knots(Knots);
1539                       //                  BSplCLib::Reparametrize(f3d,l3d,Knots);
1540                       BSplCLib::Reparametrize(fC0,lC0,Knots);
1541                       bs2d->SetKnots(Knots);
1542                       GAC2d.Load(bs2d,f3d,l3d);
1543                       curPC = bs2d;
1544
1545                   }
1546                 }
1547               }
1548
1549
1550             }
1551         }
1552
1553
1554         if(goodpc){
1555           //      Approx_SameParameter SameP(HC,HC2d,HS,Tolerance);
1556           Standard_Real aTol = (isANA && isBSP) ? 1.e-7 : theTolerance;
1557           const Handle(Adaptor3d_HCurve)& aHCurv = HC; // to avoid ambiguity
1558           const Handle(Adaptor2d_HCurve2d)& aHCurv2d = HC2d; // to avoid ambiguity
1559           Approx_SameParameter SameP(aHCurv,aHCurv2d,HS,aTol);
1560
1561           if (SameP.IsSameParameter()) {
1562             maxdist = Max(maxdist,SameP.TolReached());
1563             if(updatepc){
1564               if (i == 0) GCurve->PCurve(curPC);
1565               else GCurve->PCurve2(curPC);
1566             }
1567           }
1568           else if (SameP.IsDone()) {
1569             Standard_Real tolreached = SameP.TolReached();
1570             if(tolreached <= error) {
1571               curPC = SameP.Curve2d();
1572               updatepc = Standard_True;
1573               maxdist = Max(maxdist,tolreached);
1574             }
1575             else {
1576               maxdist = Max(maxdist, error);
1577             }
1578             if(updatepc){
1579               if (i == 0) GCurve->PCurve(curPC);
1580               else GCurve->PCurve2(curPC);
1581             }
1582           }
1583           else
1584           {
1585             //Approx_SameParameter has failed.
1586             //Consequently, the situation might be,
1587             //when 3D and 2D-curve do not have same-range.
1588             GeomLib::SameRange( Tol2d, PC[i], 
1589                                 GCurve->First(), GCurve->Last(),
1590                                 f3d,l3d,curPC);
1591
1592             if (i == 0) GCurve->PCurve(curPC);
1593             else GCurve->PCurve2(curPC);
1594
1595             IsSameP = 0;
1596           }
1597
1598         }
1599         else IsSameP = 0;
1600
1601         //  Modified by skv - Thu Jun  3 12:39:19 2004 OCC5898 Begin
1602         if (!IsSameP) {
1603           Standard_Real Prec_Surf = BRepCheck::PrecSurface(HS);
1604           Standard_Real CurTol = anEdgeTol + Max(Prec_C3d, Prec_Surf);
1605           if (CurTol >= error) {
1606             maxdist = Max(maxdist, anEdgeTol);
1607             IsSameP = Standard_True;
1608           }
1609         }
1610         //  Modified by skv - Thu Jun  3 12:39:20 2004 OCC5898 End
1611       }
1612     }
1613     It.Next() ;
1614   }
1615   B.Range(aNE,f3d,l3d);
1616   B.SameRange(aNE,Standard_True);
1617   if ( IsSameP) {
1618     // Reduce eventually the tolerance of the edge, as
1619     // all its representations are processed (except for some associated
1620     // to planes and not stored in the edge !) 
1621     // The same cannot be done with vertices that cannot be enlarged 
1622     // or left as is.
1623     if (YaPCu) {
1624       // Avoid setting too small tolerances.
1625       maxdist = Max(maxdist,Precision::Confusion());
1626       theNewTol = maxdist;
1627       aNTE->Modified(Standard_True);
1628       aNTE->Tolerance(maxdist);
1629     }
1630     B.SameParameter(aNE,Standard_True);
1631   }
1632   
1633   return aNE;
1634 }
1635
1636 //=======================================================================
1637 //function : InternalUpdateTolerances
1638 //purpose  : 
1639 //=======================================================================
1640 static void InternalUpdateTolerances(const TopoDS_Shape& theOldShape,
1641   const Standard_Boolean IsVerifyTolerance, const Standard_Boolean IsMutableInput, BRepTools_ReShape& theReshaper)
1642 {
1643   TopTools_DataMapOfShapeReal aShToTol;
1644   // Harmonize tolerances
1645   // with rule Tolerance(VERTEX)>=Tolerance(EDGE)>=Tolerance(FACE)
1646   Standard_Real tol=0;
1647   if (IsVerifyTolerance) {
1648     // Set tolerance to its minimum value
1649     Handle(Geom_Surface) S;
1650     TopLoc_Location l;
1651     TopExp_Explorer ex;
1652     Bnd_Box aB;
1653     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, dMax;
1654     for (ex.Init(theOldShape, TopAbs_FACE); ex.More(); ex.Next()) {
1655       const TopoDS_Face& curf=TopoDS::Face(ex.Current());
1656       S = BRep_Tool::Surface(curf, l);
1657       if (!S.IsNull()) {
1658         aB.SetVoid();
1659         BRepBndLib::Add(curf,aB);
1660         if (S->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
1661           S = Handle(Geom_RectangularTrimmedSurface)::DownCast (S)->BasisSurface();
1662         }
1663         GeomAdaptor_Surface AS(S);
1664         switch (AS.GetType()) {
1665         case GeomAbs_Plane: 
1666         case GeomAbs_Cylinder: 
1667         case GeomAbs_Cone: 
1668           {
1669             tol=Precision::Confusion();
1670             break;
1671           }
1672         case GeomAbs_Sphere: 
1673         case GeomAbs_Torus: 
1674           {
1675             tol=Precision::Confusion()*2;
1676             break;
1677           }
1678         default:
1679           tol=Precision::Confusion()*4;
1680         }
1681         if (!aB.IsWhole()) {
1682           aB.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
1683           dMax=1.;
1684           if (!aB.IsOpenXmin() && !aB.IsOpenXmax()) dMax=aXmax-aXmin;
1685           if (!aB.IsOpenYmin() && !aB.IsOpenYmax()) aYmin=aYmax-aYmin;
1686           if (!aB.IsOpenZmin() && !aB.IsOpenZmax()) aZmin=aZmax-aZmin;
1687           if (aYmin>dMax) dMax=aYmin;
1688           if (aZmin>dMax) dMax=aZmin;
1689           tol=tol*dMax;
1690           // Do not process tolerances > 1.
1691           if (tol>1.) tol=0.99;
1692         }
1693         aShToTol.Bind(curf, tol);
1694       }
1695     }
1696   }
1697
1698   //Process edges
1699   TopTools_IndexedDataMapOfShapeListOfShape parents;
1700   TopExp::MapShapesAndAncestors(theOldShape, TopAbs_EDGE, TopAbs_FACE, parents);
1701   TopTools_ListIteratorOfListOfShape lConx;
1702   Standard_Integer iCur;
1703   for (iCur=1; iCur<=parents.Extent(); iCur++) {
1704     tol=0;
1705     for (lConx.Initialize(parents(iCur)); lConx.More(); lConx.Next()) 
1706     {
1707       const TopoDS_Face& FF = TopoDS::Face(lConx.Value());
1708       Standard_Real Ftol;
1709       if (IsVerifyTolerance && aShToTol.IsBound(FF)) //first condition for speed-up
1710         Ftol = aShToTol(FF);
1711       else
1712         Ftol = BRep_Tool::Tolerance(FF); //tolerance have not been updated
1713       tol=Max(tol, Ftol);
1714     }
1715     // Update can only increase tolerance, so if the edge has a greater
1716     //  tolerance than its faces it is not concerned
1717     const TopoDS_Edge& EK = TopoDS::Edge(parents.FindKey(iCur));
1718     if (tol > BRep_Tool::Tolerance(EK))
1719       aShToTol.Bind(EK, tol);
1720   }
1721
1722   //Vertices are processed
1723   const Standard_Real BigTol = 1.e10;
1724   parents.Clear();
1725
1726   TopExp::MapShapesAndUniqueAncestors(theOldShape, TopAbs_VERTEX, TopAbs_EDGE, parents);
1727   TColStd_MapOfTransient Initialized;
1728   Standard_Integer nbV = parents.Extent();
1729   for (iCur=1; iCur<=nbV; iCur++) {
1730     tol=0;
1731     const TopoDS_Vertex& V = TopoDS::Vertex(parents.FindKey(iCur));
1732     Bnd_Box box;
1733     box.Add(BRep_Tool::Pnt(V));
1734     gp_Pnt p3d;
1735     for (lConx.Initialize(parents(iCur)); lConx.More(); lConx.Next()) {
1736       const TopoDS_Edge& E = TopoDS::Edge(lConx.Value());
1737       const Standard_Real* aNtol = aShToTol.Seek(E);
1738       tol=Max(tol, aNtol ? *aNtol : BRep_Tool::Tolerance(E));
1739       if(tol > BigTol) continue;
1740       if(!BRep_Tool::SameRange(E)) continue;
1741       Standard_Real par = BRep_Tool::Parameter(V,E);
1742       Handle(BRep_TEdge)& TE = *((Handle(BRep_TEdge)*)&E.TShape());
1743       BRep_ListIteratorOfListOfCurveRepresentation itcr(TE->Curves());
1744       const TopLoc_Location& Eloc = E.Location();
1745       while (itcr.More()) {
1746         // For each CurveRepresentation, check the provided parameter
1747         const Handle(BRep_CurveRepresentation)& cr = itcr.Value();
1748         const TopLoc_Location& loc = cr->Location();
1749         TopLoc_Location L = (Eloc * loc);
1750         if (cr->IsCurve3D()) {
1751           const Handle(Geom_Curve)& C = cr->Curve3D();
1752           if (!C.IsNull()) { // edge non degenerated
1753             p3d = C->Value(par);
1754             p3d.Transform(L.Transformation());
1755             box.Add(p3d);
1756           }
1757         }
1758         else if (cr->IsCurveOnSurface()) {
1759           const Handle(Geom_Surface)& Su = cr->Surface();
1760           const Handle(Geom2d_Curve)& PC = cr->PCurve();
1761           Handle(Geom2d_Curve) PC2;
1762           if (cr->IsCurveOnClosedSurface()) {
1763             PC2 = cr->PCurve2();
1764           }
1765           gp_Pnt2d p2d = PC->Value(par);
1766           p3d = Su->Value(p2d.X(),p2d.Y());
1767           p3d.Transform(L.Transformation());
1768           box.Add(p3d);
1769           if (!PC2.IsNull()) {
1770             p2d = PC2->Value(par);
1771             p3d = Su->Value(p2d.X(),p2d.Y());
1772             p3d.Transform(L.Transformation());
1773             box.Add(p3d);
1774           }
1775         }
1776         itcr.Next();
1777       }
1778     }
1779     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
1780     box.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
1781     aXmax -= aXmin; aYmax -= aYmin; aZmax -= aZmin;
1782     tol = Max(tol,sqrt(aXmax*aXmax+aYmax*aYmax+aZmax*aZmax));
1783     tol += 2.*Epsilon(tol);
1784     //
1785     Standard_Real aVTol = BRep_Tool::Tolerance(V);
1786     Standard_Boolean anUpdTol = tol > aVTol;
1787     const Handle(BRep_TVertex)& aTV = *((Handle(BRep_TVertex)*)&V.TShape());
1788     Standard_Boolean toAdd = Standard_False;
1789     if (IsVerifyTolerance) 
1790     {
1791       // ASet minimum value of the tolerance 
1792       // Attention to sharing of the vertex by other shapes      
1793       toAdd = Initialized.Add(aTV) && aVTol != tol; //if Vtol == tol => no need to update toler
1794     }
1795     //'Initialized' map is not used anywhere outside this block
1796     if (anUpdTol || toAdd)
1797       aShToTol.Bind(V, tol);
1798   }
1799
1800   UpdShTol(aShToTol, IsMutableInput, theReshaper, Standard_True);
1801 }
1802
1803 //=======================================================================
1804 //function : UpdateTolerances
1805 //purpose  : 
1806 //=======================================================================
1807 void BRepLib::UpdateTolerances (const TopoDS_Shape& S, const Standard_Boolean verifyFaceTolerance)
1808 {
1809   BRepTools_ReShape aReshaper;
1810   InternalUpdateTolerances(S, verifyFaceTolerance, Standard_True, aReshaper);
1811 }
1812
1813 //=======================================================================
1814 //function : UpdateTolerances
1815 //purpose  : 
1816 //=======================================================================
1817 void BRepLib::UpdateTolerances (const TopoDS_Shape& S, BRepTools_ReShape& theReshaper, const Standard_Boolean verifyFaceTolerance )
1818 {
1819   InternalUpdateTolerances(S, verifyFaceTolerance, Standard_False, theReshaper);
1820 }
1821
1822 //=======================================================================
1823 //function : UpdateInnerTolerances
1824 //purpose  : 
1825 //=======================================================================
1826 void  BRepLib::UpdateInnerTolerances(const TopoDS_Shape& aShape)
1827 {
1828   TopTools_IndexedDataMapOfShapeListOfShape EFmap;
1829   TopExp::MapShapesAndAncestors(aShape, TopAbs_EDGE, TopAbs_FACE, EFmap);
1830   BRep_Builder BB;
1831   for (Standard_Integer i = 1; i <= EFmap.Extent(); i++)
1832   {
1833     TopoDS_Edge anEdge = TopoDS::Edge(EFmap.FindKey(i));
1834
1835     if (!BRep_Tool::IsGeometric(anEdge))
1836     {
1837       continue;
1838     }
1839
1840     TopoDS_Vertex V1, V2;
1841     TopExp::Vertices(anEdge, V1, V2);
1842     Standard_Real fpar, lpar;
1843     BRep_Tool::Range(anEdge, fpar, lpar);
1844     Standard_Real TolEdge = BRep_Tool::Tolerance(anEdge);
1845     gp_Pnt Pnt1, Pnt2;
1846     Handle(BRepAdaptor_HCurve) anHCurve = new BRepAdaptor_HCurve();
1847     anHCurve->ChangeCurve().Initialize(anEdge);
1848     if (!V1.IsNull())
1849       Pnt1 = BRep_Tool::Pnt(V1);
1850     if (!V2.IsNull())
1851       Pnt2 = BRep_Tool::Pnt(V2);
1852     
1853     if (!BRep_Tool::Degenerated(anEdge) &&
1854         EFmap(i).Extent() > 0)
1855     {
1856       NCollection_Sequence<Handle(Adaptor3d_HCurve)> theRep;
1857       theRep.Append(anHCurve);
1858       TopTools_ListIteratorOfListOfShape itl(EFmap(i));
1859       for (; itl.More(); itl.Next())
1860       {
1861         const TopoDS_Face& aFace = TopoDS::Face(itl.Value());
1862         Handle(BRepAdaptor_HCurve) anHCurvOnSurf = new BRepAdaptor_HCurve();
1863         anHCurvOnSurf->ChangeCurve().Initialize(anEdge, aFace);
1864         theRep.Append(anHCurvOnSurf);
1865       }
1866       
1867       const Standard_Integer NbSamples = (BRep_Tool::SameParameter(anEdge))? 23 : 2;
1868       Standard_Real delta = (lpar - fpar)/(NbSamples-1);
1869       Standard_Real MaxDist = 0.;
1870       for (Standard_Integer j = 2; j <= theRep.Length(); j++)
1871       {
1872         for (Standard_Integer k = 0; k <= NbSamples; k++)
1873         {
1874           Standard_Real ParamOnCenter = (k == NbSamples)? lpar :
1875             fpar + k*delta;
1876           gp_Pnt Center = theRep(1)->Value(ParamOnCenter);
1877           Standard_Real ParamOnCurve = (BRep_Tool::SameParameter(anEdge))? ParamOnCenter
1878             : ((k == 0)? theRep(j)->FirstParameter() : theRep(j)->LastParameter());
1879           gp_Pnt aPoint = theRep(j)->Value(ParamOnCurve);
1880           Standard_Real aDist = Center.Distance(aPoint);
1881           //aDist *= 1.1;
1882           aDist += 2.*Epsilon(aDist);
1883           if (aDist > MaxDist)
1884             MaxDist = aDist;
1885
1886           //Update tolerances of vertices
1887           if (k == 0 && !V1.IsNull())
1888           {
1889             Standard_Real aDist1 = Pnt1.Distance(aPoint);
1890             aDist1 += 2.*Epsilon(aDist1);
1891             BB.UpdateVertex(V1, aDist1);
1892           }
1893           if (k == NbSamples && !V2.IsNull())
1894           {
1895             Standard_Real aDist2 = Pnt2.Distance(aPoint);
1896             aDist2 += 2.*Epsilon(aDist2);
1897             BB.UpdateVertex(V2, aDist2);
1898           }
1899         }
1900       }
1901       BB.UpdateEdge(anEdge, MaxDist);
1902     }
1903     TolEdge = BRep_Tool::Tolerance(anEdge);
1904     if (!V1.IsNull())
1905     {
1906       gp_Pnt End1 = anHCurve->Value(fpar);
1907       Standard_Real dist1 = Pnt1.Distance(End1);
1908       dist1 += 2.*Epsilon(dist1);
1909       BB.UpdateVertex(V1, Max(dist1, TolEdge));
1910     }
1911     if (!V2.IsNull())
1912     {
1913       gp_Pnt End2 = anHCurve->Value(lpar);
1914       Standard_Real dist2 = Pnt2.Distance(End2);
1915       dist2 += 2.*Epsilon(dist2);
1916       BB.UpdateVertex(V2, Max(dist2, TolEdge));
1917     }
1918   }
1919 }
1920
1921 //=======================================================================
1922 //function : OrientClosedSolid
1923 //purpose  : 
1924 //=======================================================================
1925 Standard_Boolean BRepLib::OrientClosedSolid(TopoDS_Solid& solid) 
1926 {
1927   // Set material inside the solid
1928   BRepClass3d_SolidClassifier where(solid);
1929   where.PerformInfinitePoint(Precision::Confusion());
1930   if (where.State()==TopAbs_IN) {
1931     solid.Reverse();
1932   }
1933   else if (where.State()==TopAbs_ON || where.State()==TopAbs_UNKNOWN) 
1934     return Standard_False;
1935
1936   return Standard_True;
1937 }
1938
1939 // Structure for calculation of properties, necessary for decision about continuity
1940 class SurfaceProperties
1941 {
1942 public:
1943   SurfaceProperties(const Handle(Geom_Surface)& theSurface,
1944                     const gp_Trsf&              theSurfaceTrsf,
1945                     const Handle(Geom2d_Curve)& theCurve2D,
1946                     const Standard_Boolean      theReversed)
1947     : mySurfaceProps(theSurface, 2, Precision::Confusion()),
1948       mySurfaceTrsf(theSurfaceTrsf),
1949       myCurve2d(theCurve2D),
1950       myIsReversed(theReversed)
1951   {}
1952
1953   // Calculate derivatives on surface related to the point on curve
1954   void Calculate(const Standard_Real theParamOnCurve)
1955   {
1956     gp_Pnt2d aUV;
1957     myCurve2d->D1(theParamOnCurve, aUV, myCurveTangent);
1958     mySurfaceProps.SetParameters(aUV.X(), aUV.Y());
1959   }
1960
1961   // Returns point just calculated
1962   gp_Pnt Value() 
1963   { return mySurfaceProps.Value().Transformed(mySurfaceTrsf); }
1964
1965   // Calculate a derivative orthogonal to curve's tangent vector
1966   gp_Vec Derivative()
1967   {
1968     gp_Vec aDeriv;
1969     // direction orthogonal to tangent vector of the curve
1970     gp_Vec2d anOrtho(-myCurveTangent.Y(), myCurveTangent.X());
1971     Standard_Real aLen = anOrtho.Magnitude();
1972     if (aLen < Precision::Confusion())
1973       return aDeriv;
1974     anOrtho /= aLen;
1975     if (myIsReversed)
1976       anOrtho.Reverse();
1977
1978     aDeriv.SetLinearForm(anOrtho.X(), mySurfaceProps.D1U(),
1979                          anOrtho.Y(), mySurfaceProps.D1V());
1980     return aDeriv.Transformed(mySurfaceTrsf);
1981   }
1982
1983   // Calculate principal curvatures, which consist of minimal and maximal normal curvatures and
1984   // the directions on the tangent plane (principal direction) where the extremums are reached
1985   void Curvature(gp_Dir& thePrincipalDir1, Standard_Real& theCurvature1,
1986                  gp_Dir& thePrincipalDir2, Standard_Real& theCurvature2)
1987   {
1988     mySurfaceProps.CurvatureDirections(thePrincipalDir1, thePrincipalDir2);
1989     theCurvature1 = mySurfaceProps.MaxCurvature();
1990     theCurvature2 = mySurfaceProps.MinCurvature();
1991     if (myIsReversed)
1992     {
1993       theCurvature1 = -theCurvature1;
1994       theCurvature2 = -theCurvature2;
1995     }
1996     thePrincipalDir1.Transform(mySurfaceTrsf);
1997     thePrincipalDir2.Transform(mySurfaceTrsf);
1998   }
1999
2000 private:
2001   GeomLProp_SLProps    mySurfaceProps; // properties calculator
2002   gp_Trsf              mySurfaceTrsf;
2003   Handle(Geom2d_Curve) myCurve2d;
2004   Standard_Boolean     myIsReversed; // the face based on the surface is reversed
2005
2006   // tangent vector to Pcurve in UV
2007   gp_Vec2d myCurveTangent;
2008 };
2009
2010 //=======================================================================
2011 //function : tgtfaces
2012 //purpose  : check the angle at the border between two squares.
2013 //           Two shares should have a shared front edge.
2014 //=======================================================================
2015 static GeomAbs_Shape tgtfaces(const TopoDS_Edge& Ed,
2016                               const TopoDS_Face& F1,
2017                               const TopoDS_Face& F2,
2018                               const Standard_Real theAngleTol)
2019 {
2020   Standard_Boolean isSeam = F1.IsEqual(F2);
2021
2022   TopoDS_Edge E = Ed;
2023
2024   // Check if pcurves exist on both faces of edge
2025   Standard_Real aFirst,aLast;
2026   E.Orientation(TopAbs_FORWARD);
2027   Handle(Geom2d_Curve) aCurve1 = BRep_Tool::CurveOnSurface(E, F1, aFirst, aLast);
2028   if(aCurve1.IsNull())
2029     return GeomAbs_C0;
2030
2031   if (isSeam)
2032     E.Orientation(TopAbs_REVERSED);
2033   Handle(Geom2d_Curve) aCurve2 = BRep_Tool::CurveOnSurface(E, F2, aFirst, aLast);
2034   if(aCurve2.IsNull())
2035     return GeomAbs_C0;
2036
2037   TopLoc_Location aLoc1, aLoc2;
2038   Handle(Geom_Surface) aSurface1 = BRep_Tool::Surface(F1, aLoc1);
2039   const gp_Trsf& aSurf1Trsf = aLoc1.Transformation();
2040   Handle(Geom_Surface) aSurface2 = BRep_Tool::Surface(F2, aLoc2);
2041   const gp_Trsf& aSurf2Trsf = aLoc2.Transformation();
2042
2043   if (aSurface1->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
2044     aSurface1 = Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface1)->BasisSurface();
2045   if (aSurface2->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
2046     aSurface2 = Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface2)->BasisSurface();
2047
2048   // seam edge on elementary surface is always CN
2049   Standard_Boolean isElementary =
2050     (aSurface1->IsKind(STANDARD_TYPE(Geom_ElementarySurface)) &&
2051      aSurface2->IsKind(STANDARD_TYPE(Geom_ElementarySurface)));
2052   if (isSeam && isElementary)
2053   {
2054     return GeomAbs_CN;
2055   }
2056
2057   SurfaceProperties aSP1(aSurface1, aSurf1Trsf, aCurve1, F1.Orientation() == TopAbs_REVERSED);
2058   SurfaceProperties aSP2(aSurface2, aSurf2Trsf, aCurve2, F2.Orientation() == TopAbs_REVERSED);
2059
2060   Standard_Real f, l, eps;
2061   BRep_Tool::Range(E,f,l);
2062   Extrema_LocateExtPC ext;
2063   Handle(BRepAdaptor_HCurve) aHC2;
2064
2065   eps = (l - f)/100.;
2066   f += eps; // to avoid calculations on  
2067   l -= eps; // points of pointed squares.
2068
2069   const Standard_Real anAngleTol2 = theAngleTol * theAngleTol;
2070
2071   gp_Vec aDer1, aDer2;
2072   gp_Vec aNorm1;
2073   Standard_Real aSqLen1, aSqLen2;
2074   gp_Dir aCrvDir1[2], aCrvDir2[2];
2075   Standard_Real aCrvLen1[2], aCrvLen2[2];
2076
2077   GeomAbs_Shape aCont = (isElementary ? GeomAbs_CN : GeomAbs_C2);
2078   GeomAbs_Shape aCurCont;
2079   Standard_Real u;
2080   for (Standard_Integer i = 0; i <= 20 && aCont > GeomAbs_C0; i++)
2081   {
2082     // First suppose that this is sameParameter
2083     u = f + (l-f)*i/20;
2084
2085     // Check conditions for G1 and C1 continuity:
2086     // * calculate a derivative in tangent plane of each surface
2087     //   orthogonal to curve's tangent vector
2088     // * continuity is C1 if the vectors are equal
2089     // * continuity is G1 if the vectors are just parallel
2090     aCurCont = GeomAbs_C0;
2091
2092     aSP1.Calculate(u);
2093     aSP2.Calculate(u);
2094
2095     aDer1 = aSP1.Derivative();
2096     aSqLen1 = aDer1.SquareMagnitude();
2097     aDer2 = aSP2.Derivative();
2098     aSqLen2 = aDer2.SquareMagnitude();
2099     Standard_Boolean isSmoothSuspect = (aDer1.CrossSquareMagnitude(aDer2) <= anAngleTol2 * aSqLen1 * aSqLen2);
2100     if (!isSmoothSuspect)
2101     {
2102       // Refine by projection
2103       if (aHC2.IsNull())
2104       {
2105         // adaptor for pcurve on the second surface
2106         aHC2 = new BRepAdaptor_HCurve(BRepAdaptor_Curve(E, F2));
2107         ext.Initialize(aHC2->Curve(), f, l, Precision::PConfusion());
2108       }
2109       ext.Perform(aSP1.Value(), u);
2110       if (ext.IsDone() && ext.IsMin())
2111       {
2112         const Extrema_POnCurv& poc = ext.Point();
2113         aSP2.Calculate(poc.Parameter());
2114         aDer2 = aSP2.Derivative();
2115         aSqLen2 = aDer2.SquareMagnitude();
2116       }
2117       isSmoothSuspect = (aDer1.CrossSquareMagnitude(aDer2) <= anAngleTol2 * aSqLen1 * aSqLen2);
2118     }
2119     if (isSmoothSuspect)
2120     {
2121       aCurCont = GeomAbs_G1;
2122       if (Abs(Sqrt(aSqLen1) - Sqrt(aSqLen2)) < Precision::Confusion() &&
2123           aDer1.Dot(aDer2) > Precision::SquareConfusion()) // <= check vectors are codirectional
2124         aCurCont = GeomAbs_C1;
2125     }
2126     else
2127       return GeomAbs_C0;
2128
2129     if (aCont < GeomAbs_G2)
2130       continue; // no need further processing, because maximal continuity is less than G2
2131
2132     // Check conditions for G2 and C2 continuity:
2133     // * calculate principal curvatures on each surface
2134     // * continuity is C2 if directions of principal curvatures are equal on differenct surfaces
2135     // * continuity is G2 if directions of principal curvatures are just parallel
2136     //   and values of curvatures are the same
2137     aSP1.Curvature(aCrvDir1[0], aCrvLen1[0], aCrvDir1[1], aCrvLen1[1]);
2138     aSP2.Curvature(aCrvDir2[0], aCrvLen2[0], aCrvDir2[1], aCrvLen2[1]);
2139     for (Standard_Integer aStep = 0; aStep <= 1; ++aStep)
2140     {
2141       if (aCrvDir1[0].XYZ().CrossSquareMagnitude(aCrvDir2[aStep].XYZ()) <= Precision::SquareConfusion() &&
2142           Abs(aCrvLen1[0] - aCrvLen2[aStep]) < Precision::Confusion() &&
2143           aCrvDir1[1].XYZ().CrossSquareMagnitude(aCrvDir2[1 - aStep].XYZ()) <= Precision::SquareConfusion() &&
2144           Abs(aCrvLen1[1] - aCrvLen2[1 - aStep]) < Precision::Confusion())
2145       {
2146         if (aCurCont == GeomAbs_C1 &&
2147             aCrvDir1[0].Dot(aCrvDir2[aStep]) > Precision::Confusion() &&
2148             aCrvDir1[1].Dot(aCrvDir2[1 - aStep]) > Precision::Confusion())
2149           aCurCont = GeomAbs_C2;
2150         else
2151           aCurCont = GeomAbs_G2;
2152         break;
2153       }
2154     }
2155
2156     if (aCurCont < aCont)
2157       aCont = aCurCont;
2158   }
2159
2160   // according to the list of supported elementary surfaces,
2161   // if the continuity is C2, than it is totally CN
2162   if (isElementary && aCont == GeomAbs_C2)
2163     aCont = GeomAbs_CN;
2164   return aCont;
2165 }
2166
2167 //=======================================================================
2168 // function : EncodeRegularity
2169 // purpose  : Code the regularities on all edges of the shape, boundary of 
2170 //            two faces that do not have it.
2171 //            Takes into account that compound may consists of same solid
2172 //            placed with different transformations
2173 //=======================================================================
2174 static void EncodeRegularity(const TopoDS_Shape&        theShape,
2175                              const Standard_Real        theTolAng,
2176                              TopTools_MapOfShape&       theMap,
2177                              const TopTools_MapOfShape& theEdgesToEncode = TopTools_MapOfShape())
2178 {
2179   TopoDS_Shape aShape = theShape;
2180   TopLoc_Location aNullLoc;
2181   aShape.Location(aNullLoc); // nullify location
2182   if (!theMap.Add(aShape))
2183     return; // do not need to process shape twice
2184
2185   if (aShape.ShapeType() == TopAbs_COMPOUND ||
2186       aShape.ShapeType() == TopAbs_COMPSOLID)
2187   {
2188     for (TopoDS_Iterator it(aShape); it.More(); it.Next())
2189       EncodeRegularity(it.Value(), theTolAng, theMap, theEdgesToEncode);
2190     return;
2191   }
2192
2193   try {
2194     OCC_CATCH_SIGNALS
2195
2196     TopTools_IndexedDataMapOfShapeListOfShape M;
2197     TopExp::MapShapesAndAncestors(aShape, TopAbs_EDGE, TopAbs_FACE, M);
2198     TopTools_ListIteratorOfListOfShape It;
2199     TopExp_Explorer Ex;
2200     TopoDS_Face F1,F2;
2201     Standard_Boolean found;
2202     for (Standard_Integer i = 1; i <= M.Extent(); i++){
2203       TopoDS_Edge E = TopoDS::Edge(M.FindKey(i));
2204       if (!theEdgesToEncode.IsEmpty())
2205       {
2206         // process only the edges from the list to update their regularity
2207         TopoDS_Shape aPureEdge = E.Located(aNullLoc);
2208         aPureEdge.Orientation(TopAbs_FORWARD);
2209         if (!theEdgesToEncode.Contains(aPureEdge))
2210           continue;
2211       }
2212
2213       found = Standard_False;                                     
2214       F1.Nullify();
2215       for (It.Initialize(M.FindFromIndex(i)); It.More() && !found; It.Next()){
2216         if (F1.IsNull()) { F1 = TopoDS::Face(It.Value()); }
2217         else {
2218           const TopoDS_Face& aTmpF2 = TopoDS::Face(It.Value());
2219           if (!F1.IsSame(aTmpF2)){
2220             found = Standard_True;
2221             F2 = aTmpF2;
2222           }
2223         }
2224       }
2225       if (!found && !F1.IsNull()){//is it a sewing edge?
2226         TopAbs_Orientation orE = E.Orientation();
2227         TopoDS_Edge curE;
2228         for (Ex.Init(F1, TopAbs_EDGE); Ex.More() && !found; Ex.Next()){
2229           curE = TopoDS::Edge(Ex.Current());
2230           if (E.IsSame(curE) && orE != curE.Orientation()) {
2231             found = Standard_True;
2232             F2 = F1;
2233           }
2234         }
2235       }
2236       if (found)
2237         BRepLib::EncodeRegularity(E, F1, F2, theTolAng);
2238     }
2239   }
2240   catch (Standard_Failure const& anException) {
2241 #ifdef OCCT_DEBUG
2242     std::cout << "Warning: Exception in BRepLib::EncodeRegularity(): ";
2243     anException.Print(std::cout);
2244     std::cout << std::endl;
2245 #endif
2246     (void)anException;
2247   }
2248 }
2249
2250 //=======================================================================
2251 // function : EncodeRegularity
2252 // purpose  : code the regularities on all edges of the shape, boundary of 
2253 //            two faces that do not have it.
2254 //=======================================================================
2255
2256 void BRepLib::EncodeRegularity(const TopoDS_Shape& S,
2257   const Standard_Real TolAng)
2258 {
2259   TopTools_MapOfShape aMap;
2260   ::EncodeRegularity(S, TolAng, aMap);
2261 }
2262
2263 //=======================================================================
2264 // function : EncodeRegularity
2265 // purpose  : code the regularities on all edges in the list that do not 
2266 //            have it, and which are boundary of two faces on the shape.
2267 //=======================================================================
2268
2269 void BRepLib::EncodeRegularity(const TopoDS_Shape& S,
2270   const TopTools_ListOfShape& LE,
2271   const Standard_Real TolAng)
2272 {
2273   // Collect edges without location and orientation
2274   TopTools_MapOfShape aPureEdges;
2275   TopLoc_Location aNullLoc;
2276   TopTools_ListIteratorOfListOfShape anEdgeIt(LE);
2277   for (; anEdgeIt.More(); anEdgeIt.Next())
2278   {
2279     TopoDS_Shape anEdge = anEdgeIt.Value();
2280     anEdge.Location(aNullLoc);
2281     anEdge.Orientation(TopAbs_FORWARD);
2282     aPureEdges.Add(anEdge);
2283   }
2284
2285   TopTools_MapOfShape aMap;
2286   ::EncodeRegularity(S, TolAng, aMap, aPureEdges);
2287 }
2288
2289 //=======================================================================
2290 // function : EncodeRegularity
2291 // purpose  : code the regularity between 2 faces connected by edge 
2292 //=======================================================================
2293
2294 void BRepLib::EncodeRegularity(TopoDS_Edge& E,
2295   const TopoDS_Face& F1,
2296   const TopoDS_Face& F2,
2297   const Standard_Real TolAng)
2298 {
2299   BRep_Builder B;
2300   if(BRep_Tool::Continuity(E,F1,F2)<=GeomAbs_C0){
2301     try {
2302       GeomAbs_Shape aCont = tgtfaces(E, F1, F2, TolAng);
2303       B.Continuity(E,F1,F2,aCont);
2304       
2305     }
2306     catch(Standard_Failure const&)
2307     {
2308 #ifdef OCCT_DEBUG
2309       std::cout << "Failure: Exception in BRepLib::EncodeRegularity" << std::endl;
2310 #endif
2311     }
2312   }
2313 }
2314
2315 //=======================================================================
2316 // function : EnsureNormalConsistency
2317 // purpose  : Corrects the normals in Poly_Triangulation of faces.
2318 //              Returns TRUE if any correction is done.
2319 //=======================================================================
2320 Standard_Boolean BRepLib::
2321             EnsureNormalConsistency(const TopoDS_Shape& theShape,
2322                                     const Standard_Real theAngTol,
2323                                     const Standard_Boolean theForceComputeNormals)
2324 {
2325   const Standard_Real aThresDot = cos(theAngTol);
2326
2327   Standard_Boolean aRetVal = Standard_False, isNormalsFound = Standard_False;
2328
2329   // compute normals if they are absent
2330   TopExp_Explorer anExpFace(theShape,TopAbs_FACE);
2331   for (; anExpFace.More(); anExpFace.Next())
2332   {
2333     const TopoDS_Face& aFace = TopoDS::Face(anExpFace.Current());
2334     const Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
2335     if(aSurf.IsNull())
2336       continue;
2337     TopLoc_Location aLoc;
2338     const Handle(Poly_Triangulation)& aPT = BRep_Tool::Triangulation(aFace, aLoc);
2339     if(aPT.IsNull())
2340       continue;
2341     if (!theForceComputeNormals && aPT->HasNormals())
2342     {
2343       isNormalsFound = Standard_True;
2344       continue;
2345     }
2346
2347     GeomLProp_SLProps aSLP(aSurf, 2, Precision::Confusion());
2348     const Standard_Integer anArrDim = 3*aPT->NbNodes();
2349     Handle(TShort_HArray1OfShortReal) aNormArr = new TShort_HArray1OfShortReal(1, anArrDim);
2350     Standard_Integer anNormInd = aNormArr->Lower();
2351     for(Standard_Integer i = aPT->UVNodes().Lower(); i <= aPT->UVNodes().Upper(); i++)
2352     {
2353       const gp_Pnt2d &aP2d = aPT->UVNodes().Value(i);
2354       aSLP.SetParameters(aP2d.X(), aP2d.Y());
2355
2356       gp_XYZ aNorm(0.,0.,0.);
2357       if(!aSLP.IsNormalDefined())
2358       {
2359 #ifdef OCCT_DEBUG
2360         std::cout << "BRepLib::EnsureNormalConsistency(): Cannot find normal!" << std::endl;
2361 #endif
2362       }
2363       else
2364       {
2365         aNorm = aSLP.Normal().XYZ();
2366         if (aFace.Orientation() == TopAbs_REVERSED)
2367           aNorm.Reverse();
2368       }
2369       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.X());
2370       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.Y());
2371       aNormArr->ChangeValue(anNormInd++) = static_cast<Standard_ShortReal>(aNorm.Z());
2372     }
2373
2374     aRetVal = Standard_True;
2375     isNormalsFound = Standard_True;
2376     aPT->SetNormals(aNormArr);
2377   }
2378
2379   if(!isNormalsFound)
2380   {
2381     return aRetVal;
2382   }
2383
2384   // loop by edges
2385   TopTools_IndexedDataMapOfShapeListOfShape aMapEF;
2386   TopExp::MapShapesAndAncestors(theShape,TopAbs_EDGE,TopAbs_FACE,aMapEF);
2387   for(Standard_Integer anInd = 1; anInd <= aMapEF.Extent(); anInd++)
2388   {
2389     const TopoDS_Edge& anEdg = TopoDS::Edge(aMapEF.FindKey(anInd));
2390     const TopTools_ListOfShape& anEdgList = aMapEF.FindFromIndex(anInd);
2391     if (anEdgList.Extent() != 2)
2392       continue;
2393     TopTools_ListIteratorOfListOfShape anItF(anEdgList);
2394     const TopoDS_Face aFace1 = TopoDS::Face(anItF.Value());
2395     anItF.Next();
2396     const TopoDS_Face aFace2 = TopoDS::Face(anItF.Value());
2397     TopLoc_Location aLoc1, aLoc2;
2398     const Handle(Poly_Triangulation)& aPT1 = BRep_Tool::Triangulation(aFace1, aLoc1);
2399     const Handle(Poly_Triangulation)& aPT2 = BRep_Tool::Triangulation(aFace2, aLoc2);
2400
2401     if(aPT1.IsNull() || aPT2.IsNull())
2402       continue;
2403
2404     if(!aPT1->HasNormals() || !aPT2->HasNormals())
2405       continue;
2406
2407     const Handle(Poly_PolygonOnTriangulation)& aPTEF1 = 
2408                                 BRep_Tool::PolygonOnTriangulation(anEdg, aPT1, aLoc1);
2409     const Handle(Poly_PolygonOnTriangulation)& aPTEF2 = 
2410                                 BRep_Tool::PolygonOnTriangulation(anEdg, aPT2, aLoc2);
2411
2412     TShort_Array1OfShortReal& aNormArr1 = aPT1->ChangeNormals();
2413     TShort_Array1OfShortReal& aNormArr2 = aPT2->ChangeNormals();
2414
2415     if (aPTEF1->Nodes().Lower() != aPTEF2->Nodes().Lower() || 
2416         aPTEF1->Nodes().Upper() != aPTEF2->Nodes().Upper()) 
2417       continue; 
2418
2419     for(Standard_Integer anEdgNode = aPTEF1->Nodes().Lower();
2420                          anEdgNode <= aPTEF1->Nodes().Upper(); anEdgNode++)
2421     {
2422       //Number of node
2423       const Standard_Integer aFNodF1 = aPTEF1->Nodes().Value(anEdgNode);
2424       const Standard_Integer aFNodF2 = aPTEF2->Nodes().Value(anEdgNode);
2425
2426       const Standard_Integer aFNorm1FirstIndex = aNormArr1.Lower() + 3*
2427                                                     (aFNodF1 - aPT1->Nodes().Lower());
2428       const Standard_Integer aFNorm2FirstIndex = aNormArr2.Lower() + 3*
2429                                                     (aFNodF2 - aPT2->Nodes().Lower());
2430
2431       gp_XYZ aNorm1(aNormArr1.Value(aFNorm1FirstIndex),
2432                     aNormArr1.Value(aFNorm1FirstIndex+1),
2433                     aNormArr1.Value(aFNorm1FirstIndex+2));
2434       gp_XYZ aNorm2(aNormArr2.Value(aFNorm2FirstIndex),
2435                     aNormArr2.Value(aFNorm2FirstIndex+1),
2436                     aNormArr2.Value(aFNorm2FirstIndex+2));
2437       const Standard_Real aDot = aNorm1 * aNorm2;
2438
2439       if(aDot > aThresDot)
2440       {
2441         gp_XYZ aNewNorm = (aNorm1 + aNorm2).Normalized();
2442         aNormArr1.ChangeValue(aFNorm1FirstIndex) =
2443           aNormArr2.ChangeValue(aFNorm2FirstIndex) =
2444           static_cast<Standard_ShortReal>(aNewNorm.X());
2445         aNormArr1.ChangeValue(aFNorm1FirstIndex+1) =
2446           aNormArr2.ChangeValue(aFNorm2FirstIndex+1) =
2447           static_cast<Standard_ShortReal>(aNewNorm.Y());
2448         aNormArr1.ChangeValue(aFNorm1FirstIndex+2) =
2449           aNormArr2.ChangeValue(aFNorm2FirstIndex+2) =
2450           static_cast<Standard_ShortReal>(aNewNorm.Z());
2451         aRetVal = Standard_True;
2452       }
2453     }
2454   }
2455
2456   return aRetVal;
2457 }
2458
2459 //=======================================================================
2460 //function : SortFaces
2461 //purpose  : 
2462 //=======================================================================
2463
2464 void  BRepLib::SortFaces (const TopoDS_Shape& Sh,
2465   TopTools_ListOfShape& LF)
2466 {
2467   LF.Clear();
2468   TopTools_ListOfShape LTri,LPlan,LCyl,LCon,LSphere,LTor,LOther;
2469   TopExp_Explorer exp(Sh,TopAbs_FACE);
2470   TopLoc_Location l;
2471   Handle(Geom_Surface) S;
2472
2473   for (; exp.More(); exp.Next()) {
2474     const TopoDS_Face&   F = TopoDS::Face(exp.Current());
2475     S = BRep_Tool::Surface(F, l);
2476     if (!S.IsNull()) {
2477       if (S->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
2478         S = Handle(Geom_RectangularTrimmedSurface)::DownCast (S)->BasisSurface();
2479       }
2480       GeomAdaptor_Surface AS(S);
2481       switch (AS.GetType()) {
2482       case GeomAbs_Plane: 
2483         {
2484           LPlan.Append(F);
2485           break;
2486         }
2487       case GeomAbs_Cylinder: 
2488         {
2489           LCyl.Append(F);
2490           break;
2491         }
2492       case GeomAbs_Cone: 
2493         {
2494           LCon.Append(F);
2495           break;
2496         }
2497       case GeomAbs_Sphere: 
2498         {
2499           LSphere.Append(F);
2500           break;
2501         }
2502       case GeomAbs_Torus: 
2503         {
2504           LTor.Append(F);
2505           break;
2506         }
2507       default:
2508         LOther.Append(F);
2509       }
2510     }
2511     else LTri.Append(F);
2512   }
2513   LF.Append(LPlan); LF.Append(LCyl  ); LF.Append(LCon); LF.Append(LSphere);
2514   LF.Append(LTor ); LF.Append(LOther); LF.Append(LTri); 
2515 }
2516
2517 //=======================================================================
2518 //function : ReverseSortFaces
2519 //purpose  : 
2520 //=======================================================================
2521
2522 void  BRepLib::ReverseSortFaces (const TopoDS_Shape& Sh,
2523   TopTools_ListOfShape& LF)
2524 {
2525   LF.Clear();
2526   // Use the allocator of the result LF for intermediate results
2527   TopTools_ListOfShape LTri(LF.Allocator()), LPlan(LF.Allocator()),
2528     LCyl(LF.Allocator()), LCon(LF.Allocator()), LSphere(LF.Allocator()),
2529     LTor(LF.Allocator()), LOther(LF.Allocator());
2530   TopExp_Explorer exp(Sh,TopAbs_FACE);
2531   TopLoc_Location l;
2532
2533   for (; exp.More(); exp.Next()) {
2534     const TopoDS_Face&   F = TopoDS::Face(exp.Current());
2535     const Handle(Geom_Surface)& S = BRep_Tool::Surface(F, l);
2536     if (!S.IsNull()) {
2537       GeomAdaptor_Surface AS(S);
2538       switch (AS.GetType()) {
2539       case GeomAbs_Plane: 
2540         {
2541           LPlan.Append(F);
2542           break;
2543         }
2544       case GeomAbs_Cylinder: 
2545         {
2546           LCyl.Append(F);
2547           break;
2548         }
2549       case GeomAbs_Cone: 
2550         {
2551           LCon.Append(F);
2552           break;
2553         }
2554       case GeomAbs_Sphere: 
2555         {
2556           LSphere.Append(F);
2557           break;
2558         }
2559       case GeomAbs_Torus: 
2560         {
2561           LTor.Append(F);
2562           break;
2563         }
2564       default:
2565         LOther.Append(F);
2566       }
2567     }
2568     else LTri.Append(F);
2569   }
2570   LF.Append(LTri); LF.Append(LOther); LF.Append(LTor ); LF.Append(LSphere);
2571   LF.Append(LCon); LF.Append(LCyl  ); LF.Append(LPlan);
2572
2573 }
2574
2575 //=======================================================================
2576 // function: BoundingVertex
2577 // purpose : 
2578 //=======================================================================
2579 void BRepLib::BoundingVertex(const NCollection_List<TopoDS_Shape>& theLV,
2580                              gp_Pnt& theNewCenter, Standard_Real& theNewTol)
2581 {
2582   Standard_Integer aNb;
2583   //
2584   aNb=theLV.Extent();
2585   if (aNb < 2) {
2586     return;
2587   }
2588   //
2589   else if (aNb==2) {
2590     Standard_Integer m, n;
2591     Standard_Real aR[2], dR, aD, aEps;
2592     TopoDS_Vertex aV[2];
2593     gp_Pnt aP[2];
2594     //
2595     aEps=RealEpsilon();
2596     for (m=0; m<aNb; ++m) {
2597       aV[m]=(!m)? 
2598         *((TopoDS_Vertex*)(&theLV.First())):
2599         *((TopoDS_Vertex*)(&theLV.Last()));
2600       aP[m]=BRep_Tool::Pnt(aV[m]);
2601       aR[m]=BRep_Tool::Tolerance(aV[m]);
2602     }  
2603     //
2604     m=0; // max R
2605     n=1; // min R
2606     if (aR[0]<aR[1]) {
2607       m=1;
2608       n=0;
2609     }
2610     //
2611     dR=aR[m]-aR[n]; // dR >= 0.
2612     gp_Vec aVD(aP[m], aP[n]);
2613     aD=aVD.Magnitude();
2614     //
2615     if (aD<=dR || aD<aEps) { 
2616       theNewCenter = aP[m];
2617       theNewTol = aR[m];
2618     }
2619     else {
2620       Standard_Real aRr;
2621       gp_XYZ aXYZr;
2622       gp_Pnt aPr;
2623       //
2624       aRr=0.5*(aR[m]+aR[n]+aD);
2625       aXYZr=0.5*(aP[m].XYZ()+aP[n].XYZ()-aVD.XYZ()*(dR/aD));
2626       aPr.SetXYZ(aXYZr);
2627       //
2628       theNewCenter = aPr;
2629       theNewTol = aRr;
2630       //aBB.MakeVertex (aVnew, aPr, aRr);
2631     }
2632     return;
2633   }// else if (aNb==2) {
2634   //
2635   else { // if (aNb>2)
2636     // compute the point
2637     //
2638     // issue 0027540 - sum of doubles may depend on the order
2639     // of addition, thus sort the coordinates for stable result
2640     Standard_Integer i;
2641     NCollection_Array1<gp_Pnt> aPoints(0, aNb-1);
2642     NCollection_List<TopoDS_Edge>::Iterator aIt(theLV);
2643     for (i = 0; aIt.More(); aIt.Next(), ++i) {
2644       const TopoDS_Vertex& aVi = *((TopoDS_Vertex*)(&aIt.Value()));
2645       gp_Pnt aPi = BRep_Tool::Pnt(aVi);
2646       aPoints(i) = aPi;
2647     }
2648     //
2649     std::sort(aPoints.begin(), aPoints.end(), BRepLib_ComparePoints());
2650     //
2651     gp_XYZ aXYZ(0., 0., 0.);
2652     for (i = 0; i < aNb; ++i) {
2653       aXYZ += aPoints(i).XYZ();
2654     }
2655     aXYZ.Divide((Standard_Real)aNb);
2656     //
2657     gp_Pnt aP(aXYZ);
2658     //
2659     // compute the tolerance for the new vertex
2660     Standard_Real aTi, aDi, aDmax;
2661     //
2662     aDmax=-1.;
2663     aIt.Initialize(theLV);
2664     for (; aIt.More(); aIt.Next()) {
2665       TopoDS_Vertex& aVi=*((TopoDS_Vertex*)(&aIt.Value()));
2666       gp_Pnt aPi=BRep_Tool::Pnt(aVi);
2667       aTi=BRep_Tool::Tolerance(aVi);
2668       aDi=aP.SquareDistance(aPi);
2669       aDi=sqrt(aDi);
2670       aDi=aDi+aTi;
2671       if (aDi > aDmax) {
2672         aDmax=aDi;
2673       }
2674     }
2675     //
2676     theNewCenter = aP;
2677     theNewTol = aDmax;
2678   }
2679 }
2680
2681 //=======================================================================
2682 //function : ExtendFace
2683 //purpose  :
2684 //=======================================================================
2685 void BRepLib::ExtendFace(const TopoDS_Face& theF,
2686                          const Standard_Real theExtVal,
2687                          const Standard_Boolean theExtUMin,
2688                          const Standard_Boolean theExtUMax,
2689                          const Standard_Boolean theExtVMin,
2690                          const Standard_Boolean theExtVMax,
2691                          TopoDS_Face& theFExtended)
2692 {
2693   // Get face bounds
2694   BRepAdaptor_Surface aBAS(theF);
2695   Standard_Real aFUMin = aBAS.FirstUParameter(),
2696                 aFUMax = aBAS.LastUParameter(),
2697                 aFVMin = aBAS.FirstVParameter(),
2698                 aFVMax = aBAS.LastVParameter();
2699   const Standard_Real aTol = BRep_Tool::Tolerance(theF);
2700
2701   // Surface to build the face
2702   Handle(Geom_Surface) aS;
2703
2704   const GeomAbs_SurfaceType aType = aBAS.GetType();
2705   // treat analytical surfaces first
2706   if (aType == GeomAbs_Plane ||
2707       aType == GeomAbs_Sphere ||
2708       aType == GeomAbs_Cylinder ||
2709       aType == GeomAbs_Torus ||
2710       aType == GeomAbs_Cone)
2711   {
2712     // Get basis transformed basis surface
2713     Handle(Geom_Surface) aSurf = Handle(Geom_Surface)::
2714       DownCast(aBAS.Surface().Surface()->Transformed(aBAS.Trsf()));
2715
2716     // Get bounds of the basis surface
2717     Standard_Real aSUMin, aSUMax, aSVMin, aSVMax;
2718     aSurf->Bounds(aSUMin, aSUMax, aSVMin, aSVMax);
2719
2720     Standard_Boolean isUPeriodic = aBAS.IsUPeriodic();
2721     Standard_Real anUPeriod = isUPeriodic ? aBAS.UPeriod() : 0.0;
2722     if (isUPeriodic)
2723     {
2724       // Adjust face bounds to first period
2725       Standard_Real aDelta = aFUMax - aFUMin;
2726       aFUMin = Max(aSUMin, aFUMin + anUPeriod*Ceiling((aSUMin - aFUMin) / anUPeriod));
2727       aFUMax = aFUMin + aDelta;
2728     }
2729
2730     Standard_Boolean isVPeriodic = aBAS.IsVPeriodic();
2731     Standard_Real aVPeriod = isVPeriodic ? aBAS.VPeriod() : 0.0;
2732     if (isVPeriodic)
2733     {
2734       // Adjust face bounds to first period
2735       Standard_Real aDelta = aFVMax - aFVMin;
2736       aFVMin = Max(aSVMin, aFVMin + aVPeriod*Ceiling((aSVMin - aFVMin) / aVPeriod));
2737       aFVMax = aFVMin + aDelta;
2738     }
2739
2740     // Enlarge the face
2741     Standard_Real anURes = 0., aVRes = 0.;
2742     if (theExtUMin || theExtUMax)
2743       anURes = aBAS.UResolution(theExtVal);
2744     if (theExtVMin || theExtVMax)
2745       aVRes = aBAS.VResolution(theExtVal);
2746
2747     if (theExtUMin) aFUMin = Max(aSUMin, aFUMin - anURes);
2748     if (theExtUMax) aFUMax = Min(isUPeriodic ? aFUMin + anUPeriod : aSUMax, aFUMax + anURes);
2749     if (theExtVMin) aFVMin = Max(aSVMin, aFVMin - aVRes);
2750     if (theExtVMax) aFVMax = Min(isVPeriodic ? aFVMin + aVPeriod : aSVMax, aFVMax + aVRes);
2751
2752     // Check if the periodic surface should become closed.
2753     // In this case, use the basis surface with basis bounds.
2754     const Standard_Real anEps = Precision::PConfusion();
2755     if (isUPeriodic && Abs(aFUMax - aFUMin - anUPeriod) < anEps)
2756     {
2757       aFUMin = aSUMin;
2758       aFUMax = aSUMax;
2759     }
2760     if (isVPeriodic && Abs(aFVMax - aFVMin - aVPeriod) < anEps)
2761     {
2762       aFVMin = aSVMin;
2763       aFVMax = aSVMax;
2764     }
2765
2766     aS = aSurf;
2767   }
2768   else
2769   {
2770     // General case
2771
2772     Handle(Geom_BoundedSurface) aSB =
2773       Handle(Geom_BoundedSurface)::DownCast(BRep_Tool::Surface(theF));
2774     if (aSB.IsNull())
2775     {
2776       theFExtended = theF;
2777       return;
2778     }
2779
2780     // Get surfaces bounds
2781     Standard_Real aSUMin, aSUMax, aSVMin, aSVMax;
2782     aSB->Bounds(aSUMin, aSUMax, aSVMin, aSVMax);
2783
2784     Standard_Boolean isUClosed = aSB->IsUClosed();
2785     Standard_Boolean isVClosed = aSB->IsVClosed();
2786
2787     // Check if the extension in necessary directions is done
2788     Standard_Boolean isExtUMin = Standard_False,
2789                      isExtUMax = Standard_False,
2790                      isExtVMin = Standard_False,
2791                      isExtVMax = Standard_False;
2792
2793     // UMin
2794     if (theExtUMin && !isUClosed && !Precision::IsInfinite(aSUMin)) {
2795       GeomLib::ExtendSurfByLength(aSB, theExtVal, 1, Standard_True, Standard_False);
2796       isExtUMin = Standard_True;
2797     }
2798     // UMax
2799     if (theExtUMax && !isUClosed && !Precision::IsInfinite(aSUMax)) {
2800       GeomLib::ExtendSurfByLength(aSB, theExtVal, 1, Standard_True, Standard_True);
2801       isExtUMax = Standard_True;
2802     }
2803     // VMin
2804     if (theExtVMin && !isVClosed && !Precision::IsInfinite(aSVMax)) {
2805       GeomLib::ExtendSurfByLength(aSB, theExtVal, 1, Standard_False, Standard_False);
2806       isExtVMin = Standard_True;
2807     }
2808     // VMax
2809     if (theExtVMax && !isVClosed && !Precision::IsInfinite(aSVMax)) {
2810       GeomLib::ExtendSurfByLength(aSB, theExtVal, 1, Standard_False, Standard_True);
2811       isExtVMax = Standard_True;
2812     }
2813
2814     aS = aSB;
2815
2816     // Get new bounds
2817     aS->Bounds(aSUMin, aSUMax, aSVMin, aSVMax);
2818     if (isExtUMin) aFUMin = aSUMin;
2819     if (isExtUMax) aFUMax = aSUMax;
2820     if (isExtVMin) aFVMin = aSVMin;
2821     if (isExtVMax) aFVMax = aSVMax;
2822   }
2823
2824   BRepLib_MakeFace aMF(aS, aFUMin, aFUMax, aFVMin, aFVMax, aTol);
2825   theFExtended = *(TopoDS_Face*)&aMF.Shape();
2826   if (theF.Orientation() == TopAbs_REVERSED)
2827     theFExtended.Reverse();
2828 }