1 // Created on: 2000-11-23
2 // Created by: Michael KLOKOV
3 // Copyright (c) 2000-2014 OPEN CASCADE SAS
4 //
5 // This file is part of Open CASCADE Technology software library.
6 //
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
12 //
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
15 
16 #include <IntTools_FaceFace.hxx>
17 
18 #include <BRepAdaptor_Surface.hxx>
19 #include <BRepTools.hxx>
20 #include <BRep_Tool.hxx>
21 #include <ElCLib.hxx>
22 #include <ElSLib.hxx>
23 #include <Geom2dAdaptor_Curve.hxx>
24 #include <Geom2dInt_GInter.hxx>
25 #include <Geom2d_BSplineCurve.hxx>
26 #include <Geom2d_Line.hxx>
27 #include <Geom2d_TrimmedCurve.hxx>
28 #include <GeomAPI_ProjectPointOnSurf.hxx>
29 #include <GeomAdaptor_Surface.hxx>
30 #include <GeomInt_IntSS.hxx>
31 #include <GeomInt_WLApprox.hxx>
32 #include <GeomLib_Check2dBSplineCurve.hxx>
33 #include <GeomLib_CheckBSplineCurve.hxx>
34 #include <Geom_BSplineCurve.hxx>
35 #include <Geom_Circle.hxx>
36 #include <Geom_Ellipse.hxx>
37 #include <Geom_Hyperbola.hxx>
38 #include <Geom_Line.hxx>
39 #include <Geom_OffsetSurface.hxx>
40 #include <Geom_Parabola.hxx>
41 #include <Geom_RectangularTrimmedSurface.hxx>
42 #include <Geom_TrimmedCurve.hxx>
43 #include <IntAna_QuadQuadGeo.hxx>
44 #include <IntPatch_GLine.hxx>
45 #include <IntPatch_RLine.hxx>
46 #include <IntPatch_WLine.hxx>
47 #include <IntRes2d_Domain.hxx>
48 #include <IntSurf_Quadric.hxx>
49 #include <IntTools_Context.hxx>
50 #include <IntTools_Tools.hxx>
51 #include <IntTools_TopolTool.hxx>
52 #include <IntTools_WLineTool.hxx>
53 #include <ProjLib_Plane.hxx>
54 #include <TopExp_Explorer.hxx>
55 #include <TopoDS.hxx>
56 #include <TopoDS_Edge.hxx>
57 #include <gp_Elips.hxx>
58 
59 static
60   void Parameters(const Handle(GeomAdaptor_Surface)&,
61                   const Handle(GeomAdaptor_Surface)&,
62                   const gp_Pnt&,
63                   Standard_Real&,
64                   Standard_Real&,
65                   Standard_Real&,
66                   Standard_Real&);
67 
68 static
69   void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
70                                 const Standard_Real theTolerance,
71                                 Standard_Real&      theumin,
72                                 Standard_Real&      theumax,
73                                 Standard_Real&      thevmin,
74                                 Standard_Real&      thevmax);
75 
76 static
77   Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
78                                           const Handle(Geom_Curve)& theCurve,
79                                           const TopoDS_Face&        theFace1,
80                                           const TopoDS_Face&        theFace2,
81                                           const Standard_Real       theOtherParameter,
82                                           const Standard_Boolean    bIncreasePar,
83                                           const Standard_Real       theTol,
84                                           Standard_Real&            theNewParameter,
85                                           const Handle(IntTools_Context)& );
86 
87 static
88   Standard_Boolean IsCurveValid(const Handle(Geom2d_Curve)& thePCurve);
89 
90 static
91   Standard_Boolean  ApproxWithPCurves(const gp_Cylinder& theCyl,
92                                       const gp_Sphere& theSph);
93 
94 static void  PerformPlanes(const Handle(GeomAdaptor_Surface)& theS1,
95                            const Handle(GeomAdaptor_Surface)& theS2,
96                            const Standard_Real TolF1,
97                            const Standard_Real TolF2,
98                            const Standard_Real TolAng,
99                            const Standard_Real TolTang,
100                            const Standard_Boolean theApprox1,
101                            const Standard_Boolean theApprox2,
102                            IntTools_SequenceOfCurves& theSeqOfCurve,
103                            Standard_Boolean& theTangentFaces);
104 
105 static Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_Surface)& theS,
106                                       const gp_Lin2d& theLin2d,
107                                       const Standard_Real theTol,
108                                       Standard_Real& theP1,
109                                       Standard_Real& theP2);
110 //
111 static
112   void ApproxParameters(const Handle(GeomAdaptor_Surface)& aHS1,
113                         const Handle(GeomAdaptor_Surface)& aHS2,
114                         Standard_Integer& iDegMin,
115                         Standard_Integer& iNbIter,
116                         Standard_Integer& iDegMax);
117 
118 static
119   void Tolerances(const Handle(GeomAdaptor_Surface)& aHS1,
120                   const Handle(GeomAdaptor_Surface)& aHS2,
121                   Standard_Real& aTolTang);
122 
123 static
124   Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
125                              const GeomAbs_SurfaceType aType2);
126 static
127   Standard_Integer IndexType(const GeomAbs_SurfaceType aType);
128 
129 //
130 static
131   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
132                                const TopoDS_Face& aFace,
133                                const Handle(IntTools_Context)& theCtx);
134 
135 static
136   Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
137                             const Standard_Real aT,
138                             GeomAPI_ProjectPointOnSurf& theProjPS);
139 
140 static
141   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
142                                 const Standard_Real theFirst,
143                                 const Standard_Real theLast,
144                                 GeomAPI_ProjectPointOnSurf& theProjPS,
145                                 const Standard_Real theEps);
146 
147 static
148   Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
149                                 const Standard_Real theFirst,
150                                 const Standard_Real theLast,
151                                 const TopoDS_Face& theFace,
152                                 const Handle(IntTools_Context)& theContext);
153 
154 static
155   void CorrectPlaneBoundaries(Standard_Real& aUmin,
156                               Standard_Real& aUmax,
157                               Standard_Real& aVmin,
158                               Standard_Real& aVmax);
159 
160 //=======================================================================
161 //function :
162 //purpose  :
163 //=======================================================================
IntTools_FaceFace()164 IntTools_FaceFace::IntTools_FaceFace()
165 {
166   myIsDone=Standard_False;
167   myTangentFaces=Standard_False;
168   //
169   myHS1 = new GeomAdaptor_Surface ();
170   myHS2 = new GeomAdaptor_Surface ();
171   myTolF1 = 0.;
172   myTolF2 = 0.;
173   myTol = 0.;
174   myFuzzyValue = Precision::Confusion();
175   SetParameters(Standard_True, Standard_True, Standard_True, 1.e-07);
176 }
177 //=======================================================================
178 //function : SetContext
179 //purpose  :
180 //=======================================================================
SetContext(const Handle (IntTools_Context)& aContext)181 void IntTools_FaceFace::SetContext(const Handle(IntTools_Context)& aContext)
182 {
183   myContext=aContext;
184 }
185 //=======================================================================
186 //function : Context
187 //purpose  :
188 //=======================================================================
Handle(IntTools_Context)189 const Handle(IntTools_Context)& IntTools_FaceFace::Context()const
190 {
191   return myContext;
192 }
193 //=======================================================================
194 //function : Face1
195 //purpose  :
196 //=======================================================================
Face1() const197 const TopoDS_Face& IntTools_FaceFace::Face1() const
198 {
199   return myFace1;
200 }
201 //=======================================================================
202 //function : Face2
203 //purpose  :
204 //=======================================================================
Face2() const205 const TopoDS_Face& IntTools_FaceFace::Face2() const
206 {
207   return myFace2;
208 }
209 //=======================================================================
210 //function : TangentFaces
211 //purpose  :
212 //=======================================================================
TangentFaces() const213 Standard_Boolean IntTools_FaceFace::TangentFaces() const
214 {
215   return myTangentFaces;
216 }
217 //=======================================================================
218 //function : Points
219 //purpose  :
220 //=======================================================================
Points() const221 const IntTools_SequenceOfPntOn2Faces& IntTools_FaceFace::Points() const
222 {
223   return myPnts;
224 }
225 //=======================================================================
226 //function : IsDone
227 //purpose  :
228 //=======================================================================
IsDone() const229 Standard_Boolean IntTools_FaceFace::IsDone() const
230 {
231   return myIsDone;
232 }
233 //=======================================================================
234 //function : Lines
235 //purpose  : return lines of intersection
236 //=======================================================================
Lines() const237 const IntTools_SequenceOfCurves& IntTools_FaceFace::Lines() const
238 {
239   StdFail_NotDone_Raise_if
240     (!myIsDone,
241      "IntTools_FaceFace::Lines() => myIntersector NOT DONE");
242   return mySeqOfCurve;
243 }
244 // =======================================================================
245 // function: SetParameters
246 //
247 // =======================================================================
SetParameters(const Standard_Boolean ToApproxC3d,const Standard_Boolean ToApproxC2dOnS1,const Standard_Boolean ToApproxC2dOnS2,const Standard_Real ApproximationTolerance)248 void IntTools_FaceFace::SetParameters(const Standard_Boolean ToApproxC3d,
249                                       const Standard_Boolean ToApproxC2dOnS1,
250                                       const Standard_Boolean ToApproxC2dOnS2,
251                                       const Standard_Real ApproximationTolerance)
252 {
253   myApprox = ToApproxC3d;
254   myApprox1 = ToApproxC2dOnS1;
255   myApprox2 = ToApproxC2dOnS2;
256   myTolApprox = ApproximationTolerance;
257 }
258 //=======================================================================
259 //function : SetFuzzyValue
260 //purpose  :
261 //=======================================================================
SetFuzzyValue(const Standard_Real theFuzz)262 void IntTools_FaceFace::SetFuzzyValue(const Standard_Real theFuzz)
263 {
264   myFuzzyValue = Max(theFuzz, Precision::Confusion());
265 }
266 //=======================================================================
267 //function : FuzzyValue
268 //purpose  :
269 //=======================================================================
FuzzyValue() const270 Standard_Real IntTools_FaceFace::FuzzyValue() const
271 {
272   return myFuzzyValue;
273 }
274 
275 //=======================================================================
276 //function : SetList
277 //purpose  :
278 //=======================================================================
SetList(IntSurf_ListOfPntOn2S & aListOfPnts)279 void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
280 {
281   myListOfPnts = aListOfPnts;
282 }
283 
284 
isTreatAnalityc(const BRepAdaptor_Surface & theBAS1,const BRepAdaptor_Surface & theBAS2,const Standard_Real theTol)285 static Standard_Boolean isTreatAnalityc(const BRepAdaptor_Surface& theBAS1,
286                                         const BRepAdaptor_Surface& theBAS2,
287                                         const Standard_Real theTol)
288 {
289   const Standard_Real Tolang = 1.e-8;
290   Standard_Real aHigh = 0.0;
291 
292   const GeomAbs_SurfaceType aType1=theBAS1.GetType();
293   const GeomAbs_SurfaceType aType2=theBAS2.GetType();
294 
295   gp_Pln aS1;
296   gp_Cylinder aS2;
297   if(aType1 == GeomAbs_Plane)
298   {
299     aS1=theBAS1.Plane();
300   }
301   else if(aType2 == GeomAbs_Plane)
302   {
303     aS1=theBAS2.Plane();
304   }
305   else
306   {
307     return Standard_True;
308   }
309 
310   if(aType1 == GeomAbs_Cylinder)
311   {
312     aS2=theBAS1.Cylinder();
313     const Standard_Real VMin = theBAS1.FirstVParameter();
314     const Standard_Real VMax = theBAS1.LastVParameter();
315 
316     if( Precision::IsNegativeInfinite(VMin) ||
317         Precision::IsPositiveInfinite(VMax))
318           return Standard_True;
319     else
320       aHigh = VMax - VMin;
321   }
322   else if(aType2 == GeomAbs_Cylinder)
323   {
324     aS2=theBAS2.Cylinder();
325 
326     const Standard_Real VMin = theBAS2.FirstVParameter();
327     const Standard_Real VMax = theBAS2.LastVParameter();
328 
329     if( Precision::IsNegativeInfinite(VMin) ||
330         Precision::IsPositiveInfinite(VMax))
331           return Standard_True;
332     else
333       aHigh = VMax - VMin;
334   }
335   else
336   {
337     return Standard_True;
338   }
339 
340   IntAna_QuadQuadGeo inter;
341   inter.Perform(aS1,aS2,Tolang,theTol, aHigh);
342   if(inter.TypeInter() == IntAna_Ellipse)
343   {
344     const gp_Elips anEl = inter.Ellipse(1);
345     const Standard_Real aMajorR = anEl.MajorRadius();
346     const Standard_Real aMinorR = anEl.MinorRadius();
347 
348     return (aMajorR < 100000.0 * aMinorR);
349   }
350   else
351   {
352     return inter.IsDone();
353   }
354 }
355 //=======================================================================
356 //function : Perform
357 //purpose  : intersect surfaces of the faces
358 //=======================================================================
Perform(const TopoDS_Face & aF1,const TopoDS_Face & aF2,const Standard_Boolean theToRunParallel)359 void IntTools_FaceFace::Perform (const TopoDS_Face& aF1,
360                                  const TopoDS_Face& aF2,
361                                  const Standard_Boolean theToRunParallel)
362 {
363   if (myContext.IsNull()) {
364     myContext=new IntTools_Context;
365   }
366 
367   mySeqOfCurve.Clear();
368   myIsDone = Standard_False;
369   myNbrestr=0;//?
370 
371   myFace1=aF1;
372   myFace2=aF2;
373 
374   const BRepAdaptor_Surface& aBAS1 = myContext->SurfaceAdaptor(myFace1);
375   const BRepAdaptor_Surface& aBAS2 = myContext->SurfaceAdaptor(myFace2);
376   GeomAbs_SurfaceType aType1=aBAS1.GetType();
377   GeomAbs_SurfaceType aType2=aBAS2.GetType();
378 
379   const Standard_Boolean bReverse=SortTypes(aType1, aType2);
380   if (bReverse)
381   {
382     myFace1=aF2;
383     myFace2=aF1;
384     aType1=aBAS2.GetType();
385     aType2=aBAS1.GetType();
386 
387     if (myListOfPnts.Extent())
388     {
389       Standard_Real aU1,aV1,aU2,aV2;
390       IntSurf_ListIteratorOfListOfPntOn2S aItP2S;
391       //
392       aItP2S.Initialize(myListOfPnts);
393       for (; aItP2S.More(); aItP2S.Next())
394       {
395         IntSurf_PntOn2S& aP2S=aItP2S.Value();
396         aP2S.Parameters(aU1,aV1,aU2,aV2);
397         aP2S.SetValue(aU2,aV2,aU1,aV1);
398       }
399     }
400     //
401     Standard_Boolean anAproxTmp = myApprox1;
402     myApprox1 = myApprox2;
403     myApprox2 = anAproxTmp;
404   }
405 
406 
407   const Handle(Geom_Surface) S1=BRep_Tool::Surface(myFace1);
408   const Handle(Geom_Surface) S2=BRep_Tool::Surface(myFace2);
409 
410   Standard_Real aFuzz = myFuzzyValue / 2.;
411   myTolF1 = BRep_Tool::Tolerance(myFace1) + aFuzz;
412   myTolF2 = BRep_Tool::Tolerance(myFace2) + aFuzz;
413   myTol = myTolF1 + myTolF2;
414 
415   Standard_Real TolArc = myTol;
416   Standard_Real TolTang = TolArc;
417 
418   const Standard_Boolean isFace1Quad = (aType1 == GeomAbs_Cylinder ||
419                                         aType1 == GeomAbs_Cone ||
420                                         aType1 == GeomAbs_Torus);
421 
422   const Standard_Boolean isFace2Quad = (aType2 == GeomAbs_Cylinder ||
423                                         aType2 == GeomAbs_Cone ||
424                                         aType2 == GeomAbs_Torus);
425 
426   if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane)  {
427     Standard_Real umin, umax, vmin, vmax;
428     //
429     myContext->UVBounds(myFace1, umin, umax, vmin, vmax);
430     myHS1->Load(S1, umin, umax, vmin, vmax);
431     //
432     myContext->UVBounds(myFace2, umin, umax, vmin, vmax);
433     myHS2->Load(S2, umin, umax, vmin, vmax);
434     //
435     Standard_Real TolAng = 1.e-8;
436     //
437     PerformPlanes(myHS1, myHS2,
438                   myTolF1, myTolF2, TolAng, TolTang,
439                   myApprox1, myApprox2,
440                   mySeqOfCurve, myTangentFaces);
441     //
442     myIsDone = Standard_True;
443     //
444     if (!myTangentFaces) {
445       const Standard_Integer NbLinPP = mySeqOfCurve.Length();
446       if (NbLinPP && bReverse) {
447         Handle(Geom2d_Curve) aC2D1, aC2D2;
448         const Standard_Integer aNbLin = mySeqOfCurve.Length();
449         for (Standard_Integer i = 1; i <= aNbLin; ++i) {
450           IntTools_Curve& aIC = mySeqOfCurve(i);
451           aC2D1 = aIC.FirstCurve2d();
452           aC2D2 = aIC.SecondCurve2d();
453           aIC.SetFirstCurve2d(aC2D2);
454           aIC.SetSecondCurve2d(aC2D1);
455         }
456       }
457     }
458     return;
459   }//if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane){
460 
461   if ((aType1==GeomAbs_Plane) && isFace2Quad)
462   {
463     Standard_Real umin, umax, vmin, vmax;
464     // F1
465     myContext->UVBounds(myFace1, umin, umax, vmin, vmax);
466     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
467     myHS1->Load(S1, umin, umax, vmin, vmax);
468     // F2
469     myContext->UVBounds(myFace2, umin, umax, vmin, vmax);
470     CorrectSurfaceBoundaries(myFace2, myTol * 2., umin, umax, vmin, vmax);
471     myHS2->Load(S2, umin, umax, vmin, vmax);
472   }
473   else if ((aType2==GeomAbs_Plane) && isFace1Quad)
474   {
475     Standard_Real umin, umax, vmin, vmax;
476     //F1
477     myContext->UVBounds(myFace1, umin, umax, vmin, vmax);
478     CorrectSurfaceBoundaries(myFace1, myTol * 2., umin, umax, vmin, vmax);
479     myHS1->Load(S1, umin, umax, vmin, vmax);
480     // F2
481     myContext->UVBounds(myFace2, umin, umax, vmin, vmax);
482     CorrectPlaneBoundaries(umin, umax, vmin, vmax);
483     myHS2->Load(S2, umin, umax, vmin, vmax);
484   }
485   else
486   {
487     Standard_Real umin, umax, vmin, vmax;
488     myContext->UVBounds(myFace1, umin, umax, vmin, vmax);
489     CorrectSurfaceBoundaries(myFace1, myTol * 2., umin, umax, vmin, vmax);
490     myHS1->Load(S1, umin, umax, vmin, vmax);
491     myContext->UVBounds(myFace2, umin, umax, vmin, vmax);
492     CorrectSurfaceBoundaries(myFace2, myTol * 2., umin, umax, vmin, vmax);
493     myHS2->Load(S2, umin, umax, vmin, vmax);
494   }
495 
496   const Handle(IntTools_TopolTool) dom1 = new IntTools_TopolTool(myHS1);
497   const Handle(IntTools_TopolTool) dom2 = new IntTools_TopolTool(myHS2);
498 
499   myLConstruct.Load(dom1, dom2, myHS1, myHS2);
500 
501 
502   Tolerances(myHS1, myHS2, TolTang);
503 
504   {
505     const Standard_Real UVMaxStep = 0.001;
506     const Standard_Real Deflection = 0.1;
507     myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection);
508   }
509 
510   if((aType1 != GeomAbs_BSplineSurface) &&
511       (aType1 != GeomAbs_BezierSurface)  &&
512      (aType1 != GeomAbs_OtherSurface)  &&
513      (aType2 != GeomAbs_BSplineSurface) &&
514       (aType2 != GeomAbs_BezierSurface)  &&
515      (aType2 != GeomAbs_OtherSurface))
516   {
517     if ((aType1 == GeomAbs_Torus) ||
518         (aType2 == GeomAbs_Torus))
519     {
520       myListOfPnts.Clear();
521     }
522   }
523 
524 #ifdef INTTOOLS_FACEFACE_DEBUG
525     if(!myListOfPnts.IsEmpty()) {
526       char aBuff[10000];
527 
528       Sprintf(aBuff,"bopcurves <face1 face2> -2d");
529       IntSurf_ListIteratorOfListOfPntOn2S IterLOP1(myListOfPnts);
530       for(;IterLOP1.More(); IterLOP1.Next())
531       {
532         const IntSurf_PntOn2S& aPt = IterLOP1.Value();
533         Standard_Real u1, v1, u2, v2;
534         aPt.Parameters(u1, v1, u2, v2);
535 
536         Sprintf(aBuff, "%s -p %+10.20f %+10.20f %+10.20f %+10.20f", aBuff, u1, v1, u2, v2);
537       }
538 
539       std::cout << aBuff << std::endl;
540     }
541 #endif
542 
543   const Standard_Boolean isGeomInt = isTreatAnalityc(aBAS1, aBAS2, myTol);
544   if (aF1.IsSame(aF2))
545     myIntersector.Perform(myHS1, dom1, TolArc, TolTang);
546   else
547     myIntersector.Perform(myHS1, dom1, myHS2, dom2, TolArc, TolTang,
548                           myListOfPnts, isGeomInt);
549 
550   myIsDone = myIntersector.IsDone();
551 
552   if (myIsDone)
553   {
554     myTangentFaces=myIntersector.TangentFaces();
555     if (myTangentFaces) {
556       return;
557     }
558     //
559     const Standard_Integer aNbLinIntersector = myIntersector.NbLines();
560     for (Standard_Integer i=1; i <= aNbLinIntersector; ++i) {
561       MakeCurve(i, dom1, dom2, TolArc);
562     }
563     //
564     ComputeTolReached3d (theToRunParallel);
565     //
566     if (bReverse) {
567       Handle(Geom2d_Curve) aC2D1, aC2D2;
568       //
569       const Standard_Integer aNbLinSeqOfCurve =mySeqOfCurve.Length();
570       for (Standard_Integer i=1; i<=aNbLinSeqOfCurve; ++i)
571       {
572         IntTools_Curve& aIC=mySeqOfCurve(i);
573         aC2D1=aIC.FirstCurve2d();
574         aC2D2=aIC.SecondCurve2d();
575         aIC.SetFirstCurve2d(aC2D2);
576         aIC.SetSecondCurve2d(aC2D1);
577       }
578     }
579 
580     // Points
581     Standard_Boolean bValid2D1, bValid2D2;
582     Standard_Real U1,V1,U2,V2;
583     IntTools_PntOnFace aPntOnF1, aPntOnF2;
584     IntTools_PntOn2Faces aPntOn2Faces;
585     //
586     const Standard_Integer aNbPnts = myIntersector.NbPnts();
587     for (Standard_Integer i=1; i <= aNbPnts; ++i)
588     {
589       const IntSurf_PntOn2S& aISPnt=myIntersector.Point(i).PntOn2S();
590       const gp_Pnt& aPnt=aISPnt.Value();
591       aISPnt.Parameters(U1,V1,U2,V2);
592       //
593       // check the validity of the intersection point for the faces
594       bValid2D1 = myContext->IsPointInOnFace(myFace1, gp_Pnt2d(U1, V1));
595       if (!bValid2D1) {
596         continue;
597       }
598       //
599       bValid2D2 = myContext->IsPointInOnFace(myFace2, gp_Pnt2d(U2, V2));
600       if (!bValid2D2) {
601         continue;
602       }
603       //
604       // add the intersection point
605       aPntOnF1.Init(myFace1, aPnt, U1, V1);
606       aPntOnF2.Init(myFace2, aPnt, U2, V2);
607       //
608       if (!bReverse)
609       {
610         aPntOn2Faces.SetP1(aPntOnF1);
611         aPntOn2Faces.SetP2(aPntOnF2);
612       }
613       else
614       {
615         aPntOn2Faces.SetP2(aPntOnF1);
616         aPntOn2Faces.SetP1(aPntOnF2);
617       }
618 
619       myPnts.Append(aPntOn2Faces);
620     }
621   }
622 }
623 
624 //=======================================================================
625 //function :ComputeTolReached3d
626 //purpose  :
627 //=======================================================================
ComputeTolReached3d(const Standard_Boolean theToRunParallel)628 void IntTools_FaceFace::ComputeTolReached3d (const Standard_Boolean theToRunParallel)
629 {
630   Standard_Integer i, j, aNbLin = mySeqOfCurve.Length();
631   if (!aNbLin) {
632     return;
633   }
634   //
635   // Minimal tangential tolerance for the curve
636   Standard_Real aTolFMax = Max(myTolF1, myTolF2);
637   //
638   const Handle(Geom_Surface)& aS1 = myHS1->Surface();
639   const Handle(Geom_Surface)& aS2 = myHS2->Surface();
640   //
641   for (i = 1; i <= aNbLin; ++i)
642   {
643     IntTools_Curve& aIC = mySeqOfCurve(i);
644     const Handle(Geom_Curve)& aC3D = aIC.Curve();
645     if (aC3D.IsNull())
646     {
647       continue;
648     }
649     //
650     Standard_Real aTolC = aIC.Tolerance();
651     Standard_Real aFirst = aC3D->FirstParameter();
652     Standard_Real aLast  = aC3D->LastParameter();
653     //
654     // Compute the tolerance for the curve
655     const Handle(Geom2d_Curve)& aC2D1 = aIC.FirstCurve2d();
656     const Handle(Geom2d_Curve)& aC2D2 = aIC.SecondCurve2d();
657     //
658     for (j = 0; j < 2; ++j)
659     {
660       const Handle(Geom2d_Curve)& aC2D = !j ? aC2D1 : aC2D2;
661       if (!aC2D.IsNull())
662       {
663         // Look for the maximal deviation between 3D and 2D curves
664         Standard_Real aD, aT;
665         const Handle(Geom_Surface)& aS = !j ? aS1 : aS2;
666         if (IntTools_Tools::ComputeTolerance (aC3D, aC2D, aS, aFirst, aLast, aD, aT, Precision::PConfusion(), theToRunParallel))
667         {
668           if (aD > aTolC)
669           {
670             aTolC = aD;
671           }
672         }
673       }
674       else
675       {
676         // Look for the maximal deviation between 3D curve and surface
677         const TopoDS_Face& aF = !j ? myFace1 : myFace2;
678         Standard_Real aD = FindMaxDistance(aC3D, aFirst, aLast, aF, myContext);
679         if (aD > aTolC)
680         {
681           aTolC = aD;
682         }
683       }
684     }
685     // Set the valid tolerance for the curve
686     aIC.SetTolerance(aTolC);
687     //
688     // Set the tangential tolerance for the curve.
689     // Note, that, currently, computation of the tangential tolerance is
690     // implemented for the Plane/Plane case only.
691     // Thus, set the tangential tolerance equal to maximal tolerance of faces.
692     if (aIC.TangentialTolerance() < aTolFMax) {
693       aIC.SetTangentialTolerance(aTolFMax);
694     }
695   }
696 }
697 
698 //=======================================================================
699 //function : MakeCurve
700 //purpose  :
701 //=======================================================================
MakeCurve(const Standard_Integer Index,const Handle (Adaptor3d_TopolTool)& dom1,const Handle (Adaptor3d_TopolTool)& dom2,const Standard_Real theToler)702 void IntTools_FaceFace::MakeCurve(const Standard_Integer Index,
703                                   const Handle(Adaptor3d_TopolTool)& dom1,
704                                   const Handle(Adaptor3d_TopolTool)& dom2,
705                                   const Standard_Real theToler)
706 {
707   Standard_Boolean bDone, rejectSurface, reApprox, bAvoidLineConstructor;
708   Standard_Boolean ok, bPCurvesOk;
709   Standard_Integer i, j, aNbParts;
710   Standard_Real fprm, lprm;
711   Standard_Real Tolpc;
712   Handle(IntPatch_Line) L;
713   IntPatch_IType typl;
714   Handle(Geom_Curve) newc;
715   //
716   const Standard_Real TOLCHECK    = 1.e-7;
717   const Standard_Real TOLANGCHECK = 1.e-6;
718   //
719   rejectSurface = Standard_False;
720   reApprox = Standard_False;
721   //
722   bPCurvesOk = Standard_True;
723 
724  reapprox:;
725 
726   Tolpc = myTolApprox;
727   bAvoidLineConstructor = Standard_False;
728   L = myIntersector.Line(Index);
729   typl = L->ArcType();
730   //
731   if(typl==IntPatch_Walking) {
732     Handle(IntPatch_WLine) aWLine (Handle(IntPatch_WLine)::DownCast(L));
733     if(aWLine.IsNull()) {
734       return;
735     }
736     L = aWLine;
737 
738     Standard_Integer nbp = aWLine->NbPnts();
739     const IntSurf_PntOn2S& p1 = aWLine->Point(1);
740     const IntSurf_PntOn2S& p2 = aWLine->Point(nbp);
741 
742     const gp_Pnt& P1 = p1.Value();
743     const gp_Pnt& P2 = p2.Value();
744 
745     if(P1.SquareDistance(P2) < 1.e-14) {
746       bAvoidLineConstructor = Standard_False;
747     }
748   }
749 
750   typl=L->ArcType();
751 
752   if(typl == IntPatch_Restriction)
753     bAvoidLineConstructor = Standard_True;
754 
755   //
756   // Line Constructor
757   if(!bAvoidLineConstructor) {
758     myLConstruct.Perform(L);
759     //
760     bDone=myLConstruct.IsDone();
761     if(!bDone)
762     {
763       return;
764     }
765 
766     if(typl != IntPatch_Restriction)
767     {
768       aNbParts=myLConstruct.NbParts();
769       if (aNbParts <= 0)
770       {
771         return;
772       }
773     }
774   }
775   // Do the Curve
776 
777 
778   switch (typl) {
779   //########################################
780   // Line, Parabola, Hyperbola
781   //########################################
782   case IntPatch_Lin:
783   case IntPatch_Parabola:
784   case IntPatch_Hyperbola: {
785     if (typl == IntPatch_Lin) {
786       newc =
787         new Geom_Line (Handle(IntPatch_GLine)::DownCast(L)->Line());
788     }
789 
790     else if (typl == IntPatch_Parabola) {
791       newc =
792         new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
793     }
794 
795     else if (typl == IntPatch_Hyperbola) {
796       newc =
797         new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
798     }
799     //
800     aNbParts=myLConstruct.NbParts();
801     for (i=1; i<=aNbParts; i++) {
802       Standard_Boolean bFNIt, bLPIt;
803       //
804       myLConstruct.Part(i, fprm, lprm);
805         //
806       bFNIt=Precision::IsNegativeInfinite(fprm);
807       bLPIt=Precision::IsPositiveInfinite(lprm);
808       //
809       if (!bFNIt && !bLPIt) {
810         //
811         IntTools_Curve aCurve;
812         //
813         Handle(Geom_TrimmedCurve) aCT3D=new Geom_TrimmedCurve(newc, fprm, lprm);
814         aCurve.SetCurve(aCT3D);
815         if (typl == IntPatch_Parabola) {
816           Standard_Real aTolC = IntTools_Tools::CurveTolerance(aCT3D, myTol);
817           aCurve.SetTolerance(aTolC);
818         }
819         //
820         if(myApprox1) {
821           Handle (Geom2d_Curve) C2d;
822           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
823                 myHS1->Surface(), newc, C2d);
824 
825           if (C2d.IsNull())
826             continue;
827 
828           aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d, fprm, lprm));
829         }
830         //
831         if(myApprox2) {
832           Handle (Geom2d_Curve) C2d;
833           GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
834                     myHS2->Surface(), newc, C2d);
835 
836           if (C2d.IsNull())
837             continue;
838 
839           aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d, fprm, lprm));
840         }
841         //
842         mySeqOfCurve.Append(aCurve);
843       } //if (!bFNIt && !bLPIt) {
844       else {
845         //  on regarde si on garde
846         //
847         Standard_Real aTestPrm, dT=100.;
848         //
849         aTestPrm=0.;
850         if (bFNIt && !bLPIt) {
851           aTestPrm=lprm-dT;
852         }
853         else if (!bFNIt && bLPIt) {
854           aTestPrm=fprm+dT;
855         }
856         else {
857           // i.e, if (bFNIt && bLPIt)
858           aTestPrm=IntTools_Tools::IntermediatePoint(-dT, dT);
859         }
860         //
861         gp_Pnt ptref(newc->Value(aTestPrm));
862         //
863         GeomAbs_SurfaceType typS1 = myHS1->GetType();
864         GeomAbs_SurfaceType typS2 = myHS2->GetType();
865         if( typS1 == GeomAbs_SurfaceOfExtrusion ||
866             typS1 == GeomAbs_OffsetSurface ||
867             typS1 == GeomAbs_SurfaceOfRevolution ||
868             typS2 == GeomAbs_SurfaceOfExtrusion ||
869             typS2 == GeomAbs_OffsetSurface ||
870             typS2 == GeomAbs_SurfaceOfRevolution) {
871           Handle(Geom2d_BSplineCurve) H1;
872           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
873           continue;
874         }
875 
876         Standard_Real u1, v1, u2, v2, Tol;
877 
878         Tol = Precision::Confusion();
879         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
880         ok = (dom1->Classify(gp_Pnt2d(u1, v1), Tol) != TopAbs_OUT);
881         if(ok) {
882           ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
883         }
884         if (ok) {
885           Handle(Geom2d_BSplineCurve) H1;
886           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
887         }
888       }
889     }// for (i=1; i<=aNbParts; i++) {
890   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
891     break;
892 
893   //########################################
894   // Circle and Ellipse
895   //########################################
896   case IntPatch_Circle:
897   case IntPatch_Ellipse: {
898 
899     if (typl == IntPatch_Circle) {
900       newc = new Geom_Circle
901         (Handle(IntPatch_GLine)::DownCast(L)->Circle());
902     }
903     else { //IntPatch_Ellipse
904       newc = new Geom_Ellipse
905         (Handle(IntPatch_GLine)::DownCast(L)->Ellipse());
906     }
907     //
908     aNbParts=myLConstruct.NbParts();
909     //
910     Standard_Real aPeriod, aNul;
911     TColStd_SequenceOfReal aSeqFprm,  aSeqLprm;
912 
913     aNul=0.;
914     aPeriod=M_PI+M_PI;
915 
916     for (i=1; i<=aNbParts; i++) {
917       myLConstruct.Part(i, fprm, lprm);
918 
919       if (fprm < aNul && lprm > aNul) {
920         // interval that goes through 0. is divided on two intervals;
921         while (fprm<aNul || fprm>aPeriod) fprm=fprm+aPeriod;
922         while (lprm<aNul || lprm>aPeriod) lprm=lprm+aPeriod;
923         //
924         if((aPeriod - fprm) > Tolpc) {
925           aSeqFprm.Append(fprm);
926           aSeqLprm.Append(aPeriod);
927         }
928         else {
929           gp_Pnt P1 = newc->Value(fprm);
930           gp_Pnt P2 = newc->Value(aPeriod);
931 
932           if(P1.Distance(P2) > myTol) {
933             Standard_Real anewpar = fprm;
934 
935             if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2,
936                                       lprm, Standard_False, myTol, anewpar, myContext)) {
937               fprm = anewpar;
938             }
939             aSeqFprm.Append(fprm);
940             aSeqLprm.Append(aPeriod);
941           }
942         }
943 
944         //
945         if((lprm - aNul) > Tolpc) {
946           aSeqFprm.Append(aNul);
947           aSeqLprm.Append(lprm);
948         }
949         else {
950           gp_Pnt P1 = newc->Value(aNul);
951           gp_Pnt P2 = newc->Value(lprm);
952 
953           if(P1.Distance(P2) > myTol) {
954             Standard_Real anewpar = lprm;
955 
956             if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2,
957                                       fprm, Standard_True, myTol, anewpar, myContext)) {
958               lprm = anewpar;
959             }
960             aSeqFprm.Append(aNul);
961             aSeqLprm.Append(lprm);
962           }
963         }
964       }
965       else {
966         // usual interval
967         aSeqFprm.Append(fprm);
968         aSeqLprm.Append(lprm);
969       }
970     }
971     //
972     aNbParts=aSeqFprm.Length();
973     for (i=1; i<=aNbParts; i++) {
974       fprm=aSeqFprm(i);
975       lprm=aSeqLprm(i);
976       //
977       Standard_Real aRealEpsilon=RealEpsilon();
978       if (Abs(fprm) > aRealEpsilon || Abs(lprm-2.*M_PI) > aRealEpsilon) {
979         //==============================================
980         ////
981         IntTools_Curve aCurve;
982         Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
983         aCurve.SetCurve(aTC3D);
984         fprm=aTC3D->FirstParameter();
985         lprm=aTC3D->LastParameter ();
986         ////
987         if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {////
988           if(myApprox1) {
989             Handle (Geom2d_Curve) C2d;
990             GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
991                         myHS1->Surface(), newc, C2d);
992             aCurve.SetFirstCurve2d(C2d);
993           }
994 
995           if(myApprox2) {
996             Handle (Geom2d_Curve) C2d;
997             GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
998                     myHS2->Surface(),newc,C2d);
999             aCurve.SetSecondCurve2d(C2d);
1000           }
1001         }
1002         //
1003         mySeqOfCurve.Append(aCurve);
1004           //==============================================
1005       } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
1006 
1007       else {
1008         //  on regarde si on garde
1009         //
1010         if (aNbParts==1) {
1011 //           if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
1012           if (Abs(fprm) <= aRealEpsilon && Abs(lprm-2.*M_PI) <= aRealEpsilon) {
1013             IntTools_Curve aCurve;
1014             Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
1015             aCurve.SetCurve(aTC3D);
1016             fprm=aTC3D->FirstParameter();
1017             lprm=aTC3D->LastParameter ();
1018 
1019             if(myApprox1) {
1020               Handle (Geom2d_Curve) C2d;
1021               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1022                     myHS1->Surface(),newc,C2d);
1023               aCurve.SetFirstCurve2d(C2d);
1024             }
1025 
1026             if(myApprox2) {
1027               Handle (Geom2d_Curve) C2d;
1028               GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
1029                     myHS2->Surface(),newc,C2d);
1030               aCurve.SetSecondCurve2d(C2d);
1031             }
1032             //
1033             mySeqOfCurve.Append(aCurve);
1034             break;
1035           }
1036         }
1037         //
1038         Standard_Real aTwoPIdiv17, u1, v1, u2, v2, Tol;
1039 
1040         aTwoPIdiv17=2.*M_PI/17.;
1041 
1042         for (j=0; j<=17; j++) {
1043           gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
1044           Tol = Precision::Confusion();
1045 
1046           Parameters(myHS1, myHS2, ptref, u1, v1, u2, v2);
1047           ok = (dom1->Classify(gp_Pnt2d(u1,v1),Tol) != TopAbs_OUT);
1048           if(ok) {
1049             ok = (dom2->Classify(gp_Pnt2d(u2,v2),Tol) != TopAbs_OUT);
1050           }
1051           if (ok) {
1052             IntTools_Curve aCurve;
1053             aCurve.SetCurve(newc);
1054             //==============================================
1055             if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {
1056 
1057               if(myApprox1) {
1058                 Handle (Geom2d_Curve) C2d;
1059                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1060                         myHS1->Surface(), newc, C2d);
1061                 aCurve.SetFirstCurve2d(C2d);
1062               }
1063 
1064               if(myApprox2) {
1065                 Handle (Geom2d_Curve) C2d;
1066                 GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
1067                         myHS2->Surface(), newc, C2d);
1068                 aCurve.SetSecondCurve2d(C2d);
1069               }
1070             }//  end of if (typl == IntPatch_Circle || typl == IntPatch_Ellipse)
1071             //==============================================
1072             //
1073             mySeqOfCurve.Append(aCurve);
1074             break;
1075 
1076             }//  end of if (ok) {
1077           }//  end of for (Standard_Integer j=0; j<=17; j++)
1078         }//  end of else { on regarde si on garde
1079       }// for (i=1; i<=myLConstruct.NbParts(); i++)
1080     }// IntPatch_Circle: IntPatch_Ellipse:
1081     break;
1082 
1083   case IntPatch_Analytic:
1084     //This case was processed earlier (in IntPatch_Intersection)
1085     break;
1086 
1087   case IntPatch_Walking:{
1088     Handle(IntPatch_WLine) WL =
1089       Handle(IntPatch_WLine)::DownCast(L);
1090 
1091 #ifdef INTTOOLS_FACEFACE_DEBUG
1092     WL->Dump(0);
1093 #endif
1094 
1095     //
1096     Standard_Integer ifprm, ilprm;
1097     //
1098     if (!myApprox) {
1099       aNbParts = 1;
1100       if(!bAvoidLineConstructor){
1101         aNbParts=myLConstruct.NbParts();
1102       }
1103       for (i=1; i<=aNbParts; ++i) {
1104         Handle(Geom2d_BSplineCurve) H1, H2;
1105         Handle(Geom_Curve) aBSp;
1106         //
1107         if(bAvoidLineConstructor) {
1108           ifprm = 1;
1109           ilprm = WL->NbPnts();
1110         }
1111         else {
1112           myLConstruct.Part(i, fprm, lprm);
1113           ifprm=(Standard_Integer)fprm;
1114           ilprm=(Standard_Integer)lprm;
1115         }
1116         //
1117         if(myApprox1) {
1118           H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1119         }
1120         //
1121         if(myApprox2) {
1122           H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1123         }
1124         //
1125         aBSp=GeomInt_IntSS::MakeBSpline(WL, ifprm, ilprm);
1126         IntTools_Curve aIC(aBSp, H1, H2);
1127         mySeqOfCurve.Append(aIC);
1128       }// for (i=1; i<=aNbParts; ++i) {
1129     }// if (!myApprox) {
1130     //
1131     else { // X
1132       Standard_Boolean bIsDecomposited;
1133       Standard_Integer nbiter, aNbSeqOfL;
1134       Standard_Real tol2d, aTolApproxImp;
1135       IntPatch_SequenceOfLine aSeqOfL;
1136       GeomInt_WLApprox theapp3d;
1137       Approx_ParametrizationType aParType = Approx_ChordLength;
1138       //
1139       Standard_Boolean anApprox1 = myApprox1;
1140       Standard_Boolean anApprox2 = myApprox2;
1141       //
1142       aTolApproxImp=1.e-5;
1143       tol2d = myTolApprox;
1144 
1145       GeomAbs_SurfaceType typs1, typs2;
1146       typs1 = myHS1->GetType();
1147       typs2 = myHS2->GetType();
1148       Standard_Boolean anWithPC = Standard_True;
1149 
1150       if(typs1 == GeomAbs_Cylinder && typs2 == GeomAbs_Sphere) {
1151         anWithPC =
1152           ApproxWithPCurves(myHS1->Cylinder(), myHS2->Sphere());
1153       }
1154       else if (typs1 == GeomAbs_Sphere && typs2 == GeomAbs_Cylinder) {
1155         anWithPC =
1156           ApproxWithPCurves(myHS2->Cylinder(), myHS1->Sphere());
1157       }
1158       //
1159       if(!anWithPC) {
1160         myTolApprox = aTolApproxImp;//1.e-5;
1161         anApprox1 = Standard_False;
1162         anApprox2 = Standard_False;
1163         //
1164         tol2d = myTolApprox;
1165       }
1166 
1167       if(myHS1 == myHS2) {
1168         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30, Standard_False, aParType);
1169         rejectSurface = Standard_True;
1170       }
1171       else {
1172         if(reApprox && !rejectSurface)
1173           theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30, Standard_False, aParType);
1174         else {
1175           Standard_Integer iDegMax, iDegMin, iNbIter;
1176           //
1177           ApproxParameters(myHS1, myHS2, iDegMin, iDegMax, iNbIter);
1178           theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax,
1179                                                   iNbIter, 30, Standard_True, aParType);
1180         }
1181       }
1182       //
1183       Standard_Real aReachedTol = Precision::Confusion();
1184       bIsDecomposited = IntTools_WLineTool::
1185         DecompositionOfWLine(WL,
1186                              myHS1,
1187                              myHS2,
1188                              myFace1,
1189                              myFace2,
1190                              myLConstruct,
1191                              bAvoidLineConstructor,
1192                              myTol,
1193                              aSeqOfL,
1194                              aReachedTol,
1195                              myContext);
1196       //
1197       aNbSeqOfL=aSeqOfL.Length();
1198       //
1199       Standard_Real aTolC = 0.;
1200       if (bIsDecomposited) {
1201         nbiter=aNbSeqOfL;
1202         aTolC = aReachedTol;
1203       }
1204       else {
1205         nbiter=1;
1206         aNbParts=1;
1207         if (!bAvoidLineConstructor) {
1208           aNbParts=myLConstruct.NbParts();
1209           nbiter=aNbParts;
1210         }
1211       }
1212       //
1213       for(i = 1; i <= nbiter; ++i) {
1214         if(bIsDecomposited) {
1215           WL = Handle(IntPatch_WLine)::DownCast(aSeqOfL.Value(i));
1216           ifprm = 1;
1217           ilprm = WL->NbPnts();
1218         }
1219         else {
1220           if(bAvoidLineConstructor) {
1221             ifprm = 1;
1222             ilprm = WL->NbPnts();
1223           }
1224           else {
1225             myLConstruct.Part(i, fprm, lprm);
1226             ifprm = (Standard_Integer)fprm;
1227             ilprm = (Standard_Integer)lprm;
1228           }
1229         }
1230         //-- lbr :
1231         //-- Si une des surfaces est un plan , on approxime en 2d
1232         //-- sur cette surface et on remonte les points 2d en 3d.
1233         if(typs1 == GeomAbs_Plane) {
1234           theapp3d.Perform(myHS1, myHS2, WL, Standard_False,Standard_True, myApprox2,ifprm,ilprm);
1235         }
1236         else if(typs2 == GeomAbs_Plane) {
1237           theapp3d.Perform(myHS1,myHS2,WL,Standard_False,myApprox1,Standard_True,ifprm,ilprm);
1238         }
1239         else {
1240           //
1241           if (myHS1 != myHS2){
1242             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
1243                 (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
1244 
1245               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30,
1246                                                                 Standard_True, aParType);
1247 
1248               Standard_Boolean bUseSurfaces;
1249               bUseSurfaces = IntTools_WLineTool::NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
1250               if (bUseSurfaces) {
1251                 // ######
1252                 rejectSurface = Standard_True;
1253                 // ######
1254                 theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, 30,
1255                                                                 Standard_False, aParType);
1256               }
1257             }
1258           }
1259           //
1260           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,anApprox1,anApprox2,ifprm,ilprm);
1261         }
1262           //
1263         if (!theapp3d.IsDone()) {
1264           Handle(Geom2d_BSplineCurve) H1;
1265           Handle(Geom2d_BSplineCurve) H2;
1266           //
1267           Handle(Geom_Curve) aBSp=GeomInt_IntSS::MakeBSpline(WL,ifprm, ilprm);
1268           //
1269           if(myApprox1) {
1270             H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1271           }
1272           //
1273           if(myApprox2) {
1274             H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1275           }
1276           //
1277           IntTools_Curve aIC(aBSp, H1, H2);
1278           mySeqOfCurve.Append(aIC);
1279         }
1280         else {
1281           if (typs1 == GeomAbs_Plane || typs2 == GeomAbs_Plane) {
1282             //
1283             if (typs1 == GeomAbs_Torus || typs2 == GeomAbs_Torus) {
1284               if (aTolC < 1.e-6) {
1285                 aTolC = 1.e-6;
1286               }
1287             }
1288           }
1289           //
1290           Standard_Integer aNbMultiCurves, nbpoles;
1291           aNbMultiCurves=theapp3d.NbMultiCurves();
1292           for (j=1; j<=aNbMultiCurves; j++) {
1293             if(typs1 == GeomAbs_Plane) {
1294               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1295               nbpoles = mbspc.NbPoles();
1296 
1297               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1298               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1299 
1300               mbspc.Curve(1,tpoles2d);
1301               const gp_Pln&  Pln = myHS1->Plane();
1302               //
1303               Standard_Integer ik;
1304               for(ik = 1; ik<= nbpoles; ik++) {
1305                 tpoles.SetValue(ik,
1306                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1307                                               tpoles2d.Value(ik).Y(),
1308                                               Pln));
1309               }
1310               //
1311               Handle(Geom_BSplineCurve) BS =
1312                 new Geom_BSplineCurve(tpoles,
1313                                       mbspc.Knots(),
1314                                       mbspc.Multiplicities(),
1315                                       mbspc.Degree());
1316               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1317               Check.FixTangent(Standard_True, Standard_True);
1318               //
1319               IntTools_Curve aCurve;
1320               aCurve.SetCurve(BS);
1321 
1322               if(myApprox1) {
1323                 Handle(Geom2d_BSplineCurve) BS1 =
1324                   new Geom2d_BSplineCurve(tpoles2d,
1325                                           mbspc.Knots(),
1326                                           mbspc.Multiplicities(),
1327                                           mbspc.Degree());
1328                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1329                 Check1.FixTangent(Standard_True,Standard_True);
1330                 //
1331                 // ############################################
1332                 if(!rejectSurface && !reApprox) {
1333                   Standard_Boolean isValid = IsCurveValid(BS1);
1334                   if(!isValid) {
1335                     reApprox = Standard_True;
1336                     goto reapprox;
1337                   }
1338                 }
1339                 // ############################################
1340                 aCurve.SetFirstCurve2d(BS1);
1341               }
1342 
1343               if(myApprox2) {
1344                 mbspc.Curve(2, tpoles2d);
1345 
1346                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
1347                                                                           mbspc.Knots(),
1348                                                                           mbspc.Multiplicities(),
1349                                                                           mbspc.Degree());
1350                 GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1351                 newCheck.FixTangent(Standard_True,Standard_True);
1352 
1353                 // ###########################################
1354                 if(!rejectSurface && !reApprox) {
1355                   Standard_Boolean isValid = IsCurveValid(BS2);
1356                   if(!isValid) {
1357                     reApprox = Standard_True;
1358                     goto reapprox;
1359                   }
1360                 }
1361                 // ###########################################
1362                 //
1363                 aCurve.SetSecondCurve2d(BS2);
1364               }
1365               //
1366               aCurve.SetTolerance(aTolC);
1367               //
1368               mySeqOfCurve.Append(aCurve);
1369 
1370             }//if(typs1 == GeomAbs_Plane) {
1371 
1372             else if(typs2 == GeomAbs_Plane)
1373             {
1374               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1375               nbpoles = mbspc.NbPoles();
1376 
1377               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1378               TColgp_Array1OfPnt   tpoles(1,nbpoles);
1379               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
1380               const gp_Pln&  Pln = myHS2->Plane();
1381               //
1382               Standard_Integer ik;
1383               for(ik = 1; ik<= nbpoles; ik++) {
1384                 tpoles.SetValue(ik,
1385                                 ElSLib::Value(tpoles2d.Value(ik).X(),
1386                                               tpoles2d.Value(ik).Y(),
1387                                               Pln));
1388 
1389               }
1390               //
1391               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
1392                                                                  mbspc.Knots(),
1393                                                                  mbspc.Multiplicities(),
1394                                                                  mbspc.Degree());
1395               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1396               Check.FixTangent(Standard_True,Standard_True);
1397               //
1398               IntTools_Curve aCurve;
1399               aCurve.SetCurve(BS);
1400               aCurve.SetTolerance(aTolC);
1401 
1402               if(myApprox2) {
1403                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
1404                                                                         mbspc.Knots(),
1405                                                                         mbspc.Multiplicities(),
1406                                                                         mbspc.Degree());
1407                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
1408                 Check1.FixTangent(Standard_True,Standard_True);
1409                 //
1410                 // ###########################################
1411                 if(!rejectSurface && !reApprox) {
1412                   Standard_Boolean isValid = IsCurveValid(BS1);
1413                   if(!isValid) {
1414                     reApprox = Standard_True;
1415                     goto reapprox;
1416                   }
1417                 }
1418                 // ###########################################
1419                 bPCurvesOk = CheckPCurve(BS1, myFace2, myContext);
1420                 aCurve.SetSecondCurve2d(BS1);
1421               }
1422 
1423               if(myApprox1) {
1424                 mbspc.Curve(1,tpoles2d);
1425                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
1426                                                                         mbspc.Knots(),
1427                                                                         mbspc.Multiplicities(),
1428                                                                         mbspc.Degree());
1429                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
1430                 Check2.FixTangent(Standard_True,Standard_True);
1431                 //
1432                 // ###########################################
1433                 if(!rejectSurface && !reApprox) {
1434                   Standard_Boolean isValid = IsCurveValid(BS2);
1435                   if(!isValid) {
1436                     reApprox = Standard_True;
1437                     goto reapprox;
1438                   }
1439                 }
1440                 // ###########################################
1441                 bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1, myContext);
1442                 aCurve.SetFirstCurve2d(BS2);
1443               }
1444               //
1445               //if points of the pcurves are out of the faces bounds
1446               //create 3d and 2d curves without approximation
1447               if (!bPCurvesOk) {
1448                 Handle(Geom2d_BSplineCurve) H1, H2;
1449                 bPCurvesOk = Standard_True;
1450                 //
1451                 Handle(Geom_Curve) aBSp=GeomInt_IntSS::MakeBSpline(WL,ifprm, ilprm);
1452 
1453                 if(myApprox1) {
1454                   H1 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
1455                   bPCurvesOk = CheckPCurve(H1, myFace1, myContext);
1456                 }
1457 
1458                 if(myApprox2) {
1459                   H2 = GeomInt_IntSS::MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
1460                   bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2, myContext);
1461                 }
1462                 //
1463                 //if pcurves created without approximation are out of the
1464                 //faces bounds, use approximated 3d and 2d curves
1465                 if (bPCurvesOk) {
1466                   IntTools_Curve aIC(aBSp, H1, H2, aTolC);
1467                   mySeqOfCurve.Append(aIC);
1468                 } else {
1469                   mySeqOfCurve.Append(aCurve);
1470                 }
1471               } else {
1472                 mySeqOfCurve.Append(aCurve);
1473               }
1474 
1475             }// else if(typs2 == GeomAbs_Plane)
1476             //
1477             else { //typs2 != GeomAbs_Plane && typs1 != GeomAbs_Plane
1478               Standard_Boolean bIsValid1, bIsValid2;
1479               Handle(Geom_BSplineCurve) BS;
1480               IntTools_Curve aCurve;
1481               //
1482               bIsValid1=Standard_True;
1483               bIsValid2=Standard_True;
1484               //
1485               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
1486               nbpoles = mbspc.NbPoles();
1487               TColgp_Array1OfPnt tpoles(1,nbpoles);
1488               mbspc.Curve(1,tpoles);
1489               BS=new Geom_BSplineCurve(tpoles,
1490                                                                  mbspc.Knots(),
1491                                                                  mbspc.Multiplicities(),
1492                                                                  mbspc.Degree());
1493               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
1494               Check.FixTangent(Standard_True,Standard_True);
1495               //
1496               aCurve.SetCurve(BS);
1497               aCurve.SetTolerance(aTolC);
1498               //
1499               if(myApprox1) {
1500                 if(anApprox1) {
1501                   Handle(Geom2d_BSplineCurve) BS1;
1502                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1503                   mbspc.Curve(2,tpoles2d);
1504                   //
1505                   BS1=new Geom2d_BSplineCurve(tpoles2d,
1506                                               mbspc.Knots(),
1507                                               mbspc.Multiplicities(),
1508                                               mbspc.Degree());
1509                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
1510                   newCheck.FixTangent(Standard_True,Standard_True);
1511                   //
1512                   if (!reApprox) {
1513                     bIsValid1=CheckPCurve(BS1, myFace1, myContext);
1514                   }
1515                   //
1516                   aCurve.SetFirstCurve2d(BS1);
1517                 }
1518                 else {
1519                   Handle(Geom2d_BSplineCurve) BS1;
1520                   fprm = BS->FirstParameter();
1521                   lprm = BS->LastParameter();
1522 
1523                   Handle(Geom2d_Curve) C2d;
1524                   Standard_Real aTol = myTolApprox;
1525                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
1526                             myHS1->Surface(), BS, C2d);
1527                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
1528                   aCurve.SetFirstCurve2d(BS1);
1529                 }
1530               } // if(myApprox1) {
1531                 //
1532               if(myApprox2) {
1533                 if(anApprox2) {
1534                   Handle(Geom2d_BSplineCurve) BS2;
1535                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
1536                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
1537                   BS2=new Geom2d_BSplineCurve(tpoles2d,
1538                                                                         mbspc.Knots(),
1539                                                                         mbspc.Multiplicities(),
1540                                                                         mbspc.Degree());
1541                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
1542                   newCheck.FixTangent(Standard_True,Standard_True);
1543                 //
1544                   if (!reApprox) {
1545                     bIsValid2=CheckPCurve(BS2, myFace2, myContext);
1546                   }
1547                   aCurve.SetSecondCurve2d(BS2);
1548                 }
1549                 else {
1550                   Handle(Geom2d_BSplineCurve) BS2;
1551                   fprm = BS->FirstParameter();
1552                   lprm = BS->LastParameter();
1553 
1554                   Handle(Geom2d_Curve) C2d;
1555                   Standard_Real aTol = myTolApprox;
1556                   GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
1557                             myHS2->Surface(), BS, C2d);
1558                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
1559                   aCurve.SetSecondCurve2d(BS2);
1560                 }
1561               } //if(myApprox2) {
1562               if (!bIsValid1 || !bIsValid2) {
1563                 myTolApprox=aTolApproxImp;//1.e-5;
1564                 tol2d = myTolApprox;
1565                 reApprox = Standard_True;
1566                 goto reapprox;
1567               }
1568                 //
1569               mySeqOfCurve.Append(aCurve);
1570             }
1571           }
1572         }
1573       }
1574     }// else { // X
1575   }// case IntPatch_Walking:{
1576     break;
1577 
1578   case IntPatch_Restriction:
1579     {
1580       Handle(IntPatch_RLine) RL =
1581         Handle(IntPatch_RLine)::DownCast(L);
1582 
1583 #ifdef INTTOOLS_FACEFACE_DEBUG
1584     RL->Dump(0);
1585 #endif
1586 
1587       Handle(Geom_Curve) aC3d;
1588       Handle(Geom2d_Curve) aC2d1, aC2d2;
1589       Standard_Real aTolReached;
1590       GeomInt_IntSS::TreatRLine(RL, myHS1, myHS2, aC3d,
1591                                   aC2d1, aC2d2, aTolReached);
1592 
1593       if(aC3d.IsNull())
1594         break;
1595 
1596       Bnd_Box2d aBox1, aBox2;
1597 
1598       const Standard_Real aU1f = myHS1->FirstUParameter(),
1599                           aV1f = myHS1->FirstVParameter(),
1600                           aU1l = myHS1->LastUParameter(),
1601                           aV1l = myHS1->LastVParameter();
1602       const Standard_Real aU2f = myHS2->FirstUParameter(),
1603                           aV2f = myHS2->FirstVParameter(),
1604                           aU2l = myHS2->LastUParameter(),
1605                           aV2l = myHS2->LastVParameter();
1606 
1607       aBox1.Add(gp_Pnt2d(aU1f, aV1f));
1608       aBox1.Add(gp_Pnt2d(aU1l, aV1l));
1609       aBox2.Add(gp_Pnt2d(aU2f, aV2f));
1610       aBox2.Add(gp_Pnt2d(aU2l, aV2l));
1611 
1612       GeomInt_VectorOfReal anArrayOfParameters;
1613 
1614       //We consider here that the intersection line is same-parameter-line
1615       anArrayOfParameters.Append(aC3d->FirstParameter());
1616       anArrayOfParameters.Append(aC3d->LastParameter());
1617 
1618       GeomInt_IntSS::
1619         TrimILineOnSurfBoundaries(aC2d1, aC2d2, aBox1, aBox2, anArrayOfParameters);
1620 
1621       //Intersect with true boundaries. After that, enlarge bounding-boxes in order to
1622       //correct definition, if point on curve is inscribed in the box.
1623       aBox1.Enlarge(theToler);
1624       aBox2.Enlarge(theToler);
1625 
1626       const Standard_Integer aNbIntersSolutionsm1 = anArrayOfParameters.Length() - 1;
1627 
1628       //Trim RLine found.
1629       for(Standard_Integer anInd = 0; anInd < aNbIntersSolutionsm1; anInd++)
1630       {
1631         Standard_Real &aParF = anArrayOfParameters(anInd),
1632                       &aParL = anArrayOfParameters(anInd+1);
1633 
1634         if((aParL - aParF) <= Precision::PConfusion())
1635         {
1636           //In order to more precise extending to the boundaries of source curves.
1637           if(anInd < aNbIntersSolutionsm1-1)
1638             aParL = aParF;
1639 
1640           continue;
1641         }
1642 
1643         const Standard_Real aPar = 0.5*(aParF + aParL);
1644         gp_Pnt2d aPt;
1645 
1646         Handle(Geom2d_Curve) aCurv2d1, aCurv2d2;
1647         if(!aC2d1.IsNull())
1648         {
1649           aC2d1->D0(aPar, aPt);
1650 
1651           if(aBox1.IsOut(aPt))
1652             continue;
1653 
1654           if(myApprox1)
1655             aCurv2d1 = new Geom2d_TrimmedCurve(aC2d1, aParF, aParL);
1656         }
1657 
1658         if(!aC2d2.IsNull())
1659         {
1660           aC2d2->D0(aPar, aPt);
1661 
1662           if(aBox2.IsOut(aPt))
1663             continue;
1664 
1665           if(myApprox2)
1666             aCurv2d2 = new Geom2d_TrimmedCurve(aC2d2, aParF, aParL);
1667         }
1668 
1669         Handle(Geom_Curve) aCurv3d = new Geom_TrimmedCurve(aC3d, aParF, aParL);
1670 
1671         IntTools_Curve aIC(aCurv3d, aCurv2d1, aCurv2d2);
1672         mySeqOfCurve.Append(aIC);
1673       }
1674     }
1675     break;
1676   default:
1677     break;
1678 
1679   }
1680 }
1681 
1682 //=======================================================================
1683 //function : Parameters
1684 //purpose  :
1685 //=======================================================================
Parameters(const Handle (GeomAdaptor_Surface)& HS1,const Handle (GeomAdaptor_Surface)& HS2,const gp_Pnt & Ptref,Standard_Real & U1,Standard_Real & V1,Standard_Real & U2,Standard_Real & V2)1686  void Parameters(const Handle(GeomAdaptor_Surface)& HS1,
1687                  const Handle(GeomAdaptor_Surface)& HS2,
1688                  const gp_Pnt& Ptref,
1689                  Standard_Real& U1,
1690                  Standard_Real& V1,
1691                  Standard_Real& U2,
1692                  Standard_Real& V2)
1693 {
1694 
1695   IntSurf_Quadric quad1,quad2;
1696   GeomAbs_SurfaceType typs = HS1->GetType();
1697 
1698   switch (typs) {
1699   case GeomAbs_Plane:
1700     quad1.SetValue(HS1->Plane());
1701     break;
1702   case GeomAbs_Cylinder:
1703     quad1.SetValue(HS1->Cylinder());
1704     break;
1705   case GeomAbs_Cone:
1706     quad1.SetValue(HS1->Cone());
1707     break;
1708   case GeomAbs_Sphere:
1709     quad1.SetValue(HS1->Sphere());
1710     break;
1711   case GeomAbs_Torus:
1712     quad1.SetValue(HS1->Torus());
1713     break;
1714   default:
1715     throw Standard_ConstructionError("GeomInt_IntSS::MakeCurve");
1716   }
1717 
1718   typs = HS2->GetType();
1719   switch (typs) {
1720   case GeomAbs_Plane:
1721     quad2.SetValue(HS2->Plane());
1722     break;
1723   case GeomAbs_Cylinder:
1724     quad2.SetValue(HS2->Cylinder());
1725     break;
1726   case GeomAbs_Cone:
1727     quad2.SetValue(HS2->Cone());
1728     break;
1729   case GeomAbs_Sphere:
1730     quad2.SetValue(HS2->Sphere());
1731     break;
1732   case GeomAbs_Torus:
1733     quad2.SetValue(HS2->Torus());
1734     break;
1735   default:
1736     throw Standard_ConstructionError("GeomInt_IntSS::MakeCurve");
1737   }
1738 
1739   quad1.Parameters(Ptref,U1,V1);
1740   quad2.Parameters(Ptref,U2,V2);
1741 }
1742 
1743 //=======================================================================
1744 //function : MakeBSpline
1745 //purpose  :
1746 //=======================================================================
MakeBSpline(const Handle (IntPatch_WLine)& WL,const Standard_Integer ideb,const Standard_Integer ifin)1747 Handle(Geom_Curve) MakeBSpline  (const Handle(IntPatch_WLine)& WL,
1748                                  const Standard_Integer ideb,
1749                                  const Standard_Integer ifin)
1750 {
1751   Standard_Integer i,nbpnt = ifin-ideb+1;
1752   TColgp_Array1OfPnt poles(1,nbpnt);
1753   TColStd_Array1OfReal knots(1,nbpnt);
1754   TColStd_Array1OfInteger mults(1,nbpnt);
1755   Standard_Integer ipidebm1;
1756   for(i=1,ipidebm1=i+ideb-1; i<=nbpnt;ipidebm1++, i++) {
1757     poles(i) = WL->Point(ipidebm1).Value();
1758     mults(i) = 1;
1759     knots(i) = i-1;
1760   }
1761   mults(1) = mults(nbpnt) = 2;
1762   return
1763     new Geom_BSplineCurve(poles,knots,mults,1);
1764 }
1765 
1766 //=======================================================================
1767 //function : PrepareLines3D
1768 //purpose  :
1769 //=======================================================================
PrepareLines3D(const Standard_Boolean bToSplit)1770   void IntTools_FaceFace::PrepareLines3D(const Standard_Boolean bToSplit)
1771 {
1772   Standard_Integer i, aNbCurves;
1773   GeomAbs_SurfaceType aType1, aType2;
1774   IntTools_SequenceOfCurves aNewCvs;
1775   //
1776   // 1. Treatment closed  curves
1777   aNbCurves=mySeqOfCurve.Length();
1778   for (i=1; i<=aNbCurves; ++i) {
1779     const IntTools_Curve& aIC=mySeqOfCurve(i);
1780     //
1781     if (bToSplit) {
1782       Standard_Integer j, aNbC;
1783       IntTools_SequenceOfCurves aSeqCvs;
1784       //
1785       aNbC=IntTools_Tools::SplitCurve(aIC, aSeqCvs);
1786       if (aNbC) {
1787         for (j=1; j<=aNbC; ++j) {
1788           const IntTools_Curve& aICNew=aSeqCvs(j);
1789           aNewCvs.Append(aICNew);
1790         }
1791       }
1792       else {
1793         aNewCvs.Append(aIC);
1794       }
1795     }
1796     else {
1797       aNewCvs.Append(aIC);
1798     }
1799   }
1800   //
1801   // 2. Plane\Cone intersection when we had 4 curves
1802   aType1=myHS1->GetType();
1803   aType2=myHS2->GetType();
1804   aNbCurves=aNewCvs.Length();
1805   //
1806   if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone) ||
1807       (aType2==GeomAbs_Plane && aType1==GeomAbs_Cone)) {
1808     if (aNbCurves==4) {
1809       GeomAbs_CurveType aCType1;
1810       //
1811       aCType1=aNewCvs(1).Type();
1812       if (aCType1==GeomAbs_Line) {
1813         IntTools_SequenceOfCurves aSeqIn, aSeqOut;
1814         //
1815         for (i=1; i<=aNbCurves; ++i) {
1816           const IntTools_Curve& aIC=aNewCvs(i);
1817           aSeqIn.Append(aIC);
1818         }
1819         //
1820         IntTools_Tools::RejectLines(aSeqIn, aSeqOut);
1821         //
1822         aNewCvs.Clear();
1823         aNbCurves=aSeqOut.Length();
1824         for (i=1; i<=aNbCurves; ++i) {
1825           const IntTools_Curve& aIC=aSeqOut(i);
1826           aNewCvs.Append(aIC);
1827         }
1828       }
1829     }
1830   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Cone)...
1831   //
1832   // 3. Fill  mySeqOfCurve
1833   mySeqOfCurve.Clear();
1834   aNbCurves=aNewCvs.Length();
1835   for (i=1; i<=aNbCurves; ++i) {
1836     const IntTools_Curve& aIC=aNewCvs(i);
1837     mySeqOfCurve.Append(aIC);
1838   }
1839 }
1840 //=======================================================================
1841 //function : CorrectSurfaceBoundaries
1842 //purpose  :
1843 //=======================================================================
CorrectSurfaceBoundaries(const TopoDS_Face & theFace,const Standard_Real theTolerance,Standard_Real & theumin,Standard_Real & theumax,Standard_Real & thevmin,Standard_Real & thevmax)1844  void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
1845                               const Standard_Real theTolerance,
1846                               Standard_Real&      theumin,
1847                               Standard_Real&      theumax,
1848                               Standard_Real&      thevmin,
1849                               Standard_Real&      thevmax)
1850 {
1851   Standard_Boolean enlarge, isuperiodic, isvperiodic;
1852   Standard_Real uinf, usup, vinf, vsup, delta;
1853   GeomAbs_SurfaceType aType;
1854   Handle(Geom_Surface) aSurface;
1855   //
1856   aSurface = BRep_Tool::Surface(theFace);
1857   aSurface->Bounds(uinf, usup, vinf, vsup);
1858   delta = theTolerance;
1859   enlarge = Standard_False;
1860   //
1861   GeomAdaptor_Surface anAdaptorSurface(aSurface);
1862   //
1863   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1864     Handle(Geom_Surface) aBasisSurface =
1865       (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface))->BasisSurface();
1866 
1867     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
1868        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
1869       return;
1870     }
1871   }
1872   //
1873   if(aSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
1874     Handle(Geom_Surface) aBasisSurface =
1875       (Handle(Geom_OffsetSurface)::DownCast(aSurface))->BasisSurface();
1876 
1877     if(aBasisSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)) ||
1878        aBasisSurface->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
1879       return;
1880     }
1881   }
1882   //
1883   isuperiodic = anAdaptorSurface.IsUPeriodic();
1884   isvperiodic = anAdaptorSurface.IsVPeriodic();
1885   //
1886   aType=anAdaptorSurface.GetType();
1887   if((aType==GeomAbs_BezierSurface) ||
1888      (aType==GeomAbs_BSplineSurface) ||
1889      (aType==GeomAbs_SurfaceOfExtrusion) ||
1890      (aType==GeomAbs_SurfaceOfRevolution) ||
1891      (aType==GeomAbs_Cylinder)) {
1892     enlarge=Standard_True;
1893   }
1894   //
1895   if(!isuperiodic && enlarge) {
1896 
1897     if(!Precision::IsInfinite(theumin) &&
1898         ((theumin - uinf) > delta))
1899       theumin -= delta;
1900     else {
1901       theumin = uinf;
1902     }
1903 
1904     if(!Precision::IsInfinite(theumax) &&
1905         ((usup - theumax) > delta))
1906       theumax += delta;
1907     else
1908       theumax = usup;
1909   }
1910   //
1911   if(!isvperiodic && enlarge) {
1912     if(!Precision::IsInfinite(thevmin) &&
1913         ((thevmin - vinf) > delta)) {
1914       thevmin -= delta;
1915     }
1916     else {
1917       thevmin = vinf;
1918     }
1919     if(!Precision::IsInfinite(thevmax) &&
1920         ((vsup - thevmax) > delta)) {
1921       thevmax += delta;
1922     }
1923     else {
1924       thevmax = vsup;
1925     }
1926   }
1927   //
1928   if(isuperiodic || isvperiodic) {
1929     Standard_Boolean correct = Standard_False;
1930     Standard_Boolean correctU = Standard_False;
1931     Standard_Boolean correctV = Standard_False;
1932     Bnd_Box2d aBox;
1933     TopExp_Explorer anExp;
1934 
1935     for(anExp.Init(theFace, TopAbs_EDGE); anExp.More(); anExp.Next()) {
1936       if(BRep_Tool::IsClosed(TopoDS::Edge(anExp.Current()), theFace)) {
1937         correct = Standard_True;
1938         Standard_Real f, l;
1939         TopoDS_Edge anEdge = TopoDS::Edge(anExp.Current());
1940 
1941         for(Standard_Integer i = 0; i < 2; i++) {
1942           if(i==0) {
1943             anEdge.Orientation(TopAbs_FORWARD);
1944           }
1945           else {
1946             anEdge.Orientation(TopAbs_REVERSED);
1947           }
1948           Handle(Geom2d_Curve) aCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, f, l);
1949 
1950           if(aCurve.IsNull()) {
1951             correct = Standard_False;
1952             break;
1953           }
1954           Handle(Geom2d_Line) aLine = Handle(Geom2d_Line)::DownCast(aCurve);
1955 
1956           if(aLine.IsNull()) {
1957             correct = Standard_False;
1958             break;
1959           }
1960           gp_Dir2d anUDir(1., 0.);
1961           gp_Dir2d aVDir(0., 1.);
1962           Standard_Real anAngularTolerance = Precision::Angular();
1963 
1964           correctU = correctU || aLine->Position().Direction().IsParallel(aVDir, anAngularTolerance);
1965           correctV = correctV || aLine->Position().Direction().IsParallel(anUDir, anAngularTolerance);
1966 
1967           gp_Pnt2d pp1 = aCurve->Value(f);
1968           aBox.Add(pp1);
1969           gp_Pnt2d pp2 = aCurve->Value(l);
1970           aBox.Add(pp2);
1971         }
1972         if(!correct)
1973           break;
1974       }
1975     }
1976 
1977     if(correct) {
1978       Standard_Real umin, vmin, umax, vmax;
1979       aBox.Get(umin, vmin, umax, vmax);
1980 
1981       if(isuperiodic && correctU) {
1982         if(theumin < umin)
1983           theumin = umin;
1984         if(theumax > umax) {
1985           theumax = umax;
1986         }
1987       }
1988       if(isvperiodic && correctV) {
1989         if(thevmin < vmin)
1990           thevmin = vmin;
1991         if(thevmax > vmax)
1992           thevmax = vmax;
1993       }
1994     }
1995   }
1996 }
1997 
1998 // ------------------------------------------------------------------------------------------------
1999 // static function: ParameterOutOfBoundary
2000 // purpose:         Computes a new parameter for given curve. The corresponding 2d points
2001 //                  does not lay on any boundary of given faces
2002 // ------------------------------------------------------------------------------------------------
ParameterOutOfBoundary(const Standard_Real theParameter,const Handle (Geom_Curve)& theCurve,const TopoDS_Face & theFace1,const TopoDS_Face & theFace2,const Standard_Real theOtherParameter,const Standard_Boolean bIncreasePar,const Standard_Real theTol,Standard_Real & theNewParameter,const Handle (IntTools_Context)& aContext)2003 Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
2004                                         const Handle(Geom_Curve)& theCurve,
2005                                         const TopoDS_Face&        theFace1,
2006                                         const TopoDS_Face&        theFace2,
2007                                         const Standard_Real       theOtherParameter,
2008                                         const Standard_Boolean    bIncreasePar,
2009                                         const Standard_Real       theTol,
2010                                         Standard_Real&            theNewParameter,
2011                                         const Handle(IntTools_Context)& aContext)
2012 {
2013   Standard_Boolean bIsComputed = Standard_False;
2014   theNewParameter = theParameter;
2015 
2016   Standard_Real acurpar = theParameter;
2017   TopAbs_State aState = TopAbs_ON;
2018   Standard_Integer iter = 0;
2019   Standard_Real asumtol = theTol;
2020   Standard_Real adelta = asumtol * 0.1;
2021   adelta = (adelta < Precision::Confusion()) ? Precision::Confusion() : adelta;
2022   Handle(Geom_Surface) aSurf1 = BRep_Tool::Surface(theFace1);
2023   Handle(Geom_Surface) aSurf2 = BRep_Tool::Surface(theFace2);
2024 
2025   Standard_Real u1, u2, v1, v2;
2026 
2027   GeomAPI_ProjectPointOnSurf aPrj1;
2028   aSurf1->Bounds(u1, u2, v1, v2);
2029   aPrj1.Init(aSurf1, u1, u2, v1, v2);
2030 
2031   GeomAPI_ProjectPointOnSurf aPrj2;
2032   aSurf2->Bounds(u1, u2, v1, v2);
2033   aPrj2.Init(aSurf2, u1, u2, v1, v2);
2034 
2035   while(aState == TopAbs_ON) {
2036     if(bIncreasePar)
2037       acurpar += adelta;
2038     else
2039       acurpar -= adelta;
2040     gp_Pnt aPCurrent = theCurve->Value(acurpar);
2041     aPrj1.Perform(aPCurrent);
2042     Standard_Real U=0., V=0.;
2043 
2044     if(aPrj1.IsDone()) {
2045       aPrj1.LowerDistanceParameters(U, V);
2046       aState = aContext->StatePointFace(theFace1, gp_Pnt2d(U, V));
2047     }
2048 
2049     if(aState != TopAbs_ON) {
2050       aPrj2.Perform(aPCurrent);
2051 
2052       if(aPrj2.IsDone()) {
2053         aPrj2.LowerDistanceParameters(U, V);
2054         aState = aContext->StatePointFace(theFace2, gp_Pnt2d(U, V));
2055       }
2056     }
2057 
2058     if(iter > 11) {
2059       break;
2060     }
2061     iter++;
2062   }
2063 
2064   if(iter <= 11) {
2065     theNewParameter = acurpar;
2066     bIsComputed = Standard_True;
2067 
2068     if(bIncreasePar) {
2069       if(acurpar >= theOtherParameter)
2070         theNewParameter = theOtherParameter;
2071     }
2072     else {
2073       if(acurpar <= theOtherParameter)
2074         theNewParameter = theOtherParameter;
2075     }
2076   }
2077   return bIsComputed;
2078 }
2079 
2080 //=======================================================================
2081 //function : IsCurveValid
2082 //purpose  :
2083 //=======================================================================
IsCurveValid(const Handle (Geom2d_Curve)& thePCurve)2084 Standard_Boolean IsCurveValid (const Handle(Geom2d_Curve)& thePCurve)
2085 {
2086   if(thePCurve.IsNull())
2087     return Standard_False;
2088 
2089   Standard_Real tolint = 1.e-10;
2090   Geom2dAdaptor_Curve PCA;
2091   IntRes2d_Domain PCD;
2092   Geom2dInt_GInter PCI;
2093 
2094   Standard_Real pf = 0., pl = 0.;
2095   gp_Pnt2d pntf, pntl;
2096 
2097   if(!thePCurve->IsClosed() && !thePCurve->IsPeriodic()) {
2098     pf = thePCurve->FirstParameter();
2099     pl = thePCurve->LastParameter();
2100     pntf = thePCurve->Value(pf);
2101     pntl = thePCurve->Value(pl);
2102     PCA.Load(thePCurve);
2103     if(!PCA.IsPeriodic()) {
2104       if(PCA.FirstParameter() > pf) pf = PCA.FirstParameter();
2105       if(PCA.LastParameter()  < pl) pl = PCA.LastParameter();
2106     }
2107     PCD.SetValues(pntf,pf,tolint,pntl,pl,tolint);
2108     PCI.Perform(PCA,PCD,tolint,tolint);
2109     if(PCI.IsDone())
2110       if(PCI.NbPoints() > 0) {
2111         return Standard_False;
2112       }
2113   }
2114 
2115   return Standard_True;
2116 }
2117 
2118 //=======================================================================
2119 //static function : ApproxWithPCurves
2120 //purpose  : for bug 20964 only
2121 //=======================================================================
ApproxWithPCurves(const gp_Cylinder & theCyl,const gp_Sphere & theSph)2122 Standard_Boolean ApproxWithPCurves(const gp_Cylinder& theCyl,
2123                                    const gp_Sphere& theSph)
2124 {
2125   Standard_Boolean bRes = Standard_True;
2126   Standard_Real R1 = theCyl.Radius(), R2 = theSph.Radius();
2127   //
2128   {
2129     Standard_Real aD2, aRc2, aEps;
2130     gp_Pnt aApexSph;
2131     //
2132     aEps=1.E-7;
2133     aRc2=R1*R1;
2134     //
2135     const gp_Ax3& aAx3Sph=theSph.Position();
2136     const gp_Pnt& aLocSph=aAx3Sph.Location();
2137     const gp_Dir& aDirSph=aAx3Sph.Direction();
2138     //
2139     const gp_Ax1& aAx1Cyl=theCyl.Axis();
2140     gp_Lin aLinCyl(aAx1Cyl);
2141     //
2142     aApexSph.SetXYZ(aLocSph.XYZ()+R2*aDirSph.XYZ());
2143     aD2=aLinCyl.SquareDistance(aApexSph);
2144     if (fabs(aD2-aRc2)<aEps) {
2145       return !bRes;
2146     }
2147     //
2148     aApexSph.SetXYZ(aLocSph.XYZ()-R2*aDirSph.XYZ());
2149     aD2=aLinCyl.SquareDistance(aApexSph);
2150     if (fabs(aD2-aRc2)<aEps) {
2151       return !bRes;
2152     }
2153   }
2154   //
2155 
2156   if(R1 < 2.*R2) {
2157     return bRes;
2158   }
2159   gp_Lin anCylAx(theCyl.Axis());
2160 
2161   Standard_Real aDist = anCylAx.Distance(theSph.Location());
2162   Standard_Real aDRel = Abs(aDist - R1)/R2;
2163 
2164   if(aDRel > .2) return bRes;
2165 
2166   Standard_Real par = ElCLib::Parameter(anCylAx, theSph.Location());
2167   gp_Pnt aP = ElCLib::Value(par, anCylAx);
2168   gp_Vec aV(aP, theSph.Location());
2169 
2170   Standard_Real dd = aV.Dot(theSph.Position().XDirection());
2171 
2172   if(aDist < R1 && dd > 0.) return Standard_False;
2173   if(aDist > R1 && dd < 0.) return Standard_False;
2174 
2175 
2176   return bRes;
2177 }
2178 //=======================================================================
2179 //function : PerformPlanes
2180 //purpose  :
2181 //=======================================================================
PerformPlanes(const Handle (GeomAdaptor_Surface)& theS1,const Handle (GeomAdaptor_Surface)& theS2,const Standard_Real TolF1,const Standard_Real TolF2,const Standard_Real TolAng,const Standard_Real TolTang,const Standard_Boolean theApprox1,const Standard_Boolean theApprox2,IntTools_SequenceOfCurves & theSeqOfCurve,Standard_Boolean & theTangentFaces)2182 void  PerformPlanes(const Handle(GeomAdaptor_Surface)& theS1,
2183                     const Handle(GeomAdaptor_Surface)& theS2,
2184                     const Standard_Real TolF1,
2185                     const Standard_Real TolF2,
2186                     const Standard_Real TolAng,
2187                     const Standard_Real TolTang,
2188                     const Standard_Boolean theApprox1,
2189                     const Standard_Boolean theApprox2,
2190                     IntTools_SequenceOfCurves& theSeqOfCurve,
2191                     Standard_Boolean& theTangentFaces)
2192 {
2193 
2194   gp_Pln aPln1 = theS1->Plane();
2195   gp_Pln aPln2 = theS2->Plane();
2196 
2197   IntAna_QuadQuadGeo aPlnInter(aPln1, aPln2, TolAng, TolTang);
2198 
2199   if(!aPlnInter.IsDone()) {
2200     theTangentFaces = Standard_False;
2201     return;
2202   }
2203 
2204   IntAna_ResultType aResType = aPlnInter.TypeInter();
2205 
2206   if(aResType == IntAna_Same) {
2207     theTangentFaces = Standard_True;
2208     return;
2209   }
2210 
2211   theTangentFaces = Standard_False;
2212 
2213   if(aResType == IntAna_Empty) {
2214     return;
2215   }
2216 
2217   gp_Lin aLin = aPlnInter.Line(1);
2218 
2219   ProjLib_Plane aProj;
2220 
2221   aProj.Init(aPln1);
2222   aProj.Project(aLin);
2223   gp_Lin2d aLin2d1 = aProj.Line();
2224   //
2225   aProj.Init(aPln2);
2226   aProj.Project(aLin);
2227   gp_Lin2d aLin2d2 = aProj.Line();
2228   //
2229   //classify line2d1 relatively first plane
2230   Standard_Real P11, P12;
2231   Standard_Boolean IsCrossed = ClassifyLin2d(theS1, aLin2d1, TolTang, P11, P12);
2232   if(!IsCrossed) return;
2233   //classify line2d2 relatively second plane
2234   Standard_Real P21, P22;
2235   IsCrossed = ClassifyLin2d(theS2, aLin2d2, TolTang, P21, P22);
2236   if(!IsCrossed) return;
2237 
2238   //Analysis of parametric intervals: must have common part
2239 
2240   if(P21 >= P12) return;
2241   if(P22 <= P11) return;
2242 
2243   Standard_Real pmin, pmax;
2244   pmin = Max(P11, P21);
2245   pmax = Min(P12, P22);
2246 
2247   if(pmax - pmin <= TolTang) return;
2248 
2249   Handle(Geom_Line) aGLin = new Geom_Line(aLin);
2250 
2251   IntTools_Curve aCurve;
2252   Handle(Geom_TrimmedCurve) aGTLin = new Geom_TrimmedCurve(aGLin, pmin, pmax);
2253 
2254   aCurve.SetCurve(aGTLin);
2255 
2256   if(theApprox1) {
2257     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d1);
2258     aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
2259   }
2260   else {
2261     Handle(Geom2d_Curve) H1;
2262     aCurve.SetFirstCurve2d(H1);
2263   }
2264   if(theApprox2) {
2265     Handle(Geom2d_Line) C2d = new Geom2d_Line(aLin2d2);
2266     aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d, pmin, pmax));
2267   }
2268   else {
2269     Handle(Geom2d_Curve) H1;
2270     aCurve.SetFirstCurve2d(H1);
2271   }
2272   //
2273   // Valid tolerance for the intersection curve between planar faces
2274   // is the maximal tolerance between tolerances of faces
2275   Standard_Real aTolC = Max(TolF1, TolF2);
2276   aCurve.SetTolerance(aTolC);
2277   //
2278   // Computation of the tangential tolerance
2279   Standard_Real anAngle, aDt;
2280   gp_Dir aD1, aD2;
2281   //
2282   aD1 = aPln1.Position().Direction();
2283   aD2 = aPln2.Position().Direction();
2284   anAngle = aD1.Angle(aD2);
2285   //
2286   aDt = IntTools_Tools::ComputeIntRange(TolF1, TolF2, anAngle);
2287   Standard_Real aTangTol = sqrt(aDt*aDt + TolF1*TolF1);
2288   //
2289   aCurve.SetTangentialTolerance(aTangTol);
2290   //
2291   theSeqOfCurve.Append(aCurve);
2292 }
2293 
2294 //=======================================================================
2295 //function : ClassifyLin2d
2296 //purpose  :
2297 //=======================================================================
INTER(const Standard_Real d1,const Standard_Real d2,const Standard_Real tol)2298 static inline Standard_Boolean INTER(const Standard_Real d1,
2299                                      const Standard_Real d2,
2300                                      const Standard_Real tol)
2301 {
2302   return (d1 > tol && d2 < -tol) ||
2303          (d1 < -tol && d2 > tol) ||
2304          ((d1 <= tol && d1 >= -tol) && (d2 > tol || d2 < -tol)) ||
2305          ((d2 <= tol && d2 >= -tol) && (d1 > tol || d1 < -tol));
2306 }
COINC(const Standard_Real d1,const Standard_Real d2,const Standard_Real tol)2307 static inline Standard_Boolean COINC(const Standard_Real d1,
2308                                      const Standard_Real d2,
2309                                      const Standard_Real tol)
2310 {
2311   return (d1 <= tol && d1 >= -tol) && (d2 <= tol && d2 >= -tol);
2312 }
ClassifyLin2d(const Handle (GeomAdaptor_Surface)& theS,const gp_Lin2d & theLin2d,const Standard_Real theTol,Standard_Real & theP1,Standard_Real & theP2)2313 Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_Surface)& theS,
2314                                const gp_Lin2d& theLin2d,
2315                                const Standard_Real theTol,
2316                                Standard_Real& theP1,
2317                                Standard_Real& theP2)
2318 
2319 {
2320   Standard_Real xmin, xmax, ymin, ymax, d1, d2, A, B, C;
2321   Standard_Real par[2];
2322   Standard_Integer nbi = 0;
2323 
2324   xmin = theS->FirstUParameter();
2325   xmax = theS->LastUParameter();
2326   ymin = theS->FirstVParameter();
2327   ymax = theS->LastVParameter();
2328 
2329   theLin2d.Coefficients(A, B, C);
2330 
2331   //xmin, ymin <-> xmin, ymax
2332   d1 = A*xmin + B*ymin + C;
2333   d2 = A*xmin + B*ymax + C;
2334 
2335   if(INTER(d1, d2, theTol)) {
2336     //Intersection with boundary
2337     Standard_Real y = -(C + A*xmin)/B;
2338     par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, y));
2339     nbi++;
2340   }
2341   else if (COINC(d1, d2, theTol)) {
2342     //Coincidence with boundary
2343     par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
2344     par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
2345     nbi = 2;
2346   }
2347 
2348   if(nbi == 2) {
2349 
2350     if(fabs(par[0]-par[1]) > theTol) {
2351       theP1 = Min(par[0], par[1]);
2352       theP2 = Max(par[0], par[1]);
2353       return Standard_True;
2354     }
2355     else return Standard_False;
2356 
2357   }
2358 
2359   //xmin, ymax <-> xmax, ymax
2360   d1 = d2;
2361   d2 = A*xmax + B*ymax + C;
2362 
2363   if(d1 > theTol || d1 < -theTol) {//to avoid checking of
2364                                    //coincidence with the same point
2365     if(INTER(d1, d2, theTol)) {
2366       Standard_Real x = -(C + B*ymax)/A;
2367       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymax));
2368       nbi++;
2369     }
2370     else if (COINC(d1, d2, theTol)) {
2371       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymax));
2372       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
2373       nbi = 2;
2374     }
2375   }
2376 
2377   if(nbi == 2) {
2378 
2379     if(fabs(par[0]-par[1]) > theTol) {
2380       theP1 = Min(par[0], par[1]);
2381       theP2 = Max(par[0], par[1]);
2382       return Standard_True;
2383     }
2384     else return Standard_False;
2385 
2386   }
2387 
2388   //xmax, ymax <-> xmax, ymin
2389   d1 = d2;
2390   d2 = A*xmax + B*ymin + C;
2391 
2392   if(d1 > theTol || d1 < -theTol) {
2393     if(INTER(d1, d2, theTol)) {
2394       Standard_Real y = -(C + A*xmax)/B;
2395       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, y));
2396       nbi++;
2397     }
2398     else if (COINC(d1, d2, theTol)) {
2399       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymax));
2400       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
2401       nbi = 2;
2402     }
2403   }
2404 
2405   if(nbi == 2) {
2406     if(fabs(par[0]-par[1]) > theTol) {
2407       theP1 = Min(par[0], par[1]);
2408       theP2 = Max(par[0], par[1]);
2409       return Standard_True;
2410     }
2411     else return Standard_False;
2412   }
2413 
2414   //xmax, ymin <-> xmin, ymin
2415   d1 = d2;
2416   d2 = A*xmin + B*ymin + C;
2417 
2418   if(d1 > theTol || d1 < -theTol) {
2419     if(INTER(d1, d2, theTol)) {
2420       Standard_Real x = -(C + B*ymin)/A;
2421       par[nbi] = ElCLib::Parameter(theLin2d, gp_Pnt2d(x, ymin));
2422       nbi++;
2423     }
2424     else if (COINC(d1, d2, theTol)) {
2425       par[0] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmax, ymin));
2426       par[1] = ElCLib::Parameter(theLin2d, gp_Pnt2d(xmin, ymin));
2427       nbi = 2;
2428     }
2429   }
2430 
2431   if(nbi == 2) {
2432     if(fabs(par[0]-par[1]) > theTol) {
2433       theP1 = Min(par[0], par[1]);
2434       theP2 = Max(par[0], par[1]);
2435       return Standard_True;
2436     }
2437     else return Standard_False;
2438   }
2439 
2440   return Standard_False;
2441 
2442 }
2443 //
2444 //=======================================================================
2445 //function : ApproxParameters
2446 //purpose  :
2447 //=======================================================================
ApproxParameters(const Handle (GeomAdaptor_Surface)& aHS1,const Handle (GeomAdaptor_Surface)& aHS2,Standard_Integer & iDegMin,Standard_Integer & iDegMax,Standard_Integer & iNbIter)2448 void ApproxParameters(const Handle(GeomAdaptor_Surface)& aHS1,
2449                       const Handle(GeomAdaptor_Surface)& aHS2,
2450                       Standard_Integer& iDegMin,
2451                       Standard_Integer& iDegMax,
2452                       Standard_Integer& iNbIter)
2453 
2454 {
2455   GeomAbs_SurfaceType aTS1, aTS2;
2456 
2457   //
2458   iNbIter=0;
2459   iDegMin=4;
2460   iDegMax=8;
2461   //
2462   aTS1=aHS1->GetType();
2463   aTS2=aHS2->GetType();
2464   //
2465   // Cylinder/Torus
2466   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
2467       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
2468     Standard_Real aRC, aRT, dR, aPC;
2469     gp_Cylinder aCylinder;
2470     gp_Torus aTorus;
2471     //
2472     aPC=Precision::Confusion();
2473     //
2474     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Cylinder() : aHS2->Cylinder();
2475     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Torus() : aHS2->Torus();
2476     //
2477     aRC=aCylinder.Radius();
2478     aRT=aTorus.MinorRadius();
2479     dR=aRC-aRT;
2480     if (dR<0.) {
2481       dR=-dR;
2482     }
2483     //
2484     if (dR<aPC) {
2485       iDegMax=6;
2486     }
2487   }
2488   if (aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Cylinder) {
2489     iNbIter=1;
2490   }
2491 }
2492 //=======================================================================
2493 //function : Tolerances
2494 //purpose  :
2495 //=======================================================================
Tolerances(const Handle (GeomAdaptor_Surface)& aHS1,const Handle (GeomAdaptor_Surface)& aHS2,Standard_Real & aTolTang)2496 void Tolerances(const Handle(GeomAdaptor_Surface)& aHS1,
2497                 const Handle(GeomAdaptor_Surface)& aHS2,
2498                 Standard_Real& aTolTang)
2499 {
2500   GeomAbs_SurfaceType aTS1, aTS2;
2501   //
2502   aTS1=aHS1->GetType();
2503   aTS2=aHS2->GetType();
2504   //
2505   // Cylinder/Torus
2506   if ((aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Torus) ||
2507       (aTS2==GeomAbs_Cylinder && aTS1==GeomAbs_Torus)) {
2508     Standard_Real aRC, aRT, dR, aPC;
2509     gp_Cylinder aCylinder;
2510     gp_Torus aTorus;
2511     //
2512     aPC=Precision::Confusion();
2513     //
2514     aCylinder=(aTS1==GeomAbs_Cylinder)? aHS1->Cylinder() : aHS2->Cylinder();
2515     aTorus=(aTS1==GeomAbs_Torus)? aHS1->Torus() : aHS2->Torus();
2516     //
2517     aRC=aCylinder.Radius();
2518     aRT=aTorus.MinorRadius();
2519     dR=aRC-aRT;
2520     if (dR<0.) {
2521       dR=-dR;
2522     }
2523     //
2524     if (dR<aPC) {
2525       aTolTang=0.1*aTolTang;
2526     }
2527   }
2528 }
2529 //=======================================================================
2530 //function : SortTypes
2531 //purpose  :
2532 //=======================================================================
SortTypes(const GeomAbs_SurfaceType aType1,const GeomAbs_SurfaceType aType2)2533 Standard_Boolean SortTypes(const GeomAbs_SurfaceType aType1,
2534                            const GeomAbs_SurfaceType aType2)
2535 {
2536   Standard_Boolean bRet;
2537   Standard_Integer aI1, aI2;
2538   //
2539   bRet=Standard_False;
2540   //
2541   aI1=IndexType(aType1);
2542   aI2=IndexType(aType2);
2543   if (aI1<aI2){
2544     bRet=!bRet;
2545   }
2546   return bRet;
2547 }
2548 //=======================================================================
2549 //function : IndexType
2550 //purpose  :
2551 //=======================================================================
IndexType(const GeomAbs_SurfaceType aType)2552 Standard_Integer IndexType(const GeomAbs_SurfaceType aType)
2553 {
2554   Standard_Integer aIndex;
2555   //
2556   aIndex=11;
2557   //
2558   if (aType==GeomAbs_Plane) {
2559     aIndex=0;
2560   }
2561   else if (aType==GeomAbs_Cylinder) {
2562     aIndex=1;
2563   }
2564   else if (aType==GeomAbs_Cone) {
2565     aIndex=2;
2566   }
2567   else if (aType==GeomAbs_Sphere) {
2568     aIndex=3;
2569   }
2570   else if (aType==GeomAbs_Torus) {
2571     aIndex=4;
2572   }
2573   else if (aType==GeomAbs_BezierSurface) {
2574     aIndex=5;
2575   }
2576   else if (aType==GeomAbs_BSplineSurface) {
2577     aIndex=6;
2578   }
2579   else if (aType==GeomAbs_SurfaceOfRevolution) {
2580     aIndex=7;
2581   }
2582   else if (aType==GeomAbs_SurfaceOfExtrusion) {
2583     aIndex=8;
2584   }
2585   else if (aType==GeomAbs_OffsetSurface) {
2586     aIndex=9;
2587   }
2588   else if (aType==GeomAbs_OtherSurface) {
2589     aIndex=10;
2590   }
2591   return aIndex;
2592 }
2593 
2594 //=======================================================================
2595 // Function : FindMaxDistance
2596 // purpose :
2597 //=======================================================================
FindMaxDistance(const Handle (Geom_Curve)& theCurve,const Standard_Real theFirst,const Standard_Real theLast,const TopoDS_Face & theFace,const Handle (IntTools_Context)& theContext)2598 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theCurve,
2599                               const Standard_Real theFirst,
2600                               const Standard_Real theLast,
2601                               const TopoDS_Face& theFace,
2602                               const Handle(IntTools_Context)& theContext)
2603 {
2604   Standard_Integer aNbS;
2605   Standard_Real aT1, aT2, aDt, aD, aDMax, anEps;
2606   //
2607   aNbS = 11;
2608   aDt = (theLast - theFirst) / aNbS;
2609   aDMax = 0.;
2610   anEps = 1.e-4 * aDt;
2611   //
2612   GeomAPI_ProjectPointOnSurf& aProjPS = theContext->ProjPS(theFace);
2613   aT2 = theFirst;
2614   for (;;) {
2615     aT1 = aT2;
2616     aT2 += aDt;
2617     //
2618     if (aT2 > theLast) {
2619       break;
2620     }
2621     //
2622     aD = FindMaxDistance(theCurve, aT1, aT2, aProjPS, anEps);
2623     if (aD > aDMax) {
2624       aDMax = aD;
2625     }
2626   }
2627   //
2628   return aDMax;
2629 }
2630 
2631 //=======================================================================
2632 // Function : FindMaxDistance
2633 // purpose :
2634 //=======================================================================
FindMaxDistance(const Handle (Geom_Curve)& theC,const Standard_Real theFirst,const Standard_Real theLast,GeomAPI_ProjectPointOnSurf & theProjPS,const Standard_Real theEps)2635 Standard_Real FindMaxDistance(const Handle(Geom_Curve)& theC,
2636                               const Standard_Real theFirst,
2637                               const Standard_Real theLast,
2638                               GeomAPI_ProjectPointOnSurf& theProjPS,
2639                               const Standard_Real theEps)
2640 {
2641   Standard_Real aA, aB, aCf, aX, aX1, aX2, aF1, aF2, aF;
2642   //
2643   aCf = 0.61803398874989484820458683436564;//(sqrt(5.)-1)/2.;
2644   aA = theFirst;
2645   aB = theLast;
2646   //
2647   aX1=aB - aCf*(aB-aA);
2648   aF1 = MaxDistance(theC, aX1, theProjPS);
2649   aX2 = aA + aCf * (aB - aA);
2650   aF2 = MaxDistance(theC, aX2, theProjPS);
2651 
2652   while (Abs(aX1-aX2) > theEps)
2653   {
2654     if (aF1 > aF2) {
2655       aB = aX2;
2656       aX2 = aX1;
2657       aF2 = aF1;
2658       aX1 = aB-aCf*(aB-aA);
2659       aF1 = MaxDistance(theC, aX1, theProjPS);
2660     }
2661     else {
2662       aA = aX1;
2663       aX1 = aX2;
2664       aF1 = aF2;
2665       aX2=aA+aCf*(aB-aA);
2666       aF2 = MaxDistance(theC, aX2, theProjPS);
2667     }
2668   }
2669   //
2670   aX = 0.5 * (aA + aB);
2671   aF = MaxDistance(theC, aX, theProjPS);
2672   //
2673   if (aF1 > aF) {
2674     aF = aF1;
2675   }
2676   //
2677   if (aF2 > aF) {
2678     aF = aF2;
2679   }
2680   //
2681   return aF;
2682 }
2683 
2684 //=======================================================================
2685 // Function : MaxDistance
2686 // purpose :
2687 //=======================================================================
MaxDistance(const Handle (Geom_Curve)& theC,const Standard_Real aT,GeomAPI_ProjectPointOnSurf & theProjPS)2688 Standard_Real MaxDistance(const Handle(Geom_Curve)& theC,
2689                           const Standard_Real aT,
2690                           GeomAPI_ProjectPointOnSurf& theProjPS)
2691 {
2692   Standard_Real aD;
2693   gp_Pnt aP;
2694   //
2695   theC->D0(aT, aP);
2696   theProjPS.Perform(aP);
2697   aD = theProjPS.NbPoints() ? theProjPS.LowerDistance() : 0.;
2698   //
2699   return aD;
2700 }
2701 
2702 //=======================================================================
2703 //function : CheckPCurve
2704 //purpose  : Checks if points of the pcurve are out of the face bounds.
2705 //=======================================================================
CheckPCurve(const Handle (Geom2d_Curve)& aPC,const TopoDS_Face & aFace,const Handle (IntTools_Context)& theCtx)2706   Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
2707                                const TopoDS_Face& aFace,
2708                                const Handle(IntTools_Context)& theCtx)
2709 {
2710   const Standard_Integer NPoints = 23;
2711   Standard_Integer i;
2712   Standard_Real umin,umax,vmin,vmax;
2713 
2714   theCtx->UVBounds(aFace, umin, umax, vmin, vmax);
2715   Standard_Real tolU = Max ((umax-umin)*0.01, Precision::Confusion());
2716   Standard_Real tolV = Max ((vmax-vmin)*0.01, Precision::Confusion());
2717   Standard_Real fp = aPC->FirstParameter();
2718   Standard_Real lp = aPC->LastParameter();
2719 
2720 
2721   // adjust domain for periodic surfaces
2722   TopLoc_Location aLoc;
2723   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace, aLoc);
2724   if (aSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2725     aSurf = (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf))->BasisSurface();
2726   }
2727   gp_Pnt2d pnt = aPC->Value((fp+lp)/2);
2728   Standard_Real u,v;
2729   pnt.Coord(u,v);
2730   //
2731   if (aSurf->IsUPeriodic()) {
2732     Standard_Real aPer = aSurf->UPeriod();
2733     Standard_Integer nshift = (Standard_Integer) ((u-umin)/aPer);
2734     if (u < umin+aPer*nshift) nshift--;
2735     umin += aPer*nshift;
2736     umax += aPer*nshift;
2737   }
2738   if (aSurf->IsVPeriodic()) {
2739     Standard_Real aPer = aSurf->VPeriod();
2740     Standard_Integer nshift = (Standard_Integer) ((v-vmin)/aPer);
2741     if (v < vmin+aPer*nshift) nshift--;
2742     vmin += aPer*nshift;
2743     vmax += aPer*nshift;
2744   }
2745   //
2746   //--------------------------------------------------------
2747   Standard_Boolean bRet;
2748   Standard_Integer j, aNbIntervals;
2749   Standard_Real aT, dT;
2750   gp_Pnt2d aP2D;
2751   //
2752   Geom2dAdaptor_Curve aGAC(aPC);
2753   aNbIntervals=aGAC.NbIntervals(GeomAbs_CN);
2754   //
2755   TColStd_Array1OfReal aTI(1, aNbIntervals+1);
2756   aGAC.Intervals(aTI,GeomAbs_CN);
2757   //
2758   bRet=Standard_False;
2759   //
2760   aT=aGAC.FirstParameter();
2761   for (j=1; j<=aNbIntervals; ++j) {
2762     dT=(aTI(j+1)-aTI(j))/NPoints;
2763     //
2764     for (i=1; i<NPoints; i++) {
2765       aT=aT+dT;
2766       aGAC.D0(aT, aP2D);
2767       aP2D.Coord(u,v);
2768     if (umin-u > tolU || u-umax > tolU ||
2769           vmin-v > tolV || v-vmax > tolV) {
2770         return bRet;
2771   }
2772 }
2773   }
2774   return !bRet;
2775 }
2776 //=======================================================================
2777 //function : CorrectPlaneBoundaries
2778 //purpose  :
2779 //=======================================================================
CorrectPlaneBoundaries(Standard_Real & aUmin,Standard_Real & aUmax,Standard_Real & aVmin,Standard_Real & aVmax)2780  void CorrectPlaneBoundaries(Standard_Real& aUmin,
2781                              Standard_Real& aUmax,
2782                              Standard_Real& aVmin,
2783                              Standard_Real& aVmax)
2784 {
2785   if (!(Precision::IsInfinite(aUmin) ||
2786         Precision::IsInfinite(aUmax))) {
2787     Standard_Real dU;
2788     //
2789     dU=0.1*(aUmax-aUmin);
2790     aUmin=aUmin-dU;
2791     aUmax=aUmax+dU;
2792   }
2793   if (!(Precision::IsInfinite(aVmin) ||
2794         Precision::IsInfinite(aVmax))) {
2795     Standard_Real dV;
2796     //
2797     dV=0.1*(aVmax-aVmin);
2798     aVmin=aVmin-dV;
2799     aVmax=aVmax+dV;
2800   }
2801 }
2802