0024911: Avoid using virtual functions in NCollection classes
[occt.git] / src / GeomInt / GeomInt_IntSS_1.cxx
1 // Created on: 1995-01-27
2 // Created by: Jacques GOUSSARD
3 // Copyright (c) 1995-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #include <GeomInt_IntSS.ixx>
18
19 #include <Precision.hxx>
20 #include <gp_Pnt.hxx>
21 #include <gp_Pnt2d.hxx>
22 #include <gp_Pln.hxx>
23
24 #include <TColStd_Array1OfReal.hxx>
25 #include <TColStd_Array1OfInteger.hxx>
26 #include <TColStd_Array1OfListOfInteger.hxx>
27 #include <TColStd_SequenceOfReal.hxx>
28 #include <TColStd_ListOfInteger.hxx>
29 #include <TColStd_ListIteratorOfListOfInteger.hxx>
30
31 #include <TColgp_Array1OfPnt.hxx>
32 #include <TColgp_Array1OfPnt2d.hxx>
33
34 #include <Adaptor3d_TopolTool.hxx>
35 #include <IntPatch_Line.hxx>
36 #include <IntPatch_WLine.hxx>
37 #include <IntPatch_GLine.hxx>
38 #include <IntPatch_RLine.hxx>
39 #include <IntPatch_ALineToWLine.hxx>
40 #include <IntPatch_IType.hxx>
41 #include <NCollection_IncAllocator.hxx>
42 #include <NCollection_List.hxx>
43 #include <NCollection_LocalArray.hxx>
44 #include <NCollection_StdAllocator.hxx>
45 #include <vector>
46
47 #include <AppParCurves_MultiBSpCurve.hxx>
48
49 #include <GeomAbs_SurfaceType.hxx>
50 #include <GeomAbs_CurveType.hxx>
51
52 #include <GeomAdaptor.hxx>
53 #include <Geom_Curve.hxx>
54 #include <Geom_Line.hxx>
55 #include <Geom_Parabola.hxx>
56 #include <Geom_Hyperbola.hxx>
57 #include <Geom_TrimmedCurve.hxx>
58 #include <Geom_Circle.hxx>
59 #include <Geom_Ellipse.hxx>
60 #include <Geom_BSplineCurve.hxx>
61
62 #include <Geom2dAdaptor.hxx>
63 #include <Adaptor2d_HCurve2d.hxx>
64 #include <Geom2d_Curve.hxx>
65 #include <Geom2d_BSplineCurve.hxx>
66 #include <Geom2d_TrimmedCurve.hxx>
67
68 #include <GeomLib_CheckBSplineCurve.hxx>
69 #include <GeomLib_Check2dBSplineCurve.hxx>
70 #include <GeomProjLib.hxx>
71 #include <Approx_CurveOnSurface.hxx>
72
73 #include <ElSLib.hxx>
74
75 #include <GeomInt_WLApprox.hxx>
76 #include <Extrema_ExtPS.hxx>
77
78 static
79   void BuildPCurves (Standard_Real f,
80                      Standard_Real l,
81                      Standard_Real& Tol,
82                      const Handle (Geom_Surface)& S,
83                      const Handle (Geom_Curve)&   C,
84                      Handle (Geom2d_Curve)& C2d);
85
86 static
87   void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
88                   const Handle(GeomAdaptor_HSurface)& HS2,
89                   const gp_Pnt& Ptref,
90                   Standard_Real& U1,
91                   Standard_Real& V1,
92                   Standard_Real& U2,
93                   Standard_Real& V2);
94
95 static 
96   Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
97                                    const Standard_Integer ideb,
98                                    const Standard_Integer ifin);
99
100 static
101   Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
102                                             const Standard_Integer                         ideb,
103                                             const Standard_Integer                         ifin,
104                                             const Standard_Boolean                         onFirst);
105
106 static 
107   Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine);
108
109
110 static
111   Standard_Boolean DecompositionOfWLine  (const Handle(IntPatch_WLine)& theWLine,
112                                           const Handle(GeomAdaptor_HSurface)&            theSurface1, 
113                                           const Handle(GeomAdaptor_HSurface)&            theSurface2,
114                                           const Standard_Real aTolSum, 
115                                           const GeomInt_LineConstructor&                 theLConstructor,
116                                           IntPatch_SequenceOfLine&                       theNewLines);
117
118 static
119   Standard_Real AdjustPeriodic(const Standard_Real theParameter,
120                                const Standard_Real parmin,
121                                const Standard_Real parmax,
122                                const Standard_Real thePeriod,
123                                Standard_Real&      theOffset) ;
124 static 
125   Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
126                                      const Standard_Real theFirstBoundary,
127                                      const Standard_Real theSecondBoundary,
128                                      const Standard_Real theResolution,
129                                      Standard_Boolean&   IsOnFirstBoundary);
130
131 static 
132   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
133                              const gp_Pnt2d&     theLastPoint,
134                              const Standard_Real theUmin, 
135                              const Standard_Real theUmax,
136                              const Standard_Real theVmin,
137                              const Standard_Real theVmax,
138                              gp_Pnt2d&           theNewPoint);
139
140
141 static
142   void AdjustUPeriodic (const Handle(Geom_Surface)& aS,
143                         Handle(Geom2d_Curve)& aC2D);
144 static
145   void GetQuadric(const Handle(GeomAdaptor_HSurface)& HS1,
146                   IntSurf_Quadric& quad1);
147
148 static
149     void TreatRLine(const Handle(IntPatch_RLine)& theRL, 
150                     const Handle(GeomAdaptor_HSurface)& theHS1, 
151                     const Handle(GeomAdaptor_HSurface)& theHS2, 
152                     const Standard_Boolean theApprox1, 
153                     const Standard_Boolean theApprox2, 
154                     Handle(Geom_Curve)& theC3d,
155                     Handle(Geom2d_Curve)& theC2d1, 
156                     Handle(Geom2d_Curve)& theC2d2, 
157                     Standard_Real& theTolReached);
158
159 //
160
161
162 class ProjectPointOnSurf
163 {
164  public:
165   ProjectPointOnSurf() : myIsDone (Standard_False),myIndex(0) {}
166   void Init(const Handle(Geom_Surface)& Surface,
167             const Standard_Real Umin,
168             const Standard_Real Usup,
169             const Standard_Real Vmin,
170             const Standard_Real Vsup);
171   void Init ();
172   void Perform(const gp_Pnt& P);
173   Standard_Boolean IsDone () const { return myIsDone; }
174   void LowerDistanceParameters(Standard_Real& U, Standard_Real& V ) const;
175   Standard_Real LowerDistance() const;
176  protected:
177   Standard_Boolean myIsDone;
178   Standard_Integer myIndex;
179   Extrema_ExtPS myExtPS;
180   GeomAdaptor_Surface myGeomAdaptor;
181 };
182
183 //=======================================================================
184 //function : Init
185 //purpose  : 
186 //=======================================================================
187   void ProjectPointOnSurf::Init (const Handle(Geom_Surface)& Surface,
188                                  const Standard_Real       Umin,
189                                  const Standard_Real       Usup,
190                                  const Standard_Real       Vmin,
191                                  const Standard_Real       Vsup )
192 {
193   const Standard_Real Tolerance = Precision::PConfusion();
194   //
195   myGeomAdaptor.Load(Surface, Umin,Usup,Vmin,Vsup);
196   myExtPS.Initialize(myGeomAdaptor, Umin, Usup, Vmin, Vsup, Tolerance, Tolerance);
197   myIsDone = Standard_False;
198 }
199 //=======================================================================
200 //function : Init
201 //purpose  : 
202 //=======================================================================
203   void ProjectPointOnSurf::Init ()
204 {
205   myIsDone = myExtPS.IsDone() && (myExtPS.NbExt() > 0);
206   if (myIsDone) {
207     // evaluate the lower distance and its index;
208     Standard_Real Dist2Min = myExtPS.SquareDistance(1);
209     myIndex = 1;
210     for (Standard_Integer i = 2; i <= myExtPS.NbExt(); i++)
211     {
212       const Standard_Real Dist2 = myExtPS.SquareDistance(i);
213       if (Dist2 < Dist2Min) {
214         Dist2Min = Dist2;
215         myIndex = i;
216       }
217     }
218   }
219 }
220 //=======================================================================
221 //function : Perform
222 //purpose  : 
223 //=======================================================================
224 void ProjectPointOnSurf::Perform(const gp_Pnt& P)
225 {
226   myExtPS.Perform(P);
227   Init ();
228 }
229 //=======================================================================
230 //function : LowerDistanceParameters
231 //purpose  : 
232 //=======================================================================
233 void ProjectPointOnSurf::LowerDistanceParameters (Standard_Real& U,
234                                                   Standard_Real& V ) const
235 {
236   StdFail_NotDone_Raise_if(!myIsDone, "GeomInt_IntSS::ProjectPointOnSurf::LowerDistanceParameters");
237   (myExtPS.Point(myIndex)).Parameter(U,V);
238 }
239 //=======================================================================
240 //function : LowerDistance
241 //purpose  : 
242 //=======================================================================
243 Standard_Real ProjectPointOnSurf::LowerDistance() const
244 {
245   StdFail_NotDone_Raise_if(!myIsDone, "GeomInt_IntSS::ProjectPointOnSurf::LowerDistance");
246   return sqrt(myExtPS.SquareDistance(myIndex));
247 }
248 //
249
250 //=======================================================================
251 //function : MakeCurve
252 //purpose  : 
253 //=======================================================================
254   void GeomInt_IntSS::MakeCurve(const Standard_Integer Index,
255                                 const Handle(Adaptor3d_TopolTool) & dom1,
256                                 const Handle(Adaptor3d_TopolTool) & dom2,
257                                 const Standard_Real Tol,
258                                 const Standard_Boolean Approx,
259                                 const Standard_Boolean ApproxS1,
260                                 const Standard_Boolean ApproxS2)
261
262 {
263   Standard_Boolean myApprox1, myApprox2, myApprox;
264   Standard_Real Tolpc, myTolApprox;
265   IntPatch_IType typl;
266   Handle(Geom2d_BSplineCurve) H1;
267   Handle(Geom_Surface) aS1, aS2;
268   //
269   Tolpc = Tol;
270   myApprox=Approx;
271   myApprox1=ApproxS1;
272   myApprox2=ApproxS2;
273   myTolApprox=0.0000001;
274   //
275   aS1=myHS1->ChangeSurface().Surface();
276   aS2=myHS2->ChangeSurface().Surface();
277   //
278   Handle(IntPatch_Line) L = myIntersector.Line(Index);
279   typl = L->ArcType();
280   //
281   if(typl==IntPatch_Walking) {
282     Handle(IntPatch_Line) anewL = 
283       ComputePurgedWLine(Handle(IntPatch_WLine)::DownCast(L));
284     //
285     if(anewL.IsNull()) {
286       return;
287     }
288     L = anewL;
289   }
290   //
291   // Line Constructor
292   myLConstruct.Perform(L);
293   if (!myLConstruct.IsDone() || myLConstruct.NbParts() <= 0) {
294     return;
295   }
296   // Do the Curve
297   Standard_Boolean ok;
298   Standard_Integer i, j,  aNbParts;
299   Standard_Real fprm, lprm;
300   Handle(Geom_Curve) newc;
301
302   switch (typl) {
303     //########################################  
304     // Line, Parabola, Hyperbola
305     //########################################  
306   case IntPatch_Lin:
307   case IntPatch_Parabola: 
308   case IntPatch_Hyperbola: {
309     if (typl == IntPatch_Lin) {
310       newc=new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
311     }
312     else if (typl == IntPatch_Parabola) {
313       newc=new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
314     }
315     else if (typl == IntPatch_Hyperbola) {
316       newc=new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
317     }
318     //
319     aNbParts=myLConstruct.NbParts();
320     for (i=1; i<=aNbParts; i++) {
321       myLConstruct.Part(i, fprm, lprm);
322
323       if (!Precision::IsNegativeInfinite(fprm) && 
324         !Precision::IsPositiveInfinite(lprm)) {
325           Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
326           sline.Append(aCT3D);
327           //
328           if(myApprox1) { 
329             Handle (Geom2d_Curve) C2d;
330             BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
331             if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
332               myTolReached2d=Tolpc;
333             }
334             slineS1.Append(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
335           }
336           else { 
337             slineS1.Append(H1);
338           }
339           //
340           if(myApprox2) { 
341             Handle (Geom2d_Curve) C2d;
342             BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
343             if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
344               myTolReached2d=Tolpc;
345             }
346             //
347             slineS2.Append(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
348           }
349           else { 
350             slineS2.Append(H1);
351           }
352       } // if (!Precision::IsNegativeInfinite(fprm) &&  !Precision::IsPositiveInfinite(lprm))
353
354       else {
355         GeomAbs_SurfaceType typS1 = myHS1->Surface().GetType();
356         GeomAbs_SurfaceType typS2 = myHS2->Surface().GetType();
357         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
358           typS1 == GeomAbs_OffsetSurface ||
359           typS1 == GeomAbs_SurfaceOfRevolution ||
360           typS2 == GeomAbs_SurfaceOfExtrusion ||
361           typS2 == GeomAbs_OffsetSurface ||
362           typS2 == GeomAbs_SurfaceOfRevolution) {
363             sline.Append(newc);
364             slineS1.Append(H1);
365             slineS2.Append(H1);
366             continue;
367         }
368         Standard_Boolean bFNIt, bLPIt;
369         Standard_Real aTestPrm, dT=100.;
370         Standard_Real u1, v1, u2, v2, TolX;
371         //
372         bFNIt=Precision::IsNegativeInfinite(fprm);
373         bLPIt=Precision::IsPositiveInfinite(lprm);
374
375         aTestPrm=0.;
376
377         if (bFNIt && !bLPIt) {
378           aTestPrm=lprm-dT;
379         }
380         else if (!bFNIt && bLPIt) {
381           aTestPrm=fprm+dT;
382         }
383         //
384         gp_Pnt ptref(newc->Value(aTestPrm));
385         //
386         TolX = Precision::Confusion();
387         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
388         ok = (dom1->Classify(gp_Pnt2d(u1, v1), TolX) != TopAbs_OUT);
389         if(ok) { 
390           ok = (dom2->Classify(gp_Pnt2d(u2,v2),TolX) != TopAbs_OUT); 
391         }
392         if (ok) {
393           sline.Append(newc);
394           slineS1.Append(H1);
395           slineS2.Append(H1);
396         }
397       }
398     }// end of for (i=1; i<=myLConstruct.NbParts(); i++)
399   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
400   break;
401
402                            //########################################  
403                            // Circle and Ellipse
404                            //########################################  
405   case IntPatch_Circle: 
406   case IntPatch_Ellipse: {
407
408     if (typl == IntPatch_Circle) {
409       newc = new Geom_Circle
410         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
411     }
412     else { 
413       newc = new Geom_Ellipse
414         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
415     }
416     //
417     Standard_Real aPeriod, aRealEpsilon;
418     //
419     aRealEpsilon=RealEpsilon();
420     aPeriod=M_PI+M_PI;
421     //
422     aNbParts=myLConstruct.NbParts();
423     //
424     for (i=1; i<=aNbParts; i++) {
425       myLConstruct.Part(i, fprm, lprm);
426       //
427       if (Abs(fprm) > aRealEpsilon || Abs(lprm-aPeriod) > aRealEpsilon) {
428         //==============================================
429         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
430         //
431         sline.Append(aTC3D);
432         //
433         fprm=aTC3D->FirstParameter();
434         lprm=aTC3D->LastParameter ();
435         ////    
436         if(myApprox1) { 
437           Handle (Geom2d_Curve) C2d;
438           BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
439           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
440             myTolReached2d=Tolpc;
441           }
442           slineS1.Append(C2d);
443         }
444         else { //// 
445           slineS1.Append(H1);
446         }
447         //
448         if(myApprox2) { 
449           Handle (Geom2d_Curve) C2d;
450           BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
451           if(Tolpc>myTolReached2d || myTolReached2d==0) { 
452             myTolReached2d=Tolpc;
453           }
454           slineS2.Append(C2d);
455         }
456         else { 
457           slineS2.Append(H1);  
458         }
459         //==============================================        
460       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
461       //
462       else {//  on regarde si on garde
463         //
464         if (aNbParts==1) {
465           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
466             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
467             //
468             sline.Append(aTC3D);
469             fprm=aTC3D->FirstParameter();
470             lprm=aTC3D->LastParameter ();
471
472             if(myApprox1) { 
473               Handle (Geom2d_Curve) C2d;
474               BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
475               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
476                 myTolReached2d=Tolpc;
477               }
478               slineS1.Append(C2d);
479             }
480             else { //// 
481               slineS1.Append(H1);
482             }
483
484             if(myApprox2) { 
485               Handle (Geom2d_Curve) C2d;
486               BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
487               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
488                 myTolReached2d=Tolpc;
489               }
490               slineS2.Append(C2d);
491             }
492             else { 
493               slineS2.Append(H1);
494             }
495             break;
496           }
497         }
498         //
499         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, TolX;
500         //
501         aTwoPIdiv17=2.*M_PI/17.;
502         //
503         for (j=0; j<=17; j++) {
504           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
505           TolX = Precision::Confusion();
506
507           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
508           ok = (dom1->Classify(gp_Pnt2d(u1,v1),TolX) != TopAbs_OUT);
509           if(ok) { 
510             ok = (dom2->Classify(gp_Pnt2d(u2,v2),TolX) != TopAbs_OUT);
511           }
512           if (ok) {
513             sline.Append(newc);
514             //==============================================
515             if(myApprox1) { 
516               Handle (Geom2d_Curve) C2d;
517               BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
518               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
519                 myTolReached2d=Tolpc;
520               }
521               slineS1.Append(C2d);
522             }
523             else { 
524               slineS1.Append(H1);  
525             }
526
527             if(myApprox2) { 
528               Handle (Geom2d_Curve) C2d;
529               BuildPCurves(fprm, lprm, Tolpc,myHS2->ChangeSurface().Surface(), newc, C2d);
530               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
531                 myTolReached2d=Tolpc;
532               } 
533               slineS2.Append(C2d);
534             }
535             else { 
536               slineS2.Append(H1);  
537             }
538             break;
539           }//  end of if (ok) {
540         }//  end of for (Standard_Integer j=0; j<=17; j++)
541       }//  end of else { on regarde si on garde
542     }// for (i=1; i<=myLConstruct.NbParts(); i++)
543   }// IntPatch_Circle: IntPatch_Ellipse
544   break;
545
546                          //########################################  
547                          // Analytic
548                          //######################################## 
549   case IntPatch_Analytic: {
550     IntSurf_Quadric quad1,quad2;
551     //
552     GetQuadric(myHS1, quad1);
553     GetQuadric(myHS2, quad2);
554     //=========
555     IntPatch_ALineToWLine convert (quad1, quad2);
556
557     if (!myApprox) {
558       Handle(Geom2d_BSplineCurve) aH1, aH2;
559       //
560       aNbParts=myLConstruct.NbParts();
561       for (i=1; i<=aNbParts; i++) {
562         myLConstruct.Part(i, fprm, lprm);
563         Handle(IntPatch_WLine) WL = 
564           convert.MakeWLine(Handle(IntPatch_ALine)::DownCast(L), fprm, lprm);
565         //
566         if(myApprox1) {
567           aH1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
568         }
569
570         if(myApprox2) {
571           aH2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
572         }
573         sline.Append(MakeBSpline(WL,1,WL->NbPnts()));
574         slineS1.Append(aH1);
575         slineS2.Append(aH2);
576       }
577     } // if (!myApprox)
578
579     else { // myApprox=TRUE
580       GeomInt_WLApprox theapp3d;
581       Standard_Real tol2d = myTolApprox;
582       //        
583       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
584
585       aNbParts=myLConstruct.NbParts();
586       for (i=1; i<=aNbParts; i++) {
587         myLConstruct.Part(i, fprm, lprm);
588         Handle(IntPatch_WLine) WL = 
589           convert.MakeWLine(Handle(IntPatch_ALine):: DownCast(L),fprm,lprm);
590
591         theapp3d.Perform(myHS1,myHS2,WL,Standard_True,myApprox1,myApprox2, 1, WL->NbPnts());
592         if (!theapp3d.IsDone()) {
593           //
594           Handle(Geom2d_BSplineCurve) aH1, aH2;
595
596           if(myApprox1) {
597             aH1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
598           }
599
600           if(myApprox2) {
601             aH2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
602           }
603           sline.Append(MakeBSpline(WL,1,WL->NbPnts()));
604           slineS1.Append(aH1);
605           slineS2.Append(aH1);
606         }
607         //
608         else {
609           if(myApprox1 || myApprox2) { 
610             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0) { 
611               myTolReached2d = theapp3d.TolReached2d();
612             }
613           }
614
615           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) { 
616             myTolReached3d = theapp3d.TolReached3d();
617           }
618
619           Standard_Integer aNbMultiCurves, nbpoles;
620           aNbMultiCurves=theapp3d.NbMultiCurves();
621           for (j=1; j<=aNbMultiCurves; j++) {
622             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
623
624             nbpoles = mbspc.NbPoles();
625             TColgp_Array1OfPnt tpoles(1, nbpoles);
626             mbspc.Curve(1, tpoles);
627             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
628               mbspc.Knots(),
629               mbspc.Multiplicities(),
630               mbspc.Degree());
631
632             GeomLib_CheckBSplineCurve Check(BS, myTolCheck, myTolAngCheck);
633             Check.FixTangent(Standard_True,Standard_True);
634             // 
635             sline.Append(BS);
636             //
637             if(myApprox1) { 
638               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
639               mbspc.Curve(2,tpoles2d);
640               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
641                 mbspc.Knots(),
642                 mbspc.Multiplicities(),
643                 mbspc.Degree());
644
645               GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
646               newCheck.FixTangent(Standard_True,Standard_True);
647               slineS1.Append(BS2);              
648             }
649             else {
650               slineS1.Append(H1);
651             }
652
653             if(myApprox2) { 
654               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
655               Standard_Integer TwoOrThree;
656               TwoOrThree=myApprox1 ? 3 : 2;
657               mbspc.Curve(TwoOrThree, tpoles2d);
658               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
659                 mbspc.Knots(),
660                 mbspc.Multiplicities(),
661                 mbspc.Degree());
662
663               GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
664               newCheck.FixTangent(Standard_True,Standard_True);
665               //        
666               slineS2.Append(BS2);
667             }
668             else { 
669               slineS2.Append(H1);
670             }
671             // 
672           }// for (j=1; j<=aNbMultiCurves; j++) {
673         }// else from if (!theapp3d.IsDone())
674       }// for (i=1; i<=aNbParts; i++) {
675     }// else { // myApprox=TRUE
676   }// case IntPatch_Analytic:
677   break;
678
679                           //########################################  
680                           // Walking
681                           //######################################## 
682   case IntPatch_Walking:{
683     Handle(IntPatch_WLine) WL = 
684       Handle(IntPatch_WLine)::DownCast(L);
685     //
686     Standard_Integer ifprm, ilprm;
687     //
688     if (!myApprox) {
689       aNbParts=myLConstruct.NbParts();
690       for (i=1; i<=aNbParts; i++) {
691         myLConstruct.Part(i, fprm, lprm);
692         ifprm=(Standard_Integer)fprm;
693         ilprm=(Standard_Integer)lprm;
694         //        
695         Handle(Geom2d_BSplineCurve) aH1, aH2;
696
697         if(myApprox1) {
698           aH1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
699         }
700         if(myApprox2) {
701           aH2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
702         }
703         //
704         Handle(Geom_Curve) aBSp=MakeBSpline(WL, ifprm, ilprm);
705         //      
706         sline.Append(aBSp);
707         slineS1.Append(aH1);
708         slineS2.Append(aH2);
709       }
710     }
711     //
712     else {
713       Standard_Boolean bIsDecomposited;
714       Standard_Integer nbiter, aNbSeqOfL;
715       GeomInt_WLApprox theapp3d;
716       IntPatch_SequenceOfLine aSeqOfL;
717       Standard_Real tol2d, aTolSS;
718       //        
719       tol2d = myTolApprox;
720       aTolSS=2.e-7;
721       if(myHS1 == myHS2) { 
722         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False);      
723       }
724       else { 
725         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
726       }
727       //
728       bIsDecomposited = 
729         DecompositionOfWLine(WL, myHS1, myHS2, aTolSS, myLConstruct, aSeqOfL);
730       //
731       aNbParts=myLConstruct.NbParts();
732       aNbSeqOfL=aSeqOfL.Length();
733       //
734       nbiter = (bIsDecomposited) ? aNbSeqOfL : aNbParts;
735       //
736       for(i = 1; i <= nbiter; i++) {
737         if(bIsDecomposited) {
738           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
739           ifprm = 1;
740           ilprm = WL->NbPnts();
741         }
742         else {
743           myLConstruct.Part(i, fprm, lprm);
744           ifprm = (Standard_Integer)fprm;
745           ilprm = (Standard_Integer)lprm;
746         }
747         //-- lbr : 
748         //-- Si une des surfaces est un plan , on approxime en 2d
749         //-- sur cette surface et on remonte les points 2d en 3d.
750         GeomAbs_SurfaceType typs1, typs2;
751         typs1 = myHS1->Surface().GetType();
752         typs2 = myHS2->Surface().GetType();
753         //
754         if(typs1 == GeomAbs_Plane) { 
755           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,
756             Standard_True, myApprox2,
757             ifprm,  ilprm);
758         }         
759         else if(typs2 == GeomAbs_Plane) { 
760           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,
761             myApprox1,Standard_True,
762             ifprm,  ilprm);
763         }
764         else { 
765           //
766           if (myHS1 != myHS2){
767             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
768               (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
769
770                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
771                 //Standard_Boolean bUseSurfaces;
772                 //bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
773                 //if (bUseSurfaces) {
774                 //theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False);
775                 //}
776             }
777           }
778           //
779           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,
780             myApprox1,myApprox2,
781             ifprm,  ilprm);
782         }
783
784         if (!theapp3d.IsDone()) {
785           //      
786           Handle(Geom2d_BSplineCurve) aH1, aH2;
787           //      
788           Handle(Geom_Curve) aBSp=MakeBSpline(WL, ifprm, ilprm);
789           if(myApprox1) {
790             aH1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
791           }
792           if(myApprox2) {
793             aH2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
794           }
795           //
796           sline.Append(aBSp);
797           slineS1.Append(aH1);
798           slineS2.Append(aH2);
799         }//if (!theapp3d.IsDone())
800
801         else {
802           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) { 
803             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) { 
804               myTolReached2d = theapp3d.TolReached2d();
805             }
806           }
807           if(typs1==GeomAbs_Plane || typs2==GeomAbs_Plane) { 
808             myTolReached3d = myTolReached2d;
809           }
810           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) { 
811             myTolReached3d = theapp3d.TolReached3d();
812           }
813
814           Standard_Integer aNbMultiCurves, nbpoles;
815           //
816           aNbMultiCurves=theapp3d.NbMultiCurves(); 
817           for (j=1; j<=aNbMultiCurves; j++) {
818             if(typs1 == GeomAbs_Plane) {
819               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
820               nbpoles = mbspc.NbPoles();
821
822               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
823               TColgp_Array1OfPnt   tpoles(1,nbpoles);
824
825               mbspc.Curve(1,tpoles2d);
826               const gp_Pln&  Pln = myHS1->Surface().Plane();
827               //
828               Standard_Integer ik; 
829               for(ik = 1; ik<= nbpoles; ik++) { 
830                 tpoles.SetValue(ik,
831                   ElSLib::Value(tpoles2d.Value(ik).X(),
832                   tpoles2d.Value(ik).Y(),
833                   Pln));
834               }
835               //
836               Handle(Geom_BSplineCurve) BS = 
837                 new Geom_BSplineCurve(tpoles,
838                 mbspc.Knots(),
839                 mbspc.Multiplicities(),
840                 mbspc.Degree());
841               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
842               Check.FixTangent(Standard_True, Standard_True);
843               //        
844               sline.Append(BS);
845               //
846               if(myApprox1) { 
847                 Handle(Geom2d_BSplineCurve) BS1 = 
848                   new Geom2d_BSplineCurve(tpoles2d,
849                   mbspc.Knots(),
850                   mbspc.Multiplicities(),
851                   mbspc.Degree());
852                 GeomLib_Check2dBSplineCurve Check1(BS1,myTolCheck,myTolAngCheck);
853                 Check1.FixTangent(Standard_True,Standard_True);
854                 //      
855                 AdjustUPeriodic (aS1, BS1);  
856                 //
857                 slineS1.Append(BS1);
858               }
859               else {
860                 slineS1.Append(H1);
861               }
862
863               if(myApprox2) { 
864                 mbspc.Curve(2, tpoles2d);
865
866                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
867                   mbspc.Knots(),
868                   mbspc.Multiplicities(),
869                   mbspc.Degree());
870                 GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
871                 newCheck.FixTangent(Standard_True,Standard_True);
872                 //
873                 AdjustUPeriodic (aS2, BS2);  
874                 //
875                 slineS2.Append(BS2); 
876               }
877               else { 
878                 slineS2.Append(H1); 
879               }
880             }//if(typs1 == GeomAbs_Plane) 
881             //
882             else if(typs2 == GeomAbs_Plane) { 
883               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
884               nbpoles = mbspc.NbPoles();
885
886               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
887               TColgp_Array1OfPnt   tpoles(1,nbpoles);
888               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
889               const gp_Pln&  Pln = myHS2->Surface().Plane();
890               //
891               Standard_Integer ik; 
892               for(ik = 1; ik<= nbpoles; ik++) { 
893                 tpoles.SetValue(ik,
894                   ElSLib::Value(tpoles2d.Value(ik).X(),
895                   tpoles2d.Value(ik).Y(),
896                   Pln));
897
898               }
899               //
900               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
901                 mbspc.Knots(),
902                 mbspc.Multiplicities(),
903                 mbspc.Degree());
904               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
905               Check.FixTangent(Standard_True,Standard_True);
906               //        
907               sline.Append(BS);
908               //
909               if(myApprox2) {
910                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
911                   mbspc.Knots(),
912                   mbspc.Multiplicities(),
913                   mbspc.Degree());
914                 GeomLib_Check2dBSplineCurve Check1(BS1,myTolCheck,myTolAngCheck);
915                 Check1.FixTangent(Standard_True,Standard_True);
916                 //
917                 //
918                 AdjustUPeriodic (aS2, BS1);  
919                 //
920                 slineS2.Append(BS1);
921               }
922               else {
923                 slineS2.Append(H1);
924               }
925
926               if(myApprox1) { 
927                 mbspc.Curve(1,tpoles2d);
928                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
929                   mbspc.Knots(),
930                   mbspc.Multiplicities(),
931                   mbspc.Degree());
932                 GeomLib_Check2dBSplineCurve Check2(BS2,myTolCheck,myTolAngCheck);
933                 Check2.FixTangent(Standard_True,Standard_True);
934                 // 
935                 //
936                 AdjustUPeriodic (aS1, BS2);  
937                 //      
938                 slineS1.Append(BS2);
939               }
940               else { 
941                 slineS1.Append(H1);
942               }
943             } // else if(typs2 == GeomAbs_Plane)
944             //
945             else { // typs1!=GeomAbs_Plane && typs2!=GeomAbs_Plane
946               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
947               nbpoles = mbspc.NbPoles();
948               TColgp_Array1OfPnt tpoles(1,nbpoles);
949               mbspc.Curve(1,tpoles);
950               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
951                 mbspc.Knots(),
952                 mbspc.Multiplicities(),
953                 mbspc.Degree());
954               GeomLib_CheckBSplineCurve Check(BS,myTolCheck,myTolAngCheck);
955               Check.FixTangent(Standard_True,Standard_True);
956               //        
957               //Check IsClosed()
958               Standard_Real aDist = Max(BS->StartPoint().XYZ().SquareModulus(),
959                 BS->EndPoint().XYZ().SquareModulus());
960               Standard_Real eps = Epsilon(aDist);
961               if(BS->StartPoint().SquareDistance(BS->EndPoint()) < 2.*eps &&
962                 !BS->IsClosed() && !BS->IsPeriodic())
963               {
964                 //force Closed()
965                 gp_Pnt aPm((BS->Pole(1).XYZ() + BS->Pole(BS->NbPoles()).XYZ()) / 2.);
966                 BS->SetPole(1, aPm);
967                 BS->SetPole(BS->NbPoles(), aPm);
968               }
969               sline.Append(BS);
970
971               if(myApprox1) { 
972                 TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
973                 mbspc.Curve(2,tpoles2d);
974                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
975                   mbspc.Knots(),
976                   mbspc.Multiplicities(),
977                   mbspc.Degree());
978                 GeomLib_Check2dBSplineCurve newCheck(BS1,myTolCheck,myTolAngCheck);
979                 newCheck.FixTangent(Standard_True,Standard_True);
980                 //
981                 AdjustUPeriodic (aS1, BS1);  
982                 //
983                 slineS1.Append(BS1);
984               }
985               else {
986                 slineS1.Append(H1);
987               }
988               if(myApprox2) { 
989                 TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
990                 mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
991                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
992                   mbspc.Knots(),
993                   mbspc.Multiplicities(),
994                   mbspc.Degree());
995                 GeomLib_Check2dBSplineCurve newCheck(BS2,myTolCheck,myTolAngCheck);
996                 newCheck.FixTangent(Standard_True,Standard_True);
997                 //
998                 AdjustUPeriodic (aS2, BS2);  
999                 //
1000                 slineS2.Append(BS2);
1001               }
1002               else { 
1003                 slineS2.Append(H1);
1004               }
1005             }// else { // typs1!=GeomAbs_Plane && typs2!=GeomAbs_Plane
1006           }// for (j=1; j<=aNbMultiCurves; j++
1007         }
1008       }
1009     }
1010   }
1011   break;
1012
1013   case IntPatch_Restriction: 
1014     {
1015       GeomAbs_SurfaceType typS1 = myHS1->Surface().GetType();
1016       GeomAbs_SurfaceType typS2 = myHS2->Surface().GetType();
1017       Standard_Boolean isAnalS1 = Standard_False;
1018       switch (typS1)
1019       {
1020       case GeomAbs_Plane:
1021       case GeomAbs_Cylinder:
1022       case GeomAbs_Sphere:
1023       case GeomAbs_Cone: 
1024       case GeomAbs_Torus: isAnalS1 = Standard_True; break;
1025       default: break;
1026       }
1027
1028       Standard_Integer isAnalS2 = Standard_False;
1029       switch (typS2)
1030       {
1031       case GeomAbs_Plane:
1032       case GeomAbs_Cylinder:
1033       case GeomAbs_Sphere:
1034       case GeomAbs_Cone: 
1035       case GeomAbs_Torus: isAnalS2 = Standard_True; break;
1036       default: break;
1037       }
1038
1039       if(!isAnalS1 || !isAnalS2)
1040       {
1041         Handle(IntPatch_RLine) RL = 
1042           Handle(IntPatch_RLine)::DownCast(L);
1043         Handle(Geom_Curve) aC3d;
1044         Handle(Geom2d_Curve) aC2d1, aC2d2;
1045         Standard_Real aTolReached;
1046         TreatRLine(RL, myHS1, myHS2, myApprox1, myApprox2, aC3d,
1047           aC2d1, aC2d2, aTolReached);
1048         sline.Append(aC3d);
1049         slineS1.Append(aC2d1);
1050         slineS2.Append(aC2d2);
1051       }
1052     }
1053     break;
1054
1055   }
1056 }
1057 //
1058 //=======================================================================
1059 //function : AdjustUPeriodic
1060 //purpose  : 
1061 //=======================================================================
1062 void TreatRLine(const Handle(IntPatch_RLine)& theRL, 
1063                 const Handle(GeomAdaptor_HSurface)& theHS1, 
1064                 const Handle(GeomAdaptor_HSurface)& theHS2, 
1065                 const Standard_Boolean theApprox1, 
1066                 const Standard_Boolean theApprox2, 
1067                 Handle(Geom_Curve)& theC3d,
1068                 Handle(Geom2d_Curve)& theC2d1, 
1069                 Handle(Geom2d_Curve)& theC2d2, 
1070                 Standard_Real& theTolReached)
1071 {
1072   Handle(GeomAdaptor_HSurface) aGAHS;
1073   Handle(Adaptor2d_HCurve2d) anAHC2d;
1074   Standard_Real tf, tl;
1075   gp_Lin2d aL;
1076   // It is assumed that 2d curve is 2d line (rectangular surface domain)
1077   if(theRL->IsArcOnS1())
1078   {
1079     aGAHS = theHS1;
1080     anAHC2d = theRL->ArcOnS1();
1081     theRL->ParamOnS1(tf, tl);
1082     theC2d1 = Geom2dAdaptor::MakeCurve(theRL->ArcOnS1()->Curve2d());
1083   }
1084   else if (theRL->IsArcOnS2())
1085   {
1086     aGAHS = theHS2;
1087     anAHC2d = theRL->ArcOnS2();
1088     theRL->ParamOnS2(tf, tl);
1089     theC2d2 = Geom2dAdaptor::MakeCurve(theRL->ArcOnS2()->Curve2d());
1090   }
1091   else
1092   {
1093     return;
1094   }
1095   //
1096   //To provide sameparameter it is necessary to get 3d curve as
1097   //approximation of curve on surface.
1098   Standard_Integer aMaxDeg = 8;
1099   Standard_Integer aMaxSeg = 1000;
1100   Approx_CurveOnSurface anApp(anAHC2d, aGAHS, tf, tl, Precision::Confusion(),
1101                               GeomAbs_C1, aMaxDeg, aMaxSeg, 
1102                               Standard_True, Standard_False);
1103   if(anApp.HasResult())
1104   {
1105     theC3d = anApp.Curve3d();
1106     theTolReached = anApp.MaxError3d();
1107     Standard_Real aTol = Precision::Confusion();
1108     if(theRL->IsArcOnS1() && theApprox2)
1109     {
1110       Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theHS1->Surface());
1111       BuildPCurves (tf, tl, aTol, 
1112                     aS, theC3d, theC2d2);
1113     }
1114     if(theRL->IsArcOnS2() && theApprox1)
1115     {
1116       Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theHS2->Surface());
1117       BuildPCurves (tf, tl, aTol, 
1118                     aS, theC3d, theC2d1);
1119     }
1120     theTolReached = Max(theTolReached, aTol);
1121   }
1122   //
1123
1124 }
1125 //
1126 //=======================================================================
1127 //function : AdjustUPeriodic
1128 //purpose  : 
1129 //=======================================================================
1130  void AdjustUPeriodic (const Handle(Geom_Surface)& aS, Handle(Geom2d_Curve)& aC2D)
1131 {
1132   if (aC2D.IsNull() || !aS->IsUPeriodic())
1133     return;
1134   //
1135   const Standard_Real aEps=Precision::PConfusion();//1.e-9
1136   const Standard_Real aEpsilon=Epsilon(10.);//1.77e-15 
1137   //
1138   Standard_Real umin,umax,vmin,vmax;
1139   aS->Bounds(umin,umax,vmin,vmax);
1140   const Standard_Real aPeriod = aS->UPeriod();
1141   
1142   const Standard_Real aT1=aC2D->FirstParameter();
1143   const Standard_Real aT2=aC2D->LastParameter();
1144   const Standard_Real aTx=aT1+0.467*(aT2-aT1);
1145   const gp_Pnt2d aPx=aC2D->Value(aTx);
1146   //
1147   Standard_Real aUx=aPx.X();
1148   if (fabs(aUx)<aEpsilon)
1149     aUx=0.;
1150   if (fabs(aUx-aPeriod)<aEpsilon)
1151     aUx=aPeriod;
1152   //
1153   Standard_Real dU=0.;
1154   while(aUx <(umin-aEps)) {
1155     aUx+=aPeriod;
1156     dU+=aPeriod;
1157   }
1158   while(aUx>(umax+aEps)) {
1159     aUx-=aPeriod;
1160     dU-=aPeriod;
1161   }
1162   // 
1163   if (dU!=0.) {
1164     gp_Vec2d aV2D(dU, 0.);
1165     aC2D->Translate(aV2D);
1166   }
1167 }
1168 //
1169 //
1170 //=======================================================================
1171 //function : BuildPCurves
1172 //purpose  : 
1173 //=======================================================================
1174  void BuildPCurves (Standard_Real f, Standard_Real l,
1175                     Standard_Real& Tol,
1176                     const Handle(Geom_Surface)& S,
1177                     const Handle(Geom_Curve)& C,
1178                     Handle(Geom2d_Curve)& C2d)
1179 {
1180   if (C2d.IsNull())
1181   {
1182     C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
1183     if (S->IsUPeriodic())
1184     {
1185       Standard_Real umin,umax,vmin,vmax;
1186       S->Bounds(umin,umax,vmin,vmax);
1187       // Recadre dans le domaine UV de la face
1188       const Standard_Real period = S->UPeriod();
1189       gp_Pnt2d Pf=C2d->Value(f);
1190       //
1191       Standard_Real U0=Pf.X();
1192       //
1193       const Standard_Real aEpsilon=Epsilon(10.);//1.77e-15 
1194       if (fabs(U0)<aEpsilon)
1195         U0=0.;
1196       if (fabs(U0-period)<aEpsilon)
1197         U0=period;
1198       //
1199       const Standard_Real aEps=Precision::PConfusion();//1.e-9
1200       Standard_Real du=0.;
1201       while(U0 <(umin-aEps)) { 
1202         U0+=period;
1203         du+=period;
1204       }
1205       while(U0>(umax+aEps)) { 
1206         U0-=period;
1207         du-=period;
1208       }
1209       if (du!=0.) {
1210         gp_Vec2d T1(du, 0.);
1211         C2d->Translate(T1);
1212       }
1213     }
1214   }
1215 }
1216 //
1217 //=======================================================================
1218 //function : GetQuadric
1219 //purpose  : 
1220 //=======================================================================
1221  void GetQuadric(const Handle(GeomAdaptor_HSurface)& HS1, IntSurf_Quadric& quad1)
1222 {
1223   switch (HS1->Surface().GetType())
1224   {
1225     case GeomAbs_Plane:    quad1.SetValue(HS1->Surface().Plane()); break;
1226     case GeomAbs_Cylinder: quad1.SetValue(HS1->Surface().Cylinder()); break;
1227     case GeomAbs_Cone:     quad1.SetValue(HS1->Surface().Cone()); break;
1228     case GeomAbs_Sphere:   quad1.SetValue(HS1->Surface().Sphere()); break;
1229     case GeomAbs_Torus:    quad1.SetValue(HS1->Surface().Torus()); break;
1230     default: Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve");
1231   }
1232 }
1233
1234 //=======================================================================
1235 //function : Parameters
1236 //purpose  : 
1237 //=======================================================================
1238  void Parameters(const Handle(GeomAdaptor_HSurface)& HS1,
1239                  const Handle(GeomAdaptor_HSurface)& HS2,
1240                  const gp_Pnt& Ptref,
1241                  Standard_Real& U1,
1242                  Standard_Real& V1,
1243                  Standard_Real& U2,
1244                  Standard_Real& V2)
1245 {
1246   IntSurf_Quadric quad1,quad2;
1247   //
1248   GetQuadric(HS1, quad1);
1249   GetQuadric(HS2, quad2);
1250   //
1251   quad1.Parameters(Ptref,U1,V1);
1252   quad2.Parameters(Ptref,U2,V2);
1253 }
1254
1255 //=======================================================================
1256 //function : MakeBSpline
1257 //purpose  : 
1258 //=======================================================================
1259 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
1260                                  const Standard_Integer ideb,
1261                                  const Standard_Integer ifin)
1262 {
1263   const Standard_Integer nbpnt = ifin-ideb+1;
1264   TColgp_Array1OfPnt poles(1,nbpnt);
1265   TColStd_Array1OfReal knots(1,nbpnt);
1266   TColStd_Array1OfInteger mults(1,nbpnt);
1267   Standard_Integer i = 1, ipidebm1 = ideb;
1268   for(; i<=nbpnt; ipidebm1++, i++)
1269   {
1270     poles(i) = WL->Point(ipidebm1).Value();
1271     mults(i) = 1;
1272     knots(i) = i-1;
1273   }
1274   mults(1) = mults(nbpnt) = 2;
1275   return new Geom_BSplineCurve(poles,knots,mults,1);
1276 }
1277
1278 //=======================================================================
1279 //function : MakeBSpline2d
1280 //purpose  : 
1281 //=======================================================================
1282 Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
1283                                           const Standard_Integer ideb,
1284                                           const Standard_Integer ifin,
1285                                           const Standard_Boolean onFirst)
1286 {
1287   const Standard_Integer nbpnt = ifin-ideb+1;
1288   TColgp_Array1OfPnt2d poles(1,nbpnt);
1289   TColStd_Array1OfReal knots(1,nbpnt);
1290   TColStd_Array1OfInteger mults(1,nbpnt);
1291   Standard_Integer i = 1, ipidebm1 = ideb;
1292   for(; i <= nbpnt; ipidebm1++, i++)
1293   {
1294     Standard_Real U, V;
1295     if(onFirst)
1296           theWLine->Point(ipidebm1).ParametersOnS1(U, V);
1297     else
1298           theWLine->Point(ipidebm1).ParametersOnS2(U, V);
1299     poles(i).SetCoord(U, V);
1300     mults(i) = 1;
1301     knots(i) = i-1;
1302   }
1303   mults(1) = mults(nbpnt) = 2;
1304   return new Geom2d_BSplineCurve(poles,knots,mults,1);
1305 }
1306
1307 //=========================================================================
1308 // static function : ComputePurgedWLine
1309 // purpose : Removes equal points (leave one of equal points) from theWLine
1310 //           and recompute vertex parameters.
1311 //           Returns new WLine or null WLine if the number
1312 //           of the points is less than 2.
1313 //=========================================================================
1314 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine) 
1315 {
1316   Handle(IntPatch_WLine) aResult;
1317   Handle(IntPatch_WLine) aLocalWLine;
1318   Handle(IntPatch_WLine) aTmpWLine = theWLine;
1319
1320   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1321   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
1322   Standard_Integer i, k, v, nb, nbvtx;
1323   nbvtx = theWLine->NbVertex();
1324   nb = theWLine->NbPnts();
1325
1326   for(i = 1; i <= nb; i++) {
1327     aLineOn2S->Add(theWLine->Point(i));
1328   }
1329
1330   for(v = 1; v <= nbvtx; v++) {
1331     aLocalWLine->AddVertex(theWLine->Vertex(v));
1332   }
1333   
1334   for(i = 1; i <= aLineOn2S->NbPoints(); i++) {
1335     Standard_Integer aStartIndex = i + 1;
1336     Standard_Integer anEndIndex = i + 5;
1337     nb = aLineOn2S->NbPoints();
1338     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
1339
1340     if((aStartIndex >= nb) || (anEndIndex <= 1)) {
1341       continue;
1342     }
1343     k = aStartIndex;
1344
1345     while(k <= anEndIndex) {
1346       
1347       if(i != k) {
1348         IntSurf_PntOn2S p1 = aLineOn2S->Value(i);
1349         IntSurf_PntOn2S p2 = aLineOn2S->Value(k);
1350         
1351         if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
1352           aTmpWLine = aLocalWLine;
1353           aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
1354
1355           for(v = 1; v <= aTmpWLine->NbVertex(); v++) {
1356             IntPatch_Point aVertex = aTmpWLine->Vertex(v);
1357             Standard_Integer avertexindex = (Standard_Integer)aVertex.ParameterOnLine();
1358
1359             if(avertexindex >= k) {
1360               aVertex.SetParameter(aVertex.ParameterOnLine() - 1.);
1361             }
1362             aLocalWLine->AddVertex(aVertex);
1363           }
1364           aLineOn2S->RemovePoint(k);
1365           anEndIndex--;
1366           continue;
1367         }
1368       }
1369       k++;
1370     }
1371   }
1372
1373   if(aLineOn2S->NbPoints() > 1) {
1374     aResult = aLocalWLine;
1375   }
1376   return aResult;
1377 }
1378 //=======================================================================
1379 //function : DecompositionOfWLine
1380 //purpose  : 
1381 //=======================================================================
1382 Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
1383                                       const Handle(GeomAdaptor_HSurface)&            theSurface1, 
1384                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
1385                                       const Standard_Real aTolSum, 
1386                                       const GeomInt_LineConstructor&                 theLConstructor,
1387                                       IntPatch_SequenceOfLine&                       theNewLines)
1388 {
1389   typedef NCollection_List<Standard_Integer> ListOfInteger;
1390   //have to use std::vector, not NCollection_Vector in order to use copy constructor of
1391   //ListOfInteger which will be created with specific allocator instance
1392   typedef std::vector<ListOfInteger, NCollection_StdAllocator<
1393       ListOfInteger> > ArrayOfListOfInteger;
1394
1395   Standard_Boolean bIsPrevPointOnBoundary, bIsCurrentPointOnBoundary;
1396   Standard_Integer nblines, aNbPnts, aNbParts, pit, i, j, aNbListOfPointIndex;
1397   Standard_Real aTol, umin, umax, vmin, vmax;
1398
1399   //an inc allocator, it will contain wasted space (upon list's Clear()) but it should
1400   //still be faster than the standard allocator, and wasted memory should not be
1401   //significant and will be limited by time span of this function;
1402   //this is a separate allocator from the anIncAlloc below what provides better data
1403   //locality in the latter (by avoiding wastes which will only be in anIdxAlloc)
1404   Handle(NCollection_IncAllocator) anIdxAlloc = new NCollection_IncAllocator();
1405   ListOfInteger aListOfPointIndex (anIdxAlloc);
1406
1407   //GeomAPI_ProjectPointOnSurf aPrj1, aPrj2;
1408   ProjectPointOnSurf aPrj1, aPrj2;
1409   Handle(Geom_Surface) aSurf1,  aSurf2;
1410   //
1411   aNbParts=theLConstructor.NbParts();
1412   aNbPnts=theWLine->NbPnts();
1413   //
1414   if((!aNbPnts) || (!aNbParts)){
1415     return Standard_False;
1416   }
1417   //
1418   Handle(NCollection_IncAllocator) anIncAlloc = new NCollection_IncAllocator();
1419   NCollection_StdAllocator<ListOfInteger> anAlloc (anIncAlloc);
1420   const ListOfInteger aDummy (anIncAlloc); //empty list to be copy constructed from
1421   ArrayOfListOfInteger anArrayOfLines (aNbPnts + 1, aDummy, anAlloc);
1422
1423   NCollection_LocalArray<Standard_Integer> anArrayOfLineTypeArr (aNbPnts + 1);
1424   Standard_Integer* anArrayOfLineType = anArrayOfLineTypeArr;
1425   //
1426   nblines = 0;
1427   aTol = Precision::Confusion();
1428   //
1429   aSurf1 = theSurface1->ChangeSurface().Surface();
1430   aSurf1->Bounds(umin, umax, vmin, vmax);
1431   aPrj1.Init(aSurf1, umin, umax, vmin, vmax);
1432   //
1433   aSurf2 = theSurface2->ChangeSurface().Surface();
1434   aSurf2->Bounds(umin, umax, vmin, vmax);
1435   aPrj2.Init(aSurf2, umin, umax, vmin, vmax);
1436   //
1437   //
1438   bIsPrevPointOnBoundary=Standard_False;
1439   for(pit=1; pit<=aNbPnts; pit++) {
1440     const IntSurf_PntOn2S& aPoint = theWLine->Point(pit);
1441     bIsCurrentPointOnBoundary=Standard_False;
1442     //
1443     // whether aPoint is on boundary or not
1444     //
1445     for(i=0; i<2; i++) {// exploration Surface 1,2 
1446       Handle(GeomAdaptor_HSurface) aGASurface = (!i) ? theSurface1 : theSurface2;
1447       aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
1448       //
1449       for(j=0; j<2; j++) {// exploration of coordinate U,V
1450         Standard_Boolean isperiodic;
1451         //
1452         isperiodic = (!j) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1453         if(!isperiodic) {
1454           continue;
1455         }
1456         //
1457         Standard_Real aResolution, aPeriod, alowerboundary, aupperboundary, U, V;
1458         Standard_Real aParameter, anoffset, anAdjustPar;
1459         Standard_Boolean bIsOnFirstBoundary, bIsPointOnBoundary;
1460         //
1461         aResolution = (!j) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
1462         aPeriod     = (!j) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1463         alowerboundary = (!j) ? umin : vmin;
1464         aupperboundary = (!j) ? umax : vmax;
1465         U=0.;V=0.;//?
1466         //
1467         if(!i){
1468           aPoint.ParametersOnS1(U, V);
1469         }
1470         else{
1471           aPoint.ParametersOnS2(U, V);
1472         }
1473         //
1474         aParameter = (!j) ? U : V;
1475         anoffset=0.;
1476         anAdjustPar=AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1477         //
1478         bIsOnFirstBoundary=Standard_True;
1479         //
1480         bIsPointOnBoundary=
1481           IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
1482         
1483         if(bIsPointOnBoundary) {
1484           bIsCurrentPointOnBoundary = Standard_True;
1485           break;
1486         }
1487       }// for(j=0; j<2; j++)
1488       
1489       if(bIsCurrentPointOnBoundary){
1490         break;
1491       }
1492     }// for(i=0; i<2; i++) 
1493     //
1494     if((bIsCurrentPointOnBoundary != bIsPrevPointOnBoundary)) {
1495
1496       if(!aListOfPointIndex.IsEmpty()) {
1497         nblines++;
1498         anArrayOfLines[nblines] = aListOfPointIndex;
1499         anArrayOfLineType[nblines] = bIsPrevPointOnBoundary;
1500         aListOfPointIndex.Clear();
1501       }
1502       bIsPrevPointOnBoundary = bIsCurrentPointOnBoundary;
1503     }
1504     aListOfPointIndex.Append(pit);
1505   } // for(pit=1; pit<=aNbPnts; pit++)
1506   //
1507   aNbListOfPointIndex=aListOfPointIndex.Extent();
1508   if(aNbListOfPointIndex) {
1509     nblines++;
1510     anArrayOfLines[nblines].Assign (aListOfPointIndex);
1511     anArrayOfLineType[nblines] = bIsPrevPointOnBoundary;
1512     aListOfPointIndex.Clear();
1513   }
1514   //
1515   if(nblines <= 1){
1516     return Standard_False;
1517   }
1518   //
1519   // Correct wlines.begin
1520   Standard_Integer aLineType;
1521   TColStd_Array1OfListOfInteger anArrayOfLineEnds(1, nblines);
1522   Handle(IntSurf_LineOn2S) aSeqOfPntOn2S = new IntSurf_LineOn2S (new NCollection_IncAllocator());
1523   //
1524   for(i = 1; i <= nblines; i++) {
1525     aLineType=anArrayOfLineType[i];
1526     if(aLineType) {
1527       continue;
1528     }
1529     //
1530     const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1531     if(aListOfIndex.Extent() < 2) {
1532       continue;
1533     }
1534     //
1535     TColStd_ListOfInteger aListOfFLIndex;
1536     Standard_Integer aneighbourindex, aLineTypeNeib;
1537     //
1538     for(j = 0; j < 2; j++) {// neighbour line choice 
1539       aneighbourindex = (!j) ? (i-1) : (i+1);
1540       if((aneighbourindex < 1) || (aneighbourindex > nblines)){
1541         continue;
1542       }
1543       //
1544       aLineTypeNeib=anArrayOfLineType[aneighbourindex];
1545       if(!aLineTypeNeib){
1546         continue;
1547       }
1548       //
1549       const ListOfInteger& aNeighbour = anArrayOfLines[aneighbourindex];
1550       Standard_Integer anIndex = (!j) ? aNeighbour.Last() : aNeighbour.First();
1551       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
1552       // check if need use derivative.begin .end [absence]
1553       //
1554       IntSurf_PntOn2S aNewP = aPoint;
1555       Standard_Integer surfit, parit;
1556       //
1557       for(surfit = 0; surfit < 2; ++surfit) {
1558
1559         Handle(GeomAdaptor_HSurface) aGASurface = (!surfit) ? theSurface1 : theSurface2;
1560         
1561         aGASurface->ChangeSurface().Surface()->Bounds(umin, umax, vmin, vmax);
1562         Standard_Real U=0., V=0.;
1563
1564         if(!surfit) {
1565           aNewP.ParametersOnS1(U, V);
1566         }
1567         else {
1568           aNewP.ParametersOnS2(U, V);
1569         }
1570         //
1571         Standard_Integer nbboundaries = 0;
1572         Standard_Integer bIsUBoundary = Standard_False; // use if nbboundaries == 1
1573         Standard_Integer bIsFirstBoundary = Standard_False; // use if nbboundaries == 1
1574         //
1575         for(parit = 0; parit < 2; parit++) {
1576           Standard_Boolean isperiodic = (!parit) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1577
1578           Standard_Real aResolution = (!parit) ? aGASurface->UResolution(aTol) : aGASurface->VResolution(aTol);
1579           Standard_Real alowerboundary = (!parit) ? umin : vmin;
1580           Standard_Real aupperboundary = (!parit) ? umax : vmax;
1581
1582           Standard_Real aParameter = (!parit) ? U : V;
1583           Standard_Boolean bIsOnFirstBoundary = Standard_True;
1584   
1585           if(!isperiodic) {
1586             if(IsPointOnBoundary(aParameter, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary)) {
1587               bIsUBoundary = (!parit);
1588               bIsFirstBoundary = bIsOnFirstBoundary;
1589               nbboundaries++;
1590             }
1591           }
1592           else {
1593             Standard_Real aPeriod     = (!parit) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1594             Standard_Real anoffset = 0.;
1595             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1596
1597             if(IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary)) {
1598               bIsUBoundary = (parit == 0);
1599               bIsFirstBoundary = bIsOnFirstBoundary;
1600               nbboundaries++;
1601             }
1602           }
1603         }
1604         //
1605         Standard_Boolean bComputeLineEnd = Standard_False;
1606         
1607         if(nbboundaries == 2) {
1608           bComputeLineEnd = Standard_True;
1609         }
1610         else if(nbboundaries == 1) {
1611           Standard_Boolean isperiodic = (bIsUBoundary) ? aGASurface->IsUPeriodic() : aGASurface->IsVPeriodic();
1612
1613           if(isperiodic) {
1614             Standard_Real alowerboundary = (bIsUBoundary) ? umin : vmin;
1615             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
1616             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
1617             Standard_Real aParameter = (bIsUBoundary) ? U : V;
1618             Standard_Real anoffset = 0.;
1619             Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
1620
1621             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
1622             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
1623             anotherPar += anoffset;
1624             Standard_Integer aneighbourpointindex = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1625             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex);
1626             Standard_Real nU1, nV1;
1627
1628             if(surfit == 0)
1629               aNeighbourPoint.ParametersOnS1(nU1, nV1);
1630             else
1631               aNeighbourPoint.ParametersOnS2(nU1, nV1);
1632             
1633             Standard_Real adist1 = (bIsUBoundary) ? fabs(nU1 - U) : fabs(nV1 - V);
1634             Standard_Real adist2 = (bIsUBoundary) ? fabs(nU1 - anotherPar) : fabs(nV1 - anotherPar);
1635             bComputeLineEnd = Standard_True;
1636             Standard_Boolean bCheckAngle1 = Standard_False;
1637             Standard_Boolean bCheckAngle2 = Standard_False;
1638             gp_Vec2d aNewVec;
1639             Standard_Real anewU = (bIsUBoundary) ? anotherPar : U;
1640             Standard_Real anewV = (bIsUBoundary) ? V : anotherPar;
1641             //
1642             if(((adist1 - adist2) > Precision::PConfusion()) && 
1643                (adist2 < (aPeriod / 4.))) {
1644               bCheckAngle1 = Standard_True;
1645               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
1646
1647               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
1648                 aNewP.SetValue((surfit == 0), anewU, anewV);
1649                 bCheckAngle1 = Standard_False;
1650               }
1651             }
1652             else if(adist1 < (aPeriod / 4.)) {
1653               bCheckAngle2 = Standard_True;
1654               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
1655
1656               if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
1657                 bCheckAngle2 = Standard_False;
1658               }
1659             }
1660             //
1661             if(bCheckAngle1 || bCheckAngle2) {
1662               // assume there are at least two points in line (see "if" above)
1663               Standard_Integer anindexother = aneighbourpointindex;
1664               
1665               while((anindexother <= aListOfIndex.Last()) && (anindexother >= aListOfIndex.First())) {
1666                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
1667                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
1668                 Standard_Real nU2, nV2;
1669                 
1670                 if(surfit == 0)
1671                   aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1672                 else
1673                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1674                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
1675
1676                 if(aVecOld.SquareMagnitude() <= (gp::Resolution() * gp::Resolution())) {
1677                   continue;
1678                 }
1679                 else {
1680                   Standard_Real anAngle = aNewVec.Angle(aVecOld);
1681
1682                   if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
1683
1684                     if(bCheckAngle1) {
1685                       Standard_Real U1, U2, V1, V2;
1686                       IntSurf_PntOn2S atmppoint = aNewP;
1687                       atmppoint.SetValue((surfit == 0), anewU, anewV);
1688                       atmppoint.Parameters(U1, V1, U2, V2);
1689                       gp_Pnt P1 = theSurface1->Value(U1, V1);
1690                       gp_Pnt P2 = theSurface2->Value(U2, V2);
1691                       gp_Pnt P0 = aPoint.Value();
1692
1693                       if(P0.IsEqual(P1, aTol) &&
1694                          P0.IsEqual(P2, aTol) &&
1695                          P1.IsEqual(P2, aTol)) {
1696                         bComputeLineEnd = Standard_False;
1697                         aNewP.SetValue((surfit == 0), anewU, anewV);
1698                       }
1699                     }
1700                     
1701                     if(bCheckAngle2) {
1702                       bComputeLineEnd = Standard_False;
1703                     }
1704                   }
1705                   break;
1706                 }
1707               } // end while(anindexother...)
1708             }
1709           }
1710         }
1711         //
1712         if(bComputeLineEnd) {
1713           Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
1714           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
1715           Standard_Real nU1, nV1;
1716
1717           if(surfit == 0)
1718             aNeighbourPoint.ParametersOnS1(nU1, nV1);
1719           else
1720             aNeighbourPoint.ParametersOnS2(nU1, nV1);
1721           gp_Pnt2d ap1(nU1, nV1);
1722           gp_Pnt2d ap2(nU1, nV1);
1723           Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
1724
1725           while((aneighbourpointindex2 <= aListOfIndex.Last()) && (aneighbourpointindex2 >= aListOfIndex.First())) {
1726             aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
1727             const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
1728             Standard_Real nU2, nV2;
1729
1730             if(surfit == 0)
1731               aPrevNeighbourPoint.ParametersOnS1(nU2, nV2);
1732             else
1733               aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
1734             ap2.SetX(nU2);
1735             ap2.SetY(nV2);
1736
1737             if(ap1.SquareDistance(ap2) > (gp::Resolution() * gp::Resolution())) {
1738               break;
1739             }
1740           }       
1741           gp_Pnt2d anewpoint;
1742           Standard_Boolean found = FindPoint(ap2, ap1, umin, umax, vmin, vmax, anewpoint);
1743
1744           if(found) {
1745             // check point
1746             Standard_Real aCriteria =aTolSum;// BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
1747             //GeomAPI_ProjectPointOnSurf& aProjector = (surfit == 0) ? aPrj2 : aPrj1;
1748             ProjectPointOnSurf& aProjector = (surfit == 0) ? aPrj2 : aPrj1;
1749             Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
1750
1751             gp_Pnt aP3d = aSurface->Value(anewpoint.X(), anewpoint.Y());
1752             aProjector.Perform(aP3d);
1753
1754             if(aProjector.IsDone()) {
1755               if(aProjector.LowerDistance() < aCriteria) {
1756                 Standard_Real foundU = U, foundV = V;
1757                 aProjector.LowerDistanceParameters(foundU, foundV);
1758
1759                 if(surfit == 0)
1760                   aNewP.SetValue(aP3d, anewpoint.X(), anewpoint.Y(), foundU, foundV);
1761                 else
1762                   aNewP.SetValue(aP3d, foundU, foundV, anewpoint.X(), anewpoint.Y());
1763               }
1764             }
1765           }
1766         }
1767       }
1768       aSeqOfPntOn2S->Add(aNewP);
1769       aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
1770     }
1771     anArrayOfLineEnds.SetValue(i, aListOfFLIndex);
1772   }
1773   // Correct wlines.end
1774
1775   // Split wlines.begin
1776   for(j = 1; j <= theLConstructor.NbParts(); j++) {
1777     Standard_Real fprm = 0., lprm = 0.;
1778     theLConstructor.Part(j, fprm, lprm);
1779     Standard_Integer ifprm = (Standard_Integer)fprm;
1780     Standard_Integer ilprm = (Standard_Integer)lprm;
1781     //
1782     Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1783     //
1784     for(i = 1; i <= nblines; i++) {
1785       if(anArrayOfLineType[i] != 0) {
1786         continue;
1787       }
1788       const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1789
1790       if(aListOfIndex.Extent() < 2) {
1791         continue;
1792       }
1793       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
1794       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
1795       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
1796
1797       if(!bhasfirstpoint && !aListOfFLIndex.IsEmpty()) {
1798         bhasfirstpoint = (i != 1);
1799       }
1800
1801       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
1802         bhaslastpoint = (i != nblines);
1803       }
1804       Standard_Boolean bIsFirstInside = ((ifprm >= aListOfIndex.First()) && (ifprm <= aListOfIndex.Last()));
1805       Standard_Boolean bIsLastInside =  ((ilprm >= aListOfIndex.First()) && (ilprm <= aListOfIndex.Last()));
1806
1807       if(!bIsFirstInside && !bIsLastInside) {
1808         if((ifprm < aListOfIndex.First()) && (ilprm > aListOfIndex.Last())) {
1809           // append whole line, and boundaries if neccesary
1810           if(bhasfirstpoint) {
1811             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
1812             aLineOn2S->Add(aP);
1813           }
1814           ListOfInteger::Iterator anIt(aListOfIndex);
1815
1816           for(; anIt.More(); anIt.Next()) {
1817             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1818             aLineOn2S->Add(aP);
1819           }
1820
1821           if(bhaslastpoint) {
1822             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
1823             aLineOn2S->Add(aP);
1824           }
1825
1826           // check end of split line (end is almost always)
1827           Standard_Integer aneighbour = i + 1;
1828           Standard_Boolean bIsEndOfLine = Standard_True;
1829
1830           if(aneighbour <= nblines) {
1831             const ListOfInteger& aListOfNeighbourIndex = anArrayOfLines[aneighbour];
1832
1833             if((anArrayOfLineType[aneighbour] != 0) &&
1834                (aListOfNeighbourIndex.IsEmpty())) {
1835               bIsEndOfLine = Standard_False;
1836             }
1837           }
1838
1839           if(bIsEndOfLine) {
1840             if(aLineOn2S->NbPoints() > 1) {
1841               Handle(IntPatch_WLine) aNewWLine = 
1842                 new IntPatch_WLine(aLineOn2S, Standard_False);
1843               theNewLines.Append(aNewWLine);
1844             }
1845             aLineOn2S = new IntSurf_LineOn2S();
1846           }
1847         }
1848         continue;
1849       }
1850       // end if(!bIsFirstInside && !bIsLastInside)
1851
1852       if(bIsFirstInside && bIsLastInside) {
1853         // append inside points between ifprm and ilprm
1854         ListOfInteger::Iterator anIt(aListOfIndex);
1855
1856         for(; anIt.More(); anIt.Next()) {
1857           if((anIt.Value() < ifprm) || (anIt.Value() > ilprm))
1858             continue;
1859           const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1860           aLineOn2S->Add(aP);
1861         }
1862       }
1863       else {
1864
1865         if(bIsFirstInside) {
1866           // append points from ifprm to last point + boundary point
1867           ListOfInteger::Iterator anIt(aListOfIndex);
1868
1869           for(; anIt.More(); anIt.Next()) {
1870             if(anIt.Value() < ifprm)
1871               continue;
1872             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1873             aLineOn2S->Add(aP);
1874           }
1875
1876           if(bhaslastpoint) {
1877             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
1878             aLineOn2S->Add(aP);
1879           }
1880           // check end of split line (end is almost always)
1881           Standard_Integer aneighbour = i + 1;
1882           Standard_Boolean bIsEndOfLine = Standard_True;
1883
1884           if(aneighbour <= nblines) {
1885             const ListOfInteger& aListOfNeighbourIndex = anArrayOfLines[aneighbour];
1886
1887             if((anArrayOfLineType[aneighbour] != 0) &&
1888                (aListOfNeighbourIndex.IsEmpty())) {
1889               bIsEndOfLine = Standard_False;
1890             }
1891           }
1892
1893           if(bIsEndOfLine) {
1894             if(aLineOn2S->NbPoints() > 1) {
1895               Handle(IntPatch_WLine) aNewWLine = 
1896                 new IntPatch_WLine(aLineOn2S, Standard_False);
1897               theNewLines.Append(aNewWLine);
1898             }
1899             aLineOn2S = new IntSurf_LineOn2S();
1900           }
1901         }
1902         // end if(bIsFirstInside)
1903
1904         if(bIsLastInside) {
1905           // append points from first boundary point to ilprm
1906           if(bhasfirstpoint) {
1907             const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
1908             aLineOn2S->Add(aP);
1909           }
1910           ListOfInteger::Iterator anIt(aListOfIndex);
1911
1912           for(; anIt.More(); anIt.Next()) {
1913             if(anIt.Value() > ilprm)
1914               continue;
1915             const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
1916             aLineOn2S->Add(aP);
1917           }
1918         }
1919         //end if(bIsLastInside)
1920       }
1921     }
1922
1923     if(aLineOn2S->NbPoints() > 1) {
1924       Handle(IntPatch_WLine) aNewWLine = 
1925         new IntPatch_WLine(aLineOn2S, Standard_False);
1926       theNewLines.Append(aNewWLine);
1927     }
1928   }
1929   // Split wlines.end
1930   //
1931   // cda002/I3
1932   Standard_Real fprm, lprm;
1933   Standard_Integer ifprm, ilprm, aNbPoints, aIndex;
1934   //
1935   aNbParts=theLConstructor.NbParts();
1936   //
1937   for(j = 1; j <= aNbParts; j++) {
1938     theLConstructor.Part(j, fprm, lprm);
1939     ifprm=(Standard_Integer)fprm;
1940     ilprm=(Standard_Integer)lprm;
1941     //
1942     if ((ilprm-ifprm)==1) {
1943       for(i = 1; i <= nblines; i++) {
1944         aLineType=anArrayOfLineType[i];
1945         if(aLineType) {
1946           continue;
1947         }
1948         //
1949         const ListOfInteger& aListOfIndex = anArrayOfLines[i];
1950         aNbPoints=aListOfIndex.Extent();
1951         if(aNbPoints==1) {
1952           aIndex=aListOfIndex.First();
1953           if (aIndex==ifprm || aIndex==ilprm) {
1954             Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
1955             const IntSurf_PntOn2S& aP1 = theWLine->Point(ifprm);
1956             const IntSurf_PntOn2S& aP2 = theWLine->Point(ilprm);
1957             aLineOn2S->Add(aP1);
1958             aLineOn2S->Add(aP2);
1959             Handle(IntPatch_WLine) aNewWLine = 
1960               new IntPatch_WLine(aLineOn2S, Standard_False);
1961             theNewLines.Append(aNewWLine);
1962           }
1963         }
1964       }
1965     }
1966   }
1967   //
1968   return Standard_True;
1969 }
1970
1971 //=======================================================================
1972 //function : AdjustPeriodic
1973 //purpose  : 
1974 //=======================================================================
1975 Standard_Real AdjustPeriodic(const Standard_Real theParameter,
1976                              const Standard_Real parmin,
1977                              const Standard_Real parmax,
1978                              const Standard_Real thePeriod,
1979                              Standard_Real&      theOffset)
1980 {
1981   Standard_Real aresult = theParameter;
1982   theOffset = 0.;
1983   while(aresult < parmin) {
1984     aresult += thePeriod;
1985     theOffset += thePeriod;
1986   }
1987   while(aresult > parmax) {
1988     aresult -= thePeriod;
1989     theOffset -= thePeriod;
1990   }
1991   return aresult;
1992 }
1993
1994 //=======================================================================
1995 //function : IsPointOnBoundary
1996 //purpose  : 
1997 //=======================================================================
1998 Standard_Boolean IsPointOnBoundary(const Standard_Real theParameter,
1999                                    const Standard_Real theFirstBoundary,
2000                                    const Standard_Real theSecondBoundary,
2001                                    const Standard_Real theResolution,
2002                                    Standard_Boolean&   IsOnFirstBoundary)
2003 {
2004   IsOnFirstBoundary = Standard_True;
2005   if(fabs(theParameter - theFirstBoundary) < theResolution)
2006     return Standard_True;
2007   if(fabs(theParameter - theSecondBoundary) < theResolution)
2008   {
2009     IsOnFirstBoundary = Standard_False;
2010     return Standard_True;
2011   }
2012   return Standard_False;
2013 }
2014 //=======================================================================
2015 //function : FindPoint
2016 //purpose  : 
2017 //=======================================================================
2018 Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
2019                            const gp_Pnt2d&     theLastPoint,
2020                            const Standard_Real theUmin,
2021                            const Standard_Real theUmax,
2022                            const Standard_Real theVmin,
2023                            const Standard_Real theVmax,
2024                            gp_Pnt2d&           theNewPoint)
2025 {
2026   gp_Vec2d aVec(theFirstPoint, theLastPoint);
2027   Standard_Integer i = 0, j = 0;
2028
2029   for(i = 0; i < 4; i++) {
2030     gp_Vec2d anOtherVec;
2031     gp_Vec2d anOtherVecNormal;
2032     gp_Pnt2d aprojpoint = theLastPoint;    
2033
2034     if((i % 2) == 0) {
2035       anOtherVec.SetX(0.);
2036       anOtherVec.SetY(1.);
2037       anOtherVecNormal.SetX(1.);
2038       anOtherVecNormal.SetY(0.);
2039
2040       if(i < 2)
2041         aprojpoint.SetX(theUmin);
2042       else
2043         aprojpoint.SetX(theUmax);
2044     }
2045     else {
2046       anOtherVec.SetX(1.);
2047       anOtherVec.SetY(0.);
2048       anOtherVecNormal.SetX(0.);
2049       anOtherVecNormal.SetY(1.);
2050
2051       if(i < 2)
2052         aprojpoint.SetY(theVmin);
2053       else
2054         aprojpoint.SetY(theVmax);
2055     }
2056     gp_Vec2d anormvec = aVec;
2057     anormvec.Normalize();
2058     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
2059
2060     if(fabs(adot1) < Precision::Angular())
2061       continue;
2062     Standard_Real adist = 0.;
2063
2064     if((i % 2) == 0) {
2065       adist = (i < 2) ? fabs(theLastPoint.X() - theUmin) : fabs(theLastPoint.X() - theUmax);
2066     }
2067     else {
2068       adist = (i < 2) ? fabs(theLastPoint.Y() - theVmin) : fabs(theLastPoint.Y() - theVmax);
2069     }
2070     Standard_Real anoffset = adist * anOtherVec.Dot(anormvec) / adot1;
2071
2072     for(j = 0; j < 2; j++) {
2073       anoffset = (j == 0) ? anoffset : -anoffset;
2074       gp_Pnt2d acurpoint(aprojpoint.XY() + (anOtherVec.XY()*anoffset));
2075       gp_Vec2d acurvec(theLastPoint, acurpoint);
2076
2077       //
2078       Standard_Real aDotX, anAngleX, aPC;
2079       //
2080       aDotX=aVec.Dot(acurvec);
2081       anAngleX=aVec.Angle(acurvec);
2082       aPC=Precision::PConfusion();
2083       //
2084       if(aDotX > 0. && fabs(anAngleX) < aPC) {
2085       //
2086         if((i % 2) == 0) {
2087           if((acurpoint.Y() >= theVmin) &&
2088              (acurpoint.Y() <= theVmax)) {
2089             theNewPoint = acurpoint;
2090             return Standard_True;
2091           }
2092         }
2093         else {
2094           if((acurpoint.X() >= theUmin) &&
2095              (acurpoint.X() <= theUmax)) {
2096             theNewPoint = acurpoint;
2097             return Standard_True;
2098           }
2099         }
2100       }
2101     }
2102   }
2103   return Standard_False;
2104 }
2105
2106 /*
2107 static
2108   void DumpWLine (const Handle(IntPatch_WLine)& theWLine);
2109 //=======================================================================
2110 //function : DumpWLine
2111 //purpose  : 
2112 //=======================================================================
2113   void DumpWLine (const Handle(IntPatch_WLine)& theWLine)
2114 {
2115   Standard_Integer i, nbp;
2116   Standard_Real u1,v1,u2,v2;
2117   //
2118   nbp=theWLine->NbPnts();
2119   printf(" nbp=%d\n", nbp);
2120   for(i=1;i<=nbp;i++) {
2121     const IntSurf_PntOn2S& aPoint = theWLine->Point(i);
2122     const gp_Pnt& aP=aPoint.Value();
2123     aPoint.Parameters(u1,v1,u2,v2);
2124     printf("point i%d %lf %lf %lf [ %lf %lf ] [ %lf %lf ]\n", 
2125            i, aP.X(), aP.Y(), aP.Z(), u1, v1, u2, v2);
2126   }
2127 }
2128 */