0024005: Intersecting a slightly off angle plane with a cylinder takes 7+ seconds
[occt.git] / src / IntPatch / IntPatch_Intersection.cxx
1 // Created by: Modelization
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
3 //
4 // The content of this file is subject to the Open CASCADE Technology Public
5 // License Version 6.5 (the "License"). You may not use the content of this file
6 // except in compliance with the License. Please obtain a copy of the License
7 // at http://www.opencascade.org and read it completely before using this file.
8 //
9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
11 //
12 // The Original Code and all software distributed under the License is
13 // distributed on an "AS IS" basis, without warranty of any kind, and the
14 // Initial Developer hereby disclaims all such warranties, including without
15 // limitation, any warranties of merchantability, fitness for a particular
16 // purpose or non-infringement. Please see the License for the specific terms
17 // and conditions governing the rights and limitations under the License.
18
19
20 #include <IntPatch_Intersection.ixx>
21
22 #include <IntPatch_ALineToWLine.hxx>
23 #include <IntPatch_GLine.hxx>
24 #include <IntPatch_ALine.hxx>
25 #include <IntPatch_WLine.hxx>
26 #include <IntPatch_RLine.hxx>
27 #include <IntPatch_PrmPrmIntersection.hxx>
28 #include <IntPatch_ImpPrmIntersection.hxx>
29 #include <IntPatch_ImpImpIntersection.hxx>
30 #include <IntSurf_Quadric.hxx>
31
32 #include <stdio.h>
33
34 #define DEBUG 0 
35 static const Standard_Integer aNbPointsInALine = 200;
36
37 //======================================================================
38 // function: SequenceOfLine
39 //======================================================================
40 const IntPatch_SequenceOfLine& IntPatch_Intersection::SequenceOfLine() const { return(slin); }
41
42 //======================================================================
43 // function: IntPatch_Intersection
44 //======================================================================
45 IntPatch_Intersection::IntPatch_Intersection ()
46  : done(Standard_False),
47    //empt, tgte, oppo,
48    myTolArc(0.0), myTolTang(0.0),
49    myUVMaxStep(0.0), myFleche(0.0),
50    myIsStartPnt(Standard_False)
51    //myU1Start, myV1Start, myU2Start, myV2Start
52 {
53 }
54
55 //======================================================================
56 // function: IntPatch_Intersection
57 //======================================================================
58 IntPatch_Intersection::IntPatch_Intersection(const Handle(Adaptor3d_HSurface)&  S1,
59                                              const Handle(Adaptor3d_TopolTool)& D1,
60                                              const Handle(Adaptor3d_HSurface)&  S2,
61                                              const Handle(Adaptor3d_TopolTool)& D2,
62                                              const Standard_Real TolArc,
63                                              const Standard_Real TolTang)
64  : done(Standard_False),
65    //empt, tgte, oppo,
66    myTolArc(TolArc), myTolTang(TolTang),
67    myUVMaxStep(0.0), myFleche(0.0),
68    myIsStartPnt(Standard_False)
69    //myU1Start, myV1Start, myU2Start, myV2Start
70 {
71   if(myTolArc<1e-8) myTolArc=1e-8;
72   if(myTolTang<1e-8) myTolTang=1e-8;
73   if(myTolArc>0.5) myTolArc=0.5;
74   if(myTolTang>0.5) myTolTang=0.5;
75   Perform(S1,D1,S2,D2,TolArc,TolTang);
76 }
77
78 //======================================================================
79 // function: IntPatch_Intersection
80 //======================================================================
81 IntPatch_Intersection::IntPatch_Intersection(const Handle(Adaptor3d_HSurface)&  S1,
82                                              const Handle(Adaptor3d_TopolTool)& D1,
83                                              const Standard_Real TolArc,
84                                              const Standard_Real TolTang)
85  : done(Standard_False),
86    //empt, tgte, oppo,
87    myTolArc(TolArc), myTolTang(TolTang),
88    myUVMaxStep(0.0), myFleche(0.0),
89    myIsStartPnt(Standard_False)
90    //myU1Start, myV1Start, myU2Start, myV2Start
91 {
92   Perform(S1,D1,TolArc,TolTang);
93 }
94
95 //======================================================================
96 // function: SetTolerances
97 //======================================================================
98 void IntPatch_Intersection::SetTolerances(const Standard_Real TolArc,
99                                           const Standard_Real TolTang,
100                                           const Standard_Real UVMaxStep,
101                                           const Standard_Real Fleche)
102
103   myTolArc     = TolArc;
104   myTolTang    = TolTang;
105   myUVMaxStep  = UVMaxStep;
106   myFleche     = Fleche;
107   if(myTolArc<1e-8) myTolArc=1e-8;
108   if(myTolTang<1e-8) myTolTang=1e-8;
109   if(myTolArc>0.5) myTolArc=0.5;
110   if(myTolTang>0.5) myTolTang=0.5;  
111   if(myFleche<1.0e-3) myFleche=1e-3;
112   if(myUVMaxStep<1.0e-3) myUVMaxStep=1e-3;
113   if(myFleche>10) myFleche=10;
114   if(myUVMaxStep>0.5) myUVMaxStep=0.5;
115 }
116
117 //======================================================================
118 // function: Perform
119 //======================================================================
120 void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
121                                     const Handle(Adaptor3d_TopolTool)& D1,
122                                     const Standard_Real TolArc,
123                                     const Standard_Real TolTang)
124 {
125   myTolArc = TolArc;
126   myTolTang = TolTang;
127   if(myFleche == 0.0)  myFleche = 0.01;
128   if(myUVMaxStep==0.0) myUVMaxStep = 0.01;
129
130   done = Standard_True;
131   spnt.Clear();
132   slin.Clear();
133   
134   empt = Standard_True;
135   tgte = Standard_False;
136   oppo = Standard_False;
137
138   switch (S1->GetType())
139   { 
140     case GeomAbs_Plane:
141     case GeomAbs_Cylinder:
142     case GeomAbs_Sphere:
143     case GeomAbs_Cone:
144     case GeomAbs_Torus: break;
145     default:
146     {
147           IntPatch_PrmPrmIntersection interpp;
148           interpp.Perform(S1,D1,TolArc,TolTang,myFleche,myUVMaxStep);
149           if (interpp.IsDone())
150           {
151             done = Standard_True;
152             tgte = Standard_False;
153             empt = interpp.IsEmpty();
154             const Standard_Integer nblm = interpp.NbLines();
155             for (Standard_Integer i=1; i<=nblm; i++) slin.Append(interpp.Line(i));
156           }
157     }
158     break;
159   }
160 }
161
162 /////////////////////////////////////////////////////////////////////////////
163 //  These several support functions provide methods which can help basic   //
164 //  algorithm to intersect infinite surfaces of the following types:       //
165 //                                                                         //
166 //  a.) SurfaceOfExtrusion;                                                //
167 //  b.) SurfaceOfRevolution;                                               //
168 //  c.) OffsetSurface.                                                     //
169 //                                                                         //
170 /////////////////////////////////////////////////////////////////////////////
171 #include <TColgp_Array1OfXYZ.hxx>
172 #include <TColgp_Array1OfPnt.hxx>
173 #include <TColgp_SequenceOfPnt.hxx>
174 #include <Extrema_ExtPS.hxx>
175 #include <Extrema_POnSurf.hxx>
176 #include <Geom2d_Curve.hxx>
177 #include <Geom2dAPI_InterCurveCurve.hxx>
178 #include <GeomAdaptor.hxx>
179 #include <GeomAdaptor_HCurve.hxx>
180 #include <GeomAdaptor_Curve.hxx>
181 #include <GeomAdaptor_Surface.hxx>
182 #include <Handle_GeomAdaptor_HSurface.hxx>
183 #include <Geom_Plane.hxx>
184 #include <ProjLib_ProjectOnPlane.hxx>
185 #include <GeomProjLib.hxx>
186 #include <ElCLib.hxx>
187 #include <Geom_TrimmedCurve.hxx>
188 #include <Geom_Surface.hxx>
189 #include <Geom_SurfaceOfLinearExtrusion.hxx>
190 #include <Geom_OffsetSurface.hxx>
191 #include <Geom_SurfaceOfRevolution.hxx>
192 #include <Geom_RectangularTrimmedSurface.hxx>
193
194 //===============================================================
195 //function: FUN_GetMinMaxXYZPnt
196 //===============================================================
197 static void FUN_GetMinMaxXYZPnt( const Handle(Adaptor3d_HSurface)& S,
198                                  gp_Pnt& pMin, gp_Pnt& pMax )
199 {
200   const Standard_Real DU = 0.25 * Abs(S->LastUParameter() - S->FirstUParameter());
201   const Standard_Real DV = 0.25 * Abs(S->LastVParameter() - S->FirstVParameter());
202   Standard_Real tMinXYZ = RealLast();
203   Standard_Real tMaxXYZ = -tMinXYZ;
204   gp_Pnt PUV, ptMax, ptMin;
205   for(Standard_Real U = S->FirstUParameter(); U <= S->LastUParameter(); U += DU)
206   {
207     for(Standard_Real V = S->FirstVParameter(); V <= S->LastVParameter(); V += DV)
208         {
209           S->D0(U,V,PUV);
210       const Standard_Real cXYZ = PUV.XYZ().Modulus();
211           if(cXYZ > tMaxXYZ) { tMaxXYZ = cXYZ; ptMax = PUV; }
212           if(cXYZ < tMinXYZ) { tMinXYZ = cXYZ; ptMin = PUV; }
213         }
214   }
215   pMin = ptMin;
216   pMax = ptMax;
217 }
218 //==========================================================================
219 //function: FUN_TrimInfSurf
220 //==========================================================================
221 static void FUN_TrimInfSurf(const gp_Pnt& Pmin,
222                             const gp_Pnt& Pmax,
223                             const Handle(Adaptor3d_HSurface)& InfSurf,
224                             const Standard_Real& AlternativeTrimPrm,
225                             Handle(Adaptor3d_HSurface)& TrimS)
226 {
227   Standard_Real TP = AlternativeTrimPrm;
228   Extrema_ExtPS ext1(Pmin, InfSurf->Surface(), 1.e-7, 1.e-7);
229   Extrema_ExtPS ext2(Pmax, InfSurf->Surface(), 1.e-7, 1.e-7);
230   if(ext1.IsDone() || ext2.IsDone())
231   {
232     Standard_Real Umax = -1.e+100, Umin = 1.e+100, Vmax = -1.e+100, Vmin = 1.e+100, cU, cV;
233     if(ext1.IsDone())
234         {
235       for(Standard_Integer i = 1; i <= ext1.NbExt(); i++)
236       {
237             const Extrema_POnSurf & pons = ext1.Point(i);
238             pons.Parameter(cU,cV);
239             if(cU > Umax) Umax = cU;
240             if(cU < Umin) Umin = cU;
241             if(cV > Vmax) Vmax = cV;
242             if(cV < Vmin) Vmin = cV;
243       }
244         }
245     if(ext2.IsDone())
246         {
247       for(Standard_Integer i = 1; i <= ext2.NbExt(); i++)
248       {
249         const Extrema_POnSurf & pons = ext2.Point(i);
250         pons.Parameter(cU,cV);
251         if(cU > Umax) Umax = cU;
252         if(cU < Umin) Umin = cU;
253         if(cV > Vmax) Vmax = cV;
254         if(cV < Vmin) Vmin = cV;
255       }
256         }
257     TP = Max(Abs(Umin),Max(Abs(Umax),Max(Abs(Vmin),Abs(Vmax))));
258   }
259   if(TP == 0.) { TrimS = InfSurf; return; }
260   else
261   {
262     const Standard_Boolean Uinf = Precision::IsNegativeInfinite(InfSurf->FirstUParameter()); 
263     const Standard_Boolean Usup = Precision::IsPositiveInfinite(InfSurf->LastUParameter());
264     const Standard_Boolean Vinf = Precision::IsNegativeInfinite(InfSurf->FirstVParameter()); 
265     const Standard_Boolean Vsup = Precision::IsPositiveInfinite(InfSurf->LastVParameter());
266     Handle(Adaptor3d_HSurface) TmpSS;
267     Standard_Integer IsTrimed = 0;
268     const Standard_Real tp = 1000.0 * TP;
269     if(Vinf && Vsup) { TrimS = InfSurf->VTrim(-tp, tp, 1.0e-7); IsTrimed = 1; }
270     if(Vinf && !Vsup){ TrimS = InfSurf->VTrim(-tp, InfSurf->LastVParameter(), 1.0e-7); IsTrimed = 1; }
271     if(!Vinf && Vsup){ TrimS = InfSurf->VTrim(InfSurf->FirstVParameter(), tp, 1.0e-7); IsTrimed = 1; }
272     if(IsTrimed)
273         {
274           TmpSS = TrimS;
275           if(Uinf && Usup)  TrimS = TmpSS->UTrim(-tp, tp, 1.0e-7);
276           if(Uinf && !Usup) TrimS = TmpSS->UTrim(-tp, InfSurf->LastUParameter(), 1.0e-7);
277           if(!Uinf && Usup) TrimS = TmpSS->UTrim(InfSurf->FirstUParameter(), tp, 1.0e-7);
278         }
279     else
280         {
281           if(Uinf && Usup)  TrimS = InfSurf->UTrim(-tp, tp, 1.0e-7);
282           if(Uinf && !Usup) TrimS = InfSurf->UTrim(-tp, InfSurf->LastUParameter(), 1.0e-7);
283           if(!Uinf && Usup) TrimS = InfSurf->UTrim(InfSurf->FirstUParameter(), tp, 1.0e-7);
284         }
285   }
286 }
287 //================================================================================
288 //function: FUN_GetUiso
289 //================================================================================
290 static void FUN_GetUiso(const Handle(Geom_Surface)& GS,
291                         const GeomAbs_SurfaceType&  T,
292                         const Standard_Real&        FirstV,
293                         const Standard_Real&        LastV,
294                         const Standard_Boolean&     IsVC,
295                         const Standard_Boolean&     IsVP,
296                         const Standard_Real&        U,
297                         Handle(Geom_Curve)&         I)
298 {
299   if(T !=  GeomAbs_OffsetSurface)
300   {
301     Handle(Geom_Curve) gc = GS->UIso(U);
302     if(IsVP && (FirstV == 0.0 && LastV == (2.*M_PI))) I = gc;
303     else
304         {
305           Handle(Geom_TrimmedCurve) gtc = new Geom_TrimmedCurve(gc,FirstV,LastV);
306           //szv:I = Handle(Geom_Curve)::DownCast(gtc);
307           I = gtc;
308         }
309   }
310   else//OffsetSurface
311   {
312     const Handle(Geom_OffsetSurface) gos = *(Handle_Geom_OffsetSurface*)&GS;
313     const Handle(Geom_Surface) bs = gos->BasisSurface();
314     Handle(Geom_Curve) gcbs = bs->UIso(U);
315     GeomAdaptor_Curve gac(gcbs);
316     const GeomAbs_CurveType GACT = gac.GetType();
317     if(IsVP || IsVC || GACT == GeomAbs_BSplineCurve || GACT == GeomAbs_BezierCurve || Abs(LastV - FirstV) < 1.e+5)
318         {
319           Handle(Geom_Curve) gc = gos->UIso(U);
320           if(IsVP && (FirstV == 0.0 && LastV == (2*M_PI))) I = gc;
321           else
322       {
323             Handle(Geom_TrimmedCurve) gtc = new Geom_TrimmedCurve(gc,FirstV,LastV);
324             //szv:I = Handle(Geom_Curve)::DownCast(gtc);
325             I = gtc;
326           }
327         }
328     else//Offset Line, Parab, Hyperb
329         {
330           Standard_Real VmTr, VMTr;
331           if(GACT != GeomAbs_Hyperbola)
332       {
333             if(FirstV >= 0. && LastV >= 0.){ VmTr = FirstV; VMTr = ((LastV - FirstV) > 1.e+4) ? (FirstV + 1.e+4) : LastV; }
334             else if(FirstV < 0. && LastV < 0.){ VMTr = LastV; VmTr = ((FirstV - LastV) < -1.e+4) ? (LastV - 1.e+4) : FirstV; }
335             else { VmTr = (FirstV < -1.e+4) ? -1.e+4 : FirstV; VMTr = (LastV > 1.e+4) ? 1.e+4 : LastV; }
336           }
337           else//Hyperbola
338           {
339             if(FirstV >= 0. && LastV >= 0.)
340                 {
341                   if(FirstV > 4.) return;
342                   VmTr = FirstV; VMTr = (LastV > 4.) ? 4. : LastV;
343                 }
344             else if(FirstV < 0. && LastV < 0.)
345                 {
346                   if(LastV < -4.) return;
347                   VMTr = LastV; VmTr = (FirstV < -4.) ? -4. : FirstV;
348                 }
349             else { VmTr = (FirstV < -4.) ? -4. : FirstV; VMTr = (LastV > 4.) ? 4. : LastV; }
350           }
351           //Make trimmed surface
352           Handle(Geom_RectangularTrimmedSurface) rts = new Geom_RectangularTrimmedSurface(gos,VmTr,VMTr,Standard_True);
353           I = rts->UIso(U);
354         }
355   }
356 }
357 //================================================================================
358 //function: FUN_GetViso
359 //================================================================================
360 static void FUN_GetViso(const Handle(Geom_Surface)& GS,
361                         const GeomAbs_SurfaceType&  T,
362                         const Standard_Real&        FirstU,
363                         const Standard_Real&        LastU,
364                         const Standard_Boolean&     IsUC,
365                         const Standard_Boolean&     IsUP,
366                         const Standard_Real&        V,
367                         Handle(Geom_Curve)&         I)
368 {
369   if(T !=  GeomAbs_OffsetSurface)
370   {
371     Handle(Geom_Curve) gc = GS->VIso(V);
372     if(IsUP && (FirstU == 0.0 && LastU == (2*M_PI))) I = gc;
373     else
374         {
375           Handle(Geom_TrimmedCurve) gtc = new Geom_TrimmedCurve(gc,FirstU,LastU);
376           //szv:I = Handle(Geom_Curve)::DownCast(gtc);
377           I = gtc;
378         }
379   }
380   else//OffsetSurface
381   {
382     const Handle(Geom_OffsetSurface) gos = *(Handle_Geom_OffsetSurface*)&GS;
383     const Handle(Geom_Surface) bs = gos->BasisSurface();
384     Handle(Geom_Curve) gcbs = bs->VIso(V);
385     GeomAdaptor_Curve gac(gcbs);
386     const GeomAbs_CurveType GACT = gac.GetType();
387     if(IsUP || IsUC || GACT == GeomAbs_BSplineCurve || GACT == GeomAbs_BezierCurve || Abs(LastU - FirstU) < 1.e+5)
388         {
389           Handle(Geom_Curve) gc = gos->VIso(V);
390           if(IsUP && (FirstU == 0.0 && LastU == (2*M_PI))) I = gc;
391           else
392           {
393             Handle(Geom_TrimmedCurve) gtc = new Geom_TrimmedCurve(gc,FirstU,LastU);
394             //szv:I = Handle(Geom_Curve)::DownCast(gtc);
395             I = gtc;
396           }
397         }
398     else//Offset Line, Parab, Hyperb
399         {
400           Standard_Real UmTr, UMTr;
401           if(GACT != GeomAbs_Hyperbola)
402           {
403             if(FirstU >= 0. && LastU >= 0.){ UmTr = FirstU; UMTr = ((LastU - FirstU) > 1.e+4) ? (FirstU + 1.e+4) : LastU; }
404             else if(FirstU < 0. && LastU < 0.){ UMTr = LastU; UmTr = ((FirstU - LastU) < -1.e+4) ? (LastU - 1.e+4) : FirstU; }
405             else { UmTr = (FirstU < -1.e+4) ? -1.e+4 : FirstU; UMTr = (LastU > 1.e+4) ? 1.e+4 : LastU; }
406           }
407           else//Hyperbola
408           {
409             if(FirstU >= 0. && LastU >= 0.)
410                 {
411                   if(FirstU > 4.) return;
412                   UmTr = FirstU; UMTr = (LastU > 4.) ? 4. : LastU;
413                 }
414             else if(FirstU < 0. && LastU < 0.)
415                 {
416                   if(LastU < -4.) return;
417                   UMTr = LastU; UmTr = (FirstU < -4.) ? -4. : FirstU;
418                 }
419             else { UmTr = (FirstU < -4.) ? -4. : FirstU; UMTr = (LastU > 4.) ? 4. : LastU; }
420             }
421           //Make trimmed surface
422           Handle(Geom_RectangularTrimmedSurface) rts = new Geom_RectangularTrimmedSurface(gos,UmTr,UMTr,Standard_True);
423           I = rts->VIso(V);
424         }
425   }
426 }
427 //================================================================================
428 //function: FUN_PL_Intersection
429 //================================================================================
430 static void FUN_PL_Intersection(const Handle(Adaptor3d_HSurface)& S1,
431                                 const GeomAbs_SurfaceType&        T1,
432                                 const Handle(Adaptor3d_HSurface)& S2,
433                                 const GeomAbs_SurfaceType&        T2,
434                                 Standard_Boolean&                 IsOk,
435                                 TColgp_SequenceOfPnt&             SP,
436                                 gp_Vec&                           DV)
437 {
438   IsOk = Standard_False;
439   // 1. Check: both surfaces have U(V)isos - lines.
440   DV = gp_Vec(0.,0.,1.);
441   Standard_Boolean isoS1isLine[2] = {0, 0};
442   Standard_Boolean isoS2isLine[2] = {0, 0};
443   Handle(Geom_Curve) C1, C2;
444   const GeomAdaptor_Surface & gas1 = *(GeomAdaptor_Surface*)(&(S1->Surface()));
445   const GeomAdaptor_Surface & gas2 = *(GeomAdaptor_Surface*)(&(S2->Surface()));
446   const Handle(Geom_Surface) gs1 = gas1.Surface();
447   const Handle(Geom_Surface) gs2 = gas2.Surface();
448   Standard_Real MS1[2], MS2[2];
449   MS1[0] = 0.5 * (S1->LastUParameter() + S1->FirstUParameter());
450   MS1[1] = 0.5 * (S1->LastVParameter() + S1->FirstVParameter());
451   MS2[0] = 0.5 * (S2->LastUParameter() + S2->FirstUParameter());
452   MS2[1] = 0.5 * (S2->LastVParameter() + S2->FirstVParameter());
453   if(T1 == GeomAbs_SurfaceOfExtrusion) isoS1isLine[0] = Standard_True;
454   else if(!S1->IsVPeriodic() && !S1->IsVClosed()) {
455     if(T1 != GeomAbs_OffsetSurface) C1 = gs1->UIso(MS1[0]);
456     else {
457       const Handle(Geom_OffsetSurface) gos = *(Handle_Geom_OffsetSurface*)&gs1;
458       const Handle(Geom_Surface) bs = gos->BasisSurface();
459       C1 = bs->UIso(MS1[0]);
460     }
461     GeomAdaptor_Curve gac(C1);
462     if(gac.GetType() == GeomAbs_Line) isoS1isLine[0] = Standard_True;
463   }
464   if(!S1->IsUPeriodic() && !S1->IsUClosed()) {
465     if(T1 != GeomAbs_OffsetSurface) C1 = gs1->VIso(MS1[1]);
466     else {
467       const Handle(Geom_OffsetSurface) gos = *(Handle_Geom_OffsetSurface*)&gs1;
468       const Handle(Geom_Surface) bs = gos->BasisSurface();
469       C1 = bs->VIso(MS1[1]);
470     }
471     GeomAdaptor_Curve gac(C1);
472     if(gac.GetType() == GeomAbs_Line) isoS1isLine[1] = Standard_True;
473   }
474   if(T2 == GeomAbs_SurfaceOfExtrusion) isoS2isLine[0] = Standard_True;
475   else if(!S2->IsVPeriodic() && !S2->IsVClosed()) {
476     if(T2 != GeomAbs_OffsetSurface) C2 = gs2->UIso(MS2[0]);
477     else {
478       const Handle(Geom_OffsetSurface) gos = *(Handle_Geom_OffsetSurface*)&gs2;
479       const Handle(Geom_Surface) bs = gos->BasisSurface();
480       C2 = bs->UIso(MS2[0]);
481     }
482     GeomAdaptor_Curve gac(C2);
483     if(gac.GetType() == GeomAbs_Line) isoS2isLine[0] = Standard_True;
484   }
485   if(!S2->IsUPeriodic() && !S2->IsUClosed()) {
486     if(T2 != GeomAbs_OffsetSurface) C2 = gs2->VIso(MS2[1]);
487     else {
488       const Handle(Geom_OffsetSurface) gos = *(Handle_Geom_OffsetSurface*)&gs2;
489       const Handle(Geom_Surface) bs = gos->BasisSurface();
490       C2 = bs->VIso(MS2[1]);
491     }
492     GeomAdaptor_Curve gac(C2);
493     if(gac.GetType() == GeomAbs_Line) isoS2isLine[1] = Standard_True;
494   }
495   Standard_Boolean IsBothLines = ((isoS1isLine[0] || isoS1isLine[1]) &&
496                                   (isoS2isLine[0] || isoS2isLine[1]));
497   if(!IsBothLines){
498     return;
499   }
500   // 2. Check: Uiso lines of both surfaces are collinear.
501   gp_Pnt puvS1, puvS2;
502   gp_Vec derS1[2], derS2[2];
503   S1->D1(MS1[0], MS1[1], puvS1, derS1[0], derS1[1]);
504   S2->D1(MS2[0], MS2[1], puvS2, derS2[0], derS2[1]);
505   C1.Nullify(); C2.Nullify();
506   Standard_Integer iso = 0;
507   if(isoS1isLine[0] && isoS2isLine[0] &&
508      derS1[1].IsParallel(derS2[1],Precision::Angular())) {
509     iso = 1;
510     FUN_GetViso(gs1,T1,S1->FirstUParameter(),S1->LastUParameter(),
511                 S1->IsUClosed(),S1->IsUPeriodic(),MS1[1],C1);
512     FUN_GetViso(gs2,T2,S2->FirstUParameter(),S2->LastUParameter(),
513                 S2->IsUClosed(),S2->IsUPeriodic(),MS2[1],C2);
514   }
515   else if(isoS1isLine[0] && isoS2isLine[1] &&
516           derS1[1].IsParallel(derS2[0],Precision::Angular())) {
517     iso = 1;
518     FUN_GetViso(gs1,T1,S1->FirstUParameter(),S1->LastUParameter(),
519                 S1->IsUClosed(),S1->IsUPeriodic(),MS1[1],C1);
520     FUN_GetUiso(gs2,T2,S2->FirstVParameter(),S2->LastVParameter(),
521                 S2->IsVClosed(),S2->IsVPeriodic(),MS2[0],C2);
522   }
523   else if(isoS1isLine[1] && isoS2isLine[0] &&
524           derS1[0].IsParallel(derS2[1],Precision::Angular())) {
525     iso = 0;
526     FUN_GetUiso(gs1,T1,S1->FirstVParameter(),S1->LastVParameter(),
527                 S1->IsVClosed(),S1->IsVPeriodic(),MS1[0],C1);
528     FUN_GetViso(gs2,T2,S2->FirstUParameter(),S2->LastUParameter(),
529                 S2->IsUClosed(),S2->IsUPeriodic(),MS2[1],C2);
530   }
531   else if(isoS1isLine[1] && isoS2isLine[1] &&
532           derS1[0].IsParallel(derS2[0],Precision::Angular())) {
533     iso = 0;
534     FUN_GetUiso(gs1,T1,S1->FirstVParameter(),S1->LastVParameter(),
535                 S1->IsVClosed(),S1->IsVPeriodic(),MS1[0],C1);
536     FUN_GetUiso(gs2,T2,S2->FirstVParameter(),S2->LastVParameter(),
537                 S2->IsVClosed(),S2->IsVPeriodic(),MS2[0],C2);
538   }
539   else {
540     IsOk = Standard_False;
541     return;
542   }
543   IsOk = Standard_True;
544   // 3. Make intersections of V(U)isos
545   if(C1.IsNull() || C2.IsNull()) return;
546   DV = derS1[iso];
547   Handle(Geom_Plane) GPln = new Geom_Plane(gp_Pln(puvS1,gp_Dir(DV)));
548   Handle(Geom_Curve) C1Prj =
549     GeomProjLib::ProjectOnPlane(C1,GPln,gp_Dir(DV),Standard_True);
550   Handle(Geom_Curve) C2Prj =
551     GeomProjLib::ProjectOnPlane(C2,GPln,gp_Dir(DV),Standard_True);
552   if(C1Prj.IsNull() || C2Prj.IsNull()) return;
553   Handle(Geom2d_Curve) C1Prj2d =
554     GeomProjLib::Curve2d(C1Prj,*(Handle_Geom_Surface *)&GPln);
555   Handle(Geom2d_Curve) C2Prj2d =
556     GeomProjLib::Curve2d(C2Prj,*(Handle_Geom_Surface *)&GPln);
557   Geom2dAPI_InterCurveCurve ICC(C1Prj2d,C2Prj2d,1.0e-7);
558   if(ICC.NbPoints() > 0 )
559   {
560     for(Standard_Integer ip = 1; ip <= ICC.NbPoints(); ip++)
561         {
562           gp_Pnt2d P = ICC.Point(ip);
563           gp_Pnt P3d = ElCLib::To3d(gp_Ax2(puvS1,gp_Dir(DV)),P);
564           SP.Append(P3d);
565         }
566   }
567 }
568 //================================================================================
569 //function: FUN_NewFirstLast
570 //================================================================================
571 static void FUN_NewFirstLast(const GeomAbs_CurveType& ga_ct,
572                              const Standard_Real&     Fst,
573                              const Standard_Real&     Lst,
574                              const Standard_Real&     TrVal,
575                              Standard_Real&           NewFst,
576                              Standard_Real&           NewLst,
577                              Standard_Boolean&        NeedTr)
578 {
579   NewFst = Fst; NewLst = Lst; NeedTr = Standard_False;
580   switch (ga_ct)
581   {
582     case GeomAbs_Line:
583     case GeomAbs_Parabola:
584     {
585       if(Abs(Lst - Fst) > TrVal)
586           {
587             if(Fst >= 0. && Lst >= 0.)
588         {
589               NewFst = Fst;
590               NewLst = ((Fst + TrVal) < Lst) ? (Fst + TrVal) : Lst;
591             }
592             if(Fst < 0. && Lst < 0.)
593             {
594               NewLst = Lst;
595               NewFst = ((Lst - TrVal) > Fst) ? (Lst - TrVal) : Fst;
596             }
597             else
598             {
599               NewFst = (Fst < -TrVal) ? -TrVal : Fst;
600               NewLst = (Lst > TrVal) ? TrVal : Lst;
601             }
602             NeedTr = Standard_True;
603           }
604           break;
605     }
606         case GeomAbs_Hyperbola:
607     {
608       if(Abs(Lst - Fst) > 10.)
609           { 
610             if(Fst >= 0. && Lst >= 0.)
611             {
612               if(Fst > 4.) return;
613               NewFst = Fst;
614               NewLst = (Lst > 4.) ? 4. : Lst;
615             }
616             if(Fst < 0. && Lst < 0.)
617             {
618               if(Lst < -4.) return;
619               NewLst = Lst;
620               NewFst = (Fst < -4.) ? -4. : Fst;
621             }
622             else
623             {
624               NewFst = (Fst < -4.) ? -4. : Fst;
625               NewLst = (Lst > 4.) ? 4. : Lst;
626             }
627             NeedTr = Standard_True;
628           }
629       break;
630     }
631   }
632 }
633 //================================================================================
634 //function: FUN_TrimBothSurf
635 //================================================================================               
636 static void FUN_TrimBothSurf(const Handle(Adaptor3d_HSurface)& S1,
637                              const GeomAbs_SurfaceType&        T1,
638                              const Handle(Adaptor3d_HSurface)& S2,
639                              const GeomAbs_SurfaceType&        T2,
640                              const Standard_Real&              TV,
641                              Handle(Adaptor3d_HSurface)&       NS1,
642                              Handle(Adaptor3d_HSurface)&       NS2)
643 {
644   const GeomAdaptor_Surface & gas1 = *(GeomAdaptor_Surface*)(&(S1->Surface()));
645   const GeomAdaptor_Surface & gas2 = *(GeomAdaptor_Surface*)(&(S2->Surface()));
646   const Handle(Geom_Surface) gs1 = gas1.Surface();
647   const Handle(Geom_Surface) gs2 = gas2.Surface();
648   const Standard_Real UM1 = 0.5 * (S1->LastUParameter() + S1->FirstUParameter());
649   const Standard_Real UM2 = 0.5 * (S2->LastUParameter() + S2->FirstUParameter());
650   const Standard_Real VM1 = 0.5 * (S1->LastVParameter() + S1->FirstVParameter());
651   const Standard_Real VM2 = 0.5 * (S2->LastVParameter() + S2->FirstVParameter());
652   Handle(Geom_Curve) visoS1, visoS2, uisoS1, uisoS2;
653   if(T1 != GeomAbs_OffsetSurface){ visoS1 = gs1->VIso(VM1); uisoS1 = gs1->UIso(UM1); }
654   else
655   {
656     const Handle(Geom_OffsetSurface) gos = *(Handle_Geom_OffsetSurface*)&gs1;
657     const Handle(Geom_Surface) bs = gos->BasisSurface();
658     visoS1 = bs->VIso(VM1); uisoS1 = bs->UIso(UM1);
659   }
660   if(T2 != GeomAbs_OffsetSurface){ visoS2 = gs2->VIso(VM2); uisoS2 = gs2->UIso(UM2); }
661   else
662   {
663     const Handle(Geom_OffsetSurface) gos = *(Handle_Geom_OffsetSurface*)&gs2;
664     const Handle(Geom_Surface) bs = gos->BasisSurface();
665     visoS2 = bs->VIso(VM2); uisoS2 = bs->UIso(UM2);
666   }
667   if(uisoS1.IsNull() || uisoS2.IsNull() || visoS1.IsNull() || visoS2.IsNull()){ NS1 = S1; NS2 = S2; return; }
668   GeomAdaptor_Curve gau1(uisoS1);
669   GeomAdaptor_Curve gav1(visoS1);
670   GeomAdaptor_Curve gau2(uisoS2);
671   GeomAdaptor_Curve gav2(visoS2);
672   GeomAbs_CurveType GA_U1 = gau1.GetType();
673   GeomAbs_CurveType GA_V1 = gav1.GetType();
674   GeomAbs_CurveType GA_U2 = gau2.GetType();
675   GeomAbs_CurveType GA_V2 = gav2.GetType();
676   Standard_Boolean TrmU1 = Standard_False;
677   Standard_Boolean TrmV1 = Standard_False;
678   Standard_Boolean TrmU2 = Standard_False;
679   Standard_Boolean TrmV2 = Standard_False;
680   Standard_Real V1S1,V2S1,U1S1,U2S1, V1S2,V2S2,U1S2,U2S2;
681   FUN_NewFirstLast(GA_U1,S1->FirstVParameter(),S1->LastVParameter(),TV,V1S1,V2S1,TrmV1);
682   FUN_NewFirstLast(GA_V1,S1->FirstUParameter(),S1->LastUParameter(),TV,U1S1,U2S1,TrmU1);
683   FUN_NewFirstLast(GA_U2,S2->FirstVParameter(),S2->LastVParameter(),TV,V1S2,V2S2,TrmV2);
684   FUN_NewFirstLast(GA_V2,S2->FirstUParameter(),S2->LastUParameter(),TV,U1S2,U2S2,TrmU2);
685   if(TrmV1) NS1 = S1->VTrim(V1S1, V2S1, 1.0e-7);
686   if(TrmV2) NS2 = S2->VTrim(V1S2, V2S2, 1.0e-7);
687   if(TrmU1)
688   {
689     if(TrmV1)
690         {
691           Handle(Adaptor3d_HSurface) TS = NS1;
692           NS1 = TS->UTrim(U1S1, U2S1, 1.0e-7);
693         }
694     else NS1 = S1->UTrim(U1S1, U2S1, 1.0e-7);
695   }
696   if(TrmU2)
697   {
698     if(TrmV2)
699         {
700           Handle(Adaptor3d_HSurface) TS = NS2;
701           NS2 = TS->UTrim(U1S2, U2S2, 1.0e-7);
702         }
703     else NS2 = S2->UTrim(U1S2, U2S2, 1.0e-7);
704   }
705 }
706
707 //=======================================================================
708 //function : Perform
709 //purpose  : 
710 //=======================================================================
711 void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  theS1,
712                                     const Handle(Adaptor3d_TopolTool)& theD1,
713                                     const Handle(Adaptor3d_HSurface)&  theS2,
714                                     const Handle(Adaptor3d_TopolTool)& theD2,
715                                     const Standard_Real TolArc,
716                                     const Standard_Real TolTang,
717                                     const Standard_Boolean isGeomInt)
718 {
719   myTolArc = TolArc;
720   myTolTang = TolTang;
721   if(myFleche <= Precision::PConfusion())
722     myFleche = 0.01;
723   if(myUVMaxStep <= Precision::PConfusion())
724     myUVMaxStep = 0.01;
725
726   done = Standard_False;
727   spnt.Clear();
728   slin.Clear();
729   empt = Standard_True;
730   tgte = Standard_False;
731   oppo = Standard_False;
732
733   GeomAbs_SurfaceType typs1 = theS1->GetType();
734   GeomAbs_SurfaceType typs2 = theS2->GetType();
735   
736   Standard_Boolean TreatAsBiParametric = Standard_False;
737   if(typs1 == GeomAbs_Cone)
738   {
739     const gp_Cone Con1 = theS1->Cone();
740     const Standard_Real a1 = Abs(Con1.SemiAngle());
741     if((a1 < 0.02) || (a1 > 1.55))
742     {
743       if(typs2==GeomAbs_Plane)
744       {
745         if(a1 < 0.02)
746         {
747           const gp_Pln Plan2 = theS2->Plane();
748           const gp_Dir axec = Con1.Axis().Direction();
749           const gp_Dir axep = Plan2.Axis().Direction();
750           const Standard_Real ps = Abs(axec.Dot(axep));
751           if(ps < 0.015)
752           {
753             TreatAsBiParametric = Standard_True;
754           }
755         }
756       }
757       else
758         TreatAsBiParametric = Standard_True;
759     }
760   }
761
762   if(typs2 == GeomAbs_Cone)
763   {
764     const gp_Cone Con2 = theS2->Cone();
765     const Standard_Real a2 = Abs(Con2.SemiAngle());
766     if((a2 < 0.02) || (a2 > 1.55))
767     {
768       if(typs1 == GeomAbs_Plane)
769       {
770         if(a2 < 0.02)
771         {
772           const gp_Pln Plan1 = theS1->Plane();
773           const gp_Dir axec = Con2.Axis().Direction();
774           const gp_Dir axep = Plan1.Axis().Direction();
775           const Standard_Real ps = Abs(axec.Dot(axep));
776           if(ps<0.015)
777           {
778             TreatAsBiParametric = Standard_True;
779           }
780         }
781       }
782       else
783         TreatAsBiParametric = Standard_True;
784     }
785
786     //// modified by jgv, 15.12.02 for OCC565 ////
787     if (typs1 == GeomAbs_Cone)
788     {
789       const gp_Cone Con1 = theS1->Cone();
790       const Standard_Real a1 = Abs(Con1.SemiAngle());
791       if (a1 < 0.02 && a2 < 0.02) //quasi-cylinders: if same domain, treat as canonic
792       {
793         const gp_Ax1 A1 = Con1.Axis(), A2 = Con2.Axis();
794         if (A1.IsParallel(A2,Precision::Angular()))
795         {
796           const gp_Lin L1(A1);
797           if (L1.Distance(A2.Location()) <= Precision::Confusion())
798           {
799             TreatAsBiParametric = Standard_False;
800           }
801         }
802       }
803       else if (a1 > 1.55 && a2 > 1.55) //quasi-planes: if same domain, treat as canonic
804       {
805         const gp_Ax1 A1 = Con1.Axis(), A2 = Con2.Axis();
806         if (A1.IsParallel(A2,Precision::Angular()))
807         {
808           const gp_Pnt Apex1 = Con1.Apex(), Apex2 = Con2.Apex();
809           const gp_Pln Plan1( Apex1, A1.Direction() );
810           if (Plan1.Distance( Apex2 ) <= Precision::Confusion())
811           {
812             TreatAsBiParametric = Standard_False;
813           }
814         }
815       }
816       else if ((a1 > 1.55) || (a2 > 1.55))
817       {
818         const gp_Ax1 A1 = Con1.Axis(), A2 = Con2.Axis();
819         if (A1.IsCoaxial(A2,Precision::Angular(),Precision::Confusion()))
820         {
821           TreatAsBiParametric = Standard_False;
822         }
823       }
824     }// if (typs1 == GeomAbs_Cone)    {
825   }// if(typs2 == GeomAbs_Cone)  {
826
827   if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite()) {
828     TreatAsBiParametric= Standard_False;
829   }
830
831 //  Modified by skv - Mon Sep 26 14:58:30 2005 Begin
832 //   if(TreatAsBiParametric) { typs1 = typs2 = GeomAbs_BezierSurface; }
833   if(TreatAsBiParametric)
834   {
835     if (typs1 == GeomAbs_Cone && typs2 == GeomAbs_Plane)
836       typs1 = GeomAbs_BezierSurface; // Using Imp-Prm Intersector
837     else if (typs1 == GeomAbs_Plane && typs2 == GeomAbs_Cone)
838       typs2 = GeomAbs_BezierSurface; // Using Imp-Prm Intersector
839     else {
840       // Using Prm-Prm Intersector
841       typs1 = GeomAbs_BezierSurface;
842       typs2 = GeomAbs_BezierSurface;
843     }
844   }
845 //  Modified by skv - Mon Sep 26 14:58:30 2005 End
846
847   // Surface type definition
848   Standard_Integer ts1 = 0;
849   switch (typs1)
850   {
851     case GeomAbs_Plane:
852     case GeomAbs_Cylinder:
853     case GeomAbs_Sphere:
854     case GeomAbs_Cone: ts1 = 1; break;
855     default: break;
856   }
857
858   Standard_Integer ts2 = 0;
859   switch (typs2)
860   {
861     case GeomAbs_Plane:
862     case GeomAbs_Cylinder:
863     case GeomAbs_Sphere:
864     case GeomAbs_Cone: ts2 = 1; break;
865     default: break;
866   }
867
868   // Possible intersection types: 1. ts1 == ts2 == 1 <Geom-Geom>
869   //                              2. ts1 != ts2      <Geom-Param>
870   //                              3. ts1 == ts2 == 0 <Param-Param>
871
872   // Geom - Geom
873   if(ts1 == ts2 && ts1 == 1)
874   {
875     const Standard_Boolean RestrictLine = Standard_True;
876     IntSurf_ListOfPntOn2S ListOfPnts;
877     ListOfPnts.Clear();
878     if(isGeomInt)
879     {
880       GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
881     }
882     else
883     {
884       ParamParamPerfom(theS1, theD1, theS2, theD2, TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
885     }
886   }
887
888   // Geom - Param
889   if(ts1 != ts2)
890   {
891     GeomParamPerfom(theS1, theD1, theS2, theD2, ts1 == 0, typs1, typs2);
892   }
893
894   // Param - Param 
895   if(ts1 == ts2 && ts1 == 0)
896   {
897     const Standard_Boolean RestrictLine = Standard_True;
898     IntSurf_ListOfPntOn2S ListOfPnts;
899     ListOfPnts.Clear();
900
901     ParamParamPerfom(theS1, theD1, theS2, theD2, TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
902   }
903 }
904                       
905 //=======================================================================
906 //function : Perform
907 //purpose  : 
908 //=======================================================================
909 void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  theS1,
910                                     const Handle(Adaptor3d_TopolTool)& theD1,
911                                     const Handle(Adaptor3d_HSurface)&  theS2,
912                                     const Handle(Adaptor3d_TopolTool)& theD2,
913                                     const Standard_Real TolArc,
914                                     const Standard_Real TolTang,
915                                     IntSurf_ListOfPntOn2S& ListOfPnts,
916                                     const Standard_Boolean RestrictLine,
917                                     const Standard_Boolean isGeomInt)
918 {
919   myTolArc = TolArc;
920   myTolTang = TolTang;
921   if(myFleche <= Precision::PConfusion())
922     myFleche = 0.01;
923   if(myUVMaxStep <= Precision::PConfusion())
924     myUVMaxStep = 0.01;
925     
926   done = Standard_False;
927   spnt.Clear();
928   slin.Clear();
929   empt = Standard_True;
930   tgte = Standard_False;
931   oppo = Standard_False;
932
933   GeomAbs_SurfaceType typs1 = theS1->GetType();
934   GeomAbs_SurfaceType typs2 = theS2->GetType();
935   
936   Standard_Boolean TreatAsBiParametric = Standard_False;
937   if(typs1 == GeomAbs_Cone)
938   {
939     const gp_Cone Con1 = theS1->Cone();
940     const Standard_Real a1 = Abs(Con1.SemiAngle());
941     if((a1 < 0.02) || (a1 > 1.55))
942     {
943       if(typs2==GeomAbs_Plane)
944       {
945         if(a1 < 0.02)
946         {
947           const gp_Pln Plan2 = theS2->Plane();
948           const gp_Dir axec = Con1.Axis().Direction();
949           const gp_Dir axep = Plan2.Axis().Direction();
950           const Standard_Real ps = Abs(axec.Dot(axep));
951           if(ps < 0.015)
952           {
953             TreatAsBiParametric = Standard_True;
954           }
955         }
956       }
957       else
958         TreatAsBiParametric = Standard_True;
959     }
960   }
961
962   if(typs2 == GeomAbs_Cone)
963   {
964     const gp_Cone Con2 = theS2->Cone();
965     const Standard_Real a2 = Abs(Con2.SemiAngle());
966     if((a2 < 0.02) || (a2 > 1.55))
967     {
968       if(typs1 == GeomAbs_Plane)
969       {
970         if(a2 < 0.02)
971         {
972           const gp_Pln Plan1 = theS1->Plane();
973           const gp_Dir axec = Con2.Axis().Direction();
974           const gp_Dir axep = Plan1.Axis().Direction();
975           const Standard_Real ps = Abs(axec.Dot(axep));
976           if(ps<0.015)
977           {
978             TreatAsBiParametric = Standard_True;
979           }
980         }
981       }
982       else
983         TreatAsBiParametric = Standard_True;
984     }
985
986     //// modified by jgv, 15.12.02 for OCC565 ////
987     if (typs1 == GeomAbs_Cone)
988     {
989       const gp_Cone Con1 = theS1->Cone();
990       const Standard_Real a1 = Abs(Con1.SemiAngle());
991       if (a1 < 0.02 && a2 < 0.02) //quasi-cylinders: if same domain, treat as canonic
992       {
993         const gp_Ax1 A1 = Con1.Axis(), A2 = Con2.Axis();
994         if (A1.IsParallel(A2,Precision::Angular()))
995         {
996           const gp_Lin L1(A1);
997           if (L1.Distance(A2.Location()) <= Precision::Confusion())
998           {
999             TreatAsBiParametric = Standard_False;
1000           }
1001         }
1002       }
1003       else if (a1 > 1.55 && a2 > 1.55) //quasi-planes: if same domain, treat as canonic
1004       {
1005         const gp_Ax1 A1 = Con1.Axis(), A2 = Con2.Axis();
1006         if (A1.IsParallel(A2,Precision::Angular()))
1007         {
1008           const gp_Pnt Apex1 = Con1.Apex(), Apex2 = Con2.Apex();
1009           const gp_Pln Plan1( Apex1, A1.Direction() );
1010           if (Plan1.Distance( Apex2 ) <= Precision::Confusion())
1011           {
1012             TreatAsBiParametric = Standard_False;
1013           }
1014         }
1015       }
1016       else if ((a1 > 1.55) || (a2 > 1.55))
1017       {
1018         const gp_Ax1 A1 = Con1.Axis(), A2 = Con2.Axis();
1019         if (A1.IsCoaxial(A2,Precision::Angular(),Precision::Confusion()))
1020         {
1021           TreatAsBiParametric = Standard_False;
1022         }
1023       }
1024     }// if (typs1 == GeomAbs_Cone)    {
1025   }// if(typs2 == GeomAbs_Cone)  {
1026
1027   if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite()) {
1028     TreatAsBiParametric= Standard_False;
1029   }
1030
1031   if(TreatAsBiParametric)
1032   {
1033     // Using Prm-Prm Intersector
1034     typs1 = GeomAbs_BezierSurface;
1035     typs2 = GeomAbs_BezierSurface;
1036   }
1037
1038   // Surface type definition
1039   Standard_Integer ts1 = 0;
1040   switch (typs1)
1041   {
1042     case GeomAbs_Plane:
1043     case GeomAbs_Cylinder:
1044     case GeomAbs_Sphere:
1045     case GeomAbs_Cone: ts1 = 1; break;
1046     default: break;
1047   }
1048
1049   Standard_Integer ts2 = 0;
1050   switch (typs2)
1051   {
1052     case GeomAbs_Plane:
1053     case GeomAbs_Cylinder:
1054     case GeomAbs_Sphere:
1055     case GeomAbs_Cone: ts2 = 1; break;
1056     default: break;
1057   }
1058
1059   // Possible intersection types: 1. ts1 == ts2 == 1 <Geom-Geom>
1060   //                              2. ts1 != ts2      <Geom-Param>
1061   //                              3. ts1 == ts2 == 0 <Param-Param>
1062
1063   if(!isGeomInt)
1064   {
1065     ParamParamPerfom(theS1, theD1, theS2, theD2, TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1066   }
1067   else if(ts1 != ts2)
1068   {
1069     GeomParamPerfom(theS1, theD1, theS2, theD2, ts1 == 0, typs1, typs2);
1070   }
1071   else if (ts1 == 0)
1072   {
1073     ParamParamPerfom(theS1, theD1, theS2, theD2, TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1074   }
1075   else if(ts1 == 1)
1076   {
1077     GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1078   }
1079 }
1080
1081 //=======================================================================
1082 //function : ParamParamPerfom
1083 //purpose  : 
1084 //=======================================================================
1085 void IntPatch_Intersection::ParamParamPerfom(const Handle(Adaptor3d_HSurface)&  theS1,
1086                                              const Handle(Adaptor3d_TopolTool)& theD1,
1087                                              const Handle(Adaptor3d_HSurface)&  theS2,
1088                                              const Handle(Adaptor3d_TopolTool)& theD2,
1089                                              const Standard_Real TolArc,
1090                                              const Standard_Real TolTang,
1091                                              IntSurf_ListOfPntOn2S& ListOfPnts,
1092                                              const Standard_Boolean RestrictLine,
1093                                              const GeomAbs_SurfaceType typs1,
1094                                              const GeomAbs_SurfaceType typs2)
1095 {
1096   IntPatch_PrmPrmIntersection interpp;
1097   if(!theD1->DomainIsInfinite() && !theD2->DomainIsInfinite())
1098   {
1099     Standard_Boolean ClearFlag = Standard_True;
1100     if(!ListOfPnts.IsEmpty())
1101     {
1102       interpp.Perform(theS1,theD1,theS2,theD2,TolArc,TolTang,myFleche,myUVMaxStep, ListOfPnts, RestrictLine);
1103       ClearFlag = Standard_False;
1104     }
1105     interpp.Perform(theS1,theD1,theS2,theD2,TolArc,TolTang,myFleche,myUVMaxStep,ClearFlag);   //double call!!!!!!!
1106   }
1107   else if((theD1->DomainIsInfinite()) ^ (theD2->DomainIsInfinite()))
1108   {
1109     gp_Pnt pMaxXYZ, pMinXYZ;
1110     if(theD1->DomainIsInfinite())
1111     {
1112       FUN_GetMinMaxXYZPnt( theS2, pMinXYZ, pMaxXYZ );
1113       const Standard_Real MU = Max(Abs(theS2->FirstUParameter()),Abs(theS2->LastUParameter()));
1114       const Standard_Real MV = Max(Abs(theS2->FirstVParameter()),Abs(theS2->LastVParameter()));
1115       const Standard_Real AP = Max(MU, MV);
1116       Handle(Adaptor3d_HSurface) SS;
1117       FUN_TrimInfSurf(pMinXYZ, pMaxXYZ, theS1, AP, SS);
1118       interpp.Perform(SS,theD1,theS2,theD2,TolArc,TolTang,myFleche,myUVMaxStep);
1119     }
1120     else
1121     {
1122       FUN_GetMinMaxXYZPnt( theS1, pMinXYZ, pMaxXYZ );
1123       const Standard_Real MU = Max(Abs(theS1->FirstUParameter()),Abs(theS1->LastUParameter()));
1124       const Standard_Real MV = Max(Abs(theS1->FirstVParameter()),Abs(theS1->LastVParameter()));
1125       const Standard_Real AP = Max(MU, MV);
1126       Handle(Adaptor3d_HSurface) SS;
1127       FUN_TrimInfSurf(pMinXYZ, pMaxXYZ, theS2, AP, SS);
1128       interpp.Perform(theS1, theD1, SS, theD2,TolArc,TolTang,myFleche,myUVMaxStep);
1129     }
1130   }//(theD1->DomainIsInfinite()) ^ (theD2->DomainIsInfinite())
1131   else
1132   {
1133     if(typs1 == GeomAbs_OtherSurface || typs2 == GeomAbs_OtherSurface)
1134     {
1135       done = Standard_False;
1136       return;
1137     }
1138
1139     Standard_Boolean IsPLInt = Standard_False;
1140     TColgp_SequenceOfPnt sop;
1141     gp_Vec v;
1142     FUN_PL_Intersection(theS1,typs1,theS2,typs2,IsPLInt,sop,v);
1143
1144     if(IsPLInt)
1145     {
1146       if(sop.Length() > 0)
1147       {
1148         for(Standard_Integer ip = 1; ip <= sop.Length(); ip++)
1149         {
1150           gp_Lin lin(sop.Value(ip),gp_Dir(v));
1151           Handle(IntPatch_GLine) gl = new IntPatch_GLine(lin,Standard_False);
1152           slin.Append(*(Handle_IntPatch_Line *)&gl);
1153         }
1154
1155         done = Standard_True;
1156       }
1157       else
1158         done = Standard_False;
1159
1160       return;
1161     }// 'COLLINEAR LINES'
1162     else
1163     {
1164       Handle(Adaptor3d_HSurface) nS1 = theS1;
1165       Handle(Adaptor3d_HSurface) nS2 = theS2;
1166       FUN_TrimBothSurf(theS1,typs1,theS2,typs2,1.e+8,nS1,nS2);
1167       interpp.Perform(nS1,theD1,nS2,theD2,TolArc,TolTang,myFleche,myUVMaxStep);
1168     }// 'NON - COLLINEAR LINES'
1169   }// both domains are infinite
1170
1171   if (interpp.IsDone())
1172   {
1173     done = Standard_True;
1174     tgte = Standard_False;
1175     empt = interpp.IsEmpty();
1176
1177     for(Standard_Integer i = 1; i <= interpp.NbLines(); i++)
1178     {
1179       if(interpp.Line(i)->ArcType() != IntPatch_Walking)
1180         slin.Append(interpp.Line(i));
1181     }
1182
1183     for (Standard_Integer i = 1; i <= interpp.NbLines(); i++)
1184     {
1185       if(interpp.Line(i)->ArcType() == IntPatch_Walking)
1186         slin.Append(interpp.Line(i));
1187     }
1188   }
1189 }
1190
1191 //=======================================================================
1192 ////function : GeomGeomPerfom
1193 //purpose  : 
1194 //=======================================================================
1195 void IntPatch_Intersection::GeomGeomPerfom(const Handle(Adaptor3d_HSurface)& theS1,
1196                                            const Handle(Adaptor3d_TopolTool)& theD1,
1197                                            const Handle(Adaptor3d_HSurface)& theS2,
1198                                            const Handle(Adaptor3d_TopolTool)& theD2,
1199                                            const Standard_Real TolArc,
1200                                            const Standard_Real TolTang,
1201                                            IntSurf_ListOfPntOn2S& ListOfPnts,
1202                                            const Standard_Boolean RestrictLine,
1203                                            const GeomAbs_SurfaceType typs1,
1204                                            const GeomAbs_SurfaceType typs2)
1205 {
1206   IntPatch_ImpImpIntersection interii(theS1,theD1,theS2,theD2,myTolArc,myTolTang);
1207   const Standard_Boolean anIS = interii.IsDone();
1208   if (anIS)
1209   {
1210     done = anIS;
1211     empt = interii.IsEmpty();
1212     if (!empt)
1213     {
1214       tgte = interii.TangentFaces();
1215       if (tgte)
1216         oppo = interii.OppositeFaces();
1217
1218       for (Standard_Integer i = 1; i <= interii.NbLines(); i++)
1219       {
1220         const Handle_IntPatch_Line& line = interii.Line(i);
1221         if (line->ArcType() == IntPatch_Analytic)
1222         {
1223           const GeomAbs_SurfaceType typs1 = theS1->GetType();
1224           const GeomAbs_SurfaceType typs2 = theS2->GetType();
1225           IntSurf_Quadric Quad1,Quad2;
1226           
1227           switch(typs1)
1228           {
1229           case GeomAbs_Plane:
1230             Quad1.SetValue(theS1->Plane());
1231             break;
1232
1233           case GeomAbs_Cylinder:
1234             Quad1.SetValue(theS1->Cylinder());
1235             break;
1236
1237           case GeomAbs_Sphere:
1238             Quad1.SetValue(theS1->Sphere());
1239             break;
1240
1241           case GeomAbs_Cone:
1242             Quad1.SetValue(theS1->Cone());
1243             break;
1244
1245           default:
1246             break;
1247           }
1248
1249           switch(typs2)
1250           {
1251           case GeomAbs_Plane:
1252             Quad2.SetValue(theS2->Plane());
1253             break;
1254           case GeomAbs_Cylinder:
1255             Quad2.SetValue(theS2->Cylinder());
1256             break;
1257
1258           case GeomAbs_Sphere:
1259             Quad2.SetValue(theS2->Sphere());
1260             break;
1261
1262           case GeomAbs_Cone:
1263             Quad2.SetValue(theS2->Cone());
1264             break;
1265
1266           default:
1267             break;
1268           }
1269
1270           IntPatch_ALineToWLine AToW(Quad1,Quad2,0.01,0.05,aNbPointsInALine);
1271           Handle(IntPatch_Line) wlin=AToW.MakeWLine((*((Handle_IntPatch_ALine *)(&line))));
1272           slin.Append(wlin);
1273         }
1274         else
1275           slin.Append(interii.Line(i));
1276       }
1277
1278       for (Standard_Integer i = 1; i <= interii.NbPnts(); i++)
1279       {
1280         spnt.Append(interii.Point(i));
1281       }
1282     }
1283   }
1284   else
1285     ParamParamPerfom(theS1, theD1, theS2, theD2, TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2);
1286 }
1287
1288 //=======================================================================
1289 ////function : GeomParamPerfom
1290 //purpose  : 
1291 //=======================================================================
1292 void IntPatch_Intersection::GeomParamPerfom(const Handle(Adaptor3d_HSurface)&  theS1,
1293                                             const Handle(Adaptor3d_TopolTool)& theD1,
1294                                             const Handle(Adaptor3d_HSurface)&  theS2,
1295                                             const Handle(Adaptor3d_TopolTool)& theD2,
1296                                             const Standard_Boolean isNotAnalitical,
1297                                             const GeomAbs_SurfaceType typs1,
1298                                             const GeomAbs_SurfaceType typs2)
1299 {
1300   IntPatch_ImpPrmIntersection interip;
1301   if (myIsStartPnt)
1302   {
1303     if (isNotAnalitical/*ts1 == 0*/)
1304       interip.SetStartPoint(myU1Start,myV1Start);
1305     else
1306       interip.SetStartPoint(myU2Start,myV2Start);
1307   }
1308
1309   if(theD1->DomainIsInfinite() && theD2->DomainIsInfinite())
1310   {
1311     Standard_Boolean IsPLInt = Standard_False;
1312     TColgp_SequenceOfPnt sop;
1313     gp_Vec v;
1314     FUN_PL_Intersection(theS1,typs1,theS2,typs2,IsPLInt,sop,v);
1315     
1316     if(IsPLInt)
1317     {
1318       if(sop.Length() > 0)
1319       {
1320         for(Standard_Integer ip = 1; ip <= sop.Length(); ip++)
1321         {
1322           gp_Lin lin(sop.Value(ip),gp_Dir(v));
1323           Handle(IntPatch_GLine) gl = new IntPatch_GLine(lin,Standard_False);
1324           slin.Append(*(Handle_IntPatch_Line *)&gl);
1325         }
1326
1327         done = Standard_True;
1328       }
1329       else
1330         done = Standard_False;
1331
1332       return;
1333     }
1334     else
1335     {
1336       Handle(Adaptor3d_HSurface) nS1 = theS1;
1337       Handle(Adaptor3d_HSurface) nS2 = theS2;
1338       FUN_TrimBothSurf(theS1,typs1,theS2,typs2,1.e+5,nS1,nS2);
1339       interip.Perform(nS1,theD1,nS2,theD2,myTolArc,myTolTang,myFleche,myUVMaxStep);
1340     }
1341   }
1342   else
1343     interip.Perform(theS1,theD1,theS2,theD2,myTolArc,myTolTang,myFleche,myUVMaxStep);
1344
1345   if (interip.IsDone()) 
1346   {
1347     done = Standard_True;
1348     empt = interip.IsEmpty();
1349
1350     if (!empt)
1351     {
1352       for(Standard_Integer i = 1; i <= interip.NbLines(); i++)
1353       {
1354         if(interip.Line(i)->ArcType() != IntPatch_Walking)
1355           slin.Append(interip.Line(i));
1356       }
1357
1358       for(Standard_Integer i = 1; i <= interip.NbLines(); i++)
1359       {
1360         if(interip.Line(i)->ArcType() == IntPatch_Walking)
1361           slin.Append(interip.Line(i));
1362       }
1363
1364       for (Standard_Integer i = 1; i <= interip.NbPnts(); i++)
1365         spnt.Append(interip.Point(i));
1366     }
1367   }
1368 }
1369
1370
1371 void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
1372                                     const Handle(Adaptor3d_TopolTool)& D1,
1373                                     const Handle(Adaptor3d_HSurface)&  S2,
1374                                     const Handle(Adaptor3d_TopolTool)& D2,
1375                                     const Standard_Real U1,
1376                                     const Standard_Real V1,
1377                                     const Standard_Real U2,
1378                                     const Standard_Real V2,
1379                                     const Standard_Real TolArc,
1380                                     const Standard_Real TolTang)
1381 {
1382   myTolArc = TolArc;
1383   myTolTang = TolTang;
1384   if(myFleche == 0.0) {
1385 #if DEBUG
1386     //cout<<" -- IntPatch_Intersection::myFleche fixe par defaut a 0.01 --"<<endl;
1387     //cout<<" -- Utiliser la Methode SetTolerances( ... ) "<<endl;
1388 #endif
1389     myFleche = 0.01;
1390   }
1391   if(myUVMaxStep==0.0) {
1392 #if DEBUG
1393     //cout<<" -- IntPatch_Intersection::myUVMaxStep fixe par defaut a 0.01 --"<<endl;
1394     //cout<<" -- Utiliser la Methode SetTolerances( ... ) "<<endl;
1395 #endif
1396     myUVMaxStep = 0.01;
1397   }
1398
1399   done = Standard_False;
1400   spnt.Clear();
1401   slin.Clear();
1402
1403   empt = Standard_True;
1404   tgte = Standard_False;
1405   oppo = Standard_False;
1406
1407   const GeomAbs_SurfaceType typs1 = S1->GetType();
1408   const GeomAbs_SurfaceType typs2 = S2->GetType();
1409   
1410   if(   typs1==GeomAbs_Plane 
1411      || typs1==GeomAbs_Cylinder
1412      || typs1==GeomAbs_Sphere
1413      || typs1==GeomAbs_Cone
1414      || typs2==GeomAbs_Plane 
1415      || typs2==GeomAbs_Cylinder
1416      || typs2==GeomAbs_Sphere
1417      || typs2==GeomAbs_Cone)
1418   {
1419     myIsStartPnt = Standard_True;
1420     myU1Start = U1; myV1Start = V1; myU2Start = U2; myV2Start = V2;
1421     Perform(S1,D1,S2,D2,TolArc,TolTang);
1422     myIsStartPnt = Standard_False;
1423   }
1424   else
1425   {
1426     IntPatch_PrmPrmIntersection interpp;
1427     interpp.Perform(S1,D1,S2,D2,U1,V1,U2,V2,TolArc,TolTang,myFleche,myUVMaxStep);
1428     if (interpp.IsDone())
1429         {
1430       done = Standard_True;
1431       tgte = Standard_False;
1432       empt = interpp.IsEmpty();
1433       const Standard_Integer nblm = interpp.NbLines();
1434       Standard_Integer i = 1;
1435       for (; i<=nblm; i++) slin.Append(interpp.Line(i));
1436     }
1437   }
1438 }
1439 //======================================================================
1440 #include <IntPatch_IType.hxx>
1441 #include <IntPatch_LineConstructor.hxx>
1442 #include <Handle_Adaptor2d_HCurve2d.hxx>
1443 #define MAXR 200
1444
1445
1446 //void IntPatch_Intersection__MAJ_R(Handle_Adaptor2d_HCurve2d *R1,
1447 //                                   Handle_Adaptor2d_HCurve2d *R2,
1448 //                                   int *NR1,
1449 //                                   int *NR2,
1450 //                                   Standard_Integer nbR1,
1451 //                                   Standard_Integer nbR2,
1452 //                                   const IntPatch_Point& VTX)
1453 void IntPatch_Intersection__MAJ_R(Handle_Adaptor2d_HCurve2d *,
1454                                      Handle_Adaptor2d_HCurve2d *,
1455                                      int *,
1456                                      int *,
1457                                      Standard_Integer ,
1458                                      Standard_Integer ,
1459                                      const IntPatch_Point& )
1460
1461   /*
1462   if(VTX.IsOnDomS1()) { 
1463     
1464     //-- long unsigned ptr= *((long unsigned *)(((Handle_Standard_Transient *)(&(VTX.ArcOnS1())))));
1465     for(Standard_Integer i=0; i<nbR1;i++) { 
1466       if(VTX.ArcOnS1()==R1[i]) { 
1467         NR1[i]++;
1468         printf("\n ******************************");
1469         return;
1470       }
1471     }
1472     printf("\n R Pas trouvee  (IntPatch)\n");
1473     
1474   }
1475   */
1476 }
1477
1478
1479 //void IntPatch_Intersection::Dump(const Standard_Integer Mode,
1480 void IntPatch_Intersection::Dump(const Standard_Integer ,
1481                                     const Handle(Adaptor3d_HSurface)&  S1,
1482                                     const Handle(Adaptor3d_TopolTool)& D1,
1483                                     const Handle(Adaptor3d_HSurface)&  S2,
1484                                     const Handle(Adaptor3d_TopolTool)& D2) const 
1485
1486   
1487   //-- ----------------------------------------------------------------------
1488   //--  construction de la liste des restrictions & vertex 
1489   //--
1490   int NR1[MAXR],NR2[MAXR];
1491   Handle_Adaptor2d_HCurve2d R1[MAXR],R2[MAXR];
1492   Standard_Integer nbR1=0,nbR2=0;
1493   for(D1->Init();D1->More() && nbR1<MAXR; D1->Next()) { 
1494     R1[nbR1]=D1->Value(); 
1495     NR1[nbR1]=0;
1496     nbR1++;
1497   }
1498   for(D2->Init();D2->More() && nbR2<MAXR; D2->Next()) { 
1499     R2[nbR2]=D2->Value();
1500     NR2[nbR2]=0;
1501     nbR2++;
1502   }
1503   
1504   printf("\nDUMP_INT:  ----empt:%2d  tgte:%2d  oppo:%2d ---------------------------------",empt,tgte,empt);
1505   Standard_Integer i,j,nbr1,nbr2,nbgl,nbgc,nbge,nbgp,nbgh,nbl,nbr,nbg,nbw,nba;
1506   nbl=nbr=nbg=nbw=nba=nbgl=nbge=nbr1=nbr2=nbgc=nbgp=nbgh=0;
1507   nbl=NbLines();
1508   for(i=1;i<=nbl;i++) { 
1509     const Handle(IntPatch_Line)& line=Line(i);
1510     const IntPatch_IType IType=line->ArcType();
1511     if(IType == IntPatch_Walking) nbw++;
1512     else     if(IType == IntPatch_Restriction) { 
1513       nbr++;
1514       Handle(IntPatch_RLine)& rlin =
1515         *((Handle(IntPatch_RLine) *)&line);
1516       if(rlin->IsArcOnS1()) nbr1++;
1517       if(rlin->IsArcOnS2()) nbr2++;
1518     }
1519     else     if(IType == IntPatch_Analytic) nba++;
1520     else     { nbg++; 
1521                if(IType == IntPatch_Lin) nbgl++;
1522                else if(IType == IntPatch_Circle) nbgc++;
1523                else if(IType == IntPatch_Parabola) nbgp++;
1524                else if(IType == IntPatch_Hyperbola) nbgh++;
1525                else if(IType == IntPatch_Ellipse) nbge++;
1526              }
1527   }
1528   
1529   
1530   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)",
1531          nbl,nbw,nbr,nbr1,nbr2,nba,nbg,nbgl,nbgc,nbge,nbgh,nbgp);
1532   
1533   IntPatch_LineConstructor LineConstructor(2);
1534   
1535   Standard_Integer nbllc=0;
1536   nbw=nbr=nbg=nba=0;
1537   Standard_Integer nbva,nbvw,nbvr,nbvg;
1538   nbva=nbvr=nbvw=nbvg=0;
1539   for (j=1; j<=nbl; j++) {
1540     Standard_Integer v,nbvtx;
1541     const Handle(IntPatch_Line)& intersLinej = Line(j);
1542     Standard_Integer NbLines;
1543     LineConstructor.Perform(SequenceOfLine(),intersLinej,S1,D1,S2,D2,1e-7);
1544     NbLines = LineConstructor.NbLines();
1545     
1546     for(Standard_Integer k=1;k<=NbLines;k++) { 
1547       nbllc++;
1548       const Handle(IntPatch_Line)& LineK = LineConstructor.Line(k);
1549       if (LineK->ArcType() == IntPatch_Analytic) { 
1550         Handle(IntPatch_ALine)& alin =
1551           *((Handle(IntPatch_ALine) *)&LineK);
1552         nbvtx=alin->NbVertex();
1553         nbva+=nbvtx;    nba++;
1554         for(v=1;v<=nbvtx;v++) { 
1555           IntPatch_Intersection__MAJ_R(R1,R2,NR1,NR2,nbR1,nbR2,alin->Vertex(v));
1556         }
1557       }
1558       else if (LineK->ArcType() == IntPatch_Restriction) {
1559         Handle(IntPatch_RLine)& rlin =
1560           *((Handle(IntPatch_RLine) *)&LineK);
1561         nbvtx=rlin->NbVertex();
1562         nbvr+=nbvtx;    nbr++;
1563         for(v=1;v<=nbvtx;v++) { 
1564           IntPatch_Intersection__MAJ_R(R1,R2,NR1,NR2,nbR1,nbR2,rlin->Vertex(v));
1565         }
1566       }
1567       else if (LineK->ArcType() == IntPatch_Walking) {
1568         Handle(IntPatch_WLine)& wlin =
1569           *((Handle(IntPatch_WLine) *)&LineK);
1570         nbvtx=wlin->NbVertex();
1571         nbvw+=nbvtx;    nbw++;
1572         for(v=1;v<=nbvtx;v++) { 
1573           IntPatch_Intersection__MAJ_R(R1,R2,NR1,NR2,nbR1,nbR2,wlin->Vertex(v));
1574         }
1575       }
1576       else { 
1577         Handle(IntPatch_GLine)& glin =
1578           *((Handle(IntPatch_GLine) *)&LineK);
1579         nbvtx=glin->NbVertex();
1580         nbvg+=nbvtx;    nbg++;
1581         for(v=1;v<=nbvtx;v++) { 
1582           IntPatch_Intersection__MAJ_R(R1,R2,NR1,NR2,nbR1,nbR2,glin->Vertex(v));
1583         }
1584       }
1585     }
1586   }
1587   printf("\nDUMP_LC :Lines:%2d WLin:%2d Restr:%2d Ana:%2d Geom:%2d",
1588          nbllc,nbw,nbr,nba,nbg);
1589   printf("\nDUMP_LC :vtx          :%2d     r:%2d    :%2d     :%2d",
1590          nbvw,nbvr,nbva,nbvg);
1591
1592
1593
1594    printf("\n");
1595 }