0025715: Intersection between cylinders produces excess vertices
[occt.git] / src / IntPatch / IntPatch_Intersection.cxx
1 // Created by: Modelization
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 #include <IntPatch_Intersection.ixx>
16
17 #include <IntPatch_ALineToWLine.hxx>
18 #include <IntPatch_GLine.hxx>
19 #include <IntPatch_ALine.hxx>
20 #include <IntPatch_WLine.hxx>
21 #include <IntPatch_RLine.hxx>
22 #include <IntPatch_PrmPrmIntersection.hxx>
23 #include <IntPatch_ImpPrmIntersection.hxx>
24 #include <IntPatch_ImpImpIntersection.hxx>
25 #include <IntSurf_Quadric.hxx>
26
27 #include <stdio.h>
28
29 #define DEBUG 0 
30 static const Standard_Integer aNbPointsInALine = 200;
31
32 //=======================================================================
33 //function : IsSeamOrBound
34 //purpose  : Returns TRUE if point thePt1 lies in seam-edge
35 //            (if it exists) or surface boundaries;
36 //=======================================================================
37 static Standard_Boolean IsSeamOrBound(const IntSurf_PntOn2S& thePt1,
38                                       const Standard_Real theU1Period,
39                                       const Standard_Real theU2Period,
40                                       const Standard_Real theV1Period,
41                                       const Standard_Real theV2Period,
42                                       const Standard_Real theUfSurf1,
43                                       const Standard_Real theUlSurf1,
44                                       const Standard_Real theVfSurf1,
45                                       const Standard_Real theVlSurf1,
46                                       const Standard_Real theUfSurf2,
47                                       const Standard_Real theUlSurf2,
48                                       const Standard_Real theVfSurf2,
49                                       const Standard_Real theVlSurf2)
50 {
51   Standard_Real aU11 = 0.0, aU12 = 0.0, aV11 = 0.0, aV12 = 0.0;
52   thePt1.Parameters(aU11, aV11, aU12, aV12);
53
54   Standard_Boolean aCond = Standard_False;
55   aCond = aCond || (!IsEqual(theU1Period, 0.0) &&
56                     IsEqual(fmod(aU11, theU1Period), 0.0));
57
58   aCond = aCond || (!IsEqual(theU2Period, 0.0) &&
59                     IsEqual(fmod(aU12, theU2Period), 0.0));
60
61   aCond = aCond || (!IsEqual(theV1Period, 0.0) &&
62                     IsEqual(fmod(aV11, theV1Period), 0.0));
63
64   aCond = aCond || (!IsEqual(theV2Period, 0.0) &&
65                     IsEqual(fmod(aV12, theV2Period), 0.0));
66
67   return  aCond ||
68           IsEqual(aU11, theUfSurf1) || IsEqual(aU11, theUlSurf1) ||
69           IsEqual(aU12, theUfSurf2) || IsEqual(aU12, theUlSurf2) ||
70           IsEqual(aV11, theVfSurf1) || IsEqual(aV11, theVlSurf1) ||
71           IsEqual(aV12, theVfSurf2) || IsEqual(aV12, theVlSurf2);
72 }
73
74 //=======================================================================
75 //function : JoinWLines
76 //purpose  : joins all WLines from theSlin to one if it is possible and
77 //            records the result into theSlin again.
78 //            Lines will be kept to be splitted if:
79 //              a) they are separated (has no common points);
80 //              b) resulted line (after joining) go through
81 //                 seam-edges or surface boundaries.
82 //
83 //          In addition, if points in theSPnt lies at least in one of 
84 //          the line in theSlin, this point will be deleted.
85 //=======================================================================
86 static void JoinWLines(IntPatch_SequenceOfLine& theSlin,
87                 IntPatch_SequenceOfPoint& theSPnt,
88                 const Standard_Real theTol3D,
89                 const Standard_Real theU1Period,
90                 const Standard_Real theU2Period,
91                 const Standard_Real theV1Period,
92                 const Standard_Real theV2Period,
93                 const Standard_Real theUfSurf1,
94                 const Standard_Real theUlSurf1,
95                 const Standard_Real theVfSurf1,
96                 const Standard_Real theVlSurf1,
97                 const Standard_Real theUfSurf2,
98                 const Standard_Real theUlSurf2,
99                 const Standard_Real theVfSurf2,
100                 const Standard_Real theVlSurf2)
101 {
102   if(theSlin.Length() == 0)
103     return;
104
105   for(Standard_Integer aNumOfLine1 = 1; aNumOfLine1 <= theSlin.Length(); aNumOfLine1++)
106   {
107     const Handle(IntPatch_WLine)& aWLine1 = Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNumOfLine1));
108
109     if(aWLine1.IsNull())
110     {//We must have failure to join not-point-lines
111       return;
112     }
113
114     const Standard_Integer aNbPntsWL1 = aWLine1->NbPnts();
115     const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
116     const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aNbPntsWL1);
117
118     for(Standard_Integer aNPt = 1; aNPt <= theSPnt.Length(); aNPt++)
119     {
120       const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNPt).PntOn2S();
121
122       if( aPntCur.IsSame(aPntFWL1, Precision::Confusion()) ||
123         aPntCur.IsSame(aPntLWL1, Precision::Confusion()))
124       {
125         theSPnt.Remove(aNPt);
126         aNPt--;
127       }
128     }
129
130     Standard_Boolean hasBeenRemoved = Standard_False;
131     for(Standard_Integer aNumOfLine2 = aNumOfLine1 + 1; aNumOfLine2 <= theSlin.Length(); aNumOfLine2++)
132     {
133       const Handle(IntPatch_WLine)& aWLine2 = Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNumOfLine2));
134
135       const Standard_Integer aNbPntsWL1 = aWLine1->NbPnts();
136       const Standard_Integer aNbPntsWL2 = aWLine2->NbPnts();
137
138       const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
139       const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aNbPntsWL1);
140
141       const IntSurf_PntOn2S& aPntFWL2 = aWLine2->Point(1);
142       const IntSurf_PntOn2S& aPntLWL2 = aWLine2->Point(aNbPntsWL2);
143
144       if(aPntFWL1.IsSame(aPntFWL2, Precision::Confusion()))
145       {
146         if(!IsSeamOrBound(aPntFWL1, theU1Period, theU2Period,
147                           theV1Period, theV2Period, theUfSurf1, theUlSurf1,
148                           theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
149                           theVfSurf2, theVlSurf2))
150         {
151           aWLine1->ClearVertexes();
152           for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++)
153           {
154             const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt);
155             aWLine1->Curve()->InsertBefore(1, aPt);
156           }
157
158           aWLine1->ComputeVertexParameters(theTol3D);
159
160           theSlin.Remove(aNumOfLine2);
161           aNumOfLine2--;
162           hasBeenRemoved = Standard_True;
163
164           continue;
165         }
166       }
167
168       if(aPntFWL1.IsSame(aPntLWL2, Precision::Confusion()))
169       {
170         if(!IsSeamOrBound(aPntFWL1, theU1Period, theU2Period,
171                           theV1Period, theV2Period, theUfSurf1, theUlSurf1,
172                           theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
173                           theVfSurf2, theVlSurf2))
174         {
175           aWLine1->ClearVertexes();
176           for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--)
177           {
178             const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt);
179             aWLine1->Curve()->InsertBefore(1, aPt);
180           }
181
182           aWLine1->ComputeVertexParameters(theTol3D);
183
184           theSlin.Remove(aNumOfLine2);
185           aNumOfLine2--;
186           hasBeenRemoved = Standard_True;
187
188           continue;
189         }
190       }
191
192       if(aPntLWL1.IsSame(aPntFWL2, Precision::Confusion()))
193       {
194         if(!IsSeamOrBound(aPntLWL1, theU1Period, theU2Period,
195                           theV1Period, theV2Period, theUfSurf1, theUlSurf1,
196                           theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
197                           theVfSurf2, theVlSurf2))
198         {
199           aWLine1->ClearVertexes();
200           for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++)
201           {
202             const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt);
203             aWLine1->Curve()->Add(aPt);
204           }
205
206           aWLine1->ComputeVertexParameters(theTol3D);
207
208           theSlin.Remove(aNumOfLine2);
209           aNumOfLine2--;
210           hasBeenRemoved = Standard_True;
211
212           continue;
213         }
214       }
215
216       if(aPntLWL1.IsSame(aPntLWL2, Precision::Confusion()))
217       {
218         if(!IsSeamOrBound(aPntLWL1, theU1Period, theU2Period,
219                           theV1Period, theV2Period, theUfSurf1, theUlSurf1,
220                           theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
221                           theVfSurf2, theVlSurf2))
222         {
223           aWLine1->ClearVertexes();
224           for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--)
225           {
226             const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt);
227             aWLine1->Curve()->Add(aPt);
228           }
229
230           aWLine1->ComputeVertexParameters(theTol3D);
231
232           theSlin.Remove(aNumOfLine2);
233           aNumOfLine2--;
234           hasBeenRemoved = Standard_True;
235
236           continue;
237         }
238       }
239     }
240
241     if(hasBeenRemoved)
242       aNumOfLine1--;
243   }
244 }
245
246 //======================================================================
247 // function: SequenceOfLine
248 //======================================================================
249 const IntPatch_SequenceOfLine& IntPatch_Intersection::SequenceOfLine() const { return(slin); }
250
251 //======================================================================
252 // function: IntPatch_Intersection
253 //======================================================================
254 IntPatch_Intersection::IntPatch_Intersection ()
255  : done(Standard_False),
256    //empt, tgte, oppo,
257    myTolArc(0.0), myTolTang(0.0),
258    myUVMaxStep(0.0), myFleche(0.0),
259    myIsStartPnt(Standard_False)
260    //myU1Start, myV1Start, myU2Start, myV2Start
261 {
262 }
263
264 //======================================================================
265 // function: IntPatch_Intersection
266 //======================================================================
267 IntPatch_Intersection::IntPatch_Intersection(const Handle(Adaptor3d_HSurface)&  S1,
268                                              const Handle(Adaptor3d_TopolTool)& D1,
269                                              const Handle(Adaptor3d_HSurface)&  S2,
270                                              const Handle(Adaptor3d_TopolTool)& D2,
271                                              const Standard_Real TolArc,
272                                              const Standard_Real TolTang)
273  : done(Standard_False),
274    //empt, tgte, oppo,
275    myTolArc(TolArc), myTolTang(TolTang),
276    myUVMaxStep(0.0), myFleche(0.0),
277    myIsStartPnt(Standard_False)
278    //myU1Start, myV1Start, myU2Start, myV2Start
279 {
280   if(myTolArc<1e-8) myTolArc=1e-8;
281   if(myTolTang<1e-8) myTolTang=1e-8;
282   if(myTolArc>0.5) myTolArc=0.5;
283   if(myTolTang>0.5) myTolTang=0.5;
284   Perform(S1,D1,S2,D2,TolArc,TolTang);
285 }
286
287 //======================================================================
288 // function: IntPatch_Intersection
289 //======================================================================
290 IntPatch_Intersection::IntPatch_Intersection(const Handle(Adaptor3d_HSurface)&  S1,
291                                              const Handle(Adaptor3d_TopolTool)& D1,
292                                              const Standard_Real TolArc,
293                                              const Standard_Real TolTang)
294  : done(Standard_False),
295    //empt, tgte, oppo,
296    myTolArc(TolArc), myTolTang(TolTang),
297    myUVMaxStep(0.0), myFleche(0.0),
298    myIsStartPnt(Standard_False)
299    //myU1Start, myV1Start, myU2Start, myV2Start
300 {
301   Perform(S1,D1,TolArc,TolTang);
302 }
303
304 //======================================================================
305 // function: SetTolerances
306 //======================================================================
307 void IntPatch_Intersection::SetTolerances(const Standard_Real TolArc,
308                                           const Standard_Real TolTang,
309                                           const Standard_Real UVMaxStep,
310                                           const Standard_Real Fleche)
311
312   myTolArc     = TolArc;
313   myTolTang    = TolTang;
314   myUVMaxStep  = UVMaxStep;
315   myFleche     = Fleche;
316   if(myTolArc<1e-8) myTolArc=1e-8;
317   if(myTolTang<1e-8) myTolTang=1e-8;
318   if(myTolArc>0.5) myTolArc=0.5;
319   if(myTolTang>0.5) myTolTang=0.5;  
320   if(myFleche<1.0e-3) myFleche=1e-3;
321   if(myUVMaxStep<1.0e-3) myUVMaxStep=1e-3;
322   if(myFleche>10) myFleche=10;
323   if(myUVMaxStep>0.5) myUVMaxStep=0.5;
324 }
325
326 //======================================================================
327 // function: Perform
328 //======================================================================
329 void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
330                                     const Handle(Adaptor3d_TopolTool)& D1,
331                                     const Standard_Real TolArc,
332                                     const Standard_Real TolTang)
333 {
334   myTolArc = TolArc;
335   myTolTang = TolTang;
336   if(myFleche == 0.0)  myFleche = 0.01;
337   if(myUVMaxStep==0.0) myUVMaxStep = 0.01;
338
339   done = Standard_True;
340   spnt.Clear();
341   slin.Clear();
342   
343   empt = Standard_True;
344   tgte = Standard_False;
345   oppo = Standard_False;
346
347   switch (S1->GetType())
348   { 
349     case GeomAbs_Plane:
350     case GeomAbs_Cylinder:
351     case GeomAbs_Sphere:
352     case GeomAbs_Cone:
353     case GeomAbs_Torus: break;
354     default:
355     {
356       IntPatch_PrmPrmIntersection interpp;
357       interpp.Perform(S1,D1,TolArc,TolTang,myFleche,myUVMaxStep);
358       if (interpp.IsDone())
359       {
360         done = Standard_True;
361         tgte = Standard_False;
362         empt = interpp.IsEmpty();
363         const Standard_Integer nblm = interpp.NbLines();
364         for (Standard_Integer i=1; i<=nblm; i++) slin.Append(interpp.Line(i));
365       }
366     }
367     break;
368   }
369 }
370
371 /////////////////////////////////////////////////////////////////////////////
372 //  These several support functions provide methods which can help basic   //
373 //  algorithm to intersect infinite surfaces of the following types:       //
374 //                                                                         //
375 //  a.) SurfaceOfExtrusion;                                                //
376 //  b.) SurfaceOfRevolution;                                               //
377 //  c.) OffsetSurface.                                                     //
378 //                                                                         //
379 /////////////////////////////////////////////////////////////////////////////
380 #include <TColgp_Array1OfXYZ.hxx>
381 #include <TColgp_Array1OfPnt.hxx>
382 #include <TColgp_SequenceOfPnt.hxx>
383 #include <Extrema_ExtPS.hxx>
384 #include <Extrema_POnSurf.hxx>
385 #include <Geom2d_Curve.hxx>
386 #include <Geom2dAPI_InterCurveCurve.hxx>
387 #include <GeomAdaptor.hxx>
388 #include <GeomAdaptor_HCurve.hxx>
389 #include <GeomAdaptor_Curve.hxx>
390 #include <GeomAdaptor_Surface.hxx>
391 #include <GeomAdaptor_HSurface.hxx>
392 #include <Geom_Plane.hxx>
393 #include <ProjLib_ProjectOnPlane.hxx>
394 #include <GeomProjLib.hxx>
395 #include <ElCLib.hxx>
396 #include <Geom_TrimmedCurve.hxx>
397 #include <Geom_Surface.hxx>
398 #include <Geom_SurfaceOfLinearExtrusion.hxx>
399 #include <Geom_OffsetSurface.hxx>
400 #include <Geom_SurfaceOfRevolution.hxx>
401 #include <Geom_RectangularTrimmedSurface.hxx>
402
403 //===============================================================
404 //function: FUN_GetMinMaxXYZPnt
405 //===============================================================
406 static void FUN_GetMinMaxXYZPnt( const Handle(Adaptor3d_HSurface)& S,
407                                  gp_Pnt& pMin, gp_Pnt& pMax )
408 {
409   const Standard_Real DU = 0.25 * Abs(S->LastUParameter() - S->FirstUParameter());
410   const Standard_Real DV = 0.25 * Abs(S->LastVParameter() - S->FirstVParameter());
411   Standard_Real tMinXYZ = RealLast();
412   Standard_Real tMaxXYZ = -tMinXYZ;
413   gp_Pnt PUV, ptMax, ptMin;
414   for(Standard_Real U = S->FirstUParameter(); U <= S->LastUParameter(); U += DU)
415   {
416     for(Standard_Real V = S->FirstVParameter(); V <= S->LastVParameter(); V += DV)
417     {
418       S->D0(U,V,PUV);
419       const Standard_Real cXYZ = PUV.XYZ().Modulus();
420       if(cXYZ > tMaxXYZ) { tMaxXYZ = cXYZ; ptMax = PUV; }
421       if(cXYZ < tMinXYZ) { tMinXYZ = cXYZ; ptMin = PUV; }
422     }
423   }
424   pMin = ptMin;
425   pMax = ptMax;
426 }
427 //==========================================================================
428 //function: FUN_TrimInfSurf
429 //==========================================================================
430 static void FUN_TrimInfSurf(const gp_Pnt& Pmin,
431                             const gp_Pnt& Pmax,
432                             const Handle(Adaptor3d_HSurface)& InfSurf,
433                             const Standard_Real& AlternativeTrimPrm,
434                             Handle(Adaptor3d_HSurface)& TrimS)
435 {
436   Standard_Real TP = AlternativeTrimPrm;
437   Extrema_ExtPS ext1(Pmin, InfSurf->Surface(), 1.e-7, 1.e-7);
438   Extrema_ExtPS ext2(Pmax, InfSurf->Surface(), 1.e-7, 1.e-7);
439   if(ext1.IsDone() || ext2.IsDone())
440   {
441     Standard_Real Umax = -1.e+100, Umin = 1.e+100, Vmax = -1.e+100, Vmin = 1.e+100, cU, cV;
442     if(ext1.IsDone())
443     {
444       for(Standard_Integer i = 1; i <= ext1.NbExt(); i++)
445       {
446         const Extrema_POnSurf & pons = ext1.Point(i);
447         pons.Parameter(cU,cV);
448         if(cU > Umax) Umax = cU;
449         if(cU < Umin) Umin = cU;
450         if(cV > Vmax) Vmax = cV;
451         if(cV < Vmin) Vmin = cV;
452       }
453     }
454     if(ext2.IsDone())
455     {
456       for(Standard_Integer i = 1; i <= ext2.NbExt(); i++)
457       {
458         const Extrema_POnSurf & pons = ext2.Point(i);
459         pons.Parameter(cU,cV);
460         if(cU > Umax) Umax = cU;
461         if(cU < Umin) Umin = cU;
462         if(cV > Vmax) Vmax = cV;
463         if(cV < Vmin) Vmin = cV;
464       }
465     }
466     TP = Max(Abs(Umin),Max(Abs(Umax),Max(Abs(Vmin),Abs(Vmax))));
467   }
468   if(TP == 0.) { TrimS = InfSurf; return; }
469   else
470   {
471     const Standard_Boolean Uinf = Precision::IsNegativeInfinite(InfSurf->FirstUParameter()); 
472     const Standard_Boolean Usup = Precision::IsPositiveInfinite(InfSurf->LastUParameter());
473     const Standard_Boolean Vinf = Precision::IsNegativeInfinite(InfSurf->FirstVParameter()); 
474     const Standard_Boolean Vsup = Precision::IsPositiveInfinite(InfSurf->LastVParameter());
475     Handle(Adaptor3d_HSurface) TmpSS;
476     Standard_Integer IsTrimed = 0;
477     const Standard_Real tp = 1000.0 * TP;
478     if(Vinf && Vsup) { TrimS = InfSurf->VTrim(-tp, tp, 1.0e-7); IsTrimed = 1; }
479     if(Vinf && !Vsup){ TrimS = InfSurf->VTrim(-tp, InfSurf->LastVParameter(), 1.0e-7); IsTrimed = 1; }
480     if(!Vinf && Vsup){ TrimS = InfSurf->VTrim(InfSurf->FirstVParameter(), tp, 1.0e-7); IsTrimed = 1; }
481     if(IsTrimed)
482     {
483       TmpSS = TrimS;
484       if(Uinf && Usup)  TrimS = TmpSS->UTrim(-tp, tp, 1.0e-7);
485       if(Uinf && !Usup) TrimS = TmpSS->UTrim(-tp, InfSurf->LastUParameter(), 1.0e-7);
486       if(!Uinf && Usup) TrimS = TmpSS->UTrim(InfSurf->FirstUParameter(), tp, 1.0e-7);
487     }
488     else
489     {
490       if(Uinf && Usup)  TrimS = InfSurf->UTrim(-tp, tp, 1.0e-7);
491       if(Uinf && !Usup) TrimS = InfSurf->UTrim(-tp, InfSurf->LastUParameter(), 1.0e-7);
492       if(!Uinf && Usup) TrimS = InfSurf->UTrim(InfSurf->FirstUParameter(), tp, 1.0e-7);
493     }
494   }
495 }
496 //================================================================================
497 //function: FUN_GetUiso
498 //================================================================================
499 static void FUN_GetUiso(const Handle(Geom_Surface)& GS,
500                         const GeomAbs_SurfaceType&  T,
501                         const Standard_Real&        FirstV,
502                         const Standard_Real&        LastV,
503                         const Standard_Boolean&     IsVC,
504                         const Standard_Boolean&     IsVP,
505                         const Standard_Real&        U,
506                         Handle(Geom_Curve)&         I)
507 {
508   if(T !=  GeomAbs_OffsetSurface)
509   {
510     Handle(Geom_Curve) gc = GS->UIso(U);
511     if(IsVP && (FirstV == 0.0 && LastV == (2.*M_PI))) I = gc;
512     else
513     {
514       Handle(Geom_TrimmedCurve) gtc = new Geom_TrimmedCurve(gc,FirstV,LastV);
515       //szv:I = Handle(Geom_Curve)::DownCast(gtc);
516       I = gtc;
517     }
518   }
519   else//OffsetSurface
520   {
521     const Handle(Geom_OffsetSurface) gos = *(Handle(Geom_OffsetSurface)*)&GS;
522     const Handle(Geom_Surface) bs = gos->BasisSurface();
523     Handle(Geom_Curve) gcbs = bs->UIso(U);
524     GeomAdaptor_Curve gac(gcbs);
525     const GeomAbs_CurveType GACT = gac.GetType();
526     if(IsVP || IsVC || GACT == GeomAbs_BSplineCurve || GACT == GeomAbs_BezierCurve || Abs(LastV - FirstV) < 1.e+5)
527     {
528       Handle(Geom_Curve) gc = gos->UIso(U);
529       if(IsVP && (FirstV == 0.0 && LastV == (2*M_PI))) I = gc;
530       else
531       {
532         Handle(Geom_TrimmedCurve) gtc = new Geom_TrimmedCurve(gc,FirstV,LastV);
533         //szv:I = Handle(Geom_Curve)::DownCast(gtc);
534         I = gtc;
535       }
536     }
537     else//Offset Line, Parab, Hyperb
538     {
539       Standard_Real VmTr, VMTr;
540       if(GACT != GeomAbs_Hyperbola)
541       {
542         if(FirstV >= 0. && LastV >= 0.){ VmTr = FirstV; VMTr = ((LastV - FirstV) > 1.e+4) ? (FirstV + 1.e+4) : LastV; }
543         else if(FirstV < 0. && LastV < 0.){ VMTr = LastV; VmTr = ((FirstV - LastV) < -1.e+4) ? (LastV - 1.e+4) : FirstV; }
544         else { VmTr = (FirstV < -1.e+4) ? -1.e+4 : FirstV; VMTr = (LastV > 1.e+4) ? 1.e+4 : LastV; }
545       }
546       else//Hyperbola
547       {
548         if(FirstV >= 0. && LastV >= 0.)
549         {
550           if(FirstV > 4.) return;
551           VmTr = FirstV; VMTr = (LastV > 4.) ? 4. : LastV;
552         }
553         else if(FirstV < 0. && LastV < 0.)
554         {
555           if(LastV < -4.) return;
556           VMTr = LastV; VmTr = (FirstV < -4.) ? -4. : FirstV;
557         }
558         else { VmTr = (FirstV < -4.) ? -4. : FirstV; VMTr = (LastV > 4.) ? 4. : LastV; }
559       }
560       //Make trimmed surface
561       Handle(Geom_RectangularTrimmedSurface) rts = new Geom_RectangularTrimmedSurface(gos,VmTr,VMTr,Standard_True);
562       I = rts->UIso(U);
563     }
564   }
565 }
566 //================================================================================
567 //function: FUN_GetViso
568 //================================================================================
569 static void FUN_GetViso(const Handle(Geom_Surface)& GS,
570                         const GeomAbs_SurfaceType&  T,
571                         const Standard_Real&        FirstU,
572                         const Standard_Real&        LastU,
573                         const Standard_Boolean&     IsUC,
574                         const Standard_Boolean&     IsUP,
575                         const Standard_Real&        V,
576                         Handle(Geom_Curve)&         I)
577 {
578   if(T !=  GeomAbs_OffsetSurface)
579   {
580     Handle(Geom_Curve) gc = GS->VIso(V);
581     if(IsUP && (FirstU == 0.0 && LastU == (2*M_PI))) I = gc;
582     else
583     {
584       Handle(Geom_TrimmedCurve) gtc = new Geom_TrimmedCurve(gc,FirstU,LastU);
585       //szv:I = Handle(Geom_Curve)::DownCast(gtc);
586       I = gtc;
587     }
588   }
589   else//OffsetSurface
590   {
591     const Handle(Geom_OffsetSurface) gos = *(Handle(Geom_OffsetSurface)*)&GS;
592     const Handle(Geom_Surface) bs = gos->BasisSurface();
593     Handle(Geom_Curve) gcbs = bs->VIso(V);
594     GeomAdaptor_Curve gac(gcbs);
595     const GeomAbs_CurveType GACT = gac.GetType();
596     if(IsUP || IsUC || GACT == GeomAbs_BSplineCurve || GACT == GeomAbs_BezierCurve || Abs(LastU - FirstU) < 1.e+5)
597     {
598       Handle(Geom_Curve) gc = gos->VIso(V);
599       if(IsUP && (FirstU == 0.0 && LastU == (2*M_PI))) I = gc;
600       else
601       {
602         Handle(Geom_TrimmedCurve) gtc = new Geom_TrimmedCurve(gc,FirstU,LastU);
603         //szv:I = Handle(Geom_Curve)::DownCast(gtc);
604         I = gtc;
605       }
606     }
607     else//Offset Line, Parab, Hyperb
608     {
609       Standard_Real UmTr, UMTr;
610       if(GACT != GeomAbs_Hyperbola)
611       {
612         if(FirstU >= 0. && LastU >= 0.){ UmTr = FirstU; UMTr = ((LastU - FirstU) > 1.e+4) ? (FirstU + 1.e+4) : LastU; }
613         else if(FirstU < 0. && LastU < 0.){ UMTr = LastU; UmTr = ((FirstU - LastU) < -1.e+4) ? (LastU - 1.e+4) : FirstU; }
614         else { UmTr = (FirstU < -1.e+4) ? -1.e+4 : FirstU; UMTr = (LastU > 1.e+4) ? 1.e+4 : LastU; }
615       }
616       else//Hyperbola
617       {
618         if(FirstU >= 0. && LastU >= 0.)
619         {
620           if(FirstU > 4.) return;
621           UmTr = FirstU; UMTr = (LastU > 4.) ? 4. : LastU;
622         }
623         else if(FirstU < 0. && LastU < 0.)
624         {
625           if(LastU < -4.) return;
626           UMTr = LastU; UmTr = (FirstU < -4.) ? -4. : FirstU;
627         }
628         else { UmTr = (FirstU < -4.) ? -4. : FirstU; UMTr = (LastU > 4.) ? 4. : LastU; }
629       }
630       //Make trimmed surface
631       Handle(Geom_RectangularTrimmedSurface) rts = new Geom_RectangularTrimmedSurface(gos,UmTr,UMTr,Standard_True);
632       I = rts->VIso(V);
633     }
634   }
635 }
636 //================================================================================
637 //function: FUN_PL_Intersection
638 //================================================================================
639 static void FUN_PL_Intersection(const Handle(Adaptor3d_HSurface)& S1,
640                                 const GeomAbs_SurfaceType&        T1,
641                                 const Handle(Adaptor3d_HSurface)& S2,
642                                 const GeomAbs_SurfaceType&        T2,
643                                 Standard_Boolean&                 IsOk,
644                                 TColgp_SequenceOfPnt&             SP,
645                                 gp_Vec&                           DV)
646 {
647   IsOk = Standard_False;
648   // 1. Check: both surfaces have U(V)isos - lines.
649   DV = gp_Vec(0.,0.,1.);
650   Standard_Boolean isoS1isLine[2] = {0, 0};
651   Standard_Boolean isoS2isLine[2] = {0, 0};
652   Handle(Geom_Curve) C1, C2;
653   const GeomAdaptor_Surface & gas1 = *(GeomAdaptor_Surface*)(&(S1->Surface()));
654   const GeomAdaptor_Surface & gas2 = *(GeomAdaptor_Surface*)(&(S2->Surface()));
655   const Handle(Geom_Surface) gs1 = gas1.Surface();
656   const Handle(Geom_Surface) gs2 = gas2.Surface();
657   Standard_Real MS1[2], MS2[2];
658   MS1[0] = 0.5 * (S1->LastUParameter() + S1->FirstUParameter());
659   MS1[1] = 0.5 * (S1->LastVParameter() + S1->FirstVParameter());
660   MS2[0] = 0.5 * (S2->LastUParameter() + S2->FirstUParameter());
661   MS2[1] = 0.5 * (S2->LastVParameter() + S2->FirstVParameter());
662   if(T1 == GeomAbs_SurfaceOfExtrusion) isoS1isLine[0] = Standard_True;
663   else if(!S1->IsVPeriodic() && !S1->IsVClosed()) {
664     if(T1 != GeomAbs_OffsetSurface) C1 = gs1->UIso(MS1[0]);
665     else {
666       const Handle(Geom_OffsetSurface) gos = *(Handle(Geom_OffsetSurface)*)&gs1;
667       const Handle(Geom_Surface) bs = gos->BasisSurface();
668       C1 = bs->UIso(MS1[0]);
669     }
670     GeomAdaptor_Curve gac(C1);
671     if(gac.GetType() == GeomAbs_Line) isoS1isLine[0] = Standard_True;
672   }
673   if(!S1->IsUPeriodic() && !S1->IsUClosed()) {
674     if(T1 != GeomAbs_OffsetSurface) C1 = gs1->VIso(MS1[1]);
675     else {
676       const Handle(Geom_OffsetSurface) gos = *(Handle(Geom_OffsetSurface)*)&gs1;
677       const Handle(Geom_Surface) bs = gos->BasisSurface();
678       C1 = bs->VIso(MS1[1]);
679     }
680     GeomAdaptor_Curve gac(C1);
681     if(gac.GetType() == GeomAbs_Line) isoS1isLine[1] = Standard_True;
682   }
683   if(T2 == GeomAbs_SurfaceOfExtrusion) isoS2isLine[0] = Standard_True;
684   else if(!S2->IsVPeriodic() && !S2->IsVClosed()) {
685     if(T2 != GeomAbs_OffsetSurface) C2 = gs2->UIso(MS2[0]);
686     else {
687       const Handle(Geom_OffsetSurface) gos = *(Handle(Geom_OffsetSurface)*)&gs2;
688       const Handle(Geom_Surface) bs = gos->BasisSurface();
689       C2 = bs->UIso(MS2[0]);
690     }
691     GeomAdaptor_Curve gac(C2);
692     if(gac.GetType() == GeomAbs_Line) isoS2isLine[0] = Standard_True;
693   }
694   if(!S2->IsUPeriodic() && !S2->IsUClosed()) {
695     if(T2 != GeomAbs_OffsetSurface) C2 = gs2->VIso(MS2[1]);
696     else {
697       const Handle(Geom_OffsetSurface) gos = *(Handle(Geom_OffsetSurface)*)&gs2;
698       const Handle(Geom_Surface) bs = gos->BasisSurface();
699       C2 = bs->VIso(MS2[1]);
700     }
701     GeomAdaptor_Curve gac(C2);
702     if(gac.GetType() == GeomAbs_Line) isoS2isLine[1] = Standard_True;
703   }
704   Standard_Boolean IsBothLines = ((isoS1isLine[0] || isoS1isLine[1]) &&
705                                   (isoS2isLine[0] || isoS2isLine[1]));
706   if(!IsBothLines){
707     return;
708   }
709   // 2. Check: Uiso lines of both surfaces are collinear.
710   gp_Pnt puvS1, puvS2;
711   gp_Vec derS1[2], derS2[2];
712   S1->D1(MS1[0], MS1[1], puvS1, derS1[0], derS1[1]);
713   S2->D1(MS2[0], MS2[1], puvS2, derS2[0], derS2[1]);
714   C1.Nullify(); C2.Nullify();
715   Standard_Integer iso = 0;
716   if(isoS1isLine[0] && isoS2isLine[0] &&
717      derS1[1].IsParallel(derS2[1],Precision::Angular())) {
718     iso = 1;
719     FUN_GetViso(gs1,T1,S1->FirstUParameter(),S1->LastUParameter(),
720                 S1->IsUClosed(),S1->IsUPeriodic(),MS1[1],C1);
721     FUN_GetViso(gs2,T2,S2->FirstUParameter(),S2->LastUParameter(),
722                 S2->IsUClosed(),S2->IsUPeriodic(),MS2[1],C2);
723   }
724   else if(isoS1isLine[0] && isoS2isLine[1] &&
725           derS1[1].IsParallel(derS2[0],Precision::Angular())) {
726     iso = 1;
727     FUN_GetViso(gs1,T1,S1->FirstUParameter(),S1->LastUParameter(),
728                 S1->IsUClosed(),S1->IsUPeriodic(),MS1[1],C1);
729     FUN_GetUiso(gs2,T2,S2->FirstVParameter(),S2->LastVParameter(),
730                 S2->IsVClosed(),S2->IsVPeriodic(),MS2[0],C2);
731   }
732   else if(isoS1isLine[1] && isoS2isLine[0] &&
733           derS1[0].IsParallel(derS2[1],Precision::Angular())) {
734     iso = 0;
735     FUN_GetUiso(gs1,T1,S1->FirstVParameter(),S1->LastVParameter(),
736                 S1->IsVClosed(),S1->IsVPeriodic(),MS1[0],C1);
737     FUN_GetViso(gs2,T2,S2->FirstUParameter(),S2->LastUParameter(),
738                 S2->IsUClosed(),S2->IsUPeriodic(),MS2[1],C2);
739   }
740   else if(isoS1isLine[1] && isoS2isLine[1] &&
741           derS1[0].IsParallel(derS2[0],Precision::Angular())) {
742     iso = 0;
743     FUN_GetUiso(gs1,T1,S1->FirstVParameter(),S1->LastVParameter(),
744                 S1->IsVClosed(),S1->IsVPeriodic(),MS1[0],C1);
745     FUN_GetUiso(gs2,T2,S2->FirstVParameter(),S2->LastVParameter(),
746                 S2->IsVClosed(),S2->IsVPeriodic(),MS2[0],C2);
747   }
748   else {
749     IsOk = Standard_False;
750     return;
751   }
752   IsOk = Standard_True;
753   // 3. Make intersections of V(U)isos
754   if(C1.IsNull() || C2.IsNull()) return;
755   DV = derS1[iso];
756   Handle(Geom_Plane) GPln = new Geom_Plane(gp_Pln(puvS1,gp_Dir(DV)));
757   Handle(Geom_Curve) C1Prj =
758     GeomProjLib::ProjectOnPlane(C1,GPln,gp_Dir(DV),Standard_True);
759   Handle(Geom_Curve) C2Prj =
760     GeomProjLib::ProjectOnPlane(C2,GPln,gp_Dir(DV),Standard_True);
761   if(C1Prj.IsNull() || C2Prj.IsNull()) return;
762   Handle(Geom2d_Curve) C1Prj2d =
763     GeomProjLib::Curve2d(C1Prj,*(Handle(Geom_Surface) *)&GPln);
764   Handle(Geom2d_Curve) C2Prj2d =
765     GeomProjLib::Curve2d(C2Prj,*(Handle(Geom_Surface) *)&GPln);
766   Geom2dAPI_InterCurveCurve ICC(C1Prj2d,C2Prj2d,1.0e-7);
767   if(ICC.NbPoints() > 0 )
768   {
769     for(Standard_Integer ip = 1; ip <= ICC.NbPoints(); ip++)
770     {
771       gp_Pnt2d P = ICC.Point(ip);
772       gp_Pnt P3d = ElCLib::To3d(gp_Ax2(puvS1,gp_Dir(DV)),P);
773       SP.Append(P3d);
774     }
775   }
776 }
777 //================================================================================
778 //function: FUN_NewFirstLast
779 //================================================================================
780 static void FUN_NewFirstLast(const GeomAbs_CurveType& ga_ct,
781                              const Standard_Real&     Fst,
782                              const Standard_Real&     Lst,
783                              const Standard_Real&     TrVal,
784                              Standard_Real&           NewFst,
785                              Standard_Real&           NewLst,
786                              Standard_Boolean&        NeedTr)
787 {
788   NewFst = Fst; NewLst = Lst; NeedTr = Standard_False;
789   switch (ga_ct)
790   {
791     case GeomAbs_Line:
792     case GeomAbs_Parabola:
793     {
794       if(Abs(Lst - Fst) > TrVal)
795       {
796         if(Fst >= 0. && Lst >= 0.)
797         {
798           NewFst = Fst;
799           NewLst = ((Fst + TrVal) < Lst) ? (Fst + TrVal) : Lst;
800         }
801         if(Fst < 0. && Lst < 0.)
802         {
803           NewLst = Lst;
804           NewFst = ((Lst - TrVal) > Fst) ? (Lst - TrVal) : Fst;
805         }
806         else
807         {
808           NewFst = (Fst < -TrVal) ? -TrVal : Fst;
809           NewLst = (Lst > TrVal) ? TrVal : Lst;
810         }
811         NeedTr = Standard_True;
812       }
813       break;
814     }
815     case GeomAbs_Hyperbola:
816     {
817       if(Abs(Lst - Fst) > 10.)
818       { 
819         if(Fst >= 0. && Lst >= 0.)
820         {
821           if(Fst > 4.) return;
822           NewFst = Fst;
823           NewLst = (Lst > 4.) ? 4. : Lst;
824         }
825         if(Fst < 0. && Lst < 0.)
826         {
827           if(Lst < -4.) return;
828           NewLst = Lst;
829           NewFst = (Fst < -4.) ? -4. : Fst;
830         }
831         else
832         {
833           NewFst = (Fst < -4.) ? -4. : Fst;
834           NewLst = (Lst > 4.) ? 4. : Lst;
835         }
836         NeedTr = Standard_True;
837       }
838       break;
839     }
840     default:
841     break;
842   }
843 }
844 //================================================================================
845 //function: FUN_TrimBothSurf
846 //================================================================================
847 static void FUN_TrimBothSurf(const Handle(Adaptor3d_HSurface)& S1,
848                              const GeomAbs_SurfaceType&        T1,
849                              const Handle(Adaptor3d_HSurface)& S2,
850                              const GeomAbs_SurfaceType&        T2,
851                              const Standard_Real&              TV,
852                              Handle(Adaptor3d_HSurface)&       NS1,
853                              Handle(Adaptor3d_HSurface)&       NS2)
854 {
855   const GeomAdaptor_Surface & gas1 = *(GeomAdaptor_Surface*)(&(S1->Surface()));
856   const GeomAdaptor_Surface & gas2 = *(GeomAdaptor_Surface*)(&(S2->Surface()));
857   const Handle(Geom_Surface) gs1 = gas1.Surface();
858   const Handle(Geom_Surface) gs2 = gas2.Surface();
859   const Standard_Real UM1 = 0.5 * (S1->LastUParameter() + S1->FirstUParameter());
860   const Standard_Real UM2 = 0.5 * (S2->LastUParameter() + S2->FirstUParameter());
861   const Standard_Real VM1 = 0.5 * (S1->LastVParameter() + S1->FirstVParameter());
862   const Standard_Real VM2 = 0.5 * (S2->LastVParameter() + S2->FirstVParameter());
863   Handle(Geom_Curve) visoS1, visoS2, uisoS1, uisoS2;
864   if(T1 != GeomAbs_OffsetSurface){ visoS1 = gs1->VIso(VM1); uisoS1 = gs1->UIso(UM1); }
865   else
866   {
867     const Handle(Geom_OffsetSurface) gos = *(Handle(Geom_OffsetSurface)*)&gs1;
868     const Handle(Geom_Surface) bs = gos->BasisSurface();
869     visoS1 = bs->VIso(VM1); uisoS1 = bs->UIso(UM1);
870   }
871   if(T2 != GeomAbs_OffsetSurface){ visoS2 = gs2->VIso(VM2); uisoS2 = gs2->UIso(UM2); }
872   else
873   {
874     const Handle(Geom_OffsetSurface) gos = *(Handle(Geom_OffsetSurface)*)&gs2;
875     const Handle(Geom_Surface) bs = gos->BasisSurface();
876     visoS2 = bs->VIso(VM2); uisoS2 = bs->UIso(UM2);
877   }
878   if(uisoS1.IsNull() || uisoS2.IsNull() || visoS1.IsNull() || visoS2.IsNull()){ NS1 = S1; NS2 = S2; return; }
879   GeomAdaptor_Curve gau1(uisoS1);
880   GeomAdaptor_Curve gav1(visoS1);
881   GeomAdaptor_Curve gau2(uisoS2);
882   GeomAdaptor_Curve gav2(visoS2);
883   GeomAbs_CurveType GA_U1 = gau1.GetType();
884   GeomAbs_CurveType GA_V1 = gav1.GetType();
885   GeomAbs_CurveType GA_U2 = gau2.GetType();
886   GeomAbs_CurveType GA_V2 = gav2.GetType();
887   Standard_Boolean TrmU1 = Standard_False;
888   Standard_Boolean TrmV1 = Standard_False;
889   Standard_Boolean TrmU2 = Standard_False;
890   Standard_Boolean TrmV2 = Standard_False;
891   Standard_Real V1S1,V2S1,U1S1,U2S1, V1S2,V2S2,U1S2,U2S2;
892   FUN_NewFirstLast(GA_U1,S1->FirstVParameter(),S1->LastVParameter(),TV,V1S1,V2S1,TrmV1);
893   FUN_NewFirstLast(GA_V1,S1->FirstUParameter(),S1->LastUParameter(),TV,U1S1,U2S1,TrmU1);
894   FUN_NewFirstLast(GA_U2,S2->FirstVParameter(),S2->LastVParameter(),TV,V1S2,V2S2,TrmV2);
895   FUN_NewFirstLast(GA_V2,S2->FirstUParameter(),S2->LastUParameter(),TV,U1S2,U2S2,TrmU2);
896   if(TrmV1) NS1 = S1->VTrim(V1S1, V2S1, 1.0e-7);
897   if(TrmV2) NS2 = S2->VTrim(V1S2, V2S2, 1.0e-7);
898   if(TrmU1)
899   {
900     if(TrmV1)
901     {
902       Handle(Adaptor3d_HSurface) TS = NS1;
903       NS1 = TS->UTrim(U1S1, U2S1, 1.0e-7);
904     }
905     else NS1 = S1->UTrim(U1S1, U2S1, 1.0e-7);
906   }
907   if(TrmU2)
908   {
909     if(TrmV2)
910     {
911       Handle(Adaptor3d_HSurface) TS = NS2;
912       NS2 = TS->UTrim(U1S2, U2S2, 1.0e-7);
913     }
914     else NS2 = S2->UTrim(U1S2, U2S2, 1.0e-7);
915   }
916 }
917
918 //=======================================================================
919 //function : Perform
920 //purpose  : 
921 //=======================================================================
922 void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  theS1,
923                                     const Handle(Adaptor3d_TopolTool)& theD1,
924                                     const Handle(Adaptor3d_HSurface)&  theS2,
925                                     const Handle(Adaptor3d_TopolTool)& theD2,
926                                     const Standard_Real TolArc,
927                                     const Standard_Real TolTang,
928                                     const Standard_Boolean isGeomInt)
929 {
930   myTolArc = TolArc;
931   myTolTang = TolTang;
932   if(myFleche <= Precision::PConfusion())
933     myFleche = 0.01;
934   if(myUVMaxStep <= Precision::PConfusion())
935     myUVMaxStep = 0.01;
936
937   done = Standard_False;
938   spnt.Clear();
939   slin.Clear();
940   empt = Standard_True;
941   tgte = Standard_False;
942   oppo = Standard_False;
943
944   GeomAbs_SurfaceType typs1 = theS1->GetType();
945   GeomAbs_SurfaceType typs2 = theS2->GetType();
946   
947   //treatment of the cases with cone or torus
948   Standard_Boolean TreatAsBiParametric = Standard_False;
949   Standard_Integer bGeomGeom = 0;
950   //
951   if (typs1 == GeomAbs_Cone  || typs2 == GeomAbs_Cone ||
952       typs1 == GeomAbs_Torus || typs2 == GeomAbs_Torus) {
953     gp_Ax1 aCTAx, aGeomAx;
954     GeomAbs_SurfaceType aCTType;
955     Standard_Boolean bToCheck;
956     //
957     const Handle(Adaptor3d_HSurface)& aCTSurf = 
958       (typs1 == GeomAbs_Cone || typs1 == GeomAbs_Torus) ? theS1 : theS2;
959     const Handle(Adaptor3d_HSurface)& aGeomSurf = 
960       (typs1 == GeomAbs_Cone || typs1 == GeomAbs_Torus) ? theS2 : theS1;
961     //
962     aCTType = aCTSurf->GetType();
963     bToCheck = Standard_False;
964     //
965     if (typs1 == GeomAbs_Cone  || typs2 == GeomAbs_Cone) {
966       const gp_Cone aCon1 = (aCTType == GeomAbs_Cone) ? 
967         aCTSurf->Cone() : aGeomSurf->Cone();
968       Standard_Real a1 = Abs(aCon1.SemiAngle());
969       bToCheck = (a1 < 0.02) || (a1 > 1.55);
970       //
971       if (typs1 == typs2) {
972         const gp_Cone aCon2 = aGeomSurf->Cone();
973         Standard_Real a2 = Abs(aCon2.SemiAngle());
974         bToCheck = bToCheck || (a2 < 0.02) || (a2 > 1.55);
975         //
976         if (a1 > 1.55 && a2 > 1.55) {//quasi-planes: if same domain, treat as canonic
977           const gp_Ax1 A1 = aCon1.Axis(), A2 = aCon2.Axis();
978           if (A1.IsParallel(A2,Precision::Angular())) {
979             const gp_Pnt Apex1 = aCon1.Apex(), Apex2 = aCon2.Apex();
980             const gp_Pln Plan1( Apex1, A1.Direction() );
981             if (Plan1.Distance( Apex2 ) <= Precision::Confusion()) {
982               bToCheck = Standard_False;
983             }
984           }
985         }
986       }
987       //
988       TreatAsBiParametric = bToCheck;
989       if (aCTType == GeomAbs_Cone) {
990         aCTAx = aCon1.Axis();
991       }
992     }
993     //
994     if (typs1 == GeomAbs_Torus || typs2 == GeomAbs_Torus) {
995       const gp_Torus aTor1 = (aCTType == GeomAbs_Torus) ? 
996         aCTSurf->Torus() : aGeomSurf->Torus();
997       bToCheck = aTor1.MajorRadius() > aTor1.MinorRadius();
998       if (typs1 == typs2) {
999         const gp_Torus aTor2 = aGeomSurf->Torus();
1000         bToCheck = aTor2.MajorRadius() > aTor2.MinorRadius();
1001       }
1002       //
1003       if (aCTType == GeomAbs_Torus) {
1004         aCTAx = aTor1.Axis();
1005       }
1006     }
1007     //
1008     if (bToCheck) {
1009       const gp_Lin aL1(aCTAx);
1010       //
1011       switch (aGeomSurf->GetType()) {
1012       case GeomAbs_Plane: {
1013         aGeomAx = aGeomSurf->Plane().Axis();
1014         if (aCTType == GeomAbs_Cone) {
1015           bGeomGeom = 1;
1016           if (Abs(aCTSurf->Cone().SemiAngle()) < 0.02) {
1017             Standard_Real ps = Abs(aCTAx.Direction().Dot(aGeomAx.Direction()));
1018             if(ps < 0.015) {
1019               bGeomGeom = 0;
1020             }
1021           }
1022         }
1023         else {
1024           if (aCTAx.IsParallel(aGeomAx, Precision::Angular()) ||
1025               (aCTAx.IsNormal(aGeomAx, Precision::Angular()) && 
1026                (aGeomSurf->Plane().Distance(aCTAx.Location()) < Precision::Confusion()))) {
1027             bGeomGeom = 1;
1028           }
1029         }
1030         bToCheck = Standard_False;
1031         break;
1032       }
1033       case GeomAbs_Sphere: {
1034         if (aL1.Distance(aGeomSurf->Sphere().Location()) < Precision::Confusion()) {
1035           bGeomGeom = 1;
1036         }
1037         bToCheck = Standard_False;
1038         break;
1039       }
1040       case GeomAbs_Cylinder:
1041         aGeomAx = aGeomSurf->Cylinder().Axis();
1042         break;
1043       case GeomAbs_Cone: 
1044         aGeomAx = aGeomSurf->Cone().Axis();
1045         break;
1046       case GeomAbs_Torus: 
1047         aGeomAx = aGeomSurf->Torus().Axis();
1048         break;
1049       default: 
1050         bToCheck = Standard_False;
1051         break;
1052       }
1053       //
1054       if (bToCheck) {
1055         if (aCTAx.IsParallel(aGeomAx, Precision::Angular()) &&
1056             (aL1.Distance(aGeomAx.Location()) <= Precision::Confusion())) {
1057           bGeomGeom = 1;
1058         }
1059       }
1060       //
1061       if (bGeomGeom == 1) {
1062         TreatAsBiParametric = Standard_False;
1063       }
1064     }
1065   }
1066   //
1067
1068   if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite()) {
1069     TreatAsBiParametric= Standard_False;
1070   }
1071
1072 //  Modified by skv - Mon Sep 26 14:58:30 2005 Begin
1073 //   if(TreatAsBiParametric) { typs1 = typs2 = GeomAbs_BezierSurface; }
1074   if(TreatAsBiParametric)
1075   {
1076     if (typs1 == GeomAbs_Cone && typs2 == GeomAbs_Plane)
1077       typs1 = GeomAbs_BezierSurface; // Using Imp-Prm Intersector
1078     else if (typs1 == GeomAbs_Plane && typs2 == GeomAbs_Cone)
1079       typs2 = GeomAbs_BezierSurface; // Using Imp-Prm Intersector
1080     else {
1081       // Using Prm-Prm Intersector
1082       typs1 = GeomAbs_BezierSurface;
1083       typs2 = GeomAbs_BezierSurface;
1084     }
1085   }
1086 //  Modified by skv - Mon Sep 26 14:58:30 2005 End
1087
1088   // Surface type definition
1089   Standard_Integer ts1 = 0;
1090   switch (typs1)
1091   {
1092     case GeomAbs_Plane:
1093     case GeomAbs_Cylinder:
1094     case GeomAbs_Sphere:
1095     case GeomAbs_Cone: ts1 = 1; break;
1096     case GeomAbs_Torus: ts1 = bGeomGeom; break;
1097     default: break;
1098   }
1099
1100   Standard_Integer ts2 = 0;
1101   switch (typs2)
1102   {
1103     case GeomAbs_Plane:
1104     case GeomAbs_Cylinder:
1105     case GeomAbs_Sphere:
1106     case GeomAbs_Cone: ts2 = 1; break;
1107     case GeomAbs_Torus: ts2 = bGeomGeom; break;
1108     default: break;
1109   }
1110   //
1111   // treatment of the cases with torus and any other geom surface
1112   //
1113   // Possible intersection types: 1. ts1 == ts2 == 1 <Geom-Geom>
1114   //                              2. ts1 != ts2      <Geom-Param>
1115   //                              3. ts1 == ts2 == 0 <Param-Param>
1116
1117   // Geom - Geom
1118   if(ts1 == ts2 && ts1 == 1)
1119   {
1120     const Standard_Boolean RestrictLine = Standard_True;
1121     IntSurf_ListOfPntOn2S ListOfPnts;
1122     ListOfPnts.Clear();
1123     if(isGeomInt)
1124     {
1125       if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite())
1126       {
1127         GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, 
1128                       TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1129       }
1130       else
1131       {
1132         GeomGeomPerfomTrimSurf(theS1, theD1, theS2, theD2,
1133               TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1134       }
1135     }
1136     else
1137     {
1138       ParamParamPerfom(theS1, theD1, theS2, theD2, 
1139               TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1140     }
1141   }
1142
1143   // Geom - Param
1144   if(ts1 != ts2)
1145   {
1146     GeomParamPerfom(theS1, theD1, theS2, theD2, ts1 == 0, typs1, typs2);
1147   }
1148
1149   // Param - Param 
1150   if(ts1 == ts2 && ts1 == 0)
1151   {
1152     const Standard_Boolean RestrictLine = Standard_True;
1153     IntSurf_ListOfPntOn2S ListOfPnts;
1154     ListOfPnts.Clear();
1155
1156     ParamParamPerfom(theS1, theD1, theS2, theD2, TolArc,
1157                         TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1158   }
1159 }
1160
1161 //=======================================================================
1162 //function : Perform
1163 //purpose  : 
1164 //=======================================================================
1165 void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  theS1,
1166                                     const Handle(Adaptor3d_TopolTool)& theD1,
1167                                     const Handle(Adaptor3d_HSurface)&  theS2,
1168                                     const Handle(Adaptor3d_TopolTool)& theD2,
1169                                     const Standard_Real TolArc,
1170                                     const Standard_Real TolTang,
1171                                     IntSurf_ListOfPntOn2S& ListOfPnts,
1172                                     const Standard_Boolean RestrictLine,
1173                                     const Standard_Boolean isGeomInt)
1174 {
1175   myTolArc = TolArc;
1176   myTolTang = TolTang;
1177   if(myFleche <= Precision::PConfusion())
1178     myFleche = 0.01;
1179   if(myUVMaxStep <= Precision::PConfusion())
1180     myUVMaxStep = 0.01;
1181     
1182   done = Standard_False;
1183   spnt.Clear();
1184   slin.Clear();
1185   empt = Standard_True;
1186   tgte = Standard_False;
1187   oppo = Standard_False;
1188
1189   GeomAbs_SurfaceType typs1 = theS1->GetType();
1190   GeomAbs_SurfaceType typs2 = theS2->GetType();
1191   //
1192   //treatment of the cases with cone or torus
1193   Standard_Boolean TreatAsBiParametric = Standard_False;
1194   Standard_Integer bGeomGeom = 0;
1195   //
1196   if (typs1 == GeomAbs_Cone  || typs2 == GeomAbs_Cone ||
1197       typs1 == GeomAbs_Torus || typs2 == GeomAbs_Torus) {
1198     gp_Ax1 aCTAx, aGeomAx;
1199     GeomAbs_SurfaceType aCTType;
1200     Standard_Boolean bToCheck;
1201     //
1202     const Handle(Adaptor3d_HSurface)& aCTSurf = 
1203       (typs1 == GeomAbs_Cone || typs1 == GeomAbs_Torus) ? theS1 : theS2;
1204     const Handle(Adaptor3d_HSurface)& aGeomSurf = 
1205       (typs1 == GeomAbs_Cone || typs1 == GeomAbs_Torus) ? theS2 : theS1;
1206     //
1207     aCTType = aCTSurf->GetType();
1208     bToCheck = Standard_False;
1209     //
1210     if (typs1 == GeomAbs_Cone  || typs2 == GeomAbs_Cone) {
1211       const gp_Cone aCon1 = (aCTType == GeomAbs_Cone) ? 
1212         aCTSurf->Cone() : aGeomSurf->Cone();
1213       Standard_Real a1 = Abs(aCon1.SemiAngle());
1214       bToCheck = (a1 < 0.02) || (a1 > 1.55);
1215       //
1216       if (typs1 == typs2) {
1217         const gp_Cone aCon2 = aGeomSurf->Cone();
1218         Standard_Real a2 = Abs(aCon2.SemiAngle());
1219         bToCheck = bToCheck || (a2 < 0.02) || (a2 > 1.55);
1220         //
1221         if (a1 > 1.55 && a2 > 1.55) {//quasi-planes: if same domain, treat as canonic
1222           const gp_Ax1 A1 = aCon1.Axis(), A2 = aCon2.Axis();
1223           if (A1.IsParallel(A2,Precision::Angular())) {
1224             const gp_Pnt Apex1 = aCon1.Apex(), Apex2 = aCon2.Apex();
1225             const gp_Pln Plan1( Apex1, A1.Direction() );
1226             if (Plan1.Distance( Apex2 ) <= Precision::Confusion()) {
1227               bToCheck = Standard_False;
1228             }
1229           }
1230         }
1231       }
1232       //
1233       TreatAsBiParametric = bToCheck;
1234       if (aCTType == GeomAbs_Cone) {
1235         aCTAx = aCon1.Axis();
1236       }
1237     }
1238     //
1239     if (typs1 == GeomAbs_Torus || typs2 == GeomAbs_Torus) {
1240       const gp_Torus aTor1 = (aCTType == GeomAbs_Torus) ? 
1241         aCTSurf->Torus() : aGeomSurf->Torus();
1242       bToCheck = aTor1.MajorRadius() > aTor1.MinorRadius();
1243       if (typs1 == typs2) {
1244         const gp_Torus aTor2 = aGeomSurf->Torus();
1245         bToCheck = aTor2.MajorRadius() > aTor2.MinorRadius();
1246       }
1247       //
1248       if (aCTType == GeomAbs_Torus) {
1249         aCTAx = aTor1.Axis();
1250       }
1251     }
1252     //
1253     if (bToCheck) {
1254       const gp_Lin aL1(aCTAx);
1255       //
1256       switch (aGeomSurf->GetType()) {
1257       case GeomAbs_Plane: {
1258         aGeomAx = aGeomSurf->Plane().Axis();
1259         if (aCTType == GeomAbs_Cone) {
1260           bGeomGeom = 1;
1261           if (Abs(aCTSurf->Cone().SemiAngle()) < 0.02) {
1262             Standard_Real ps = Abs(aCTAx.Direction().Dot(aGeomAx.Direction()));
1263             if(ps < 0.015) {
1264               bGeomGeom = 0;
1265             }
1266           }
1267         }
1268         else {
1269           if (aCTAx.IsParallel(aGeomAx, Precision::Angular()) ||
1270               (aCTAx.IsNormal(aGeomAx, Precision::Angular()) && 
1271                (aGeomSurf->Plane().Distance(aCTAx.Location()) < Precision::Confusion()))) {
1272             bGeomGeom = 1;
1273           }
1274         }
1275         bToCheck = Standard_False;
1276         break;
1277       }
1278       case GeomAbs_Sphere: {
1279         if (aL1.Distance(aGeomSurf->Sphere().Location()) < Precision::Confusion()) {
1280           bGeomGeom = 1;
1281         }
1282         bToCheck = Standard_False;
1283         break;
1284       }
1285       case GeomAbs_Cylinder:
1286         aGeomAx = aGeomSurf->Cylinder().Axis();
1287         break;
1288       case GeomAbs_Cone: 
1289         aGeomAx = aGeomSurf->Cone().Axis();
1290         break;
1291       case GeomAbs_Torus: 
1292         aGeomAx = aGeomSurf->Torus().Axis();
1293         break;
1294       default: 
1295         bToCheck = Standard_False;
1296         break;
1297       }
1298       //
1299       if (bToCheck) {
1300         if (aCTAx.IsParallel(aGeomAx, Precision::Angular()) &&
1301             (aL1.Distance(aGeomAx.Location()) <= Precision::Confusion())) {
1302           bGeomGeom = 1;
1303         }
1304       }
1305       //
1306       if (bGeomGeom == 1) {
1307         TreatAsBiParametric = Standard_False;
1308       }
1309     }
1310   }
1311   //
1312
1313   if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite()) {
1314     TreatAsBiParametric= Standard_False;
1315   }
1316
1317   if(TreatAsBiParametric)
1318   {
1319     // Using Prm-Prm Intersector
1320     typs1 = GeomAbs_BezierSurface;
1321     typs2 = GeomAbs_BezierSurface;
1322   }
1323
1324   // Surface type definition
1325   Standard_Integer ts1 = 0;
1326   switch (typs1)
1327   {
1328     case GeomAbs_Plane:
1329     case GeomAbs_Cylinder:
1330     case GeomAbs_Sphere:
1331     case GeomAbs_Cone: ts1 = 1; break;
1332     case GeomAbs_Torus: ts1 = bGeomGeom; break;
1333     default: break;
1334   }
1335
1336   Standard_Integer ts2 = 0;
1337   switch (typs2)
1338   {
1339     case GeomAbs_Plane:
1340     case GeomAbs_Cylinder:
1341     case GeomAbs_Sphere:
1342     case GeomAbs_Cone: ts2 = 1; break;
1343     case GeomAbs_Torus: ts2 = bGeomGeom; break;
1344     default: break;
1345   }
1346   //
1347   // Possible intersection types: 1. ts1 == ts2 == 1 <Geom-Geom>
1348   //                              2. ts1 != ts2      <Geom-Param>
1349   //                              3. ts1 == ts2 == 0 <Param-Param>
1350
1351   if(!isGeomInt)
1352   {
1353     ParamParamPerfom(theS1, theD1, theS2, theD2, 
1354                 TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1355   }
1356   else if(ts1 != ts2)
1357   {
1358     GeomParamPerfom(theS1, theD1, theS2, theD2, ts1 == 0, typs1, typs2);
1359   }
1360   else if (ts1 == 0)
1361   {
1362     ParamParamPerfom(theS1, theD1, theS2, theD2,
1363                 TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1364   }
1365   else if(ts1 == 1)
1366   {
1367     if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite())
1368     {
1369       GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, 
1370                       TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1371     }
1372     else
1373     {
1374       GeomGeomPerfomTrimSurf(theS1, theD1, theS2, theD2,
1375               TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1376     }
1377   }
1378 }
1379
1380 //=======================================================================
1381 //function : ParamParamPerfom
1382 //purpose  : 
1383 //=======================================================================
1384 void IntPatch_Intersection::ParamParamPerfom(const Handle(Adaptor3d_HSurface)&  theS1,
1385                                              const Handle(Adaptor3d_TopolTool)& theD1,
1386                                              const Handle(Adaptor3d_HSurface)&  theS2,
1387                                              const Handle(Adaptor3d_TopolTool)& theD2,
1388                                              const Standard_Real TolArc,
1389                                              const Standard_Real TolTang,
1390                                              IntSurf_ListOfPntOn2S& ListOfPnts,
1391                                              const Standard_Boolean RestrictLine,
1392                                              const GeomAbs_SurfaceType typs1,
1393                                              const GeomAbs_SurfaceType typs2)
1394 {
1395   IntPatch_PrmPrmIntersection interpp;
1396   if(!theD1->DomainIsInfinite() && !theD2->DomainIsInfinite())
1397   {
1398     Standard_Boolean ClearFlag = Standard_True;
1399     if(!ListOfPnts.IsEmpty())
1400     {
1401       interpp.Perform(theS1,theD1,theS2,theD2,TolArc,TolTang,myFleche,myUVMaxStep, ListOfPnts, RestrictLine);
1402       ClearFlag = Standard_False;
1403     }
1404     interpp.Perform(theS1,theD1,theS2,theD2,TolArc,TolTang,myFleche,myUVMaxStep,ClearFlag);   //double call!!!!!!!
1405   }
1406   else if((theD1->DomainIsInfinite()) ^ (theD2->DomainIsInfinite()))
1407   {
1408     gp_Pnt pMaxXYZ, pMinXYZ;
1409     if(theD1->DomainIsInfinite())
1410     {
1411       FUN_GetMinMaxXYZPnt( theS2, pMinXYZ, pMaxXYZ );
1412       const Standard_Real MU = Max(Abs(theS2->FirstUParameter()),Abs(theS2->LastUParameter()));
1413       const Standard_Real MV = Max(Abs(theS2->FirstVParameter()),Abs(theS2->LastVParameter()));
1414       const Standard_Real AP = Max(MU, MV);
1415       Handle(Adaptor3d_HSurface) SS;
1416       FUN_TrimInfSurf(pMinXYZ, pMaxXYZ, theS1, AP, SS);
1417       interpp.Perform(SS,theD1,theS2,theD2,TolArc,TolTang,myFleche,myUVMaxStep);
1418     }
1419     else
1420     {
1421       FUN_GetMinMaxXYZPnt( theS1, pMinXYZ, pMaxXYZ );
1422       const Standard_Real MU = Max(Abs(theS1->FirstUParameter()),Abs(theS1->LastUParameter()));
1423       const Standard_Real MV = Max(Abs(theS1->FirstVParameter()),Abs(theS1->LastVParameter()));
1424       const Standard_Real AP = Max(MU, MV);
1425       Handle(Adaptor3d_HSurface) SS;
1426       FUN_TrimInfSurf(pMinXYZ, pMaxXYZ, theS2, AP, SS);
1427       interpp.Perform(theS1, theD1, SS, theD2,TolArc,TolTang,myFleche,myUVMaxStep);
1428     }
1429   }//(theD1->DomainIsInfinite()) ^ (theD2->DomainIsInfinite())
1430   else
1431   {
1432     if(typs1 == GeomAbs_OtherSurface || typs2 == GeomAbs_OtherSurface)
1433     {
1434       done = Standard_False;
1435       return;
1436     }
1437
1438     Standard_Boolean IsPLInt = Standard_False;
1439     TColgp_SequenceOfPnt sop;
1440     gp_Vec v;
1441     FUN_PL_Intersection(theS1,typs1,theS2,typs2,IsPLInt,sop,v);
1442
1443     if(IsPLInt)
1444     {
1445       if(sop.Length() > 0)
1446       {
1447         for(Standard_Integer ip = 1; ip <= sop.Length(); ip++)
1448         {
1449           gp_Lin lin(sop.Value(ip),gp_Dir(v));
1450           Handle(IntPatch_GLine) gl = new IntPatch_GLine(lin,Standard_False);
1451           slin.Append(*(Handle(IntPatch_Line) *)&gl);
1452         }
1453
1454         done = Standard_True;
1455       }
1456       else
1457         done = Standard_False;
1458
1459       return;
1460     }// 'COLLINEAR LINES'
1461     else
1462     {
1463       Handle(Adaptor3d_HSurface) nS1 = theS1;
1464       Handle(Adaptor3d_HSurface) nS2 = theS2;
1465       FUN_TrimBothSurf(theS1,typs1,theS2,typs2,1.e+8,nS1,nS2);
1466       interpp.Perform(nS1,theD1,nS2,theD2,TolArc,TolTang,myFleche,myUVMaxStep);
1467     }// 'NON - COLLINEAR LINES'
1468   }// both domains are infinite
1469
1470   if (interpp.IsDone())
1471   {
1472     done = Standard_True;
1473     tgte = Standard_False;
1474     empt = interpp.IsEmpty();
1475
1476     for(Standard_Integer i = 1; i <= interpp.NbLines(); i++)
1477     {
1478       if(interpp.Line(i)->ArcType() != IntPatch_Walking)
1479         slin.Append(interpp.Line(i));
1480     }
1481
1482     for (Standard_Integer i = 1; i <= interpp.NbLines(); i++)
1483     {
1484       if(interpp.Line(i)->ArcType() == IntPatch_Walking)
1485         slin.Append(interpp.Line(i));
1486     }
1487   }
1488 }
1489
1490 //=======================================================================
1491 ////function : GeomGeomPerfom
1492 //purpose  : 
1493 //=======================================================================
1494 void IntPatch_Intersection::GeomGeomPerfom(const Handle(Adaptor3d_HSurface)& theS1,
1495                                            const Handle(Adaptor3d_TopolTool)& theD1,
1496                                            const Handle(Adaptor3d_HSurface)& theS2,
1497                                            const Handle(Adaptor3d_TopolTool)& theD2,
1498                                            const Standard_Real TolArc,
1499                                            const Standard_Real TolTang,
1500                                            IntSurf_ListOfPntOn2S& ListOfPnts,
1501                                            const Standard_Boolean RestrictLine,
1502                                            const GeomAbs_SurfaceType typs1,
1503                                            const GeomAbs_SurfaceType typs2)
1504 {
1505   IntPatch_ImpImpIntersection interii(theS1,theD1,theS2,theD2,myTolArc,myTolTang);
1506   const Standard_Boolean anIS = interii.IsDone();
1507   if (anIS)
1508   {
1509     done = anIS;
1510     empt = interii.IsEmpty();
1511     if (!empt)
1512     {
1513       tgte = interii.TangentFaces();
1514       if (tgte)
1515         oppo = interii.OppositeFaces();
1516
1517       for (Standard_Integer i = 1; i <= interii.NbLines(); i++)
1518       {
1519         const Handle(IntPatch_Line)& line = interii.Line(i);
1520         if (line->ArcType() == IntPatch_Analytic)
1521         {
1522           const GeomAbs_SurfaceType typs1 = theS1->GetType();
1523           const GeomAbs_SurfaceType typs2 = theS2->GetType();
1524           IntSurf_Quadric Quad1,Quad2;
1525           
1526           switch(typs1)
1527           {
1528           case GeomAbs_Plane:
1529             Quad1.SetValue(theS1->Plane());
1530             break;
1531
1532           case GeomAbs_Cylinder:
1533             Quad1.SetValue(theS1->Cylinder());
1534             break;
1535
1536           case GeomAbs_Sphere:
1537             Quad1.SetValue(theS1->Sphere());
1538             break;
1539
1540           case GeomAbs_Cone:
1541             Quad1.SetValue(theS1->Cone());
1542             break;
1543
1544           case GeomAbs_Torus:
1545             Quad1.SetValue(theS1->Torus());
1546             break;
1547
1548           default:
1549             break;
1550           }
1551
1552           switch(typs2)
1553           {
1554           case GeomAbs_Plane:
1555             Quad2.SetValue(theS2->Plane());
1556             break;
1557           case GeomAbs_Cylinder:
1558             Quad2.SetValue(theS2->Cylinder());
1559             break;
1560
1561           case GeomAbs_Sphere:
1562             Quad2.SetValue(theS2->Sphere());
1563             break;
1564
1565           case GeomAbs_Cone:
1566             Quad2.SetValue(theS2->Cone());
1567             break;
1568
1569           case GeomAbs_Torus:
1570             Quad2.SetValue(theS2->Torus());
1571             break;
1572
1573           default:
1574             break;
1575           }
1576
1577           IntPatch_ALineToWLine AToW(Quad1,Quad2,0.01,0.05,aNbPointsInALine);
1578           Handle(IntPatch_Line) wlin=AToW.MakeWLine((*((Handle(IntPatch_ALine) *)(&line))));
1579           slin.Append(wlin);
1580         }
1581         else
1582           slin.Append(interii.Line(i));
1583       }
1584
1585       for (Standard_Integer i = 1; i <= interii.NbPnts(); i++)
1586       {
1587         spnt.Append(interii.Point(i));
1588       }
1589     }
1590   }
1591   else
1592     ParamParamPerfom(theS1, theD1, theS2, theD2, 
1593                 TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1594 }
1595
1596 //=======================================================================
1597 //function : GeomParamPerfom
1598 //purpose  : 
1599 //=======================================================================
1600 void IntPatch_Intersection::
1601   GeomParamPerfom(const Handle(Adaptor3d_HSurface)&  theS1,
1602                   const Handle(Adaptor3d_TopolTool)& theD1,
1603                   const Handle(Adaptor3d_HSurface)&  theS2,
1604                   const Handle(Adaptor3d_TopolTool)& theD2,
1605                   const Standard_Boolean isNotAnalitical,
1606                   const GeomAbs_SurfaceType typs1,
1607                   const GeomAbs_SurfaceType typs2)
1608 {
1609   IntPatch_ImpPrmIntersection interip;
1610   if (myIsStartPnt)
1611   {
1612     if (isNotAnalitical/*ts1 == 0*/)
1613       interip.SetStartPoint(myU1Start,myV1Start);
1614     else
1615       interip.SetStartPoint(myU2Start,myV2Start);
1616   }
1617
1618   if(theD1->DomainIsInfinite() && theD2->DomainIsInfinite())
1619   {
1620     Standard_Boolean IsPLInt = Standard_False;
1621     TColgp_SequenceOfPnt sop;
1622     gp_Vec v;
1623     FUN_PL_Intersection(theS1,typs1,theS2,typs2,IsPLInt,sop,v);
1624     
1625     if(IsPLInt)
1626     {
1627       if(sop.Length() > 0)
1628       {
1629         for(Standard_Integer ip = 1; ip <= sop.Length(); ip++)
1630         {
1631           gp_Lin lin(sop.Value(ip),gp_Dir(v));
1632           Handle(IntPatch_GLine) gl = new IntPatch_GLine(lin,Standard_False);
1633           slin.Append(*(Handle(IntPatch_Line) *)&gl);
1634         }
1635
1636         done = Standard_True;
1637       }
1638       else
1639         done = Standard_False;
1640
1641       return;
1642     }
1643     else
1644     {
1645       Handle(Adaptor3d_HSurface) nS1 = theS1;
1646       Handle(Adaptor3d_HSurface) nS2 = theS2;
1647       FUN_TrimBothSurf(theS1,typs1,theS2,typs2,1.e+5,nS1,nS2);
1648       interip.Perform(nS1,theD1,nS2,theD2,myTolArc,myTolTang,myFleche,myUVMaxStep);
1649     }
1650   }
1651   else
1652     interip.Perform(theS1,theD1,theS2,theD2,myTolArc,myTolTang,myFleche,myUVMaxStep);
1653
1654   if (interip.IsDone()) 
1655   {
1656     done = Standard_True;
1657     empt = interip.IsEmpty();
1658
1659     if (!empt)
1660     {
1661       const Standard_Integer aNbLines = interip.NbLines();
1662       for(Standard_Integer i = 1; i <= aNbLines; i++)
1663       {
1664         if(interip.Line(i)->ArcType() != IntPatch_Walking)
1665           slin.Append(interip.Line(i));
1666       }
1667
1668       for(Standard_Integer i = 1; i <= aNbLines; i++)
1669       {
1670         if(interip.Line(i)->ArcType() == IntPatch_Walking)
1671           slin.Append(interip.Line(i));
1672       }
1673
1674       for (Standard_Integer i = 1; i <= interip.NbPnts(); i++)
1675         spnt.Append(interip.Point(i));
1676     }
1677   }
1678 }
1679
1680 //=======================================================================
1681 //function : GeomGeomPerfomTrimSurf
1682 //purpose  : This function returns ready walking-line (which is not need
1683 //            in convertation) as an intersection line between two
1684 //            trimmed surfaces.
1685 //=======================================================================
1686 void IntPatch_Intersection::
1687   GeomGeomPerfomTrimSurf( const Handle(Adaptor3d_HSurface)& theS1,
1688                           const Handle(Adaptor3d_TopolTool)& theD1,
1689                           const Handle(Adaptor3d_HSurface)& theS2,
1690                           const Handle(Adaptor3d_TopolTool)& theD2,
1691                           const Standard_Real theTolArc,
1692                           const Standard_Real theTolTang,
1693                           IntSurf_ListOfPntOn2S& theListOfPnts,
1694                           const Standard_Boolean RestrictLine,
1695                           const GeomAbs_SurfaceType theTyps1,
1696                           const GeomAbs_SurfaceType theTyps2)
1697 {
1698   IntSurf_Quadric Quad1,Quad2;
1699
1700   if((theTyps1 == GeomAbs_Cylinder) && (theTyps2 == GeomAbs_Cylinder))
1701   {
1702     IntPatch_ImpImpIntersection anInt;
1703     anInt.Perform(theS1, theD1, theS2, theD2, myTolArc, myTolTang, Standard_True);
1704
1705     done = anInt.IsDone();
1706
1707     if(done)
1708     {
1709       empt = anInt.IsEmpty();
1710       if (!empt)
1711       {
1712         tgte = anInt.TangentFaces();
1713         if (tgte)
1714           oppo = anInt.OppositeFaces();
1715
1716         const Standard_Integer aNbLin = anInt.NbLines();
1717         const Standard_Integer aNbPts = anInt.NbPnts();
1718
1719         for(Standard_Integer aLID = 1; aLID <= aNbLin; aLID++)
1720         {
1721           const Handle(IntPatch_Line)& aLine = anInt.Line(aLID);
1722           slin.Append(aLine);
1723         }
1724
1725         for(Standard_Integer aPID = 1; aPID <= aNbPts; aPID++)
1726         {
1727           const IntPatch_Point& aPoint = anInt.Point(aPID);
1728           spnt.Append(aPoint);
1729         }
1730
1731         JoinWLines( slin, spnt, theTolTang,
1732                     theS1->IsUPeriodic()? theS1->UPeriod() : 0.0,
1733                     theS2->IsUPeriodic()? theS2->UPeriod() : 0.0,
1734                     theS1->IsVPeriodic()? theS1->VPeriod() : 0.0,
1735                     theS2->IsVPeriodic()? theS2->VPeriod() : 0.0,
1736                     theS1->FirstUParameter(), theS1->LastUParameter(),
1737                     theS1->FirstVParameter(), theS1->LastVParameter(),
1738                     theS2->FirstUParameter(), theS2->LastUParameter(),
1739                     theS2->FirstVParameter(), theS2->LastVParameter());
1740       }
1741     }
1742   }
1743   else
1744   {
1745     GeomGeomPerfom(theS1, theD1, theS2, theD2,
1746             theTolArc, theTolTang, theListOfPnts, RestrictLine, theTyps1, theTyps2);
1747   }
1748 }
1749
1750
1751 void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
1752                                     const Handle(Adaptor3d_TopolTool)& D1,
1753                                     const Handle(Adaptor3d_HSurface)&  S2,
1754                                     const Handle(Adaptor3d_TopolTool)& D2,
1755                                     const Standard_Real U1,
1756                                     const Standard_Real V1,
1757                                     const Standard_Real U2,
1758                                     const Standard_Real V2,
1759                                     const Standard_Real TolArc,
1760                                     const Standard_Real TolTang)
1761 {
1762   myTolArc = TolArc;
1763   myTolTang = TolTang;
1764   if(myFleche == 0.0) {
1765 #if DEBUG
1766     //cout<<" -- IntPatch_Intersection::myFleche fixe par defaut a 0.01 --"<<endl;
1767     //cout<<" -- Utiliser la Methode SetTolerances( ... ) "<<endl;
1768 #endif
1769     myFleche = 0.01;
1770   }
1771   if(myUVMaxStep==0.0) {
1772 #if DEBUG
1773     //cout<<" -- IntPatch_Intersection::myUVMaxStep fixe par defaut a 0.01 --"<<endl;
1774     //cout<<" -- Utiliser la Methode SetTolerances( ... ) "<<endl;
1775 #endif
1776     myUVMaxStep = 0.01;
1777   }
1778
1779   done = Standard_False;
1780   spnt.Clear();
1781   slin.Clear();
1782
1783   empt = Standard_True;
1784   tgte = Standard_False;
1785   oppo = Standard_False;
1786
1787   const GeomAbs_SurfaceType typs1 = S1->GetType();
1788   const GeomAbs_SurfaceType typs2 = S2->GetType();
1789   
1790   if(   typs1==GeomAbs_Plane 
1791      || typs1==GeomAbs_Cylinder
1792      || typs1==GeomAbs_Sphere
1793      || typs1==GeomAbs_Cone
1794      || typs2==GeomAbs_Plane 
1795      || typs2==GeomAbs_Cylinder
1796      || typs2==GeomAbs_Sphere
1797      || typs2==GeomAbs_Cone)
1798   {
1799     myIsStartPnt = Standard_True;
1800     myU1Start = U1; myV1Start = V1; myU2Start = U2; myV2Start = V2;
1801     Perform(S1,D1,S2,D2,TolArc,TolTang);
1802     myIsStartPnt = Standard_False;
1803   }
1804   else
1805   {
1806     IntPatch_PrmPrmIntersection interpp;
1807     interpp.Perform(S1,D1,S2,D2,U1,V1,U2,V2,TolArc,TolTang,myFleche,myUVMaxStep);
1808     if (interpp.IsDone())
1809     {
1810       done = Standard_True;
1811       tgte = Standard_False;
1812       empt = interpp.IsEmpty();
1813       const Standard_Integer nblm = interpp.NbLines();
1814       Standard_Integer i = 1;
1815       for (; i<=nblm; i++) slin.Append(interpp.Line(i));
1816     }
1817   }
1818 }
1819 //======================================================================
1820 #include <IntPatch_IType.hxx>
1821 #include <IntPatch_LineConstructor.hxx>
1822 #include <Adaptor2d_HCurve2d.hxx>
1823 #define MAXR 200
1824
1825
1826 //void IntPatch_Intersection__MAJ_R(Handle(Adaptor2d_HCurve2d) *R1,
1827 //                                     Handle(Adaptor2d_HCurve2d) *R2,
1828 //                                     int *NR1,
1829 //                                     int *NR2,
1830 //                                     Standard_Integer nbR1,
1831 //                                     Standard_Integer nbR2,
1832 //                                     const IntPatch_Point& VTX)
1833 void IntPatch_Intersection__MAJ_R(Handle(Adaptor2d_HCurve2d) *,
1834                                   Handle(Adaptor2d_HCurve2d) *,
1835                                   int *,
1836                                   int *,
1837                                   Standard_Integer ,
1838                                   Standard_Integer ,
1839                                   const IntPatch_Point& )
1840
1841   /*
1842   if(VTX.IsOnDomS1()) { 
1843     
1844     //-- long unsigned ptr= *((long unsigned *)(((Handle(Standard_Transient) *)(&(VTX.ArcOnS1())))));
1845     for(Standard_Integer i=0; i<nbR1;i++) { 
1846       if(VTX.ArcOnS1()==R1[i]) { 
1847         NR1[i]++;
1848         printf("\n ******************************");
1849         return;
1850       }
1851     }
1852     printf("\n R Pas trouvee  (IntPatch)\n");
1853     
1854   }
1855   */
1856 }
1857
1858
1859 //void IntPatch_Intersection::Dump(const Standard_Integer Mode,
1860 void IntPatch_Intersection::Dump(const Standard_Integer ,
1861                                  const Handle(Adaptor3d_HSurface)&  S1,
1862                                  const Handle(Adaptor3d_TopolTool)& D1,
1863                                  const Handle(Adaptor3d_HSurface)&  S2,
1864                                  const Handle(Adaptor3d_TopolTool)& D2) const 
1865
1866   
1867   //-- ----------------------------------------------------------------------
1868   //--  construction de la liste des restrictions & vertex 
1869   //--
1870   int NR1[MAXR],NR2[MAXR];
1871   Handle(Adaptor2d_HCurve2d) R1[MAXR],R2[MAXR];
1872   Standard_Integer nbR1=0,nbR2=0;
1873   for(D1->Init();D1->More() && nbR1<MAXR; D1->Next()) { 
1874     R1[nbR1]=D1->Value(); 
1875     NR1[nbR1]=0;
1876     nbR1++;
1877   }
1878   for(D2->Init();D2->More() && nbR2<MAXR; D2->Next()) { 
1879     R2[nbR2]=D2->Value();
1880     NR2[nbR2]=0;
1881     nbR2++;
1882   }
1883   
1884   printf("\nDUMP_INT:  ----empt:%2ud  tgte:%2ud  oppo:%2ud ---------------------------------",empt,tgte,empt);
1885   Standard_Integer i,j,nbr1,nbr2,nbgl,nbgc,nbge,nbgp,nbgh,nbl,nbr,nbg,nbw,nba;
1886   nbl=nbr=nbg=nbw=nba=nbgl=nbge=nbr1=nbr2=nbgc=nbgp=nbgh=0;
1887   nbl=NbLines();
1888   for(i=1;i<=nbl;i++) { 
1889     const Handle(IntPatch_Line)& line=Line(i);
1890     const IntPatch_IType IType=line->ArcType();
1891     if(IType == IntPatch_Walking) nbw++;
1892     else     if(IType == IntPatch_Restriction) { 
1893       nbr++;
1894       Handle(IntPatch_RLine)& rlin =
1895         *((Handle(IntPatch_RLine) *)&line);
1896       if(rlin->IsArcOnS1()) nbr1++;
1897       if(rlin->IsArcOnS2()) nbr2++;
1898     }
1899     else     if(IType == IntPatch_Analytic) nba++;
1900     else     { 
1901       nbg++; 
1902       if(IType == IntPatch_Lin) nbgl++;
1903       else if(IType == IntPatch_Circle) nbgc++;
1904       else if(IType == IntPatch_Parabola) nbgp++;
1905       else if(IType == IntPatch_Hyperbola) nbgh++;
1906       else if(IType == IntPatch_Ellipse) nbge++;
1907     }
1908   }
1909   
1910   
1911   printf("\nDUMP_INT:Lines:%2d Wlin:%2d Restr:%2d(On1:%2d On2:%2d) Ana:%2d Geom:%2d(L:%2d C:%2d E:%2d H:%2d P:%2d)",
1912          nbl,nbw,nbr,nbr1,nbr2,nba,nbg,nbgl,nbgc,nbge,nbgh,nbgp);
1913   
1914   IntPatch_LineConstructor LineConstructor(2);
1915   
1916   Standard_Integer nbllc=0;
1917   nbw=nbr=nbg=nba=0;
1918   Standard_Integer nbva,nbvw,nbvr,nbvg;
1919   nbva=nbvr=nbvw=nbvg=0;
1920   for (j=1; j<=nbl; j++) {
1921     Standard_Integer v,nbvtx;
1922     const Handle(IntPatch_Line)& intersLinej = Line(j);
1923     Standard_Integer NbLines;
1924     LineConstructor.Perform(SequenceOfLine(),intersLinej,S1,D1,S2,D2,1e-7);
1925     NbLines = LineConstructor.NbLines();
1926     
1927     for(Standard_Integer k=1;k<=NbLines;k++) { 
1928       nbllc++;
1929       const Handle(IntPatch_Line)& LineK = LineConstructor.Line(k);
1930       if (LineK->ArcType() == IntPatch_Analytic) { 
1931         Handle(IntPatch_ALine)& alin =
1932           *((Handle(IntPatch_ALine) *)&LineK);
1933         nbvtx=alin->NbVertex();
1934         nbva+=nbvtx;        nba++;
1935         for(v=1;v<=nbvtx;v++) { 
1936           IntPatch_Intersection__MAJ_R(R1,R2,NR1,NR2,nbR1,nbR2,alin->Vertex(v));
1937         }
1938       }
1939       else if (LineK->ArcType() == IntPatch_Restriction) {
1940         Handle(IntPatch_RLine)& rlin =
1941           *((Handle(IntPatch_RLine) *)&LineK);
1942         nbvtx=rlin->NbVertex();
1943         nbvr+=nbvtx;        nbr++;
1944         for(v=1;v<=nbvtx;v++) { 
1945           IntPatch_Intersection__MAJ_R(R1,R2,NR1,NR2,nbR1,nbR2,rlin->Vertex(v));
1946         }
1947       }
1948       else if (LineK->ArcType() == IntPatch_Walking) {
1949         Handle(IntPatch_WLine)& wlin =
1950           *((Handle(IntPatch_WLine) *)&LineK);
1951         nbvtx=wlin->NbVertex();
1952         nbvw+=nbvtx;        nbw++;
1953         for(v=1;v<=nbvtx;v++) { 
1954           IntPatch_Intersection__MAJ_R(R1,R2,NR1,NR2,nbR1,nbR2,wlin->Vertex(v));
1955         }
1956       }
1957       else { 
1958         Handle(IntPatch_GLine)& glin =
1959           *((Handle(IntPatch_GLine) *)&LineK);
1960         nbvtx=glin->NbVertex();
1961         nbvg+=nbvtx;        nbg++;
1962         for(v=1;v<=nbvtx;v++) { 
1963           IntPatch_Intersection__MAJ_R(R1,R2,NR1,NR2,nbR1,nbR2,glin->Vertex(v));
1964         }
1965       }
1966     }
1967   }
1968   printf("\nDUMP_LC :Lines:%2d WLin:%2d Restr:%2d Ana:%2d Geom:%2d",
1969          nbllc,nbw,nbr,nba,nbg);
1970   printf("\nDUMP_LC :vtx          :%2d     r:%2d    :%2d     :%2d",
1971          nbvw,nbvr,nbva,nbvg);
1972
1973
1974
1975    printf("\n");
1976 }