0027896: Prm-Prm intersection algo returns wrong result if it is called with start...
[occt.git] / src / Draft / Draft_Modification_1.cxx
1 // Created on: 1994-12-02
2 // Created by: Jacques GOUSSARD
3 // Copyright (c) 1994-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #include <Adaptor3d_CurveOnSurface.hxx>
18 #include <Adaptor3d_HCurveOnSurface.hxx>
19 #include <GeomAdaptor_SurfaceOfLinearExtrusion.hxx>
20 #include <Approx_CurveOnSurface.hxx>
21 #include <BRep_Builder.hxx>
22 #include <BRep_Tool.hxx>
23 #include <BRepExtrema_ExtCF.hxx>
24 #include <BRepExtrema_ExtPC.hxx>
25 #include <BRepLib_MakeFace.hxx>
26 #include <BRepTools.hxx>
27 #include <Draft_EdgeInfo.hxx>
28 #include <Draft_FaceInfo.hxx>
29 #include <Draft_Modification.hxx>
30 #include <Draft_VertexInfo.hxx>
31 #include <ElCLib.hxx>
32 #include <ElSLib.hxx>
33 #include <Extrema_ExtPC.hxx>
34 #include <Geom2d_BezierCurve.hxx>
35 #include <Geom2d_BSplineCurve.hxx>
36 #include <Geom2d_Curve.hxx>
37 #include <Geom2d_Line.hxx>
38 #include <Geom2d_TrimmedCurve.hxx>
39 #include <Geom2dAdaptor_HCurve.hxx>
40 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
41 #include <Geom2dConvert.hxx>
42 #include <Geom2dConvert_CompCurveToBSplineCurve.hxx>
43 #include <Geom_BSplineCurve.hxx>
44 #include <Geom_Circle.hxx>
45 #include <Geom_Conic.hxx>
46 #include <Geom_ConicalSurface.hxx>
47 #include <Geom_Curve.hxx>
48 #include <Geom_CylindricalSurface.hxx>
49 #include <Geom_ElementarySurface.hxx>
50 #include <Geom_Ellipse.hxx>
51 #include <Geom_Hyperbola.hxx>
52 #include <Geom_Line.hxx>
53 #include <Geom_Parabola.hxx>
54 #include <Geom_Plane.hxx>
55 #include <Geom_RectangularTrimmedSurface.hxx>
56 #include <Geom_Surface.hxx>
57 #include <Geom_SurfaceOfLinearExtrusion.hxx>
58 #include <Geom_TrimmedCurve.hxx>
59 #include <GeomAdaptor_Curve.hxx>
60 #include <GeomAdaptor_HCurve.hxx>
61 #include <GeomAdaptor_HSurface.hxx>
62 #include <GeomAdaptor_Surface.hxx>
63 #include <GeomAPI_ProjectPointOnCurve.hxx>
64 #include <GeomAPI_ProjectPointOnSurf.hxx>
65 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
66 #include <GeomInt_IntSS.hxx>
67 #include <GeomProjLib.hxx>
68 #include <gp.hxx>
69 #include <gp_Circ.hxx>
70 #include <gp_Dir.hxx>
71 #include <gp_Elips.hxx>
72 #include <gp_Hypr.hxx>
73 #include <gp_Lin.hxx>
74 #include <gp_Parab.hxx>
75 #include <gp_Pln.hxx>
76 #include <gp_Pnt.hxx>
77 #include <gp_Vec.hxx>
78 #include <IntAna_IntConicQuad.hxx>
79 #include <IntAna_QuadQuadGeo.hxx>
80 #include <IntCurveSurface_HInter.hxx>
81 #include <IntCurveSurface_IntersectionPoint.hxx>
82 #include <Precision.hxx>
83 #include <ProjLib_CompProjectedCurve.hxx>
84 #include <ProjLib_HCompProjectedCurve.hxx>
85 #include <Standard_ConstructionError.hxx>
86 #include <Standard_DomainError.hxx>
87 #include <Standard_Failure.hxx>
88 #include <Standard_NoSuchObject.hxx>
89 #include <Standard_NotImplemented.hxx>
90 #include <StdFail_NotDone.hxx>
91 #include <TColgp_Array1OfPnt2d.hxx>
92 #include <TopExp.hxx>
93 #include <TopExp_Explorer.hxx>
94 #include <TopLoc_Location.hxx>
95 #include <TopoDS.hxx>
96 #include <TopoDS_Edge.hxx>
97 #include <TopoDS_Face.hxx>
98 #include <TopoDS_Shape.hxx>
99 #include <TopoDS_Vertex.hxx>
100 #include <TopTools_ListIteratorOfListOfShape.hxx>
101 #include <TopTools_MapIteratorOfMapOfShape.hxx>
102 #include <TopTools_MapOfShape.hxx>
103
104 static Standard_Boolean Choose(const Draft_IndexedDataMapOfFaceFaceInfo&,
105   Draft_IndexedDataMapOfEdgeEdgeInfo&,
106   const TopoDS_Vertex&,
107   Draft_VertexInfo&,
108   GeomAdaptor_Curve&,
109   GeomAdaptor_Surface&);
110
111 static Standard_Real Parameter(const Handle(Geom_Curve)&,
112   const gp_Pnt&,
113   Standard_Integer&);
114
115 static Standard_Real SmartParameter(Draft_EdgeInfo&,
116   const Standard_Real EdgeTol,
117   const gp_Pnt&,
118   const Standard_Integer,
119   const Handle(Geom_Surface)&,
120   const Handle(Geom_Surface)&);
121
122 static TopAbs_Orientation Orientation(const TopoDS_Shape&,
123   const TopoDS_Face&);
124
125 static Standard_Boolean FindRotation(const gp_Pln&,
126   const TopAbs_Orientation,
127   const gp_Dir&,
128   const Standard_Real,
129   const gp_Pln&,
130   gp_Ax1&,
131   Standard_Real&);
132
133
134 //=======================================================================
135 //function : InternalAdd
136 //purpose  : 
137 //=======================================================================
138
139 Standard_Boolean Draft_Modification::InternalAdd(const TopoDS_Face& F,
140   const gp_Dir& Direction,
141   const Standard_Real Angle,
142   const gp_Pln& NeutralPlane,
143   const Standard_Boolean Flag)
144 {
145
146   if (myFMap.Contains(F)) {
147     return (badShape.IsNull());
148   }
149
150   TopAbs_Orientation oris = Orientation(myShape,F);
151   TopLoc_Location Lo;
152   //gp_Dir NewDirection = Direction;
153   //Standard_Real NewAngle = Angle;
154   Handle(Geom_Surface) S = BRep_Tool::Surface(F,Lo);
155   S = Handle(Geom_Surface)::DownCast(S->Transformed(Lo.Transformation()));
156   if (S->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
157     S = Handle(Geom_RectangularTrimmedSurface)::DownCast(S)->BasisSurface();
158   }
159   Handle(Geom_Surface) NewS;
160   Handle(Geom_Circle) theCircle;
161
162   Standard_Boolean postponed = (Flag == Standard_False);
163   if (postponed) {
164     Handle(Standard_Type) typS = S->DynamicType();
165     if (typS == STANDARD_TYPE(Geom_CylindricalSurface) ||
166         typS == STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion)) {
167       gp_Circ Cir;
168       if (typS == STANDARD_TYPE(Geom_CylindricalSurface)) {
169         gp_Cylinder cyl = 
170           Handle(Geom_CylindricalSurface)::DownCast(S)->Cylinder();
171         gp_Ax1 axcyl = cyl.Axis();
172         Cir = ElSLib::CylinderVIso( cyl.Position(), cyl.Radius(), 0.);
173         gp_Vec VV(cyl.Location(),NeutralPlane.Location());
174         Cir.Translate(VV.Dot(axcyl.Direction())*axcyl.Direction());
175       }
176       else {
177         Handle(Geom_Curve) Cbas = 
178           Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(S)->BasisCurve();
179         gp_Dir theDirextr = 
180           Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(S)->Direction();
181
182         if (Cbas->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
183           Cbas = Handle(Geom_TrimmedCurve)::DownCast(Cbas)->BasisCurve();
184         }
185         if (Cbas->IsKind(STANDARD_TYPE(Geom_Circle))) {
186           Cir = Handle(Geom_Circle)::DownCast(Cbas)->Circ();
187           gp_Dir dircir = Cir.Axis().Direction();
188           if (!Direction.IsParallel(dircir,Precision::Angular())) {
189             badShape = F;
190             errStat = Draft_FaceRecomputation;
191             return Standard_False;
192           }
193         }
194         else {
195           badShape = F;
196           errStat = Draft_FaceRecomputation;
197           return Standard_False;
198         }
199
200         gp_Ax3 Axis = NeutralPlane.Position();
201         Standard_Real L =
202           gp_Vec(Cir.Location(),Axis.Location()).
203           Dot(Axis.Direction());
204         Standard_Real Cos = theDirextr.Dot(Axis.Direction());
205         gp_Vec VV = ( L / Cos) * theDirextr;
206         Cir.Translate(VV);
207       }
208
209       theCircle = new Geom_Circle(Cir);
210
211     }
212     else {
213       postponed = Standard_False;
214     }
215   }
216
217
218   if (!postponed) {
219     NewS = NewSurface(S,oris,Direction,Angle,NeutralPlane); 
220     if (NewS.IsNull()) {
221       badShape = F;
222       errStat = Draft_FaceRecomputation;
223       return Standard_False;
224     }
225     // To avoid some problems with infinite restrictions
226     const Handle(Standard_Type)& typs = NewS->DynamicType();
227     if (typs == STANDARD_TYPE(Geom_CylindricalSurface) ||
228         typs == STANDARD_TYPE(Geom_ConicalSurface)) {
229       Standard_Real umin,umax,vmin,vmax;
230       BRepTools::UVBounds(F,umin,umax,vmin,vmax);
231       if (!Precision::IsNegativeInfinite(vmin) &&
232           !Precision::IsPositiveInfinite(vmax)) {
233         Standard_Real deltav = 10.*(vmax-vmin);
234         if(typs == STANDARD_TYPE(Geom_CylindricalSurface)) {
235           vmin = vmin - deltav;
236           vmax = vmax + deltav;
237         }
238         else {
239           gp_Cone Co = Handle(Geom_ConicalSurface)::DownCast(NewS)->Cone();
240           Standard_Real Vapex = - Co.RefRadius()/Sin(Co.SemiAngle());
241           if (vmin < Vapex) { // vmax should not exceed Vapex
242             if (vmax + deltav > Vapex) {
243               vmax = Vapex;
244               vmin = vmin - 10.*(vmax - vmin);
245               // JAG debug to avoid apex
246               vmax = vmax-Precision::Confusion();
247             }
248             else {
249               vmin = vmin - deltav;
250               vmax = vmax + deltav;
251             }
252           }
253           else { // Vapex <= vmin < vmax
254             if (vmin - deltav < Vapex) {
255               vmin = Vapex;
256               vmax = vmax + 10.*(vmax - vmin);
257               // JAG debug to avoid apex
258               vmin = vmin+Precision::Confusion();
259             }
260             else {
261               vmin = vmin - deltav;
262               vmax = vmax + deltav;
263             }
264           }
265         }
266         NewS = new Geom_RectangularTrimmedSurface(NewS,0.,2.*M_PI,vmin,vmax);
267       }
268     }
269   }
270
271   if (postponed || S != NewS) {
272     Draft_FaceInfo FI(NewS,Standard_True);
273     FI.RootFace(curFace);
274     myFMap.Add(F,FI);
275     if (postponed) {
276       myFMap.ChangeFromKey(F).ChangeCurve() = theCircle;
277     }
278   }    
279
280   TopExp_Explorer expl(F,TopAbs_EDGE);
281   TopTools_MapOfShape MapOfE;
282   while (expl.More() && badShape.IsNull()) {
283     const TopoDS_Edge& edg = TopoDS::Edge(expl.Current());
284     if (!myEMap.Contains(edg)) {
285       Standard_Boolean addedg  = Standard_False;
286       Standard_Boolean addface = Standard_False;
287       TopoDS_Face OtherF;
288       //if (BRep_Tool::IsClosed(edg,F)) {
289       if (BRepTools::IsReallyClosed(edg,F)) {
290         addedg = Standard_True;
291         addface = Standard_False;
292       }
293       else {
294         // Find the other face containing the edge.
295         TopTools_ListIteratorOfListOfShape it;
296         it.Initialize(myEFMap.FindFromKey(edg));
297         Standard_Integer nbother = 0;
298         while (it.More()) {
299           if (!it.Value().IsSame(F)) {
300             if (OtherF.IsNull()) {
301               OtherF = TopoDS::Face(it.Value());
302             }
303             nbother++;
304           }
305           it.Next();
306         }         
307         if (nbother >=2) {
308           badShape = edg;
309           errStat = Draft_EdgeRecomputation;
310         }
311         else if (! OtherF.IsNull() && 
312                  BRep_Tool::Continuity(edg,F,OtherF) >= GeomAbs_G1) {
313           addface= Standard_True;
314           addedg = Standard_True;
315         }
316         else if (nbother == 0) {
317           //        badShape = F;
318         }
319       }
320       if (addedg) {
321         if (postponed) {
322           myFMap.ChangeFromKey(F).Add(OtherF);
323         }
324         Standard_Real f,l;
325         TopLoc_Location L;
326         Handle(Geom_Curve) C = BRep_Tool::Curve(edg,L,f,l);
327         C = Handle(Geom_Curve)::DownCast(C->Transformed(L));
328         if (C->DynamicType() == STANDARD_TYPE(Geom_TrimmedCurve)) {
329           C = Handle(Geom_TrimmedCurve)::DownCast(C)->BasisCurve();
330         }
331         Handle(Geom_Curve) NewC;
332         Draft_EdgeInfo EInf(Standard_True);     
333         if(postponed) {
334           EInf.Add(F);
335           EInf.Add(OtherF);
336
337           // find fixed point 
338           Handle(Geom_Line) aLocalGeom = Handle(Geom_Line)::DownCast(C);
339           if (aLocalGeom.IsNull()) {
340             badShape = edg;
341             errStat = Draft_EdgeRecomputation;
342           }
343           else {
344             gp_Lin lin = aLocalGeom->Lin();
345             IntAna_IntConicQuad ilipl(lin,NeutralPlane,Precision::Angular());
346             if (ilipl.IsDone() && ilipl.NbPoints() != 0){
347               EInf.Tangent(ilipl.Point(1));
348             }
349             else {
350               badShape = edg;
351               errStat = Draft_EdgeRecomputation;
352             }
353           }
354         }
355         else {
356           NewC = NewCurve(C,S,oris,Direction,Angle,NeutralPlane, Flag);
357           if (NewC.IsNull()) {
358             badShape = edg;
359             errStat = Draft_EdgeRecomputation;
360           }
361         }
362
363         Handle(Geom_TrimmedCurve) T = Handle(Geom_TrimmedCurve)::DownCast(NewC);
364         if (!T.IsNull()) NewC = T->BasisCurve();
365         EInf.ChangeGeometry() = NewC;
366
367         EInf.RootFace(curFace);
368         myEMap.Add(edg,EInf);
369         MapOfE.Add(edg);
370         if (addface) {
371           Standard_Boolean Fl = Flag;
372           Handle(Geom_Surface) alocalSurface = BRep_Tool::Surface(OtherF,Lo);
373           if (alocalSurface->DynamicType() == 
374               STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
375             alocalSurface = Handle(Geom_RectangularTrimmedSurface)::
376               DownCast(alocalSurface)->BasisSurface();
377           }
378           Handle(Standard_Type) typS = alocalSurface->DynamicType();
379           if (typS == STANDARD_TYPE(Geom_CylindricalSurface) || 
380               typS == STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion)) {
381             if ( myFMap.Contains(F)) {
382               if ( Flag == Standard_False && !postponed) {
383                 Standard_Integer IndToReplace = myFMap.FindIndex(F);
384                 if (IndToReplace) {
385                   Standard_Integer LInd = myFMap.Extent();
386                   TopoDS_Face LF = myFMap.FindKey(LInd);
387                   Draft_FaceInfo LFInfo = myFMap.FindFromIndex(LInd);
388                   myFMap.RemoveLast();
389
390                   if (IndToReplace != LInd)
391                     myFMap.Substitute(IndToReplace, LF, LFInfo);
392                 }
393                 TopTools_MapIteratorOfMapOfShape itm(MapOfE);
394                 for ( ; itm.More(); itm.Next())
395                 {
396                   IndToReplace = myEMap.FindIndex(TopoDS::Edge(itm.Key()));
397                   if ( IndToReplace )
398                   {
399                     Standard_Integer LInd = myEMap.Extent();
400                     TopoDS_Edge LE = myEMap.FindKey(LInd);
401                     Draft_EdgeInfo LEInfo = myEMap.FindFromIndex(LInd);
402                     myEMap.RemoveLast();
403
404                     if (IndToReplace != LInd)
405                       myEMap.Substitute(IndToReplace, LE, LEInfo);
406                   }
407                 }
408               }
409             }
410           } 
411           InternalAdd(OtherF,Direction,Angle,NeutralPlane, Fl);
412         }
413       }
414     }
415     expl.Next();
416   }
417   return (badShape.IsNull());
418 }
419
420
421 //=======================================================================
422 //function : Propagate
423 //purpose  : 
424 //=======================================================================
425
426 Standard_Boolean Draft_Modification::Propagate ()
427 {
428
429   if (!badShape.IsNull()) return Standard_False;
430
431   // Set all edges and vertices of modified faces
432   TopoDS_Face F;
433   TopoDS_Edge E;
434   TopoDS_Vertex V;
435   TopExp_Explorer editer;
436   TopExp_Explorer vtiter;
437
438   for (Standard_Integer i = 1; i <= myFMap.Extent(); i++)
439   {
440     const TopoDS_Face& Fc = myFMap.FindKey(i);
441
442     // Exploration of the edges of the face
443     editer.Init(Fc,TopAbs_EDGE);
444     while (editer.More()) {
445       E = TopoDS::Edge(editer.Current());
446
447       if (!myEMap.Contains(E)) {
448         Draft_EdgeInfo EInf(Standard_True);     
449         myEMap.Add(E,EInf);
450       }
451       myEMap.ChangeFromKey(E).Add(Fc);
452
453       // Exploration of the vertices of the edge
454       vtiter.Init(E,TopAbs_VERTEX);
455       while (vtiter.More()) {
456         V = TopoDS::Vertex(vtiter.Current());
457         if (!myVMap.Contains(V)) {
458           Draft_VertexInfo VInf;
459           myVMap.Add(V,VInf);
460         }
461
462         myVMap.ChangeFromKey(V).Add(E);
463         myVMap.ChangeFromKey(V).ChangeParameter(E) = BRep_Tool::Parameter(V, E);
464         vtiter.Next();
465       }
466       editer.Next();
467     }
468   }
469
470
471   TopExp_Explorer anc;
472   Standard_Boolean found;
473
474   // Set edges containing modified vertices.
475
476   for (Standard_Integer i = 1; i <= myVMap.Extent(); i++)
477   {
478     const TopoDS_Vertex& Vt = myVMap.FindKey(i);
479
480     // Exploration of the ancestors of the vertex
481     anc.Init(myShape,TopAbs_EDGE);
482
483     while (anc.More()) {
484       E = TopoDS::Edge(anc.Current());
485       vtiter.Init(E,TopAbs_VERTEX);
486       found = Standard_False;
487       while (vtiter.More()) {
488         if (Vt.IsSame(TopoDS::Vertex(vtiter.Current()))) {
489           found = Standard_True;
490           break;
491         }
492         vtiter.Next();
493       }
494       if (found) {
495         if (!myEMap.Contains(E)) {
496           Draft_EdgeInfo EInf(Standard_False);
497           myEMap.Add(E,EInf);
498         }
499         myVMap.ChangeFromKey(Vt).Add(E);
500         myVMap.ChangeFromKey(Vt).ChangeParameter(E) = BRep_Tool::Parameter(Vt, E);
501       }
502       anc.Next();
503     }
504   }
505
506
507   // Set faces containing modified edges
508   for (Standard_Integer i = 1; i <= myEMap.Extent(); i++)
509   {
510     const TopoDS_Edge& Ed = myEMap.FindKey(i);
511     TopTools_ListIteratorOfListOfShape it;
512     for (it.Initialize(myEFMap.FindFromKey(Ed)); it.More(); it.Next()) {
513       F = TopoDS::Face(it.Value());
514       if (!myFMap.Contains(F)) {
515         TopLoc_Location L;
516         Handle(Geom_Surface) S = BRep_Tool::Surface(F,L);
517         Handle(Geom_Surface) NewS = 
518           Handle(Geom_Surface)::DownCast(S->Transformed(L.Transformation()));
519
520         const Handle(Standard_Type)& typs = S->DynamicType();
521         if (typs == STANDARD_TYPE(Geom_CylindricalSurface) ||
522             typs == STANDARD_TYPE(Geom_ConicalSurface)) {
523           Standard_Real umin,umax,vmin,vmax;
524           BRepTools::UVBounds(F,umin,umax,vmin,vmax);
525           if (!Precision::IsNegativeInfinite(vmin) &&
526               !Precision::IsPositiveInfinite(vmax)) {
527             Standard_Real deltav = 10.*(vmax-vmin);
528             vmin = vmin - deltav;
529             vmax = vmax + deltav;
530             NewS = new Geom_RectangularTrimmedSurface(NewS,0.,2.*M_PI,vmin,vmax);
531           }
532         }
533
534         Draft_FaceInfo FInf(NewS,Standard_False);
535         myFMap.Add(F,FInf);
536       }
537       myEMap.ChangeFromKey(Ed).Add(F);
538     }
539   }
540
541   //  Try to add faces for free borders...
542   // JAG 09.11.95
543   for (Standard_Integer i = 1; i <= myEMap.Extent(); i++)
544   {
545     Draft_EdgeInfo& Einf = myEMap.ChangeFromIndex(i);
546     if (Einf.NewGeometry() && 
547         Einf.Geometry().IsNull() &&     
548         Einf.SecondFace().IsNull()) {
549
550       TopLoc_Location Loc;
551       Handle(Geom_Surface) S1 = BRep_Tool::Surface(Einf.FirstFace(),Loc);
552       S1 = Handle(Geom_Surface)::
553         DownCast(S1->Transformed(Loc.Transformation()));
554       Handle(Geom_Surface) S2;
555
556       Standard_Real f,l;
557       const TopoDS_Edge& EK = myEMap.FindKey(i);
558       Handle(Geom_Curve) C = BRep_Tool::Curve(EK,Loc,f,l);
559       C = Handle(Geom_Curve)::DownCast(C->Transformed(Loc.Transformation()));
560       if (C->DynamicType() == STANDARD_TYPE(Geom_TrimmedCurve)) {
561         C = Handle(Geom_TrimmedCurve)::DownCast(C)->BasisCurve();
562       }
563       if (!S1->IsKind(STANDARD_TYPE(Geom_Plane))) {
564         if (C->IsKind(STANDARD_TYPE(Geom_Conic))) {
565           gp_Ax3 thePl(Handle(Geom_Conic)::DownCast(C)->Position());
566           S2 = new Geom_Plane(thePl);
567         }
568         else if (C->DynamicType() == STANDARD_TYPE(Geom_Line)) {
569           gp_Ax1 axis;
570           if (S1->DynamicType()== STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
571             axis = Handle(Geom_ElementarySurface)::DownCast
572               (Handle(Geom_RectangularTrimmedSurface)::DownCast(S1)->
573               BasisSurface())->Axis();
574           }
575           else {
576             axis = Handle(Geom_ElementarySurface)::DownCast(S1)->Axis();
577           }
578           gp_Vec they(axis.Location(), C->Value(0.));
579           gp_Dir axz(axis.Direction().Crossed(they));
580           S2=new Geom_Plane(gp_Ax3(axis.Location(),axz,axis.Direction()));
581
582         }
583         else {
584           badShape = EK;
585           errStat = Draft_EdgeRecomputation;
586           break; // leave from for
587         }
588       }
589       else { // on the plane
590         for (Standard_Integer j = 1; j <= myVMap.Extent(); j++)
591         {
592           Draft_VertexInfo& Vinf = myVMap.ChangeFromIndex(j);
593           for (Vinf.InitEdgeIterator();Vinf.MoreEdge();Vinf.NextEdge()) {
594             if (Vinf.Edge().IsSame(EK)) {
595               break;
596             }
597           }
598           if (Vinf.MoreEdge()) {
599             for (Vinf.InitEdgeIterator();Vinf.MoreEdge();Vinf.NextEdge()) {
600               const TopoDS_Edge& edg = Vinf.Edge();
601               if (!edg.IsSame(EK)) {
602                 const Draft_EdgeInfo& EI = myEMap.FindFromKey(edg);
603                 if (!EI.FirstFace().IsSame(Einf.FirstFace()) &&
604                   (EI.SecondFace().IsNull() || 
605                   !EI.SecondFace().IsSame(Einf.FirstFace()))) {
606                     break;
607                 }
608               }
609             }
610             if (Vinf.MoreEdge()) {
611               Handle(Geom_Curve) C2 = BRep_Tool::Curve(Vinf.Edge(), Loc,f,l);
612               Handle(GeomAdaptor_HCurve) HCur;
613               gp_Vec Direc;
614               C2 = Handle(Geom_Curve)::DownCast
615                 (C2->Transformed(Loc.Transformation()));
616               if (C2->DynamicType() == STANDARD_TYPE(Geom_TrimmedCurve)) {
617                 C2 = Handle(Geom_TrimmedCurve)::DownCast(C2)->BasisCurve();
618               }
619               if (C->DynamicType() == STANDARD_TYPE(Geom_Line)) {
620                 Direc = Handle(Geom_Line)::DownCast(C)->Lin().Direction();
621                 HCur = new GeomAdaptor_HCurve(C2);
622
623               }
624               else if (C2->DynamicType() == STANDARD_TYPE(Geom_Line)) {
625                 Direc = Handle(Geom_Line)::DownCast(C2)->Lin().Direction();
626                 HCur = new GeomAdaptor_HCurve(C);
627               }
628               else {
629                 badShape = EK;
630                 errStat = Draft_EdgeRecomputation;
631                 break; // leave from while
632               }
633               GeomAdaptor_SurfaceOfLinearExtrusion SLE(HCur,Direc);
634               switch(SLE.GetType()){
635
636               case GeomAbs_Plane :
637                 {
638                   S2 = new Geom_Plane(SLE.Plane());
639                 }
640                 break;
641               case GeomAbs_Cylinder :
642                 {
643                   S2 =   new Geom_CylindricalSurface(SLE.Cylinder());
644                 }
645                 break;
646               default :
647                 {
648                   S2 = new Geom_SurfaceOfLinearExtrusion(HCur->ChangeCurve().
649                     Curve(),
650                     Direc);
651                 }
652                 break;
653               }
654
655             }
656             else {
657               badShape = EK;
658               errStat = Draft_EdgeRecomputation;
659               break; // leave from while
660             }
661             break;
662           }
663           //j++;
664         }
665       }
666
667       if (badShape.IsNull()) {
668         BRep_Builder B;
669         TopoDS_Face TheNewFace;
670         B.MakeFace(TheNewFace,S2,Precision::Confusion());
671         Einf.Add(TheNewFace);
672         Draft_FaceInfo FI(S2,Standard_False);
673         myFMap.Add(TheNewFace,FI);
674       }
675       else {
676         break; // leave from for
677       }
678       // Fin JAG 09.11.95
679     }
680   }
681   return (badShape.IsNull());
682 }
683
684
685
686
687 //=======================================================================
688 //function : Perform
689 //purpose  : 
690 //=======================================================================
691
692 void Draft_Modification::Perform ()
693 {
694   if (!badShape.IsNull())  Standard_ConstructionError::Raise();
695
696   if (!myComp) {
697     myComp = Standard_True;
698     if (!Propagate()) {
699       return;
700     }
701
702     // Calculate eventual faces
703
704     for (Standard_Integer i = 1; i <= myFMap.Extent(); i++) 
705     {
706       const TopoDS_Face& FK = myFMap.FindKey(i);
707       Draft_FaceInfo& Finf = myFMap.ChangeFromIndex(i);
708       if (Finf.NewGeometry() && Finf.Geometry().IsNull()) {
709         const TopoDS_Face& F1 = Finf.FirstFace();
710         const TopoDS_Face& F2 = Finf.SecondFace();
711
712         if (F1.IsNull() || F2.IsNull()) {
713           errStat = Draft_FaceRecomputation;
714           badShape = FK;
715           return;
716         }
717         Handle(Geom_Surface) S1 = myFMap.FindFromKey(F1).Geometry();
718         Handle(Geom_Surface) S2 = myFMap.FindFromKey(F2).Geometry();
719         if (S1.IsNull() || S2.IsNull()) {
720           errStat = Draft_FaceRecomputation;
721           badShape = FK;
722           return;
723         }
724         if (S1->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
725           S1 = Handle(Geom_RectangularTrimmedSurface)::
726             DownCast(S1)->BasisSurface();
727         }
728         if (S2->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
729           S2 = Handle(Geom_RectangularTrimmedSurface)::
730             DownCast(S2)->BasisSurface();
731         }
732         Handle(Geom_Plane) P1 = Handle(Geom_Plane)::DownCast(S1);
733         Handle(Geom_Plane) P2 = Handle(Geom_Plane)::DownCast(S2);
734         if (P1.IsNull() || P2.IsNull()) {
735           errStat = Draft_FaceRecomputation;
736           badShape = FK;
737           return;
738         }
739         gp_Pln pp1 = P1->Pln();
740         gp_Pln pp2 = P2->Pln();
741         IntAna_QuadQuadGeo i2p(pp1,pp2,
742           Precision::Angular(),Precision::Confusion());
743         if (!i2p.IsDone() || i2p.TypeInter() != IntAna_Line) {
744           errStat = Draft_FaceRecomputation;
745           badShape = FK;
746           return;
747         }
748
749         gp_Dir extrdir = i2p.Line(1).Direction();
750
751         // Preserve the same direction as the base face
752         Handle(Geom_Surface) RefSurf = 
753           BRep_Tool::Surface(FK);
754         if (RefSurf->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
755           RefSurf = 
756             Handle(Geom_RectangularTrimmedSurface)::DownCast(RefSurf)
757             ->BasisSurface();
758         }
759         gp_Dir DirRef;
760
761         if ( RefSurf->DynamicType() == STANDARD_TYPE(Geom_CylindricalSurface)) {
762           gp_Ax3 AxeRef = 
763             Handle(Geom_CylindricalSurface)::DownCast(RefSurf)
764             ->Cylinder().Position();
765           DirRef = AxeRef.Direction();
766         }
767         else if (RefSurf->DynamicType() == 
768                  STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion)) {
769           DirRef = 
770             Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(RefSurf)->Direction();
771         }
772
773         if (extrdir.Dot(DirRef) < 0.) extrdir.Reverse();
774
775         // it is possible to accelerate speed by storing the info during
776         // InternalAdd --> modification of FaceInfo to preserve the circle
777
778         Handle(Geom_Circle) CCir = 
779           Handle(Geom_Circle)::DownCast(Finf.Curve());
780         Handle(Geom_Surface) NewS = 
781           new Geom_SurfaceOfLinearExtrusion(CCir, extrdir);    
782
783         Standard_Real umin, umax, vmin, vmax;
784         BRepTools::UVBounds(FK,umin,umax,vmin,vmax);
785         if (!Precision::IsNegativeInfinite(vmin) &&
786             !Precision::IsPositiveInfinite(vmax)) {
787           Standard_Real deltav = 2.*(vmax-vmin);
788           vmin = vmin - deltav;
789           vmax = vmax + deltav;
790         }
791
792         // very temporary
793         else {
794           vmax = 300;
795           vmin = -300;
796         }
797
798         NewS = new Geom_RectangularTrimmedSurface(NewS,0.,1.9*M_PI,vmin,vmax);
799         Finf.ChangeGeometry() = NewS;
800       }
801     }
802
803     // Calculate new edges.
804     for (Standard_Integer ii = 1; ii <= myEMap.Extent(); ii++) 
805     {
806       Draft_EdgeInfo& Einf = myEMap.ChangeFromIndex(ii); 
807
808       const TopoDS_Edge& theEdge = TopoDS::Edge(myEMap.FindKey(ii));
809
810       Handle(Geom_Surface) S1,S2;
811       Handle(Geom_Curve) C, newC;
812       Standard_Real f,l;
813       TopLoc_Location L;
814       C = BRep_Tool::Curve(theEdge,L,f,l);
815       C = Handle(Geom_Curve)::DownCast(C->Transformed(L.Transformation()));
816
817       if (Einf.NewGeometry() && Einf.Geometry().IsNull()) {
818         gp_Pnt ptfixe;
819         if (!Einf.IsTangent(ptfixe)) {
820           const TopoDS_Face& FirstFace = Einf.FirstFace();
821           const TopoDS_Face& SecondFace = Einf.SecondFace();
822
823           S1 = myFMap.FindFromKey(FirstFace).Geometry();
824           S2 = myFMap.FindFromKey(SecondFace).Geometry();
825
826           Standard_Integer detrompeur = 0;
827
828           // Return FirstVertex and the tangent at this point.
829           TopoDS_Vertex FV = TopExp::FirstVertex(theEdge);
830           TopoDS_Vertex LV = TopExp::LastVertex(theEdge);
831           Standard_Real pmin = 0.;
832           Standard_Real prmfv = BRep_Tool::Parameter(FV, theEdge);
833           Standard_Real prmlv = BRep_Tool::Parameter(LV, theEdge);
834           gp_Pnt pfv, plv;
835           gp_Vec d1fv,d1lv, newd1;
836           C->D1(prmfv,pfv,d1fv);
837           C->D1(prmlv,plv,d1lv);
838
839           Standard_Real TolF1 = BRep_Tool::Tolerance (FirstFace);
840           Standard_Real TolF2 = BRep_Tool::Tolerance (SecondFace);
841
842           //Pass the tolerance of the face to project
843           GeomAPI_ProjectPointOnSurf proj1 (pfv, S1, TolF1);
844           GeomAPI_ProjectPointOnSurf proj2 (plv, S1, TolF1);
845           GeomAPI_ProjectPointOnSurf proj3 (pfv, S2, TolF2);
846           GeomAPI_ProjectPointOnSurf proj4 (plv, S2, TolF2);
847
848           if (proj1.IsDone () && proj2.IsDone ()) {
849             if(proj1.LowerDistance()<= Precision::Confusion() &&
850                proj2.LowerDistance()<= Precision::Confusion())  {
851                detrompeur = 1;
852             }
853           }
854
855           if (proj3.IsDone () && proj4.IsDone ()) {
856             if(proj3.LowerDistance() <= Precision::Confusion() &&
857                proj4.LowerDistance() <= Precision::Confusion())  {
858                detrompeur = 2;
859             }
860           }
861
862           gp_Dir TheDirExtr;
863           gp_Ax3 Axis;
864           Handle(Geom_Curve) TheNewCurve;
865           Standard_Boolean KPart = Standard_False;
866
867           if ( S1->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
868             S1 = Handle(Geom_RectangularTrimmedSurface)::
869               DownCast(S1)->BasisSurface();
870           }
871           if ( S2->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
872             S2 = Handle(Geom_RectangularTrimmedSurface)::
873               DownCast(S2)->BasisSurface();
874           }
875
876           Standard_Boolean PC1 = Standard_True; // KPart on S1
877           if (S1->DynamicType() == STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion) &&
878               S2->DynamicType() == STANDARD_TYPE(Geom_Plane) ) {
879             KPart = Standard_True;
880             Axis = Handle(Geom_Plane)::DownCast(S2)->Position();
881             TheNewCurve = Handle(Geom_SurfaceOfLinearExtrusion)::
882               DownCast(S1)->BasisCurve();
883             TheDirExtr = Handle(Geom_SurfaceOfLinearExtrusion)::
884               DownCast(S1)->Direction();
885           }
886           else if (S2->DynamicType() == STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion) &&
887                    S1->DynamicType() == STANDARD_TYPE(Geom_Plane) ) {
888             KPart = Standard_True;
889             PC1 = Standard_False;
890             Axis = Handle(Geom_Plane)::DownCast(S1)->Position();
891             TheNewCurve = Handle(Geom_SurfaceOfLinearExtrusion)::
892               DownCast(S2)->BasisCurve();
893             TheDirExtr = Handle(Geom_SurfaceOfLinearExtrusion)::
894               DownCast(S2)->Direction();
895           }
896           Handle(Geom_Circle) aCirc ;
897           if ( KPart) {  // very temporary on circles !!!
898             aCirc = Handle(Geom_Circle)::DownCast(TheNewCurve);
899             if (aCirc.IsNull())
900               KPart = Standard_False;
901             else
902             {
903               gp_Dir AxofCirc = aCirc->Position().Direction();
904               if (AxofCirc.IsParallel(Axis.Direction(),Precision::Angular()))
905                 KPart = Standard_True;
906               else
907                 KPart = Standard_False;
908             }
909           }
910
911           Standard_Integer imin;
912           GeomInt_IntSS i2s;
913           if ( KPart) {
914             // direct calculation of NewC
915             Standard_Real aLocalReal =
916               gp_Vec(aCirc->Circ().Location(),Axis.Location()).
917               Dot(Axis.Direction());
918             Standard_Real Cos = TheDirExtr.Dot(Axis.Direction());
919             gp_Vec VV = ( aLocalReal / Cos) * TheDirExtr;
920             newC = Handle(Geom_Curve)::DownCast(TheNewCurve->Translated(VV));
921             // it is possible to calculate PCurve
922             Handle(Geom2d_Line) L2d 
923               = new Geom2d_Line(gp_Pnt2d(0.,aLocalReal/Cos),
924               gp::DX2d());
925
926             if ( PC1) 
927               Einf.ChangeFirstPC() = L2d;
928             else
929               Einf.ChangeSecondPC() = L2d;
930           }
931           else {
932             S1 = myFMap.FindFromKey(Einf.FirstFace()).Geometry();
933             S2 = myFMap.FindFromKey(Einf.SecondFace()).Geometry();
934
935
936             // PCurves are not calculated immediately for 2 reasons:
937             // 1 - If ProjLib should make an Approx, it is stupid to approximate the 
938             //     entire intersection curve.
939             // 2 - Additionally, if YaRev, there is a risk to not be SameRange.
940             i2s.Perform(S1,S2,Precision::Confusion(),
941               Standard_True,Standard_False,Standard_False);
942
943             if (!i2s.IsDone() || i2s.NbLines() <= 0) {
944               errStat = Draft_EdgeRecomputation;
945               badShape = theEdge;
946               return;
947             }
948
949             Standard_Real Glob2Min = RealLast();
950             GeomAdaptor_Curve TheCurve;
951
952             Standard_Integer i,j; //,jmin;
953
954             if (i2s.Line(1)->DynamicType() != STANDARD_TYPE(Geom_BSplineCurve))
955             {
956               Standard_Real Dist2Min = RealLast();
957               imin = 0;
958               for (i=1; i<= i2s.NbLines(); i++) {
959                 TheCurve.Load(i2s.Line(i));
960                 Extrema_ExtPC myExtPC(pfv,TheCurve);
961
962                 Standard_Real locpmin = 0.;
963                 if (myExtPC.IsDone()) {
964                   if(myExtPC.NbExt() >= 1) {
965                     Dist2Min = myExtPC.SquareDistance(1);
966                     locpmin = myExtPC.Point(1).Parameter();
967                   }
968                   if(myExtPC.NbExt() == 2 && Dist2Min > Precision::SquareConfusion()) {
969                     //to avoid incorrectly choosing the image 
970                     //of the first vertex of the initial edge
971                     Standard_Real d1_2 = myExtPC.SquareDistance(1);
972                     Standard_Real d2_2 = myExtPC.SquareDistance(2);
973                     if(d1_2 > 1.21*d2_2) {
974                       Dist2Min = myExtPC.SquareDistance(2);
975                       locpmin = myExtPC.Point(2).Parameter();
976                     }
977                     else if(d2_2 > 1.21*d1_2) {
978                       Dist2Min = myExtPC.SquareDistance(1);
979                       locpmin = myExtPC.Point(1).Parameter();
980                     }
981                     else {
982                       Standard_Real pfvpar = myExtPC.Point(1).Parameter();
983                       Standard_Real plvpar = myExtPC.Point(2).Parameter();
984                       newC = i2s.Line(i);
985
986                       gp_Pnt pfvprim, plvprim;
987
988                       newC->D0(pfvpar,pfvprim);
989                       newC->D0(plvpar,plvprim);
990
991                       Handle(Geom_Surface) theSurf;
992                       if(detrompeur == 1) {
993                         if(S1->DynamicType() == 
994                            STANDARD_TYPE(Geom_RectangularTrimmedSurface))  
995                           S1 = Handle(Geom_RectangularTrimmedSurface)::
996                           DownCast(S1)->BasisSurface();
997                         theSurf = S1;
998
999                       }
1000                       else if(detrompeur == 2) {
1001                         if(S2->DynamicType() == 
1002                            STANDARD_TYPE(Geom_RectangularTrimmedSurface)) 
1003                           S2 = Handle(Geom_RectangularTrimmedSurface)::
1004                           DownCast(S2)->BasisSurface();
1005                         theSurf = S2;               
1006                       }
1007                       if(detrompeur != 0 && detrompeur != 4) {
1008                         Standard_Real ul = 0., vl = 0., uf = 0., vf = 0.;
1009                         Standard_Real ufprim = 0., ulprim = 0., vfprim = 0., vlprim = 0.;
1010
1011                         if(theSurf->DynamicType() == STANDARD_TYPE(Geom_Plane)) {       
1012                           gp_Pln pl = Handle(Geom_Plane)::DownCast(S2)->Pln();
1013                           ElSLib::Parameters(pl, plv, ul, vl);
1014                           ElSLib::Parameters(pl, pfv, uf, vf);
1015                           ElSLib::Parameters(pl, plvprim, ulprim, vlprim);
1016                           ElSLib::Parameters(pl, pfvprim, ufprim, vfprim);
1017                         }
1018                         else if(theSurf->DynamicType() == 
1019                                 STANDARD_TYPE(Geom_CylindricalSurface)) {
1020                           gp_Cylinder cy = Handle(Geom_CylindricalSurface)
1021                             ::DownCast(S2)->Cylinder();
1022                           ElSLib::Parameters(cy, plv, ul, vl);
1023                           ElSLib::Parameters(cy, pfv, uf, vf);
1024                           ElSLib::Parameters(cy, plvprim, ulprim, vlprim);
1025                           ElSLib::Parameters(cy, pfvprim, ufprim, vfprim);
1026                         }
1027                         else detrompeur = 4;
1028
1029                         if(detrompeur == 1 || detrompeur == 2) {
1030                           gp_Vec2d v1((ul-ufprim), (vl-vfprim));
1031                           gp_Vec2d norm((vf-vfprim), (ufprim-uf)); 
1032                           gp_Vec2d v2((ulprim-ufprim), (vlprim-vfprim));
1033                           if( (v1.Dot(norm))*(v2.Dot(norm)) < 0) {
1034                             Dist2Min = myExtPC.SquareDistance(2);
1035                             locpmin = myExtPC.Point(2).Parameter();
1036                           }
1037                         }
1038                       }
1039                     }
1040                   }
1041                   if (myExtPC.NbExt() == 1 || myExtPC.NbExt() > 2 || detrompeur ==4) {
1042                     Dist2Min = myExtPC.SquareDistance(1);
1043                     locpmin = myExtPC.Point(1).Parameter();
1044                     for (j=2; j<=myExtPC.NbExt(); j++) {
1045                       const Standard_Real Dist2 = myExtPC.SquareDistance(j);
1046                       if (Dist2 < Dist2Min) {
1047                         Dist2Min = Dist2;
1048                         locpmin = myExtPC.Point(j).Parameter();
1049                       }
1050                     }
1051                   }
1052                   else if(myExtPC.NbExt() < 1){
1053                     Standard_Real dist1_2,dist2_2;
1054                     gp_Pnt p1b,p2b;
1055                     myExtPC.TrimmedSquareDistances(dist1_2,dist2_2,p1b,p2b);
1056                     if (dist1_2 < dist2_2) {
1057                       Dist2Min = dist1_2;
1058                       locpmin = TheCurve.FirstParameter();
1059                     }
1060                     else {
1061                       Dist2Min = dist2_2;
1062                       locpmin = TheCurve.LastParameter();
1063                     }
1064                   }
1065
1066                   if (Dist2Min  < Glob2Min) {
1067                     Glob2Min = Dist2Min;
1068                     imin = i;
1069                     pmin = locpmin;
1070                   }
1071                 }
1072               }
1073               if (imin == 0) {
1074                 errStat = Draft_EdgeRecomputation;
1075                 badShape = theEdge;
1076                 return;
1077               }
1078
1079               newC = i2s.Line(imin);
1080
1081               newC->D1(pmin,pfv,newd1);
1082               Standard_Boolean YaRev = d1fv.Dot(newd1) <0.; 
1083
1084               if (YaRev)
1085                 newC->Reverse();
1086
1087               if (i2s.HasLineOnS1(imin)) {
1088                 Einf.ChangeFirstPC() = i2s.LineOnS1(imin);
1089                 if ( YaRev) 
1090                   Einf.ChangeFirstPC()->Reverse();
1091               }
1092
1093               if (i2s.HasLineOnS2(imin)) {
1094                 Einf.ChangeSecondPC() = i2s.LineOnS2(imin);
1095                 if ( YaRev) 
1096                   Einf.ChangeSecondPC()->Reverse();
1097               }
1098             } // if (i2s.Line(1)->DynamicType() != STANDARD_TYPE(Geom_BSplineCurve))
1099             else // i2s.Line(1) is BSplineCurve
1100             {
1101               //Find the first curve to glue
1102               TColGeom_SequenceOfCurve Candidates;
1103               if (S1->DynamicType() == STANDARD_TYPE(Geom_CylindricalSurface) ||
1104                   S1->DynamicType() == STANDARD_TYPE(Geom_ConicalSurface)) 
1105               {
1106                 for (i = 1; i <= i2s.NbLines(); i++)
1107                 {
1108                   Handle( Geom_Curve ) aCurve = i2s.Line(i);
1109                   gp_Pnt Pnt = aCurve->Value( aCurve->FirstParameter() );
1110                   GeomAPI_ProjectPointOnSurf projector( Pnt, S1, Precision::Confusion() );
1111                   Standard_Real U, V;
1112                   projector.LowerDistanceParameters( U, V );
1113                   if (Abs(U) <= Precision::Confusion() || Abs(U-2.*M_PI) <= Precision::Confusion())
1114                     Candidates.Append( aCurve );
1115                   else
1116                   {
1117                     Pnt = aCurve->Value( aCurve->LastParameter() );
1118                     projector.Init( Pnt, S1, Precision::Confusion() );
1119                     projector.LowerDistanceParameters( U, V );
1120                     if (Abs(U) <= Precision::Confusion() || Abs(U-2.*M_PI) <= Precision::Confusion())
1121                     {
1122                       aCurve->Reverse();
1123                       Candidates.Append( aCurve );
1124                     }
1125                   }
1126                 }
1127
1128                 if(Candidates.Length() == 0) 
1129                 {
1130                   //errStat = Draft_EdgeRecomputation;
1131                   //badShape = TopoDS::Edge(ite.Key());
1132                   //return;
1133                   for (i = 1; i <= i2s.NbLines(); i++)
1134                     Candidates.Append( i2s.Line(i) );
1135                 }
1136               }
1137               else 
1138               {
1139                 for (i = 1; i <= i2s.NbLines(); i++)
1140                   Candidates.Append( i2s.Line(i) );
1141               }
1142
1143               Handle( Geom_Curve ) FirstCurve;
1144               if (Candidates.Length() > 1)
1145               {
1146                 Standard_Real DistMin = Precision::Infinite();
1147                 for (i = 1; i <= Candidates.Length(); i++)
1148                 {
1149                   Handle( Geom_Curve ) aCurve = Candidates(i);
1150                   gp_Pnt Pnt = aCurve->Value( aCurve->FirstParameter() );
1151                   const Standard_Real Dist = Pnt.Distance( pfv );
1152                   if (Dist - DistMin < -Precision::Confusion())
1153                   {
1154                     DistMin = Dist;
1155                     FirstCurve = aCurve;
1156                   }
1157                 }
1158               }
1159               else
1160                 FirstCurve = Candidates(1);
1161
1162               //Glueing
1163               TColGeom_SequenceOfCurve Curves;
1164               for (i = 1; i <= i2s.NbLines(); i++)
1165                 if (FirstCurve != i2s.Line(i))
1166                   Curves.Append( i2s.Line(i) );
1167
1168               TColGeom_SequenceOfCurve ToGlue;
1169               gp_Pnt EndPoint = FirstCurve->Value( FirstCurve->LastParameter() );
1170               Standard_Boolean added = Standard_True;
1171               while (added)
1172               {
1173                 added = Standard_False;
1174                 for (i = 1; i <= Curves.Length(); i++)
1175                 {
1176                   Handle( Geom_Curve ) aCurve = Curves(i);
1177                   gp_Pnt pfirst, plast;
1178                   pfirst = aCurve->Value( aCurve->FirstParameter() );
1179                   plast = aCurve->Value( aCurve->LastParameter() );
1180                   if (pfirst.Distance( EndPoint ) <= Precision::Confusion())
1181                   {
1182                     ToGlue.Append( aCurve );
1183                     EndPoint = plast;
1184                     Curves.Remove(i);
1185                     added = Standard_True;
1186                     break;
1187                   }
1188                   if (plast.Distance( EndPoint ) <= Precision::Confusion())
1189                   {
1190                     aCurve->Reverse();
1191                     ToGlue.Append( aCurve );
1192                     EndPoint = pfirst;
1193                     Curves.Remove(i);
1194                     added = Standard_True;
1195                     break;
1196                   }
1197                 }
1198               }
1199
1200               if (FirstCurve.IsNull()) {
1201                 errStat = Draft_EdgeRecomputation;
1202                 badShape = theEdge;
1203                 return;
1204               }
1205
1206               GeomConvert_CompCurveToBSplineCurve Concat( Handle(Geom_BSplineCurve)::DownCast(FirstCurve) );
1207               for (i = 1; i <= ToGlue.Length(); i++)
1208                 Concat.Add( Handle(Geom_BSplineCurve)::DownCast(ToGlue(i)), Precision::Confusion(), Standard_True );
1209
1210               newC = Concat.BSplineCurve();
1211
1212               TheCurve.Load( newC );
1213               Extrema_ExtPC myExtPC( pfv, TheCurve );
1214               Standard_Real Dist2Min = RealLast();
1215               for (i = 1; i <= myExtPC.NbExt(); i++)
1216               {
1217                 if (myExtPC.IsMin(i))
1218                 {
1219                   const Standard_Real Dist2 = myExtPC.SquareDistance(i);
1220                   if (Dist2 < Dist2Min)
1221                   {
1222                     Dist2Min = Dist2;
1223                     pmin = myExtPC.Point(i).Parameter();
1224                   }
1225                 }
1226               }
1227               newC->D1(pmin,pfv,newd1);
1228               Standard_Boolean YaRev = d1fv.Dot(newd1) < 0.; 
1229
1230               if (YaRev)
1231                 newC->Reverse();
1232               /*
1233               if (i2s.HasLineOnS1(imin)) {
1234               Einf.ChangeFirstPC() = i2s.LineOnS1(imin);
1235               if ( YaRev) 
1236               Einf.ChangeFirstPC()->Reverse();
1237               }
1238
1239               if (i2s.HasLineOnS2(imin)) {
1240               Einf.ChangeSecondPC() = i2s.LineOnS2(imin);
1241               if ( YaRev) 
1242               Einf.ChangeSecondPC()->Reverse();
1243               }
1244               */
1245             } // else: i2s.NbLines() > 2 && S1 is Cylinder or Cone
1246
1247             Einf.Tolerance(Max(Einf.Tolerance(), i2s.TolReached3d()));
1248           }  // End step KPart
1249         }
1250         else { // case of tangency
1251           const TopoDS_Face& F1 = Einf.FirstFace();
1252           const TopoDS_Face& F2 = Einf.SecondFace();
1253
1254           Handle(Geom_Surface) aLocalS1 = myFMap.FindFromKey(F1).Geometry();
1255           Handle(Geom_Surface) aLocalS2 = myFMap.FindFromKey(F2).Geometry();
1256           if (aLocalS1.IsNull() || aLocalS2.IsNull()) {
1257             errStat = Draft_EdgeRecomputation;
1258             badShape = theEdge;
1259             return;
1260           }
1261           if (aLocalS1->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
1262             aLocalS1 = Handle(Geom_RectangularTrimmedSurface)::
1263               DownCast(aLocalS1)->BasisSurface();
1264           }
1265           if (aLocalS2->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
1266             aLocalS2 = Handle(Geom_RectangularTrimmedSurface)::
1267               DownCast(aLocalS2)->BasisSurface();
1268           }
1269
1270           gp_Dir dirextr;
1271           //Standard_Boolean dirfound = Standard_False;
1272           if (aLocalS1->DynamicType() == STANDARD_TYPE(Geom_CylindricalSurface)) {
1273             gp_Cylinder cyl = 
1274               Handle(Geom_CylindricalSurface)::DownCast(aLocalS1)->Cylinder();
1275             dirextr = cyl.Axis().Direction();
1276             //dirfound = Standard_True;
1277             // see direction...
1278
1279           }
1280           else if (aLocalS1->DynamicType() == STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion)) {
1281             dirextr = Handle(Geom_SurfaceOfLinearExtrusion)::
1282               DownCast(aLocalS1)->Direction();
1283             //dirfound = Standard_True;
1284             // see direction...
1285
1286             // Here it is possible to calculate PCurve.
1287             Handle(Geom_SurfaceOfLinearExtrusion) SEL = 
1288               Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(aLocalS1);
1289             Handle(Geom_Circle) GCir = 
1290               Handle(Geom_Circle)::DownCast(SEL->BasisCurve());
1291             if ( !GCir.IsNull()) {
1292               Standard_Real U = ElCLib::Parameter(GCir->Circ(),ptfixe);
1293               Handle(Geom2d_Line) PC1 = 
1294                 new Geom2d_Line(gp_Pnt2d(U,0.),gp::DY2d());
1295               Einf.ChangeFirstPC() = PC1;
1296             }
1297           }
1298
1299           else if (aLocalS2->DynamicType() == STANDARD_TYPE(Geom_CylindricalSurface)) {
1300             gp_Cylinder cyl = 
1301               Handle(Geom_CylindricalSurface)::DownCast(aLocalS2)->Cylinder();
1302             dirextr = cyl.Axis().Direction();
1303             // dirfound = Standard_True;
1304             // see direction...
1305
1306           }
1307           else if (aLocalS2->DynamicType() == STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion)) {
1308             dirextr = Handle(Geom_SurfaceOfLinearExtrusion)::
1309               DownCast(aLocalS2)->Direction();
1310             // dirfound = Standard_True;
1311             // see direction...
1312
1313             // Here it is possible to calculate PCurve.
1314             Handle(Geom_SurfaceOfLinearExtrusion) SEL = 
1315               Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(aLocalS2);
1316             Handle(Geom_Circle) GCir = 
1317               Handle(Geom_Circle)::DownCast(SEL->BasisCurve());
1318             if ( !GCir.IsNull()) {
1319               Standard_Real U = ElCLib::Parameter(GCir->Circ(),ptfixe);
1320               Handle(Geom2d_Line) PC2 = 
1321                 new Geom2d_Line(gp_Pnt2d(U,0.),gp::DY2d());
1322               Einf.ChangeSecondPC() = PC2;
1323             }
1324           }
1325           newC = new Geom_Line(ptfixe,dirextr);
1326
1327           gp_Pnt pfv;
1328           gp_Vec d1fv,newd1;
1329           C->D1(0.,pfv,d1fv);
1330           newC->D1(0.,pfv,newd1);
1331           Standard_Boolean YaRev = d1fv.Dot(newd1) <0.; 
1332           if (YaRev) {
1333             newC->Reverse();
1334             if(!Einf.FirstPC().IsNull()) {
1335               Einf.ChangeFirstPC()->Reverse();
1336             }
1337             if(!Einf.SecondPC().IsNull()) {
1338               Einf.ChangeSecondPC()->Reverse();
1339             }
1340           }
1341         }
1342
1343         Handle(Geom_TrimmedCurve) T = Handle(Geom_TrimmedCurve)::DownCast(newC);
1344         if (!T.IsNull()) newC = T->BasisCurve();
1345         Einf.ChangeGeometry() = newC;
1346       }
1347       else if (!Einf.NewGeometry()){
1348         // set existing curve 3D
1349         Handle(Geom_TrimmedCurve) T = Handle(Geom_TrimmedCurve)::DownCast(C);
1350         if (!T.IsNull()) C = T->BasisCurve();
1351         Einf.ChangeGeometry() = C;
1352       }
1353     }
1354
1355     // Calculate new vertices.
1356
1357     Handle(GeomAdaptor_HCurve)   HAC = new GeomAdaptor_HCurve;
1358     Handle(GeomAdaptor_HSurface) HAS = new GeomAdaptor_HSurface;
1359
1360     for (Standard_Integer ii = 1; ii <= myVMap.Extent(); ii++)
1361     {
1362       GeomAdaptor_Curve    AC;
1363       GeomAdaptor_Surface  AS;
1364
1365       const TopoDS_Vertex& TVV = myVMap.FindKey(ii);
1366       Draft_VertexInfo& Vinf = myVMap.ChangeFromIndex(ii);
1367       if (!Choose(myFMap,myEMap,TVV,Vinf,AC,AS)) {
1368
1369         // no concerted edge => alignment of two consecutive edges.
1370         gp_Pnt pvt;
1371         Vinf.ChangeGeometry() = pvt;
1372         Vinf.InitEdgeIterator();
1373         if (Vinf.MoreEdge()) {
1374           const TopoDS_Edge& Edg1 = Vinf.Edge();
1375           //const Draft_EdgeInfo& Einf1 = myEMap(Edg1);
1376           Draft_EdgeInfo& Einf1 = myEMap.ChangeFromKey(Edg1);
1377           gp_Pnt vtori = BRep_Tool::Pnt(TVV);
1378           //Einf1.Geometry()->D0(Vinf.Parameter(Edg1), pvt);
1379           GeomAPI_ProjectPointOnCurve Projector( vtori, Einf1.Geometry() ); //patch
1380           pvt = Projector.NearestPoint();
1381
1382 #ifdef OCCT_DEBUG
1383           static Standard_Integer VertexRecomp = 1;
1384           if (VertexRecomp!=0) {
1385             cout << "pori :" << vtori.X() << " " << vtori.Y() << " " << vtori.Z() << endl;
1386             cout << "  Edg 1 :" << Vinf.Parameter(Edg1) << endl;
1387             cout << "pvt :" << pvt.X() << " " << pvt.Y() << " " << pvt.Z() << endl;
1388           }
1389 #endif
1390
1391           Standard_Real dion=pvt.SquareDistance(vtori);
1392           Vinf.NextEdge();
1393           if (Vinf.MoreEdge()) {
1394             const TopoDS_Edge& Edg2 = Vinf.Edge();
1395             //const Draft_EdgeInfo& Einf2 = myEMap(Edg2);
1396             Draft_EdgeInfo& Einf2 = myEMap.ChangeFromKey(Edg2);
1397             //      Standard_Real f;
1398             gp_Pnt opvt;
1399             Einf2.Geometry()->D0(Vinf.Parameter(Edg2), opvt);
1400
1401 #ifdef OCCT_DEBUG
1402             if (VertexRecomp!=0) {
1403               cout << "  Edg 2 :" << Vinf.Parameter(Vinf.Edge()) << endl;
1404               cout << "opvt " << opvt.X() << " " << opvt.Y() << " " << opvt.Z() << endl;
1405             }
1406 #endif
1407
1408             if (opvt.SquareDistance(vtori) < dion) {
1409               pvt = opvt;
1410             }
1411             //Vinf.ChangeParameter(Edg2) = Parameter(Einf2.Geometry(), pvt);
1412             Standard_Integer done;
1413             Standard_Real param = Parameter(Einf2.Geometry(), pvt, done);
1414             if (done != 0)
1415             {
1416               Handle(Geom_Surface) S1 = myFMap.FindFromKey(Einf2.FirstFace()).Geometry();
1417               Handle(Geom_Surface) S2 = myFMap.FindFromKey(Einf2.SecondFace()).Geometry();
1418               Vinf.ChangeParameter(Edg2) = SmartParameter( Einf2, BRep_Tool::Tolerance(Edg2), pvt, done, S1, S2 );
1419             }
1420             else
1421               Vinf.ChangeParameter(Edg2) = param;
1422           }
1423
1424           Vinf.ChangeGeometry() = pvt;
1425           //Vinf.ChangeParameter(Edg1) = Parameter(Einf1.Geometry(), pvt);
1426           Standard_Integer done;
1427           Standard_Real param = Parameter(Einf1.Geometry(), pvt, done);
1428           if (done != 0)
1429           {
1430             Handle(Geom_Surface) S1 = myFMap.FindFromKey(Einf1.FirstFace()).Geometry();
1431             Handle(Geom_Surface) S2 = myFMap.FindFromKey(Einf1.SecondFace()).Geometry();
1432             Vinf.ChangeParameter(Edg1) = SmartParameter( Einf1, BRep_Tool::Tolerance(Edg1), pvt, done, S1, S2 );
1433           }
1434           else
1435             Vinf.ChangeParameter(Edg1) = param;
1436           continue;
1437         }
1438
1439
1440         errStat = Draft_VertexRecomputation;
1441         badShape = TVV;
1442         return;
1443       }
1444
1445       IntCurveSurface_HInter myintcs;
1446       HAC->Set(AC);
1447       HAS->Set(AS);
1448
1449       myintcs.Perform(HAC,HAS);
1450
1451
1452       if (!myintcs.IsDone()) {
1453         errStat = Draft_VertexRecomputation;
1454         badShape = TVV;
1455         return;
1456       }
1457
1458       gp_Pnt vtori = BRep_Tool::Pnt(TVV);
1459       gp_Pnt pvt;
1460
1461       Standard_Integer nbsol = myintcs.NbPoints();
1462       if (nbsol <= 0)
1463       {
1464         Extrema_ExtCS extr( AC, AS, Precision::PConfusion(), Precision::PConfusion() );
1465
1466         if(!extr.IsDone() || extr.NbExt() == 0) {
1467           errStat = Draft_VertexRecomputation;
1468           badShape = TVV;
1469           return;
1470         }
1471
1472
1473         Standard_Real disref = RealLast();
1474         Standard_Integer iref = 0;
1475         Extrema_POnCurv Pc;
1476         Extrema_POnSurf Ps;
1477         for (Standard_Integer i = 1; i <= extr.NbExt(); i++)
1478         {
1479           extr.Points( i, Pc, Ps );
1480           Standard_Real distemp = Pc.Value().SquareDistance(vtori);
1481           if ( distemp < disref)
1482           {
1483             disref = distemp;
1484             iref = i;
1485           }
1486         } 
1487         extr.Points( iref, Pc, Ps );
1488         pvt = Pc.Value();
1489       }
1490       else
1491       {
1492         Standard_Real disref = RealLast();
1493         Standard_Integer iref = 0;
1494         for (Standard_Integer i = 1; i <= nbsol; i++)
1495         {
1496           Standard_Real distemp = myintcs.Point(i).Pnt().SquareDistance(vtori);
1497           if ( distemp < disref)
1498           {
1499             disref = distemp;
1500             iref = i;
1501           }
1502         } 
1503         pvt = myintcs.Point(iref).Pnt();
1504       }
1505
1506       Vinf.ChangeGeometry() = pvt;
1507
1508       for (Vinf.InitEdgeIterator();Vinf.MoreEdge(); Vinf.NextEdge()) {
1509         const TopoDS_Edge& Edg = Vinf.Edge();
1510         Standard_Real initpar = Vinf.Parameter(Edg);
1511         //const Draft_EdgeInfo& Einf = myEMap(Edg);
1512         Draft_EdgeInfo& Einf = myEMap.ChangeFromKey(Edg);
1513         //Vinf.ChangeParameter(Edg) = Parameter(Einf.Geometry(),pvt);
1514         Standard_Integer done;
1515         Standard_Real param = Parameter(Einf.Geometry(), pvt, done);
1516         if (done != 0)
1517         {
1518           Handle(Geom_Surface) S1 = myFMap.FindFromKey(Einf.FirstFace()).Geometry();
1519           Handle(Geom_Surface) S2 = myFMap.FindFromKey(Einf.SecondFace()).Geometry();
1520           Vinf.ChangeParameter(Edg) = SmartParameter( Einf, BRep_Tool::Tolerance(Edg), pvt, done, S1, S2 );
1521         }
1522         else
1523         {
1524           if(Abs(initpar - param) > Precision::PConfusion())
1525           {
1526             Standard_Real f, l;
1527             TopLoc_Location Loc;
1528             const Handle(Geom_Curve)& aC = BRep_Tool::Curve(Edg, Loc, f, l);
1529             if(aC->DynamicType() == STANDARD_TYPE(Geom_TrimmedCurve))
1530             {
1531               Einf.SetNewGeometry(Standard_True);
1532             }
1533           }
1534           Vinf.ChangeParameter(Edg) = param;
1535         }
1536       }
1537     }
1538   }
1539
1540   // small loop of validation/protection
1541
1542   for (Standard_Integer i = 1; i <= myEMap.Extent(); i++) 
1543   {
1544     const TopoDS_Edge& edg = TopoDS::Edge(myEMap.FindKey(i));
1545
1546     TopoDS_Vertex Vf,Vl;
1547     TopExp::Vertices(edg,Vf,Vl);
1548     if (edg.Orientation() == TopAbs_REVERSED) {
1549       Vf.Reverse();
1550       Vl.Reverse();
1551     }
1552
1553     if(myVMap.Contains(Vf) && myVMap.Contains(Vl))
1554     {
1555       //Here, we compare directions of the source edge (from input shape)
1556       //and corresponding selected part of the intersection edge.
1557       //If these directions are opposite then we reverse intersection edge
1558       //and recompute corresponding vertex-parameters.
1559
1560       Standard_Real aParF = myVMap.ChangeFromKey(Vf).Parameter(edg);
1561       Standard_Real aParL = myVMap.ChangeFromKey(Vl).Parameter(edg);
1562
1563       if(aParL < aParF)
1564       {
1565         Draft_EdgeInfo& aEinf = myEMap.ChangeFromKey(edg);
1566         TopLoc_Location aLoc;
1567         Standard_Real aF = 0.0, aL = 0.0;
1568         const Handle(Geom_Curve) aSCurve = BRep_Tool::Curve(edg, aF, aL);
1569         Handle(Geom_Curve) anIntCurv = aEinf.Geometry();
1570         gp_Pnt aPf, aPl;
1571         gp_Vec aDirNF, aDirNL, aDirOF, aDirOL;
1572         aSCurve->D1(BRep_Tool::Parameter(Vf, edg), aPf, aDirOF);
1573         aSCurve->D1(BRep_Tool::Parameter(Vl, edg), aPl, aDirOL);
1574         anIntCurv->D1(aParF, aPf, aDirNF);
1575         anIntCurv->D1(aParL, aPl, aDirNL);
1576
1577         Standard_Real aSqMagn = aDirNF.SquareMagnitude();
1578
1579         if (aSqMagn > Precision::SquareConfusion())
1580           aDirNF.Divide(sqrt(aSqMagn));
1581
1582         aSqMagn = aDirNL.SquareMagnitude();
1583         if (aSqMagn > Precision::SquareConfusion())
1584           aDirNL.Divide(sqrt(aSqMagn));
1585
1586         aSqMagn = aDirOF.SquareMagnitude();
1587         if (aSqMagn > Precision::SquareConfusion())
1588           aDirOF.Divide(sqrt(aSqMagn));
1589
1590         aSqMagn = aDirOL.SquareMagnitude();
1591         if (aSqMagn > Precision::SquareConfusion())
1592           aDirOL.Divide(sqrt(aSqMagn));
1593
1594         const Standard_Real aCosF = aDirNF.Dot(aDirOF), aCosL = aDirNL.Dot(aDirOL);
1595         const Standard_Real aCosMax = Abs(aCosF) > Abs(aCosL) ? aCosF : aCosL;
1596
1597         if(aCosMax < 0.0)
1598         {
1599           Standard_Integer anErr = 0;
1600           anIntCurv->Reverse();
1601           aEinf.ChangeGeometry() = anIntCurv;
1602           Standard_Real aPar = Parameter(aEinf.Geometry(), aPf, anErr);
1603           if(anErr == 0)
1604           {
1605             myVMap.ChangeFromKey(Vf).ChangeParameter(edg) = aPar;
1606           }
1607           aPar = Parameter(aEinf.Geometry(), aPl, anErr);
1608           if(anErr == 0)
1609           {
1610             myVMap.ChangeFromKey(Vl).ChangeParameter(edg) = aPar;
1611           }
1612         }
1613       }
1614     }
1615
1616     Standard_Real pf,pl,tolerance;
1617     if (!NewParameter(Vf,edg,pf,tolerance)) {
1618       pf = BRep_Tool::Parameter(Vf,edg);
1619     }
1620     if (!NewParameter(Vl,edg,pl,tolerance)) {
1621       pl = BRep_Tool::Parameter(Vl,edg);
1622     }
1623     if (pl <= pf) {
1624       //      const Handle(Geom_Curve) gc=ite.Value().Geometry();
1625       //      if (!gc.IsNull()) {
1626       //        pl = gc->LastParameter();
1627       //        pf = gc->FirstParameter();
1628       //      }
1629       Handle( Geom_Curve ) theCurve = myEMap.FindFromKey(edg).Geometry();
1630       if (theCurve->IsClosed())
1631       {
1632         // pf >= pl
1633         Standard_Real FirstPar = theCurve->FirstParameter(), LastPar = theCurve->LastParameter();
1634         Standard_Real pconf = Precision::PConfusion();
1635         if (Abs( pf - LastPar ) <= pconf)
1636           pf = FirstPar;
1637         else if (Abs( pl - FirstPar ) <= pconf)
1638           pl = LastPar;
1639
1640         if(pl <= pf) {
1641           pl += (LastPar-FirstPar);
1642         }
1643
1644       }
1645       if (pl <= pf) {
1646         errStat = Draft_EdgeRecomputation;
1647         badShape = edg;
1648         return;
1649       }
1650     }
1651     if (myVMap.Contains( Vf ))
1652       myVMap.ChangeFromKey(Vf).ChangeParameter(edg) = pf;
1653     if (myVMap.Contains( Vl ))
1654       myVMap.ChangeFromKey(Vl).ChangeParameter(edg) = pl;
1655   }
1656 }
1657
1658
1659
1660 //=======================================================================
1661 //function : NewSurface
1662 //purpose  : 
1663 //=======================================================================
1664
1665 Handle(Geom_Surface) Draft_Modification::NewSurface
1666   (const Handle(Geom_Surface)& S,
1667   const TopAbs_Orientation Oris,
1668   const gp_Dir& Direction,
1669   const Standard_Real Angle,
1670   const gp_Pln& NeutralPlane) 
1671 {
1672   Handle(Geom_Surface) NewS;
1673
1674   Handle(Standard_Type) TypeS = S->DynamicType();
1675
1676   if (TypeS == STANDARD_TYPE(Geom_Plane)) {
1677     gp_Pln Pl = Handle(Geom_Plane)::DownCast(S)->Pln();
1678     gp_Ax1 Axe;
1679     Standard_Real Theta;
1680     if (FindRotation(Pl,Oris,Direction,Angle,NeutralPlane,Axe,Theta)) {
1681       if ( Abs(Theta) > Precision::Angular()) {
1682         NewS = Handle(Geom_Surface)::DownCast(S->Rotated(Axe,Theta));
1683       }
1684       else {
1685         NewS = S;
1686       }
1687     }
1688   }
1689   else if (TypeS == STANDARD_TYPE(Geom_CylindricalSurface)) {
1690     Standard_Real testdir = Direction.Dot(NeutralPlane.Axis().Direction());
1691     if (Abs(testdir) <= 1.-Precision::Angular()) {      
1692 #ifdef OCCT_DEBUG
1693       cout << "NewSurfaceCyl:Draft_Direction_and_Neutral_Perpendicular" << endl;
1694 #endif
1695       return NewS;      
1696     }     
1697     gp_Cylinder Cy = Handle(Geom_CylindricalSurface)::DownCast(S)->Cylinder();     
1698     testdir = Direction.Dot(Cy.Axis().Direction());
1699     if (Abs(testdir) <= 1.-Precision::Angular()) {
1700 #ifdef OCCT_DEBUG
1701       cout << "NewSurfaceCyl:Draft_Direction_and_Cylinder_Perpendicular" << endl;
1702 #endif
1703       return NewS;
1704     }
1705     if (Abs(Angle) > Precision::Angular())
1706     {
1707       IntAna_QuadQuadGeo i2s;
1708       i2s.Perform(NeutralPlane,Cy,Precision::Angular(),Precision::Confusion());
1709       Standard_Boolean isIntDone = i2s.IsDone();
1710
1711       if(i2s.TypeInter() == IntAna_Ellipse)
1712       {
1713         const gp_Elips anEl = i2s.Ellipse(1);
1714         const Standard_Real aMajorR = anEl.MajorRadius();
1715         const Standard_Real aMinorR = anEl.MinorRadius();
1716         isIntDone = (aMajorR < 100000.0 * aMinorR);
1717       }
1718
1719       if (!isIntDone || i2s.TypeInter() != IntAna_Circle) {
1720 #ifdef OCCT_DEBUG
1721         cout << "NewSurfaceCyl:Draft_Intersection_Neutral_Cylinder_NotDone" << endl;
1722 #endif
1723         return NewS;
1724       } 
1725       gp_Ax3 axcone = Cy.Position();
1726       // Pb : Where is the material???
1727       Standard_Real alpha = Angle;
1728       Standard_Boolean direct(axcone.Direct());
1729       if ((direct && Oris == TopAbs_REVERSED) ||
1730          (!direct && Oris == TopAbs_FORWARD)) {
1731         alpha = -alpha;
1732       }
1733
1734       gp_Pnt Center = i2s.Circle(1).Location();
1735       if (testdir <0.) {
1736         alpha = -alpha;
1737       }
1738       Standard_Real Z = ElCLib::LineParameter(Cy.Axis(),Center);
1739       Standard_Real Rad = Cy.Radius()+Z*Tan(alpha);
1740       if (Rad < 0.) {
1741         Rad = -Rad;
1742       }
1743       else {
1744         alpha = -alpha;
1745       }
1746       gp_Cone co(axcone,alpha,Rad);
1747       NewS = new Geom_ConicalSurface(co);
1748     }
1749     else {
1750       NewS = S;
1751     }
1752   }
1753   else if (TypeS == STANDARD_TYPE(Geom_ConicalSurface)) {
1754
1755     Standard_Real testdir = Direction.Dot(NeutralPlane.Axis().Direction());
1756     if (Abs(testdir) <= 1.-Precision::Angular()) {      
1757 #ifdef OCCT_DEBUG
1758       cout << "NewSurfaceCone:Draft_Direction_and_Neutral_Perpendicular" << endl;
1759 #endif
1760       return NewS;      
1761     }   
1762
1763     gp_Cone Co1 = Handle(Geom_ConicalSurface)::DownCast(S)->Cone();
1764
1765     testdir = Direction.Dot(Co1.Axis().Direction());
1766     if (Abs(testdir) <= 1.-Precision::Angular()) {
1767 #ifdef OCCT_DEBUG
1768       cout << "NewSurfaceCone:Draft_Direction_and_Cone_Perpendicular" << endl;
1769 #endif
1770       return NewS;
1771     }
1772
1773
1774     IntAna_QuadQuadGeo i2s;
1775     i2s.Perform(NeutralPlane,Co1,Precision::Angular(),Precision::Confusion());
1776     if (!i2s.IsDone() || i2s.TypeInter() != IntAna_Circle) {
1777 #ifdef OCCT_DEBUG
1778       cout << "NewSurfaceCone:Draft_Intersection_Neutral_Conical_NotDone" << endl;
1779 #endif
1780       return NewS;
1781     }
1782     gp_Ax3 axcone = Co1.Position();
1783     // Pb : Where is the material???
1784     Standard_Real alpha = Angle;
1785     Standard_Boolean direct(axcone.Direct());
1786     if ((direct && Oris == TopAbs_REVERSED) ||
1787        (!direct && Oris == TopAbs_FORWARD)) {
1788       alpha = -alpha;
1789     }
1790
1791     gp_Pnt Center = i2s.Circle(1).Location();
1792     if (Abs(Angle) > Precision::Angular()) {
1793       if (testdir <0.) {
1794         alpha = -alpha;
1795       }
1796       Standard_Real Z = ElCLib::LineParameter(Co1.Axis(),Center);
1797       Standard_Real Rad = i2s.Circle(1).Radius()+Z*Tan(alpha);
1798       if (Rad < 0.) {
1799         Rad = -Rad;
1800       }
1801       else {
1802         alpha = -alpha;
1803       }
1804       if (Abs(alpha-Co1.SemiAngle()) < Precision::Angular()) {
1805         NewS = S;
1806       }
1807       else {
1808         gp_Cone co(axcone,alpha,Rad);
1809         NewS = new Geom_ConicalSurface(co);
1810       }
1811     }
1812     else {
1813       NewS = new
1814         Geom_CylindricalSurface(gp_Cylinder(axcone,i2s.Circle(1).Radius()));
1815     }
1816   }
1817   else {
1818 #ifdef OCCT_DEBUG
1819     cout << "NewSurface:Draft_SurfNotYetImplemented" << endl;
1820 #endif
1821   }
1822   return NewS;
1823
1824
1825
1826 //=======================================================================
1827 //function : NewCurve
1828 //purpose  : 
1829 //=======================================================================
1830
1831 Handle(Geom_Curve) Draft_Modification::NewCurve
1832   (const Handle(Geom_Curve)& C,
1833   const Handle(Geom_Surface)& S,
1834   const TopAbs_Orientation Oris,
1835   const gp_Dir& Direction,
1836   const Standard_Real Angle,
1837   const gp_Pln& NeutralPlane,
1838   const Standard_Boolean )
1839
1840 {
1841   Handle(Geom_Curve) NewC;
1842
1843   Handle(Standard_Type) TypeS = S->DynamicType();
1844
1845   if (TypeS == STANDARD_TYPE(Geom_Plane)) {
1846     gp_Pln Pl = Handle(Geom_Plane)::DownCast(S)->Pln();
1847     gp_Ax1 Axe;
1848     Standard_Real Theta;
1849     if (FindRotation(Pl,Oris,Direction,Angle,NeutralPlane,Axe,Theta)) {
1850       if ( Abs(Theta) > Precision::Angular()) {
1851         NewC = Handle(Geom_Curve)::DownCast(C->Rotated(Axe,Theta));
1852       }
1853       else {
1854         NewC = C;
1855       }
1856     }
1857     return NewC;
1858   }
1859
1860
1861   if (C->DynamicType() != STANDARD_TYPE(Geom_Line)) {
1862     return NewC;
1863   }
1864
1865
1866   gp_Lin lin = Handle(Geom_Line)::DownCast(C)->Lin();
1867   //  Standard_Real testdir = Direction.Dot(lin.Direction());
1868   //  if (Abs(testdir) <= 1.-Precision::Angular()) {
1869   //    return NewC;
1870   //  }
1871   gp_Dir Norm;
1872   if (TypeS == STANDARD_TYPE(Geom_CylindricalSurface)) {
1873     Standard_Real U,V;
1874     gp_Vec d1u,d1v;
1875     gp_Pnt pbid;
1876     gp_Cylinder Cy = Handle(Geom_CylindricalSurface)::DownCast(S)->Cylinder();
1877     ElSLib::Parameters(Cy,lin.Location(),U,V);
1878     ElSLib::D1(U,V,Cy,pbid,d1u,d1v);
1879     Norm = d1u.Crossed(d1v);
1880   }
1881   else if (TypeS == STANDARD_TYPE(Geom_ConicalSurface)) {
1882     Standard_Real U,V;
1883     gp_Vec d1u,d1v;
1884     gp_Pnt pbid;
1885     gp_Cone Co = Handle(Geom_ConicalSurface)::DownCast(S)->Cone();
1886     ElSLib::Parameters(Co,lin.Location(),U,V);
1887     ElSLib::D1(U,V,Co,pbid,d1u,d1v);
1888     Norm = d1u.Crossed(d1v);
1889   }
1890
1891   IntAna_IntConicQuad ilipl(lin,NeutralPlane,Precision::Angular());
1892   if (ilipl.IsDone() && ilipl.NbPoints() != 0){
1893     if (Oris == TopAbs_REVERSED) {
1894       Norm.Reverse();
1895     }
1896     gp_Ax1 axrot(ilipl.Point(1), Norm.Crossed(Direction));
1897     gp_Lin lires = gp_Lin(gp_Ax1(ilipl.Point(1),Direction)).
1898       Rotated(axrot,Angle);
1899     if (lires.Direction().Dot(lin.Direction()) < 0.) {
1900       lires.Reverse();
1901     }
1902     NewC = new Geom_Line(lires);
1903   }
1904   return NewC;
1905 }
1906
1907
1908 //=======================================================================
1909 //function : Choose
1910 //purpose  : 
1911 //=======================================================================
1912
1913 static Standard_Boolean Choose(const Draft_IndexedDataMapOfFaceFaceInfo& theFMap,
1914   Draft_IndexedDataMapOfEdgeEdgeInfo& theEMap,
1915   const TopoDS_Vertex& Vtx,
1916   Draft_VertexInfo& Vinf,
1917   GeomAdaptor_Curve& AC,
1918   GeomAdaptor_Surface& AS)
1919 {
1920   gp_Vec tgref; 
1921   Vinf.InitEdgeIterator();
1922
1923   // Find a regular edge with null SecondFace
1924   while (Vinf.MoreEdge()) {
1925     const TopoDS_Edge& E1 = Vinf.Edge();
1926     const Draft_EdgeInfo& Einf1 = theEMap.FindFromKey(E1);
1927     if (Einf1.SecondFace().IsNull()) {
1928       break;
1929     }
1930     else {
1931       GeomAbs_Shape te = BRep_Tool::Continuity(E1,Einf1.FirstFace(),
1932         Einf1.SecondFace());
1933       if (te >= GeomAbs_G1) {
1934         break;
1935       }
1936     }
1937     Vinf.NextEdge();
1938   }
1939   if (!Vinf.MoreEdge()) { // take the first edge
1940     Vinf.InitEdgeIterator();
1941   }
1942
1943   const TopoDS_Edge& Eref = Vinf.Edge();
1944   //const Draft_EdgeInfo& Einf = theEMap(Eref);
1945   Draft_EdgeInfo& Einf = theEMap.ChangeFromKey(Eref);
1946
1947   AC.Load(Einf.Geometry());
1948
1949   Standard_Real f,l,prm;
1950   TopLoc_Location Loc;
1951   Handle(Geom_Curve) C = BRep_Tool::Curve(Eref,Loc,f,l);
1952   C = Handle(Geom_Curve)::DownCast(C->Transformed(Loc.Transformation()));
1953   gp_Pnt ptbid;
1954   //prm = Parameter(C,BRep_Tool::Pnt(Vtx));
1955   Standard_Integer done;
1956   Standard_Real param = Parameter( C, BRep_Tool::Pnt(Vtx), done );
1957   if (done != 0)
1958   {
1959     Handle( Geom_Surface ) S1 = theFMap.FindFromKey(Einf.FirstFace()).Geometry();
1960     Handle( Geom_Surface ) S2 = theFMap.FindFromKey(Einf.SecondFace()).Geometry();
1961     prm = SmartParameter( Einf, BRep_Tool::Tolerance(Eref), BRep_Tool::Pnt(Vtx), done, S1, S2 );
1962   }
1963   else
1964     prm = param;
1965   C->D1(prm,ptbid,tgref);
1966
1967
1968   Vinf.InitEdgeIterator();
1969   while (Vinf.MoreEdge()) {
1970     // Find a non tangent edge
1971     const TopoDS_Edge& Edg = Vinf.Edge();
1972     if (!Edg.IsSame(Eref)) {
1973       //const Draft_EdgeInfo& Einfo = theEMap(Edg);
1974       Draft_EdgeInfo& Einfo = theEMap.ChangeFromKey(Edg);
1975       if (!Einfo.SecondFace().IsNull() &&
1976           BRep_Tool::Continuity(Edg,Einfo.FirstFace(),Einfo.SecondFace()) 
1977           <= GeomAbs_C0) {
1978         C = BRep_Tool::Curve(Edg,Loc,f,l);
1979         C = Handle(Geom_Curve)::DownCast(C->Transformed(Loc.Transformation()));
1980         //prm = Parameter(C,BRep_Tool::Pnt(Vtx));
1981         Standard_Integer anewdone;
1982         Standard_Real anewparam = Parameter( C, BRep_Tool::Pnt(Vtx), anewdone );
1983         if (anewdone != 0)
1984         {
1985           Handle( Geom_Surface ) S1 = theFMap.FindFromKey(Einfo.FirstFace()).Geometry();
1986           Handle( Geom_Surface ) S2 = theFMap.FindFromKey(Einfo.SecondFace()).Geometry();
1987           prm = SmartParameter( Einfo, BRep_Tool::Tolerance(Edg), BRep_Tool::Pnt(Vtx), anewdone, S1, S2 );
1988         }
1989         else
1990           prm = anewparam;
1991         gp_Vec tg;
1992         C->D1(prm,ptbid,tg);
1993         if (tg.CrossMagnitude(tgref) > Precision::Confusion()) {
1994           break;
1995         }
1996       }
1997     }
1998     Vinf.NextEdge();
1999   }
2000   if (!Vinf.MoreEdge()) {
2001     return Standard_False;
2002   }
2003
2004   const Draft_EdgeInfo& Einf2 = theEMap.FindFromKey(Vinf.Edge());
2005   if (!Einf.SecondFace().IsNull()) {
2006
2007     if (Einf2.FirstFace().IsSame(Einf.FirstFace()) ||
2008         Einf2.FirstFace().IsSame(Einf.SecondFace())) {
2009       AS.Load(theFMap.FindFromKey(Einf2.SecondFace()).Geometry());
2010     }
2011     else {
2012       AS.Load(theFMap.FindFromKey(Einf2.FirstFace()).Geometry());
2013     }
2014   }
2015   else {
2016     if (Einf2.FirstFace().IsSame(Einf.FirstFace())) {
2017       AS.Load(theFMap.FindFromKey(Einf2.SecondFace()).Geometry());
2018     }
2019     else {
2020       AS.Load(theFMap.FindFromKey(Einf2.FirstFace()).Geometry());
2021     }
2022   }
2023   return Standard_True;
2024 }
2025
2026
2027 //=======================================================================
2028 //function : Parameter
2029 //purpose  : 
2030 //=======================================================================
2031
2032 static Standard_Real Parameter(const Handle(Geom_Curve)& C,
2033   const gp_Pnt& P,
2034   Standard_Integer& done)
2035 {
2036   done = 0;
2037   Handle(Geom_Curve) cbase = C;
2038   Handle(Standard_Type) ctyp = C->DynamicType();
2039   if (ctyp == STANDARD_TYPE(Geom_TrimmedCurve)) {
2040     cbase = Handle(Geom_TrimmedCurve)::DownCast(C)->BasisCurve();
2041     ctyp = cbase->DynamicType();
2042   }
2043   Standard_Real param;
2044   if (ctyp == STANDARD_TYPE(Geom_Line)) {
2045     param = ElCLib::Parameter(Handle(Geom_Line)::DownCast(cbase)->Lin(),P);
2046   }
2047   else if (ctyp == STANDARD_TYPE(Geom_Circle)) {
2048     param = ElCLib::Parameter(Handle(Geom_Circle)::DownCast(cbase)->Circ(),P);
2049     if (Abs(2.*M_PI-param) <=Epsilon(2.*M_PI)) {
2050       param = 0.;
2051     }
2052   }
2053   else if (ctyp == STANDARD_TYPE(Geom_Ellipse)) {
2054     param = ElCLib::Parameter(Handle(Geom_Ellipse)::DownCast(cbase)->Elips(),P);
2055     if (Abs(2.*M_PI-param) <=Epsilon(2.*M_PI)) {
2056       param = 0.;
2057     }
2058   }
2059   else if (ctyp == STANDARD_TYPE(Geom_Parabola)) {
2060     param = ElCLib::Parameter(Handle(Geom_Parabola)::DownCast(cbase)->Parab(),P);
2061   }
2062   else if (ctyp == STANDARD_TYPE(Geom_Hyperbola)) {
2063     param = ElCLib::Parameter(Handle(Geom_Hyperbola)::DownCast(cbase)->Hypr(),P);
2064   }
2065   else {
2066     GeomAdaptor_Curve TheCurve(C);
2067     Extrema_ExtPC myExtPC(P,TheCurve);
2068     if (!myExtPC.IsDone()) {
2069       Standard_Failure::Raise("Draft_Modification_1::Parameter: ExtremaPC not done.");
2070     }
2071     if (myExtPC.NbExt() >= 1) {
2072       Standard_Real Dist2, Dist2Min = myExtPC.SquareDistance(1);
2073       Standard_Integer j, jmin = 1;
2074       for (j = 2; j <= myExtPC.NbExt(); j++) {
2075         Dist2 = myExtPC.SquareDistance(j);
2076         if (Dist2 < Dist2Min) {
2077           Dist2Min = Dist2;
2078           jmin = j;
2079         }
2080       }
2081       param = myExtPC.Point(jmin).Parameter();
2082     }
2083     else {
2084       Standard_Real dist1_2,dist2_2;
2085       gp_Pnt p1b,p2b;
2086       myExtPC.TrimmedSquareDistances(dist1_2,dist2_2,p1b,p2b);
2087       if (dist1_2 < dist2_2) {
2088         done = -1;
2089         param = TheCurve.FirstParameter();
2090       }
2091       else {
2092         done = 1;
2093         param = TheCurve.LastParameter();
2094       }
2095     }
2096
2097     if (cbase->IsPeriodic()) {
2098       Standard_Real Per  = cbase->Period();
2099       Standard_Real Tolp = Precision::Parametric(Precision::Confusion());  
2100       if (Abs(Per-param) <= Tolp) {
2101         param = 0.;
2102       }
2103     }
2104   }
2105   return param;
2106 }
2107
2108 //=======================================================================
2109 //function : SmartParameter
2110 //purpose  : 
2111 //=======================================================================
2112
2113 static Standard_Real SmartParameter(Draft_EdgeInfo& Einf,
2114   const Standard_Real EdgeTol,
2115   const gp_Pnt& Pnt,
2116   const Standard_Integer sign,
2117   const Handle(Geom_Surface)& S1,
2118   const Handle(Geom_Surface)& S2)
2119 {
2120   Handle( Geom2d_Curve ) NewC2d;
2121   Standard_Real Tol = Precision::Confusion();
2122   Standard_Real Etol = EdgeTol;
2123
2124   Handle( Geom2d_Curve ) pcu1 = Einf.FirstPC();
2125   Handle( Geom2d_Curve ) pcu2 = Einf.SecondPC();
2126
2127   if (pcu1.IsNull())
2128   {
2129     Handle( Geom_Curve ) theCurve = Einf.Geometry();
2130     pcu1 = GeomProjLib::Curve2d( theCurve, theCurve->FirstParameter(), theCurve->LastParameter(), S1, Etol );
2131     Einf.ChangeFirstPC() = pcu1;
2132   }
2133   if (pcu2.IsNull())
2134   {
2135     Handle( Geom_Curve ) theCurve = Einf.Geometry();
2136     pcu2 = GeomProjLib::Curve2d( theCurve, theCurve->FirstParameter(), theCurve->LastParameter(), S2, Etol );
2137     Einf.ChangeSecondPC() = pcu2;
2138   }
2139
2140   GeomAPI_ProjectPointOnSurf Projector( Pnt, S1 );
2141   Standard_Real U, V;
2142   Projector.LowerDistanceParameters( U, V );
2143
2144   NewC2d = Einf.FirstPC();
2145   if (NewC2d->DynamicType() == STANDARD_TYPE(Geom2d_TrimmedCurve))
2146     NewC2d = (Handle(Geom2d_TrimmedCurve)::DownCast(NewC2d))->BasisCurve();
2147
2148   gp_Pnt2d P2d( U, V );
2149   Geom2dAPI_ProjectPointOnCurve Projector2d( P2d, NewC2d );
2150   if (Projector2d.NbPoints() == 0 || Projector2d.LowerDistance() > Tol)
2151   {
2152     Handle( Geom2d_BSplineCurve ) BCurve;
2153     if (NewC2d->DynamicType() != STANDARD_TYPE(Geom2d_BSplineCurve))
2154       BCurve = Geom2dConvert::CurveToBSplineCurve( NewC2d );
2155     else
2156       BCurve = Handle( Geom2d_BSplineCurve )::DownCast( NewC2d );
2157     if (sign == -1)
2158     {
2159       TColgp_Array1OfPnt2d PntArray( 1, 2 );
2160       PntArray(1) = P2d;
2161       PntArray(2) = BCurve->Pole(1);
2162       Handle( Geom2d_BezierCurve ) Patch = new Geom2d_BezierCurve( PntArray );
2163       Geom2dConvert_CompCurveToBSplineCurve Concat( BCurve, Convert_QuasiAngular );
2164       Concat.Add( Patch, Tol, Standard_False );
2165       BCurve = Concat.BSplineCurve();
2166     }
2167     else
2168     {
2169       TColgp_Array1OfPnt2d PntArray( 1, 2 );
2170       PntArray(1) = BCurve->Pole( BCurve->NbPoles() );
2171       PntArray(2) = P2d;
2172       Handle( Geom2d_BezierCurve ) Patch = new Geom2d_BezierCurve( PntArray );
2173       Geom2dConvert_CompCurveToBSplineCurve Concat( BCurve, Convert_QuasiAngular );
2174       Concat.Add( Patch, Tol, Standard_True );
2175       BCurve = Concat.BSplineCurve();
2176     }
2177     NewC2d = BCurve;
2178   }
2179   Einf.ChangeFirstPC() = NewC2d;
2180   Handle( Geom2dAdaptor_HCurve ) hcur = new Geom2dAdaptor_HCurve( NewC2d );
2181   Handle( GeomAdaptor_HSurface ) hsur = new GeomAdaptor_HSurface( S1 );
2182   Adaptor3d_CurveOnSurface cons( hcur, hsur );
2183   Handle( Adaptor3d_HCurveOnSurface ) hcons = new Adaptor3d_HCurveOnSurface( cons );
2184   Handle( GeomAdaptor_HSurface ) hsur2 = new GeomAdaptor_HSurface( S2 );
2185   ProjLib_CompProjectedCurve ProjCurve( hsur2, hcons, Tol, Tol );
2186   Handle(ProjLib_HCompProjectedCurve) HProjector = new ProjLib_HCompProjectedCurve();
2187   HProjector->Set( ProjCurve );
2188   Standard_Real Udeb, Ufin;
2189   ProjCurve.Bounds(1, Udeb, Ufin);
2190   Standard_Integer MaxSeg = 20 + HProjector->NbIntervals(GeomAbs_C3);
2191   Approx_CurveOnSurface appr( HProjector, hsur2, Udeb, Ufin, Tol,
2192     GeomAbs_C1, 10, MaxSeg, 
2193     Standard_False, Standard_False );
2194   Einf.ChangeSecondPC() = appr.Curve2d();
2195   Einf.ChangeGeometry() = appr.Curve3d();
2196   Einf.SetNewGeometry( Standard_True );
2197
2198   if (sign == -1)
2199     return Einf.Geometry()->FirstParameter();
2200   else
2201     return Einf.Geometry()->LastParameter();
2202
2203 }
2204
2205 //=======================================================================
2206 //function : Orientation
2207 //purpose  : 
2208 //=======================================================================
2209
2210 static TopAbs_Orientation Orientation(const TopoDS_Shape& S,
2211   const TopoDS_Face& F)
2212 {
2213   //
2214   // change porting NT
2215   //
2216   TopExp_Explorer expl ;
2217   expl.Init(S,
2218     TopAbs_FACE) ;
2219   while (expl.More()) {
2220     if (TopoDS::Face(expl.Current()).IsSame(F)) {
2221       return expl.Current().Orientation();
2222     }
2223     expl.Next();
2224   }
2225   return TopAbs_FORWARD;
2226 }
2227
2228
2229 //=======================================================================
2230 //function : FindRotation
2231 //purpose  : 
2232 //=======================================================================
2233
2234 static Standard_Boolean FindRotation(const gp_Pln& Pl,
2235   const TopAbs_Orientation Oris,
2236   const gp_Dir& Direction,
2237   const Standard_Real Angle,
2238   const gp_Pln& NeutralPlane,
2239   gp_Ax1& Axe,
2240   Standard_Real& theta)
2241 {                                    
2242   IntAna_QuadQuadGeo i2pl(Pl,NeutralPlane,
2243     Precision::Angular(),Precision::Confusion());
2244
2245   if (i2pl.IsDone() && i2pl.TypeInter() == IntAna_Line) {
2246     gp_Lin li = i2pl.Line(1);
2247     // Try to turn around this line
2248     gp_Dir nx = li.Direction();
2249     gp_Dir ny = Pl.Axis().Direction().Crossed(nx);
2250     Standard_Real a = Direction.Dot(nx);
2251     if (Abs(a) <=1-Precision::Angular()) { 
2252       Standard_Real b = Direction.Dot(ny);
2253       Standard_Real c = Direction.Dot(Pl.Axis().Direction());
2254       Standard_Boolean direct(Pl.Position().Direct());
2255       if ((direct && Oris == TopAbs_REVERSED) ||
2256         (!direct && Oris == TopAbs_FORWARD)) {
2257           b = -b;
2258           c = -c;
2259       }
2260       Standard_Real denom = Sqrt(1-a*a);
2261       Standard_Real Sina = Sin(Angle);
2262       if (denom>Abs(Sina)) {
2263         Standard_Real phi = ATan2(b/denom,c/denom);
2264         Standard_Real theta0 = ACos(Sina/denom); 
2265         theta = theta0 - phi;
2266         if (Cos(theta) <0.) {
2267           theta = -theta0 -phi;
2268         }
2269         //  modified by NIZHNY-EAP Tue Nov 16 15:51:38 1999 ___BEGIN___
2270         while (Abs(theta)>M_PI) {
2271           theta = theta + M_PI*(theta<0 ? 1 : -1);
2272         }
2273         //  modified by NIZHNY-EAP Tue Nov 16 15:53:32 1999 ___END___
2274         Axe = li.Position();
2275         return Standard_True;
2276       }
2277     }
2278   }
2279   return Standard_False;
2280 }
2281