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 <Adaptor3d_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               Adaptor3d_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 Dist2, Dist2Min = 0., 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               imin = 0;
 957               for (i=1; i<= i2s.NbLines(); i++) {
 958                 TheCurve.Load(i2s.Line(i));
 959                 Extrema_ExtPC myExtPC(pfv,TheCurve);
 960
 961                 Standard_Real locpmin = 0.;
 962                 if (myExtPC.IsDone()) {
 963                   if(myExtPC.NbExt() >= 1) {
 964                     Dist2Min = myExtPC.SquareDistance(1);
 965                     locpmin = myExtPC.Point(1).Parameter();
 966                   }
 967                   if(myExtPC.NbExt() == 2 && Dist2Min > Precision::SquareConfusion()) {
 968                     //to avoid incorrectly choosing the image
 969                     //of the first vertex of the initial edge
 970                     Standard_Real d1_2 = myExtPC.SquareDistance(1);
 971                     Standard_Real d2_2 = myExtPC.SquareDistance(2);
 972                     if(d1_2 > 1.21*d2_2) {
 973                       Dist2Min = myExtPC.SquareDistance(2);
 974                       locpmin = myExtPC.Point(2).Parameter();
 975                     }
 976                     else if(d2_2 > 1.21*d1_2) {
 977                       Dist2Min = myExtPC.SquareDistance(1);
 978                       locpmin = myExtPC.Point(1).Parameter();
 979                     }
 980                     else {
 981                       Standard_Real pfvpar = myExtPC.Point(1).Parameter();
 982                       Standard_Real plvpar = myExtPC.Point(2).Parameter();
 983                       newC = i2s.Line(i);
 984
 985                       gp_Pnt pfvprim, plvprim;
 986
 987                       newC->D0(pfvpar,pfvprim);
 988                       newC->D0(plvpar,plvprim);
 989
 990                       Handle(Geom_Surface) theSurf;
 991                       if(detrompeur == 1) {
 992                         if(S1->DynamicType() ==
 993                            STANDARD_TYPE(Geom_RectangularTrimmedSurface))
 994                           S1 = Handle(Geom_RectangularTrimmedSurface)::
 995                           DownCast(S1)->BasisSurface();
 996                         theSurf = S1;
 997
 998                       }
 999                       else if(detrompeur == 2) {
1000                         if(S2->DynamicType() ==
1001                            STANDARD_TYPE(Geom_RectangularTrimmedSurface))
1002                           S2 = Handle(Geom_RectangularTrimmedSurface)::
1003                           DownCast(S2)->BasisSurface();
1004                         theSurf = S2;
1005                       }
1006                       if(detrompeur != 0 && detrompeur != 4) {
1007                         Standard_Real ul = 0., vl = 0., uf = 0., vf = 0.;
1008                         Standard_Real ufprim = 0., ulprim = 0., vfprim = 0., vlprim = 0.;
1009
1010                         if(theSurf->DynamicType() == STANDARD_TYPE(Geom_Plane)) {
1011                           gp_Pln pl = Handle(Geom_Plane)::DownCast(S2)->Pln();
1012                           ElSLib::Parameters(pl, plv, ul, vl);
1013                           ElSLib::Parameters(pl, pfv, uf, vf);
1014                           ElSLib::Parameters(pl, plvprim, ulprim, vlprim);
1015                           ElSLib::Parameters(pl, pfvprim, ufprim, vfprim);
1016                         }
1017                         else if(theSurf->DynamicType() ==
1018                                 STANDARD_TYPE(Geom_CylindricalSurface)) {
1019                           gp_Cylinder cy = Handle(Geom_CylindricalSurface)
1020                             ::DownCast(S2)->Cylinder();
1021                           ElSLib::Parameters(cy, plv, ul, vl);
1022                           ElSLib::Parameters(cy, pfv, uf, vf);
1023                           ElSLib::Parameters(cy, plvprim, ulprim, vlprim);
1024                           ElSLib::Parameters(cy, pfvprim, ufprim, vfprim);
1025                         }
1026                         else detrompeur = 4;
1027
1028                         if(detrompeur == 1 || detrompeur == 2) {
1029                           gp_Vec2d v1((ul-ufprim), (vl-vfprim));
1030                           gp_Vec2d norm((vf-vfprim), (ufprim-uf));
1031                           gp_Vec2d v2((ulprim-ufprim), (vlprim-vfprim));
1032                           if( (v1.Dot(norm))*(v2.Dot(norm)) < 0) {
1033                             Dist2Min = myExtPC.SquareDistance(2);
1034                             locpmin = myExtPC.Point(2).Parameter();
1035                           }
1036                         }
1037                       }
1038                     }
1039                   }
1040                   if (myExtPC.NbExt() == 1 || myExtPC.NbExt() > 2 || detrompeur ==4) {
1041                     Dist2Min = myExtPC.SquareDistance(1);
1042                     locpmin = myExtPC.Point(1).Parameter();
1043                     for (j=2; j<=myExtPC.NbExt(); j++) {
1044                       Dist2 = myExtPC.SquareDistance(j);
1045                       if (Dist2 < Dist2Min) {
1046                         Dist2Min = Dist2;
1047                         locpmin = myExtPC.Point(j).Parameter();
1048                       }
1049                     }
1050                   }
1051                   else if(myExtPC.NbExt() < 1){
1052                     Standard_Real dist1_2,dist2_2;
1053                     gp_Pnt p1b,p2b;
1054                     myExtPC.TrimmedSquareDistances(dist1_2,dist2_2,p1b,p2b);
1055                     if (dist1_2 < dist2_2) {
1056                       Dist2Min = dist1_2;
1057                       locpmin = TheCurve.FirstParameter();
1058                     }
1059                     else {
1060                       Dist2Min = dist2_2;
1061                       locpmin = TheCurve.LastParameter();
1062                     }
1063                   }
1064
1065                   if (Dist2Min  < Glob2Min) {
1066                     Glob2Min = Dist2Min;
1067                     imin = i;
1068                     pmin = locpmin;
1069                   }
1070                 }
1071               }
1072               if (imin == 0) {
1073                 errStat = Draft_EdgeRecomputation;
1074                 badShape = theEdge;
1075                 return;
1076               }
1077
1078               newC = i2s.Line(imin);
1079
1080               newC->D1(pmin,pfv,newd1);
1081               Standard_Boolean YaRev = d1fv.Dot(newd1) <0.;
1082
1083               if (YaRev)
1084                 newC->Reverse();
1085
1086               if (i2s.HasLineOnS1(imin)) {
1087                 Einf.ChangeFirstPC() = i2s.LineOnS1(imin);
1088                 if ( YaRev)
1089                   Einf.ChangeFirstPC()->Reverse();
1090               }
1091
1092               if (i2s.HasLineOnS2(imin)) {
1093                 Einf.ChangeSecondPC() = i2s.LineOnS2(imin);
1094                 if ( YaRev)
1095                   Einf.ChangeSecondPC()->Reverse();
1096               }
1097             } // if (i2s.Line(1)->DynamicType() != STANDARD_TYPE(Geom_BSplineCurve))
1098             else // i2s.Line(1) is BSplineCurve
1099             {
1100               //Find the first curve to glue
1101               TColGeom_SequenceOfCurve Candidates;
1102               if (S1->DynamicType() == STANDARD_TYPE(Geom_CylindricalSurface) ||
1103                   S1->DynamicType() == STANDARD_TYPE(Geom_ConicalSurface))
1104               {
1105                 for (i = 1; i <= i2s.NbLines(); i++)
1106                 {
1107                   Handle( Geom_Curve ) aCurve = i2s.Line(i);
1108                   gp_Pnt Pnt = aCurve->Value( aCurve->FirstParameter() );
1109                   GeomAPI_ProjectPointOnSurf projector( Pnt, S1, Precision::Confusion() );
1110                   Standard_Real U, V;
1111                   projector.LowerDistanceParameters( U, V );
1112                   if (Abs(U) <= Precision::Confusion() || Abs(U-2.*M_PI) <= Precision::Confusion())
1113                     Candidates.Append( aCurve );
1114                   else
1115                   {
1116                     Pnt = aCurve->Value( aCurve->LastParameter() );
1117                     projector.Init( Pnt, S1, Precision::Confusion() );
1118                     projector.LowerDistanceParameters( U, V );
1119                     if (Abs(U) <= Precision::Confusion() || Abs(U-2.*M_PI) <= Precision::Confusion())
1120                     {
1121                       aCurve->Reverse();
1122                       Candidates.Append( aCurve );
1123                     }
1124                   }
1125                 }
1126
1127                 if(Candidates.Length() == 0)
1128                 {
1129                   //errStat = Draft_EdgeRecomputation;
1130                   //badShape = TopoDS::Edge(ite.Key());
1131                   //return;
1132                   for (i = 1; i <= i2s.NbLines(); i++)
1133                     Candidates.Append( i2s.Line(i) );
1134                 }
1135               }
1136               else
1137               {
1138                 for (i = 1; i <= i2s.NbLines(); i++)
1139                   Candidates.Append( i2s.Line(i) );
1140               }
1141
1142               Handle( Geom_Curve ) FirstCurve;
1143               if (Candidates.Length() > 1)
1144               {
1145                 Dist2Min = RealLast();
1146                 for (i = 1; i <= Candidates.Length(); i++)
1147                 {
1148                   Handle( Geom_Curve ) aCurve = Candidates(i);
1149                   gp_Pnt Pnt = aCurve->Value( aCurve->FirstParameter() );
1150                   Dist2 = Pnt.SquareDistance( pfv );
1151                   if (Dist2 < Dist2Min)
1152                   {
1153                     Dist2Min = Dist2;
1154                     FirstCurve = aCurve;
1155                   }
1156                 }
1157               }
1158               else
1159                 FirstCurve = Candidates(1);
1160
1161               //Glueing
1162               TColGeom_SequenceOfCurve Curves;
1163               for (i = 1; i <= i2s.NbLines(); i++)
1164                 if (FirstCurve != i2s.Line(i))
1165                   Curves.Append( i2s.Line(i) );
1166
1167               TColGeom_SequenceOfCurve ToGlue;
1168               gp_Pnt EndPoint = FirstCurve->Value( FirstCurve->LastParameter() );
1169               Standard_Boolean added = Standard_True;
1170               while (added)
1171               {
1172                 added = Standard_False;
1173                 for (i = 1; i <= Curves.Length(); i++)
1174                 {
1175                   Handle( Geom_Curve ) aCurve = Curves(i);
1176                   gp_Pnt pfirst, plast;
1177                   pfirst = aCurve->Value( aCurve->FirstParameter() );
1178                   plast = aCurve->Value( aCurve->LastParameter() );
1179                   if (pfirst.Distance( EndPoint ) <= Precision::Confusion())
1180                   {
1181                     ToGlue.Append( aCurve );
1182                     EndPoint = plast;
1183                     Curves.Remove(i);
1184                     added = Standard_True;
1185                     break;
1186                   }
1187                   if (plast.Distance( EndPoint ) <= Precision::Confusion())
1188                   {
1189                     aCurve->Reverse();
1190                     ToGlue.Append( aCurve );
1191                     EndPoint = pfirst;
1192                     Curves.Remove(i);
1193                     added = Standard_True;
1194                     break;
1195                   }
1196                 }
1197               }
1198
1199               if (FirstCurve.IsNull()) {
1200                 errStat = Draft_EdgeRecomputation;
1201                 badShape = theEdge;
1202                 return;
1203               }
1204
1205               GeomConvert_CompCurveToBSplineCurve Concat( Handle(Geom_BSplineCurve)::DownCast(FirstCurve) );
1206               for (i = 1; i <= ToGlue.Length(); i++)
1207                 Concat.Add( Handle(Geom_BSplineCurve)::DownCast(ToGlue(i)), Precision::Confusion(), Standard_True );
1208
1209               newC = Concat.BSplineCurve();
1210
1211               TheCurve.Load( newC );
1212               Extrema_ExtPC myExtPC( pfv, TheCurve );
1213               Dist2Min = RealLast();
1214               for (i = 1; i <= myExtPC.NbExt(); i++)
1215               {
1216                 if (myExtPC.IsMin(i))
1217                 {
1218                   Dist2 = myExtPC.SquareDistance(i);
1219                   if (Dist2 < Dist2Min)
1220                   {
1221                     Dist2Min = Dist2;
1222                     pmin = myExtPC.Point(i).Parameter();
1223                   }
1224                 }
1225               }
1226               newC->D1(pmin,pfv,newd1);
1227               Standard_Boolean YaRev = d1fv.Dot(newd1) < 0.;
1228
1229               if (YaRev)
1230                 newC->Reverse();
1231               /*
1232               if (i2s.HasLineOnS1(imin)) {
1233               Einf.ChangeFirstPC() = i2s.LineOnS1(imin);
1234               if ( YaRev)
1235               Einf.ChangeFirstPC()->Reverse();
1236               }
1237
1238               if (i2s.HasLineOnS2(imin)) {
1239               Einf.ChangeSecondPC() = i2s.LineOnS2(imin);
1240               if ( YaRev)
1241               Einf.ChangeSecondPC()->Reverse();
1242               }
1243               */
1244             } // else: i2s.NbLines() > 2 && S1 is Cylinder or Cone
1245
1246             Einf.Tolerance(Max(Einf.Tolerance(), i2s.TolReached3d()));
1247           }  // End step KPart
1248         }
1249         else { // case of tangency
1250           const TopoDS_Face& F1 = Einf.FirstFace();
1251           const TopoDS_Face& F2 = Einf.SecondFace();
1252
1253           Handle(Geom_Surface) aLocalS1 = myFMap.FindFromKey(F1).Geometry();
1254           Handle(Geom_Surface) aLocalS2 = myFMap.FindFromKey(F2).Geometry();
1255           if (aLocalS1.IsNull() || aLocalS2.IsNull()) {
1256             errStat = Draft_EdgeRecomputation;
1257             badShape = theEdge;
1258             return;
1259           }
1260           if (aLocalS1->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
1261             aLocalS1 = Handle(Geom_RectangularTrimmedSurface)::
1262               DownCast(aLocalS1)->BasisSurface();
1263           }
1264           if (aLocalS2->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
1265             aLocalS2 = Handle(Geom_RectangularTrimmedSurface)::
1266               DownCast(aLocalS2)->BasisSurface();
1267           }
1268
1269           gp_Dir dirextr;
1270           //Standard_Boolean dirfound = Standard_False;
1271           if (aLocalS1->DynamicType() == STANDARD_TYPE(Geom_CylindricalSurface)) {
1272             gp_Cylinder cyl =
1273               Handle(Geom_CylindricalSurface)::DownCast(aLocalS1)->Cylinder();
1274             dirextr = cyl.Axis().Direction();
1275             //dirfound = Standard_True;
1276             // see direction...
1277
1278           }
1279           else if (aLocalS1->DynamicType() == STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion)) {
1280             dirextr = Handle(Geom_SurfaceOfLinearExtrusion)::
1281               DownCast(aLocalS1)->Direction();
1282             //dirfound = Standard_True;
1283             // see direction...
1284
1285             // Here it is possible to calculate PCurve.
1286             Handle(Geom_SurfaceOfLinearExtrusion) SEL =
1287               Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(aLocalS1);
1288             Handle(Geom_Circle) GCir =
1289               Handle(Geom_Circle)::DownCast(SEL->BasisCurve());
1290             if ( !GCir.IsNull()) {
1291               Standard_Real U = ElCLib::Parameter(GCir->Circ(),ptfixe);
1292               Handle(Geom2d_Line) PC1 =
1293                 new Geom2d_Line(gp_Pnt2d(U,0.),gp::DY2d());
1294               Einf.ChangeFirstPC() = PC1;
1295             }
1296           }
1297
1298           else if (aLocalS2->DynamicType() == STANDARD_TYPE(Geom_CylindricalSurface)) {
1299             gp_Cylinder cyl =
1300               Handle(Geom_CylindricalSurface)::DownCast(aLocalS2)->Cylinder();
1301             dirextr = cyl.Axis().Direction();
1302             // dirfound = Standard_True;
1303             // see direction...
1304
1305           }
1306           else if (aLocalS2->DynamicType() == STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion)) {
1307             dirextr = Handle(Geom_SurfaceOfLinearExtrusion)::
1308               DownCast(aLocalS2)->Direction();
1309             // dirfound = Standard_True;
1310             // see direction...
1311
1312             // Here it is possible to calculate PCurve.
1313             Handle(Geom_SurfaceOfLinearExtrusion) SEL =
1314               Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(aLocalS2);
1315             Handle(Geom_Circle) GCir =
1316               Handle(Geom_Circle)::DownCast(SEL->BasisCurve());
1317             if ( !GCir.IsNull()) {
1318               Standard_Real U = ElCLib::Parameter(GCir->Circ(),ptfixe);
1319               Handle(Geom2d_Line) PC2 =
1320                 new Geom2d_Line(gp_Pnt2d(U,0.),gp::DY2d());
1321               Einf.ChangeSecondPC() = PC2;
1322             }
1323           }
1324           newC = new Geom_Line(ptfixe,dirextr);
1325
1326           gp_Pnt pfv;
1327           gp_Vec d1fv,newd1;
1328           C->D1(0.,pfv,d1fv);
1329           newC->D1(0.,pfv,newd1);
1330           Standard_Boolean YaRev = d1fv.Dot(newd1) <0.;
1331           if (YaRev) {
1332             newC->Reverse();
1333             if(!Einf.FirstPC().IsNull()) {
1334               Einf.ChangeFirstPC()->Reverse();
1335             }
1336             if(!Einf.SecondPC().IsNull()) {
1337               Einf.ChangeSecondPC()->Reverse();
1338             }
1339           }
1340         }
1341
1342         Handle(Geom_TrimmedCurve) T = Handle(Geom_TrimmedCurve)::DownCast(newC);
1343         if (!T.IsNull()) newC = T->BasisCurve();
1344         Einf.ChangeGeometry() = newC;
1345       }
1346       else if (!Einf.NewGeometry()){
1347         // set existing curve 3D
1348         Handle(Geom_TrimmedCurve) T = Handle(Geom_TrimmedCurve)::DownCast(C);
1349         if (!T.IsNull()) C = T->BasisCurve();
1350         Einf.ChangeGeometry() = C;
1351       }
1352     }
1353
1354     // Calculate new vertices.
1355
1356     Handle(GeomAdaptor_HCurve)   HAC = new GeomAdaptor_HCurve;
1357     Handle(GeomAdaptor_HSurface) HAS = new GeomAdaptor_HSurface;
1358
1359     for (Standard_Integer ii = 1; ii <= myVMap.Extent(); ii++)
1360     {
1361       GeomAdaptor_Curve    AC;
1362       GeomAdaptor_Surface  AS;
1363
1364       const TopoDS_Vertex& TVV = myVMap.FindKey(ii);
1365       Draft_VertexInfo& Vinf = myVMap.ChangeFromIndex(ii);
1366       if (!Choose(myFMap,myEMap,TVV,Vinf,AC,AS)) {
1367
1368         // no concerted edge => alignment of two consecutive edges.
1369         gp_Pnt pvt;
1370         Vinf.ChangeGeometry() = pvt;
1371         Vinf.InitEdgeIterator();
1372         if (Vinf.MoreEdge()) {
1373           const TopoDS_Edge& Edg1 = Vinf.Edge();
1374           //const Draft_EdgeInfo& Einf1 = myEMap(Edg1);
1375           Draft_EdgeInfo& Einf1 = myEMap.ChangeFromKey(Edg1);
1376           gp_Pnt vtori = BRep_Tool::Pnt(TVV);
1377           //Einf1.Geometry()->D0(Vinf.Parameter(Edg1), pvt);
1378           GeomAPI_ProjectPointOnCurve Projector( vtori, Einf1.Geometry() ); //patch
1379           pvt = Projector.NearestPoint();
1380
1381 #ifdef OCCT_DEBUG
1382           static Standard_Integer VertexRecomp = 1;
1383           if (VertexRecomp!=0) {
1384             cout << "pori :" << vtori.X() << " " << vtori.Y() << " " << vtori.Z() << endl;
1385             cout << "  Edg 1 :" << Vinf.Parameter(Edg1) << endl;
1386             cout << "pvt :" << pvt.X() << " " << pvt.Y() << " " << pvt.Z() << endl;
1387           }
1388 #endif
1389
1390           Standard_Real dion=pvt.SquareDistance(vtori);
1391           Vinf.NextEdge();
1392           if (Vinf.MoreEdge()) {
1393             const TopoDS_Edge& Edg2 = Vinf.Edge();
1394             //const Draft_EdgeInfo& Einf2 = myEMap(Edg2);
1395             Draft_EdgeInfo& Einf2 = myEMap.ChangeFromKey(Edg2);
1396             //      Standard_Real f;
1397             gp_Pnt opvt;
1398             Einf2.Geometry()->D0(Vinf.Parameter(Edg2), opvt);
1399
1400 #ifdef OCCT_DEBUG
1401             if (VertexRecomp!=0) {
1402               cout << "  Edg 2 :" << Vinf.Parameter(Vinf.Edge()) << endl;
1403               cout << "opvt " << opvt.X() << " " << opvt.Y() << " " << opvt.Z() << endl;
1404             }
1405 #endif
1406
1407             if (opvt.SquareDistance(vtori) < dion) {
1408               pvt = opvt;
1409             }
1410             //Vinf.ChangeParameter(Edg2) = Parameter(Einf2.Geometry(), pvt);
1411             Standard_Integer done;
1412             Standard_Real param = Parameter(Einf2.Geometry(), pvt, done);
1413             if (done != 0)
1414             {
1415               Handle(Geom_Surface) S1 = myFMap.FindFromKey(Einf2.FirstFace()).Geometry();
1416               Handle(Geom_Surface) S2 = myFMap.FindFromKey(Einf2.SecondFace()).Geometry();
1417               Vinf.ChangeParameter(Edg2) = SmartParameter( Einf2, BRep_Tool::Tolerance(Edg2), pvt, done, S1, S2 );
1418             }
1419             else
1420               Vinf.ChangeParameter(Edg2) = param;
1421           }
1422
1423           Vinf.ChangeGeometry() = pvt;
1424           //Vinf.ChangeParameter(Edg1) = Parameter(Einf1.Geometry(), pvt);
1425           Standard_Integer done;
1426           Standard_Real param = Parameter(Einf1.Geometry(), pvt, done);
1427           if (done != 0)
1428           {
1429             Handle(Geom_Surface) S1 = myFMap.FindFromKey(Einf1.FirstFace()).Geometry();
1430             Handle(Geom_Surface) S2 = myFMap.FindFromKey(Einf1.SecondFace()).Geometry();
1431             Vinf.ChangeParameter(Edg1) = SmartParameter( Einf1, BRep_Tool::Tolerance(Edg1), pvt, done, S1, S2 );
1432           }
1433           else
1434             Vinf.ChangeParameter(Edg1) = param;
1435           continue;
1436         }
1437
1438
1439         errStat = Draft_VertexRecomputation;
1440         badShape = TVV;
1441         return;
1442       }
1443
1444       IntCurveSurface_HInter myintcs;
1445       HAC->Set(AC);
1446       HAS->Set(AS);
1447
1448       myintcs.Perform(HAC,HAS);
1449
1450
1451       if (!myintcs.IsDone()) {
1452         errStat = Draft_VertexRecomputation;
1453         badShape = TVV;
1454         return;
1455       }
1456
1457       gp_Pnt vtori = BRep_Tool::Pnt(TVV);
1458       gp_Pnt pvt;
1459
1460       Standard_Integer nbsol = myintcs.NbPoints();
1461       if (nbsol <= 0)
1462       {
1463         Extrema_ExtCS extr( AC, AS, Precision::PConfusion(), Precision::PConfusion() );
1464
1465         if(!extr.IsDone() || extr.NbExt() == 0) {
1466           errStat = Draft_VertexRecomputation;
1467           badShape = TVV;
1468           return;
1469         }
1470
1471
1472         Standard_Real disref = RealLast();
1473         Standard_Integer iref = 0;
1474         Extrema_POnCurv Pc;
1475         Extrema_POnSurf Ps;
1476         for (Standard_Integer i = 1; i <= extr.NbExt(); i++)
1477         {
1478           extr.Points( i, Pc, Ps );
1479           Standard_Real distemp = Pc.Value().SquareDistance(vtori);
1480           if ( distemp < disref)
1481           {
1482             disref = distemp;
1483             iref = i;
1484           }
1485         }
1486         extr.Points( iref, Pc, Ps );
1487         pvt = Pc.Value();
1488       }
1489       else
1490       {
1491         Standard_Real disref = RealLast();
1492         Standard_Integer iref = 0;
1493         for (Standard_Integer i = 1; i <= nbsol; i++)
1494         {
1495           Standard_Real distemp = myintcs.Point(i).Pnt().SquareDistance(vtori);
1496           if ( distemp < disref)
1497           {
1498             disref = distemp;
1499             iref = i;
1500           }
1501         }
1502         pvt = myintcs.Point(iref).Pnt();
1503       }
1504
1505       Vinf.ChangeGeometry() = pvt;
1506
1507       for (Vinf.InitEdgeIterator();Vinf.MoreEdge(); Vinf.NextEdge()) {
1508         const TopoDS_Edge& Edg = Vinf.Edge();
1509         Standard_Real initpar = Vinf.Parameter(Edg);
1510         //const Draft_EdgeInfo& Einf = myEMap(Edg);
1511         Draft_EdgeInfo& Einf = myEMap.ChangeFromKey(Edg);
1512         //Vinf.ChangeParameter(Edg) = Parameter(Einf.Geometry(),pvt);
1513         Standard_Integer done;
1514         Standard_Real param = Parameter(Einf.Geometry(), pvt, done);
1515         if (done != 0)
1516         {
1517           Handle(Geom_Surface) S1 = myFMap.FindFromKey(Einf.FirstFace()).Geometry();
1518           Handle(Geom_Surface) S2 = myFMap.FindFromKey(Einf.SecondFace()).Geometry();
1519           Vinf.ChangeParameter(Edg) = SmartParameter( Einf, BRep_Tool::Tolerance(Edg), pvt, done, S1, S2 );
1520         }
1521         else
1522         {
1523           if(Abs(initpar - param) > Precision::PConfusion())
1524           {
1525             Standard_Real f, l;
1526             TopLoc_Location Loc;
1527             const Handle(Geom_Curve)& aC = BRep_Tool::Curve(Edg, Loc, f, l);
1528             if(aC->DynamicType() == STANDARD_TYPE(Geom_TrimmedCurve))
1529             {
1530               Einf.SetNewGeometry(Standard_True);
1531             }
1532           }
1533           Vinf.ChangeParameter(Edg) = param;
1534         }
1535       }
1536     }
1537   }
1538
1539   // small loop of validation/protection
1540
1541   for (Standard_Integer i = 1; i <= myEMap.Extent(); i++)
1542   {
1543     const TopoDS_Edge& edg = TopoDS::Edge(myEMap.FindKey(i));
1544
1545     TopoDS_Vertex Vf,Vl;
1546     TopExp::Vertices(edg,Vf,Vl);
1547     if (edg.Orientation() == TopAbs_REVERSED) {
1548       Vf.Reverse();
1549       Vl.Reverse();
1550     }
1551     Standard_Real pf,pl,tolerance;
1552     if (!NewParameter(Vf,edg,pf,tolerance)) {
1553       pf = BRep_Tool::Parameter(Vf,edg);
1554     }
1555     if (!NewParameter(Vl,edg,pl,tolerance)) {
1556       pl = BRep_Tool::Parameter(Vl,edg);
1557     }
1558     if (pl <= pf) {
1559       //      const Handle(Geom_Curve) gc=ite.Value().Geometry();
1560       //      if (!gc.IsNull()) {
1561       //        pl = gc->LastParameter();
1562       //        pf = gc->FirstParameter();
1563       //      }
1564       Handle( Geom_Curve ) theCurve = myEMap.FindFromKey(edg).Geometry();
1565       if (theCurve->IsClosed())
1566       {
1567         // pf >= pl
1568         Standard_Real FirstPar = theCurve->FirstParameter(), LastPar = theCurve->LastParameter();
1569         Standard_Real pconf = Precision::PConfusion();
1570         if (Abs( pf - LastPar ) <= pconf)
1571           pf = FirstPar;
1572         else if (Abs( pl - FirstPar ) <= pconf)
1573           pl = LastPar;
1574
1575         if(pl <= pf) {
1576           pl += (LastPar-FirstPar);
1577         }
1578
1579       }
1580       if (pl <= pf) {
1581         errStat = Draft_EdgeRecomputation;
1582         badShape = edg;
1583         return;
1584       }
1585     }
1586     if (myVMap.Contains( Vf ))
1587       myVMap.ChangeFromKey(Vf).ChangeParameter(edg) = pf;
1588     if (myVMap.Contains( Vl ))
1589       myVMap.ChangeFromKey(Vl).ChangeParameter(edg) = pl;
1590   }
1591 }
1592
1593
1594
1595 //=======================================================================
1596 //function : NewSurface
1597 //purpose  :
1598 //=======================================================================
1599
1600 Handle(Geom_Surface) Draft_Modification::NewSurface
1601   (const Handle(Geom_Surface)& S,
1602   const TopAbs_Orientation Oris,
1603   const gp_Dir& Direction,
1604   const Standard_Real Angle,
1605   const gp_Pln& NeutralPlane)
1606 {
1607   Handle(Geom_Surface) NewS;
1608
1609   Handle(Standard_Type) TypeS = S->DynamicType();
1610
1611   if (TypeS == STANDARD_TYPE(Geom_Plane)) {
1612     gp_Pln Pl = Handle(Geom_Plane)::DownCast(S)->Pln();
1613     gp_Ax1 Axe;
1614     Standard_Real Theta;
1615     if (FindRotation(Pl,Oris,Direction,Angle,NeutralPlane,Axe,Theta)) {
1616       if ( Abs(Theta) > Precision::Angular()) {
1617         NewS = Handle(Geom_Surface)::DownCast(S->Rotated(Axe,Theta));
1618       }
1619       else {
1620         NewS = S;
1621       }
1622     }
1623   }
1624   else if (TypeS == STANDARD_TYPE(Geom_CylindricalSurface)) {
1625     Standard_Real testdir = Direction.Dot(NeutralPlane.Axis().Direction());
1626     if (Abs(testdir) <= 1.-Precision::Angular()) {
1627 #ifdef OCCT_DEBUG
1628       cout << "NewSurfaceCyl:Draft_Direction_and_Neutral_Perpendicular" << endl;
1629 #endif
1630       return NewS;
1631     }
1632     gp_Cylinder Cy = Handle(Geom_CylindricalSurface)::DownCast(S)->Cylinder();
1633     testdir = Direction.Dot(Cy.Axis().Direction());
1634     if (Abs(testdir) <= 1.-Precision::Angular()) {
1635 #ifdef OCCT_DEBUG
1636       cout << "NewSurfaceCyl:Draft_Direction_and_Cylinder_Perpendicular" << endl;
1637 #endif
1638       return NewS;
1639     }
1640     if (Abs(Angle) > Precision::Angular())
1641     {
1642       IntAna_QuadQuadGeo i2s;
1643       i2s.Perform(NeutralPlane,Cy,Precision::Angular(),Precision::Confusion());
1644       Standard_Boolean isIntDone = i2s.IsDone();
1645
1646       if(i2s.TypeInter() == IntAna_Ellipse)
1647       {
1648         const gp_Elips anEl = i2s.Ellipse(1);
1649         const Standard_Real aMajorR = anEl.MajorRadius();
1650         const Standard_Real aMinorR = anEl.MinorRadius();
1651         isIntDone = (aMajorR < 100000.0 * aMinorR);
1652       }
1653
1654       if (!isIntDone || i2s.TypeInter() != IntAna_Circle) {
1655 #ifdef OCCT_DEBUG
1656         cout << "NewSurfaceCyl:Draft_Intersection_Neutral_Cylinder_NotDone" << endl;
1657 #endif
1658         return NewS;
1659       }
1660       gp_Ax3 axcone = Cy.Position();
1661       // Pb : Where is the material???
1662       Standard_Real alpha = Angle;
1663       Standard_Boolean direct(axcone.Direct());
1664       if ((direct && Oris == TopAbs_REVERSED) ||
1665          (!direct && Oris == TopAbs_FORWARD)) {
1666         alpha = -alpha;
1667       }
1668
1669       gp_Pnt Center = i2s.Circle(1).Location();
1670       if (testdir <0.) {
1671         alpha = -alpha;
1672       }
1673       Standard_Real Z = ElCLib::LineParameter(Cy.Axis(),Center);
1674       Standard_Real Rad = Cy.Radius()+Z*Tan(alpha);
1675       if (Rad < 0.) {
1676         Rad = -Rad;
1677       }
1678       else {
1679         alpha = -alpha;
1680       }
1681       gp_Cone co(axcone,alpha,Rad);
1682       NewS = new Geom_ConicalSurface(co);
1683     }
1684     else {
1685       NewS = S;
1686     }
1687   }
1688   else if (TypeS == STANDARD_TYPE(Geom_ConicalSurface)) {
1689
1690     Standard_Real testdir = Direction.Dot(NeutralPlane.Axis().Direction());
1691     if (Abs(testdir) <= 1.-Precision::Angular()) {
1692 #ifdef OCCT_DEBUG
1693       cout << "NewSurfaceCone:Draft_Direction_and_Neutral_Perpendicular" << endl;
1694 #endif
1695       return NewS;
1696     }
1697
1698     gp_Cone Co1 = Handle(Geom_ConicalSurface)::DownCast(S)->Cone();
1699
1700     testdir = Direction.Dot(Co1.Axis().Direction());
1701     if (Abs(testdir) <= 1.-Precision::Angular()) {
1702 #ifdef OCCT_DEBUG
1703       cout << "NewSurfaceCone:Draft_Direction_and_Cone_Perpendicular" << endl;
1704 #endif
1705       return NewS;
1706     }
1707
1708
1709     IntAna_QuadQuadGeo i2s;
1710     i2s.Perform(NeutralPlane,Co1,Precision::Angular(),Precision::Confusion());
1711     if (!i2s.IsDone() || i2s.TypeInter() != IntAna_Circle) {
1712 #ifdef OCCT_DEBUG
1713       cout << "NewSurfaceCone:Draft_Intersection_Neutral_Conical_NotDone" << endl;
1714 #endif
1715       return NewS;
1716     }
1717     gp_Ax3 axcone = Co1.Position();
1718     // Pb : Where is the material???
1719     Standard_Real alpha = Angle;
1720     Standard_Boolean direct(axcone.Direct());
1721     if ((direct && Oris == TopAbs_REVERSED) ||
1722        (!direct && Oris == TopAbs_FORWARD)) {
1723       alpha = -alpha;
1724     }
1725
1726     gp_Pnt Center = i2s.Circle(1).Location();
1727     if (Abs(Angle) > Precision::Angular()) {
1728       if (testdir <0.) {
1729         alpha = -alpha;
1730       }
1731       Standard_Real Z = ElCLib::LineParameter(Co1.Axis(),Center);
1732       Standard_Real Rad = i2s.Circle(1).Radius()+Z*Tan(alpha);
1733       if (Rad < 0.) {
1734         Rad = -Rad;
1735       }
1736       else {
1737         alpha = -alpha;
1738       }
1739       if (Abs(alpha-Co1.SemiAngle()) < Precision::Angular()) {
1740         NewS = S;
1741       }
1742       else {
1743         gp_Cone co(axcone,alpha,Rad);
1744         NewS = new Geom_ConicalSurface(co);
1745       }
1746     }
1747     else {
1748       NewS = new
1749         Geom_CylindricalSurface(gp_Cylinder(axcone,i2s.Circle(1).Radius()));
1750     }
1751   }
1752   else {
1753 #ifdef OCCT_DEBUG
1754     cout << "NewSurface:Draft_SurfNotYetImplemented" << endl;
1755 #endif
1756   }
1757   return NewS;
1758 }
1759
1760
1761 //=======================================================================
1762 //function : NewCurve
1763 //purpose  :
1764 //=======================================================================
1765
1766 Handle(Geom_Curve) Draft_Modification::NewCurve
1767   (const Handle(Geom_Curve)& C,
1768   const Handle(Geom_Surface)& S,
1769   const TopAbs_Orientation Oris,
1770   const gp_Dir& Direction,
1771   const Standard_Real Angle,
1772   const gp_Pln& NeutralPlane,
1773   const Standard_Boolean )
1774
1775 {
1776   Handle(Geom_Curve) NewC;
1777
1778   Handle(Standard_Type) TypeS = S->DynamicType();
1779
1780   if (TypeS == STANDARD_TYPE(Geom_Plane)) {
1781     gp_Pln Pl = Handle(Geom_Plane)::DownCast(S)->Pln();
1782     gp_Ax1 Axe;
1783     Standard_Real Theta;
1784     if (FindRotation(Pl,Oris,Direction,Angle,NeutralPlane,Axe,Theta)) {
1785       if ( Abs(Theta) > Precision::Angular()) {
1786         NewC = Handle(Geom_Curve)::DownCast(C->Rotated(Axe,Theta));
1787       }
1788       else {
1789         NewC = C;
1790       }
1791     }
1792     return NewC;
1793   }
1794
1795
1796   if (C->DynamicType() != STANDARD_TYPE(Geom_Line)) {
1797     return NewC;
1798   }
1799
1800
1801   gp_Lin lin = Handle(Geom_Line)::DownCast(C)->Lin();
1802   //  Standard_Real testdir = Direction.Dot(lin.Direction());
1803   //  if (Abs(testdir) <= 1.-Precision::Angular()) {
1804   //    return NewC;
1805   //  }
1806   gp_Dir Norm;
1807   if (TypeS == STANDARD_TYPE(Geom_CylindricalSurface)) {
1808     Standard_Real U,V;
1809     gp_Vec d1u,d1v;
1810     gp_Pnt pbid;
1811     gp_Cylinder Cy = Handle(Geom_CylindricalSurface)::DownCast(S)->Cylinder();
1812     ElSLib::Parameters(Cy,lin.Location(),U,V);
1813     ElSLib::D1(U,V,Cy,pbid,d1u,d1v);
1814     Norm = d1u.Crossed(d1v);
1815   }
1816   else if (TypeS == STANDARD_TYPE(Geom_ConicalSurface)) {
1817     Standard_Real U,V;
1818     gp_Vec d1u,d1v;
1819     gp_Pnt pbid;
1820     gp_Cone Co = Handle(Geom_ConicalSurface)::DownCast(S)->Cone();
1821     ElSLib::Parameters(Co,lin.Location(),U,V);
1822     ElSLib::D1(U,V,Co,pbid,d1u,d1v);
1823     Norm = d1u.Crossed(d1v);
1824   }
1825
1826   IntAna_IntConicQuad ilipl(lin,NeutralPlane,Precision::Angular());
1827   if (ilipl.IsDone() && ilipl.NbPoints() != 0){
1828     if (Oris == TopAbs_REVERSED) {
1829       Norm.Reverse();
1830     }
1831     gp_Ax1 axrot(ilipl.Point(1), Norm.Crossed(Direction));
1832     gp_Lin lires = gp_Lin(gp_Ax1(ilipl.Point(1),Direction)).
1833       Rotated(axrot,Angle);
1834     if (lires.Direction().Dot(lin.Direction()) < 0.) {
1835       lires.Reverse();
1836     }
1837     NewC = new Geom_Line(lires);
1838   }
1839   return NewC;
1840 }
1841
1842
1843 //=======================================================================
1844 //function : Choose
1845 //purpose  :
1846 //=======================================================================
1847
1848 static Standard_Boolean Choose(const Draft_IndexedDataMapOfFaceFaceInfo& theFMap,
1849   Draft_IndexedDataMapOfEdgeEdgeInfo& theEMap,
1850   const TopoDS_Vertex& Vtx,
1851   Draft_VertexInfo& Vinf,
1852   GeomAdaptor_Curve& AC,
1853   GeomAdaptor_Surface& AS)
1854 {
1855   gp_Vec tgref;
1856   Vinf.InitEdgeIterator();
1857
1858   // Find a regular edge with null SecondFace
1859   while (Vinf.MoreEdge()) {
1860     const TopoDS_Edge& E1 = Vinf.Edge();
1861     const Draft_EdgeInfo& Einf1 = theEMap.FindFromKey(E1);
1862     if (Einf1.SecondFace().IsNull()) {
1863       break;
1864     }
1865     else {
1866       GeomAbs_Shape te = BRep_Tool::Continuity(E1,Einf1.FirstFace(),
1867         Einf1.SecondFace());
1868       if (te >= GeomAbs_G1) {
1869         break;
1870       }
1871     }
1872     Vinf.NextEdge();
1873   }
1874   if (!Vinf.MoreEdge()) { // take the first edge
1875     Vinf.InitEdgeIterator();
1876   }
1877
1878   const TopoDS_Edge& Eref = Vinf.Edge();
1879   //const Draft_EdgeInfo& Einf = theEMap(Eref);
1880   Draft_EdgeInfo& Einf = theEMap.ChangeFromKey(Eref);
1881
1882   AC.Load(Einf.Geometry());
1883
1884   Standard_Real f,l,prm;
1885   TopLoc_Location Loc;
1886   Handle(Geom_Curve) C = BRep_Tool::Curve(Eref,Loc,f,l);
1887   C = Handle(Geom_Curve)::DownCast(C->Transformed(Loc.Transformation()));
1888   gp_Pnt ptbid;
1889   //prm = Parameter(C,BRep_Tool::Pnt(Vtx));
1890   Standard_Integer done;
1891   Standard_Real param = Parameter( C, BRep_Tool::Pnt(Vtx), done );
1892   if (done != 0)
1893   {
1894     Handle( Geom_Surface ) S1 = theFMap.FindFromKey(Einf.FirstFace()).Geometry();
1895     Handle( Geom_Surface ) S2 = theFMap.FindFromKey(Einf.SecondFace()).Geometry();
1896     prm = SmartParameter( Einf, BRep_Tool::Tolerance(Eref), BRep_Tool::Pnt(Vtx), done, S1, S2 );
1897   }
1898   else
1899     prm = param;
1900   C->D1(prm,ptbid,tgref);
1901
1902
1903   Vinf.InitEdgeIterator();
1904   while (Vinf.MoreEdge()) {
1905     // Find a non tangent edge
1906     const TopoDS_Edge& Edg = Vinf.Edge();
1907     if (!Edg.IsSame(Eref)) {
1908       //const Draft_EdgeInfo& Einfo = theEMap(Edg);
1909       Draft_EdgeInfo& Einfo = theEMap.ChangeFromKey(Edg);
1910       if (!Einfo.SecondFace().IsNull() &&
1911           BRep_Tool::Continuity(Edg,Einfo.FirstFace(),Einfo.SecondFace())
1912           <= GeomAbs_C0) {
1913         C = BRep_Tool::Curve(Edg,Loc,f,l);
1914         C = Handle(Geom_Curve)::DownCast(C->Transformed(Loc.Transformation()));
1915         //prm = Parameter(C,BRep_Tool::Pnt(Vtx));
1916         Standard_Integer anewdone;
1917         Standard_Real anewparam = Parameter( C, BRep_Tool::Pnt(Vtx), anewdone );
1918         if (anewdone != 0)
1919         {
1920           Handle( Geom_Surface ) S1 = theFMap.FindFromKey(Einfo.FirstFace()).Geometry();
1921           Handle( Geom_Surface ) S2 = theFMap.FindFromKey(Einfo.SecondFace()).Geometry();
1922           prm = SmartParameter( Einfo, BRep_Tool::Tolerance(Edg), BRep_Tool::Pnt(Vtx), anewdone, S1, S2 );
1923         }
1924         else
1925           prm = anewparam;
1926         gp_Vec tg;
1927         C->D1(prm,ptbid,tg);
1928         if (tg.CrossMagnitude(tgref) > Precision::Confusion()) {
1929           break;
1930         }
1931       }
1932     }
1933     Vinf.NextEdge();
1934   }
1935   if (!Vinf.MoreEdge()) {
1936     return Standard_False;
1937   }
1938
1939   const Draft_EdgeInfo& Einf2 = theEMap.FindFromKey(Vinf.Edge());
1940   if (!Einf.SecondFace().IsNull()) {
1941
1942     if (Einf2.FirstFace().IsSame(Einf.FirstFace()) ||
1943         Einf2.FirstFace().IsSame(Einf.SecondFace())) {
1944       AS.Load(theFMap.FindFromKey(Einf2.SecondFace()).Geometry());
1945     }
1946     else {
1947       AS.Load(theFMap.FindFromKey(Einf2.FirstFace()).Geometry());
1948     }
1949   }
1950   else {
1951     if (Einf2.FirstFace().IsSame(Einf.FirstFace())) {
1952       AS.Load(theFMap.FindFromKey(Einf2.SecondFace()).Geometry());
1953     }
1954     else {
1955       AS.Load(theFMap.FindFromKey(Einf2.FirstFace()).Geometry());
1956     }
1957   }
1958   return Standard_True;
1959 }
1960
1961
1962 //=======================================================================
1963 //function : Parameter
1964 //purpose  :
1965 //=======================================================================
1966
1967 static Standard_Real Parameter(const Handle(Geom_Curve)& C,
1968   const gp_Pnt& P,
1969   Standard_Integer& done)
1970 {
1971   done = 0;
1972   Handle(Geom_Curve) cbase = C;
1973   Handle(Standard_Type) ctyp = C->DynamicType();
1974   if (ctyp == STANDARD_TYPE(Geom_TrimmedCurve)) {
1975     cbase = Handle(Geom_TrimmedCurve)::DownCast(C)->BasisCurve();
1976     ctyp = cbase->DynamicType();
1977   }
1978   Standard_Real param;
1979   if (ctyp == STANDARD_TYPE(Geom_Line)) {
1980     param = ElCLib::Parameter(Handle(Geom_Line)::DownCast(cbase)->Lin(),P);
1981   }
1982   else if (ctyp == STANDARD_TYPE(Geom_Circle)) {
1983     param = ElCLib::Parameter(Handle(Geom_Circle)::DownCast(cbase)->Circ(),P);
1984     if (Abs(2.*M_PI-param) <=Epsilon(2.*M_PI)) {
1985       param = 0.;
1986     }
1987   }
1988   else if (ctyp == STANDARD_TYPE(Geom_Ellipse)) {
1989     param = ElCLib::Parameter(Handle(Geom_Ellipse)::DownCast(cbase)->Elips(),P);
1990     if (Abs(2.*M_PI-param) <=Epsilon(2.*M_PI)) {
1991       param = 0.;
1992     }
1993   }
1994   else if (ctyp == STANDARD_TYPE(Geom_Parabola)) {
1995     param = ElCLib::Parameter(Handle(Geom_Parabola)::DownCast(cbase)->Parab(),P);
1996   }
1997   else if (ctyp == STANDARD_TYPE(Geom_Hyperbola)) {
1998     param = ElCLib::Parameter(Handle(Geom_Hyperbola)::DownCast(cbase)->Hypr(),P);
1999   }
2000   else {
2001     GeomAdaptor_Curve TheCurve(C);
2002     Extrema_ExtPC myExtPC(P,TheCurve);
2003     if (!myExtPC.IsDone()) {
2004       Standard_Failure::Raise("Draft_Modification_1::Parameter: ExtremaPC not done.");
2005     }
2006     if (myExtPC.NbExt() >= 1) {
2007       Standard_Real Dist2, Dist2Min = myExtPC.SquareDistance(1);
2008       Standard_Integer j, jmin = 1;
2009       for (j = 2; j <= myExtPC.NbExt(); j++) {
2010         Dist2 = myExtPC.SquareDistance(j);
2011         if (Dist2 < Dist2Min) {
2012           Dist2Min = Dist2;
2013           jmin = j;
2014         }
2015       }
2016       param = myExtPC.Point(jmin).Parameter();
2017     }
2018     else {
2019       Standard_Real dist1_2,dist2_2;
2020       gp_Pnt p1b,p2b;
2021       myExtPC.TrimmedSquareDistances(dist1_2,dist2_2,p1b,p2b);
2022       if (dist1_2 < dist2_2) {
2023         done = -1;
2024         param = TheCurve.FirstParameter();
2025       }
2026       else {
2027         done = 1;
2028         param = TheCurve.LastParameter();
2029       }
2030     }
2031
2032     if (cbase->IsPeriodic()) {
2033       Standard_Real Per  = cbase->Period();
2034       Standard_Real Tolp = Precision::Parametric(Precision::Confusion());
2035       if (Abs(Per-param) <= Tolp) {
2036         param = 0.;
2037       }
2038     }
2039   }
2040   return param;
2041 }
2042
2043 //=======================================================================
2044 //function : SmartParameter
2045 //purpose  :
2046 //=======================================================================
2047
2048 static Standard_Real SmartParameter(Draft_EdgeInfo& Einf,
2049   const Standard_Real EdgeTol,
2050   const gp_Pnt& Pnt,
2051   const Standard_Integer sign,
2052   const Handle(Geom_Surface)& S1,
2053   const Handle(Geom_Surface)& S2)
2054 {
2055   Handle( Geom2d_Curve ) NewC2d;
2056   Standard_Real Tol = Precision::Confusion();
2057   Standard_Real Etol = EdgeTol;
2058
2059   Handle( Geom2d_Curve ) pcu1 = Einf.FirstPC();
2060   Handle( Geom2d_Curve ) pcu2 = Einf.SecondPC();
2061
2062   if (pcu1.IsNull())
2063   {
2064     Handle( Geom_Curve ) theCurve = Einf.Geometry();
2065     pcu1 = GeomProjLib::Curve2d( theCurve, theCurve->FirstParameter(), theCurve->LastParameter(), S1, Etol );
2066     Einf.ChangeFirstPC() = pcu1;
2067   }
2068   if (pcu2.IsNull())
2069   {
2070     Handle( Geom_Curve ) theCurve = Einf.Geometry();
2071     pcu2 = GeomProjLib::Curve2d( theCurve, theCurve->FirstParameter(), theCurve->LastParameter(), S2, Etol );
2072     Einf.ChangeSecondPC() = pcu2;
2073   }
2074
2075   GeomAPI_ProjectPointOnSurf Projector( Pnt, S1 );
2076   Standard_Real U, V;
2077   Projector.LowerDistanceParameters( U, V );
2078
2079   NewC2d = Einf.FirstPC();
2080   if (NewC2d->DynamicType() == STANDARD_TYPE(Geom2d_TrimmedCurve))
2081     NewC2d = (Handle(Geom2d_TrimmedCurve)::DownCast(NewC2d))->BasisCurve();
2082
2083   gp_Pnt2d P2d( U, V );
2084   Geom2dAPI_ProjectPointOnCurve Projector2d( P2d, NewC2d );
2085   if (Projector2d.NbPoints() == 0 || Projector2d.LowerDistance() > Tol)
2086   {
2087     Handle( Geom2d_BSplineCurve ) BCurve;
2088     if (NewC2d->DynamicType() != STANDARD_TYPE(Geom2d_BSplineCurve))
2089       BCurve = Geom2dConvert::CurveToBSplineCurve( NewC2d );
2090     else
2091       BCurve = Handle( Geom2d_BSplineCurve )::DownCast( NewC2d );
2092     if (sign == -1)
2093     {
2094       TColgp_Array1OfPnt2d PntArray( 1, 2 );
2095       PntArray(1) = P2d;
2096       PntArray(2) = BCurve->Pole(1);
2097       Handle( Geom2d_BezierCurve ) Patch = new Geom2d_BezierCurve( PntArray );
2098       Geom2dConvert_CompCurveToBSplineCurve Concat( BCurve, Convert_QuasiAngular );
2099       Concat.Add( Patch, Tol, Standard_False );
2100       BCurve = Concat.BSplineCurve();
2101     }
2102     else
2103     {
2104       TColgp_Array1OfPnt2d PntArray( 1, 2 );
2105       PntArray(1) = BCurve->Pole( BCurve->NbPoles() );
2106       PntArray(2) = P2d;
2107       Handle( Geom2d_BezierCurve ) Patch = new Geom2d_BezierCurve( PntArray );
2108       Geom2dConvert_CompCurveToBSplineCurve Concat( BCurve, Convert_QuasiAngular );
2109       Concat.Add( Patch, Tol, Standard_True );
2110       BCurve = Concat.BSplineCurve();
2111     }
2112     NewC2d = BCurve;
2113   }
2114   Einf.ChangeFirstPC() = NewC2d;
2115   Handle( Geom2dAdaptor_HCurve ) hcur = new Geom2dAdaptor_HCurve( NewC2d );
2116   Handle( GeomAdaptor_HSurface ) hsur = new GeomAdaptor_HSurface( S1 );
2117   Adaptor3d_CurveOnSurface cons( hcur, hsur );
2118   Handle( Adaptor3d_HCurveOnSurface ) hcons = new Adaptor3d_HCurveOnSurface( cons );
2119   Handle( GeomAdaptor_HSurface ) hsur2 = new GeomAdaptor_HSurface( S2 );
2120   ProjLib_CompProjectedCurve ProjCurve( hsur2, hcons, Tol, Tol );
2121   Handle(ProjLib_HCompProjectedCurve) HProjector = new ProjLib_HCompProjectedCurve();
2122   HProjector->Set( ProjCurve );
2123   Standard_Real Udeb, Ufin;
2124   ProjCurve.Bounds(1, Udeb, Ufin);
2125   Standard_Integer MaxSeg = 20 + HProjector->NbIntervals(GeomAbs_C3);
2126   Approx_CurveOnSurface appr( HProjector, hsur2, Udeb, Ufin, Tol,
2127     GeomAbs_C1, 10, MaxSeg,
2128     Standard_False, Standard_False );
2129   Einf.ChangeSecondPC() = appr.Curve2d();
2130   Einf.ChangeGeometry() = appr.Curve3d();
2131   Einf.SetNewGeometry( Standard_True );
2132
2133   if (sign == -1)
2134     return Einf.Geometry()->FirstParameter();
2135   else
2136     return Einf.Geometry()->LastParameter();
2137
2138 }
2139
2140 //=======================================================================
2141 //function : Orientation
2142 //purpose  :
2143 //=======================================================================
2144
2145 static TopAbs_Orientation Orientation(const TopoDS_Shape& S,
2146   const TopoDS_Face& F)
2147 {
2148   //
2149   // change porting NT
2150   //
2151   TopExp_Explorer expl ;
2152   expl.Init(S,
2153     TopAbs_FACE) ;
2154   while (expl.More()) {
2155     if (TopoDS::Face(expl.Current()).IsSame(F)) {
2156       return expl.Current().Orientation();
2157     }
2158     expl.Next();
2159   }
2160   return TopAbs_FORWARD;
2161 }
2162
2163
2164 //=======================================================================
2165 //function : FindRotation
2166 //purpose  :
2167 //=======================================================================
2168
2169 static Standard_Boolean FindRotation(const gp_Pln& Pl,
2170   const TopAbs_Orientation Oris,
2171   const gp_Dir& Direction,
2172   const Standard_Real Angle,
2173   const gp_Pln& NeutralPlane,
2174   gp_Ax1& Axe,
2175   Standard_Real& theta)
2176 {
2177   IntAna_QuadQuadGeo i2pl(Pl,NeutralPlane,
2178     Precision::Angular(),Precision::Confusion());
2179
2180   if (i2pl.IsDone() && i2pl.TypeInter() == IntAna_Line) {
2181     gp_Lin li = i2pl.Line(1);
2182     // Try to turn around this line
2183     gp_Dir nx = li.Direction();
2184     gp_Dir ny = Pl.Axis().Direction().Crossed(nx);
2185     Standard_Real a = Direction.Dot(nx);
2186     if (Abs(a) <=1-Precision::Angular()) {
2187       Standard_Real b = Direction.Dot(ny);
2188       Standard_Real c = Direction.Dot(Pl.Axis().Direction());
2189       Standard_Boolean direct(Pl.Position().Direct());
2190       if ((direct && Oris == TopAbs_REVERSED) ||
2191         (!direct && Oris == TopAbs_FORWARD)) {
2192           b = -b;
2193           c = -c;
2194       }
2195       Standard_Real denom = Sqrt(1-a*a);
2196       Standard_Real Sina = Sin(Angle);
2197       if (denom>Abs(Sina)) {
2198         Standard_Real phi = ATan2(b/denom,c/denom);
2199         Standard_Real theta0 = ACos(Sina/denom);
2200         theta = theta0 - phi;
2201         if (Cos(theta) <0.) {
2202           theta = -theta0 -phi;
2203         }
2204         //  modified by NIZHNY-EAP Tue Nov 16 15:51:38 1999 ___BEGIN___
2205         while (Abs(theta)>M_PI) {
2206           theta = theta + M_PI*(theta<0 ? 1 : -1);
2207         }
2208         //  modified by NIZHNY-EAP Tue Nov 16 15:53:32 1999 ___END___
2209         Axe = li.Position();
2210         return Standard_True;
2211       }
2212     }
2213   }
2214   return Standard_False;
2215 }
2216