0028639: Improve performance of the IntPolyh_MaillageAffinage algorithm
[occt.git] / src / IntPolyh / IntPolyh_MaillageAffinage.cxx
1 // Created on: 1999-03-05
2 // Created by: Fabrice SERVANT
3 // Copyright (c) 1999-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 //  modified by Edward AGAPOV (eap) Tue Jan 22 2002 (bug occ53)
18 //  - improve SectionLine table management (avoid memory reallocation)
19 //  - some protection against arrays overflow
20 //  modified by Edward AGAPOV (eap) Thu Feb 14 2002 (occ139)
21 //  - make Section Line parts rightly connected (prepend 2nd part to the 1st)
22 //  - TriangleShape() for debugging purpose
23 //  Modified by skv - Thu Sep 25 17:42:42 2003 OCC567
24 //  modified by ofv Thu Apr  8 14:58:13 2004 fip
25
26 #include <Adaptor3d_HSurface.hxx>
27 #include <Bnd_BoundSortBox.hxx>
28 #include <Bnd_Box.hxx>
29 #include <Bnd_HArray1OfBox.hxx>
30 #include <gp.hxx>
31 #include <gp_Pnt.hxx>
32 #include <IntCurveSurface_ThePolyhedronOfHInter.hxx>
33 #include <IntPolyh_ListOfCouples.hxx>
34 #include <IntPolyh_Couple.hxx>
35 #include <IntPolyh_Edge.hxx>
36 #include <IntPolyh_MaillageAffinage.hxx>
37 #include <IntPolyh_Point.hxx>
38 #include <IntPolyh_SectionLine.hxx>
39 #include <IntPolyh_StartPoint.hxx>
40 #include <IntPolyh_Triangle.hxx>
41 #include <Precision.hxx>
42 #include <TColStd_Array1OfInteger.hxx>
43 #include <TColStd_MapOfInteger.hxx>
44 #include <TColStd_ListIteratorOfListOfInteger.hxx>
45 #include <NCollection_UBTree.hxx>
46 #include <NCollection_UBTreeFiller.hxx>
47 #include <algorithm>
48 #include <NCollection_IndexedDataMap.hxx>
49
50 typedef NCollection_Array1<Standard_Integer> IntPolyh_ArrayOfInteger;
51 typedef NCollection_IndexedDataMap
52   <Standard_Integer,
53    IntPolyh_ArrayOfInteger,
54    TColStd_MapIntegerHasher> IntPolyh_IndexedDataMapOfIntegerArrayOfInteger;
55
56
57 static Standard_Real MyTolerance=10.0e-7;
58 static Standard_Real MyConfusionPrecision=10.0e-12;
59 static Standard_Real SquareMyConfusionPrecision=10.0e-24;
60 //
61 static
62   inline Standard_Real maxSR(const Standard_Real a,
63                              const Standard_Real b,
64                              const Standard_Real c);
65
66 static
67   inline Standard_Real minSR(const Standard_Real a,
68                              const Standard_Real b,
69                              const Standard_Real c);
70 static
71   Standard_Integer project6(const IntPolyh_Point &ax, 
72                             const IntPolyh_Point &p1,
73                             const IntPolyh_Point &p2, 
74                             const IntPolyh_Point &p3, 
75                             const IntPolyh_Point &q1,
76                             const IntPolyh_Point &q2, 
77                             const IntPolyh_Point &q3);
78 static
79   void TestNbPoints(const Standard_Integer ,
80                     Standard_Integer &NbPoints,
81                     Standard_Integer &NbPointsTotal,
82                     const IntPolyh_StartPoint &Pt1,
83                     const IntPolyh_StartPoint &Pt2,
84                     IntPolyh_StartPoint &SP1,
85                     IntPolyh_StartPoint &SP2);
86 static
87   void CalculPtsInterTriEdgeCoplanaires(const Standard_Integer TriSurfID,
88                                         const IntPolyh_Point &NormaleTri,
89                                         const IntPolyh_Triangle &Tri1,
90                                         const IntPolyh_Triangle &Tri2,
91                                         const IntPolyh_Point &PE1,
92                                         const IntPolyh_Point &PE2,
93                                         const IntPolyh_Point &Edge,
94                                         const Standard_Integer EdgeIndex,
95                                         const IntPolyh_Point &PT1,
96                                         const IntPolyh_Point &PT2,
97                                         const IntPolyh_Point &Cote,
98                                         const Standard_Integer CoteIndex,
99                                         IntPolyh_StartPoint &SP1,
100                                         IntPolyh_StartPoint &SP2,
101                                         Standard_Integer &NbPoints); 
102 static
103   Standard_Boolean CheckCoupleAndGetAngle(const Standard_Integer T1, 
104                                           const Standard_Integer T2,
105                                           Standard_Real& Angle, 
106                                           IntPolyh_ListOfCouples &TTrianglesContacts);
107 static
108   Standard_Boolean CheckCoupleAndGetAngle2(const Standard_Integer T1,
109                                            const Standard_Integer T2,
110                                            const Standard_Integer T11, 
111                                            const Standard_Integer T22,
112                                            IntPolyh_ListIteratorOfListOfCouples& theItCT11,
113                                            IntPolyh_ListIteratorOfListOfCouples& theItCT22,
114                                            Standard_Real & Angle,
115                                            IntPolyh_ListOfCouples &TTrianglesContacts);
116 static
117   Standard_Integer CheckNextStartPoint(IntPolyh_SectionLine & SectionLine,
118                                        IntPolyh_ArrayOfTangentZones & TTangentZones,
119                                        IntPolyh_StartPoint & SP,
120                                        const Standard_Boolean Prepend=Standard_False); 
121
122 static
123   Standard_Boolean IsDegenerated(const Handle(Adaptor3d_HSurface)& aS,
124                                  const Standard_Integer aIndex,
125                                  const Standard_Real aTol2,
126                                  Standard_Real& aDegX);
127 static
128   void DegeneratedIndex(const TColStd_Array1OfReal& Xpars,
129                         const Standard_Integer aNbX,
130                         const Handle(Adaptor3d_HSurface)& aS,
131                         const Standard_Integer aIsoDirection,
132                         Standard_Integer& aI1,
133                         Standard_Integer& aI2);
134 static
135   void EnlargeZone(const Handle(Adaptor3d_HSurface)& MaSurface,
136                    Standard_Real &u0, 
137                    Standard_Real &u1, 
138                    Standard_Real &v0, 
139                    Standard_Real &v1);
140
141 //=======================================================================
142 //class : IntPolyh_BoxBndTreeSelector
143 //purpose  : Select interfering bounding boxes
144 //=======================================================================
145 typedef NCollection_UBTree<Standard_Integer, Bnd_Box> IntPolyh_BoxBndTree;
146 class IntPolyh_BoxBndTreeSelector : public IntPolyh_BoxBndTree::Selector {
147  public:
148   IntPolyh_BoxBndTreeSelector(const Bnd_Box& theBox) : myBox(theBox) {}
149   //
150   virtual Standard_Boolean Reject(const Bnd_Box& theOther) const
151   {
152     return myBox.IsOut(theOther);
153   }
154   //
155   virtual Standard_Boolean Accept(const Standard_Integer &theInd)
156   {
157     myIndices.Append(theInd);
158     return Standard_True;
159   }
160   //
161   const TColStd_ListOfInteger& Indices() const
162   {
163     return myIndices;
164   }
165  private:
166   Bnd_Box myBox;
167   TColStd_ListOfInteger myIndices;
168 };
169
170 //=======================================================================
171 //function : GetInterferingTriangles
172 //purpose  : Returns indices of the triangles with interfering bounding boxes
173 //=======================================================================
174 static
175   void GetInterferingTriangles(IntPolyh_ArrayOfTriangles& theTriangles1,
176                                const IntPolyh_ArrayOfPoints& thePoints1,
177                                IntPolyh_ArrayOfTriangles& theTriangles2,
178                                const IntPolyh_ArrayOfPoints& thePoints2,
179                                IntPolyh_IndexedDataMapOfIntegerArrayOfInteger& theCouples)
180 {
181   // To find the triangles with interfering bounding boxes
182   // use the algorithm of unbalanced binary tree of overlapping bounding boxes
183   IntPolyh_BoxBndTree aBBTree;
184   NCollection_UBTreeFiller <Standard_Integer, Bnd_Box> aTreeFiller(aBBTree);
185   // 1. Fill the tree with the boxes of the triangles from second surface
186   Standard_Integer i, aNbT2 = theTriangles2.NbItems();
187   for (i = 0; i < aNbT2; ++i) {
188     IntPolyh_Triangle& aT = theTriangles2[i];
189     if (!aT.IsIntersectionPossible() || aT.IsDegenerated()) {
190       continue;
191     }
192     //
193     const Bnd_Box& aBox = aT.BoundingBox(thePoints2);
194     aTreeFiller.Add(i, aBox);
195   }
196   //
197   // 2. Shake the tree filler
198   aTreeFiller.Fill();
199   //
200   // 3. Find boxes interfering with the first triangles
201   Standard_Integer aNbT1 = theTriangles1.NbItems();
202   for (i = 0; i < aNbT1; ++i) {
203     IntPolyh_Triangle& aT = theTriangles1[i];
204     if (!aT.IsIntersectionPossible() || aT.IsDegenerated()) {
205       continue;
206     }
207     //
208     const Bnd_Box& aBox = aT.BoundingBox(thePoints1);
209     //
210     IntPolyh_BoxBndTreeSelector aSelector(aBox);
211     if (!aBBTree.Select(aSelector)) {
212       continue;
213     }
214     //
215     const TColStd_ListOfInteger& aLI = aSelector.Indices();
216     // Sort the indices
217     IntPolyh_ArrayOfInteger anArr(1, aLI.Extent());
218     TColStd_ListIteratorOfListOfInteger aItLI(aLI);
219     for (Standard_Integer j = 1; aItLI.More(); aItLI.Next(), ++j) {
220       anArr(j) = aItLI.Value();
221     }
222     //
223     std::sort(anArr.begin(), anArr.end());
224     //
225     theCouples.Add(i, anArr);
226   }
227 }
228
229 //=======================================================================
230 //function : IntPolyh_MaillageAffinage
231 //purpose  : 
232 //=======================================================================
233 IntPolyh_MaillageAffinage::IntPolyh_MaillageAffinage
234   (const Handle(Adaptor3d_HSurface)& Surface1,
235    const Handle(Adaptor3d_HSurface)& Surface2,
236    const Standard_Integer )
237 :
238   MaSurface1(Surface1), 
239   MaSurface2(Surface2), 
240   NbSamplesU1(10), 
241   NbSamplesU2(10), 
242   NbSamplesV1(10), 
243   NbSamplesV2(10),
244   FlecheMax1(0.0), 
245   FlecheMax2(0.0), 
246   FlecheMin1(0.0), 
247   FlecheMin2(0.0),
248   FlecheMoy1(0.0), 
249   FlecheMoy2(0.0), 
250   myEnlargeZone(Standard_False) 
251
252 }
253 //=======================================================================
254 //function : IntPolyh_MaillageAffinage
255 //purpose  : 
256 //=======================================================================
257 IntPolyh_MaillageAffinage::IntPolyh_MaillageAffinage
258   (const Handle(Adaptor3d_HSurface)& Surface1,
259    const Standard_Integer NbSU1,
260    const Standard_Integer NbSV1,
261    const Handle(Adaptor3d_HSurface)& Surface2,
262    const Standard_Integer NbSU2,
263    const Standard_Integer NbSV2,
264    const Standard_Integer )
265 :
266   MaSurface1(Surface1), 
267   MaSurface2(Surface2), 
268   NbSamplesU1(NbSU1), 
269   NbSamplesU2(NbSU2), 
270   NbSamplesV1(NbSV1), 
271   NbSamplesV2(NbSV2),
272   FlecheMax1(0.0), 
273   FlecheMax2(0.0), 
274   FlecheMin1(0.0), 
275   FlecheMin2(0.0),
276   FlecheMoy1(0.0), 
277   FlecheMoy2(0.0), 
278   myEnlargeZone(Standard_False)
279
280 }
281 //=======================================================================
282 //function : FillArrayOfPnt
283 //purpose  : Compute points on one surface and fill an array of points
284 //=======================================================================
285 void IntPolyh_MaillageAffinage::FillArrayOfPnt
286   (const Standard_Integer SurfID)
287 {
288   Standard_Integer NbSamplesU, NbSamplesV, i, aNbSamplesU1, aNbSamplesV1;
289   Standard_Real u0, u1, v0, v1, aU, aV, dU, dV;
290   //
291   const Handle(Adaptor3d_HSurface)& MaSurface=(SurfID==1)? MaSurface1 : MaSurface2;
292   NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
293   NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
294   //
295   u0 = (MaSurface)->FirstUParameter();  
296   u1 = (MaSurface)->LastUParameter();
297   v0 = (MaSurface)->FirstVParameter();  
298   v1 = (MaSurface)->LastVParameter();  
299
300   if(myEnlargeZone) { 
301     EnlargeZone(MaSurface, u0, u1, v0, v1);
302   }
303   //
304   TColStd_Array1OfReal aUpars(1, NbSamplesU);
305   TColStd_Array1OfReal aVpars(1, NbSamplesV);
306   //
307   aNbSamplesU1=NbSamplesU-1;
308   aNbSamplesV1=NbSamplesV-1;
309   //
310   dU=(u1-u0)/Standard_Real(aNbSamplesU1);
311   dV=(v1-v0)/Standard_Real(aNbSamplesV1);
312   //
313   for (i=0; i<NbSamplesU; ++i) {
314     aU=u0+i*dU;
315     if (i==aNbSamplesU1) {
316       aU=u1;
317     }
318     aUpars.SetValue(i+1, aU);
319   }
320   //
321   for (i=0; i<NbSamplesV; ++i) {
322     aV=v0+i*dV;
323     if (i==aNbSamplesV1) {
324       aV=v1;
325     }
326     aVpars.SetValue(i+1, aV);
327   }
328   //
329   FillArrayOfPnt(SurfID, aUpars, aVpars);
330 }
331 //=======================================================================
332 //function : FillArrayOfPnt
333 //purpose  : Compute points on one surface and fill an array of points
334 //           FILL AN ARRAY OF POINTS
335 //=======================================================================
336 void IntPolyh_MaillageAffinage::FillArrayOfPnt
337   (const Standard_Integer SurfID,
338    const Standard_Boolean isShiftFwd)
339 {
340   Standard_Integer NbSamplesU, NbSamplesV, i, aNbSamplesU1, aNbSamplesV1; 
341   Standard_Real u0, u1, v0, v1, aU, aV, dU, dV;
342   const Handle(Adaptor3d_HSurface)& MaSurface=(SurfID==1)? MaSurface1 : MaSurface2;
343   NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
344   NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
345
346   u0 = (MaSurface)->FirstUParameter();  
347   u1 = (MaSurface)->LastUParameter();
348   v0 = (MaSurface)->FirstVParameter();  
349   v1 = (MaSurface)->LastVParameter();  
350
351   if(myEnlargeZone) {
352     EnlargeZone(MaSurface, u0, u1, v0, v1);
353   }
354   //
355   TColStd_Array1OfReal aUpars(1, NbSamplesU);
356   TColStd_Array1OfReal aVpars(1, NbSamplesV);
357   //
358   aNbSamplesU1=NbSamplesU-1;
359   aNbSamplesV1=NbSamplesV-1;
360   //
361   dU=(u1-u0)/Standard_Real(aNbSamplesU1);
362   dV=(v1-v0)/Standard_Real(aNbSamplesV1);
363   //
364   for (i=0; i<NbSamplesU; ++i) {
365     aU=u0+i*dU;
366     if (i==aNbSamplesU1) {
367       aU=u1;
368     }
369     aUpars.SetValue(i+1, aU);
370   }
371   //
372   for (i=0; i<NbSamplesV; ++i) {
373     aV=v0+i*dV;
374     if (i==aNbSamplesV1) {
375       aV=v1;
376     }
377     aVpars.SetValue(i+1, aV);
378   }
379   //
380   FillArrayOfPnt(SurfID, isShiftFwd, aUpars, aVpars);  
381 }
382 //=======================================================================
383 //function : FillArrayOfPnt
384 //purpose  : Compute points on one surface and fill an array of points
385 //=======================================================================
386 void IntPolyh_MaillageAffinage::FillArrayOfPnt
387   (const Standard_Integer SurfID, 
388    const TColStd_Array1OfReal& Upars,
389    const TColStd_Array1OfReal& Vpars)
390 {
391   Standard_Boolean bDegI, bDeg;
392   Standard_Integer aNbU, aNbV, iCnt, i, j;
393   Standard_Integer aID1, aID2, aJD1, aJD2;
394   Standard_Real aTol, aU, aV, aX, aY, aZ;
395   gp_Pnt aP;
396   //
397   aNbU=(SurfID==1)? NbSamplesU1 : NbSamplesU2;
398   aNbV=(SurfID==1)? NbSamplesV1 : NbSamplesV2;
399   Bnd_Box& aBox = (SurfID==1) ? MyBox1 : MyBox2;
400   Handle(Adaptor3d_HSurface)& aS=(SurfID==1)? MaSurface1:MaSurface2;
401   IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
402   //
403   aJD1=0;
404   aJD2=0;
405   aID1=0;
406   aID2=0;
407   DegeneratedIndex(Vpars, aNbV, aS, 1, aJD1, aJD2);
408   if (!(aJD1 || aJD2)) {
409     DegeneratedIndex(Upars, aNbU, aS, 2, aID1, aID2);
410   }
411   //
412   TPoints.Init(aNbU*aNbV);
413   iCnt=0;
414   for(i=1; i<=aNbU; ++i){
415     bDegI=(aID1==i || aID2==i);
416     aU=Upars(i);
417     for(j=1; j<=aNbV; ++j){
418       aV=Vpars(j);
419       aP=aS->Value(aU, aV);
420       aP.Coord(aX, aY, aZ);
421       IntPolyh_Point& aIP=TPoints[iCnt];
422       aIP.Set(aX, aY, aZ, aU, aV);
423       //
424       bDeg=bDegI || (aJD1==j || aJD2==j);
425       if (bDeg) {
426         aIP.SetDegenerated(bDeg);
427       }
428       ++iCnt;
429       aBox.Add(aP);
430     }
431   }
432   //
433   TPoints.SetNbItems(iCnt);
434   //
435   IntCurveSurface_ThePolyhedronOfHInter polyhedron(aS, Upars, Vpars);
436   //
437   aTol=polyhedron.DeflectionOverEstimation();
438   aTol*=1.2;
439
440   Standard_Real a1,a2,a3,b1,b2,b3;
441   //
442   aBox.Get(a1,a2,a3,b1,b2,b3);
443   aBox.Update(a1-aTol,a2-aTol,a3-aTol,b1+aTol,b2+aTol,b3+aTol);
444   aBox.Enlarge(MyTolerance);
445 }
446
447 //=======================================================================
448 //function : FillArrayOfPnt
449 //purpose  : Compute points on one surface and fill an array of points
450 //           REMPLISSAGE DU TABLEAU DE POINTS
451 //=======================================================================
452 void IntPolyh_MaillageAffinage::FillArrayOfPnt
453   (const Standard_Integer SurfID,
454    const Standard_Boolean isShiftFwd,
455    const TColStd_Array1OfReal& Upars,
456    const TColStd_Array1OfReal& Vpars)
457 {
458   Standard_Boolean bDegI, bDeg;
459   Standard_Integer aNbU, aNbV, iCnt, i, j;
460   Standard_Integer aID1, aID2, aJD1, aJD2;
461   Standard_Real Tol, resol, aU, aV, aMag;
462   Standard_Real aX, aY, aZ;
463   gp_Pnt aP;
464   gp_Vec aDU, aDV, aNorm;
465   //
466   aNbU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
467   aNbV=(SurfID==1)? NbSamplesV1:NbSamplesV2; 
468   Bnd_Box& aBox = (SurfID==1) ? MyBox1 : MyBox2;
469   Handle(Adaptor3d_HSurface) aS=(SurfID==1)? MaSurface1:MaSurface2;
470   IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
471   //
472   resol = gp::Resolution();
473   //
474   IntCurveSurface_ThePolyhedronOfHInter polyhedron(aS, Upars, Vpars);
475   Tol=polyhedron.DeflectionOverEstimation();
476   aJD1=0;
477   aJD2=0;
478   aID1=0;
479   aID2=0;
480   DegeneratedIndex(Vpars, aNbV, aS, 1, aJD1, aJD2);
481   if (!(aJD1 || aJD2)) {
482     DegeneratedIndex(Upars, aNbU, aS, 2, aID1, aID2);
483   }
484   //
485   TPoints.Init(aNbU*aNbV);
486   iCnt=0;
487   for(i=1; i<=aNbU; ++i){
488     bDegI=(aID1==i || aID2==i);
489     aU = Upars(i);
490     for(j=1; j<=aNbV; ++j){
491       aV = Vpars(j);
492       aS->D1(aU, aV, aP, aDU, aDV);
493       
494       aNorm = aDU.Crossed(aDV);
495       aMag = aNorm.Magnitude();
496       if (aMag > resol) {
497         aNorm /= aMag;
498         aNorm.Multiply(Tol*1.5);
499         //
500         if (isShiftFwd) {
501           aP.Translate(aNorm);
502         }
503         else{
504           aP.Translate(aNorm.Reversed());
505         }
506       }
507       
508       IntPolyh_Point& aIP=TPoints[iCnt];
509       aP.Coord(aX, aY, aZ);
510       aIP.Set(aX, aY, aZ, aU, aV);
511       //
512       bDeg=bDegI || (aJD1==j || aJD2==j);
513       if (bDeg) {
514         aIP.SetDegenerated(bDeg);
515       }
516       ++iCnt;
517       aBox.Add(aP);
518     }
519   }
520   //
521   TPoints.SetNbItems(iCnt);
522   //
523   Tol*=1.2;
524   //
525   Standard_Real a1,a2,a3,b1,b2,b3;
526   //
527   aBox.Get(a1,a2,a3,b1,b2,b3);
528   aBox.Update(a1-Tol,a2-Tol,a3-Tol,b1+Tol,b2+Tol,b3+Tol);
529   aBox.Enlarge(MyTolerance);
530 }
531 //=======================================================================
532 //function : CommonBox
533 //purpose  : Compute the common box  witch is the intersection
534 //           of the two bounding boxes,  and mark the points of
535 //           the two surfaces that are inside.
536 //           REJECTION BOUNDING BOXES
537 //           DETERMINATION OF THE COMMON BOX
538 //=======================================================================
539 void IntPolyh_MaillageAffinage::CommonBox (const Bnd_Box &,
540                                            const Bnd_Box &,
541                                            Standard_Real &XMin,
542                                            Standard_Real &YMin,
543                                            Standard_Real &ZMin,
544                                            Standard_Real &XMax,
545                                            Standard_Real &YMax,
546                                            Standard_Real &ZMax) 
547 {
548   Standard_Real x10,y10,z10,x11,y11,z11;
549   Standard_Real x20,y20,z20,x21,y21,z21;
550
551   MyBox1.Get(x10,y10,z10,x11,y11,z11);
552   MyBox2.Get(x20,y20,z20,x21,y21,z21);
553   XMin = 0.;
554   YMin = 0.;
555   ZMin = 0.;
556   XMax = 0.;
557   YMax = 0.;
558   ZMax = 0.;
559
560   if((x10>x21)||(x20>x11)||(y10>y21)||(y20>y11)||(z10>z21)||(z20>z11)) {
561   }
562   else {
563     if(x11<=x21) XMax=x11; else { if(x21<=x11) XMax=x21;}
564     if(x20<=x10) XMin=x10; else { if(x10<=x20) XMin=x20;}
565     if(y11<=y21) YMax=y11; else { if(y21<=y11) YMax=y21;}
566     if(y20<=y10) YMin=y10; else { if(y10<=y20) YMin=y20;}
567     if(z11<=z21) ZMax=z11; else { if(z21<=z11) ZMax=z21;}
568     if(z20<=z10) ZMin=z10; else { if(z10<=z20) ZMin=z20;}
569
570     if(((XMin==XMax)&&(!(YMin==YMax)&&!(ZMin==ZMax)))
571         ||((YMin==YMax)&&(!(XMin==XMax)&&!(ZMin==ZMax)))//ou exclusif ??
572         ||((ZMin==ZMax)&&(!(XMin==XMax)&&!(YMin==YMax)))) {
573     }
574   }
575   //
576   Standard_Real X,Y,Z;
577   X=XMax-XMin; 
578   Y=YMax-YMin; 
579   Z=ZMax-ZMin; 
580   //extension of the box
581   if( (X==0)&&(Y!=0) ) X=Y*0.1;
582   else if( (X==0)&&(Z!=0) ) X=Z*0.1;
583   else X*=0.1;
584
585   if( (Y==0)&&(X!=0) ) Y=X*0.1;
586   else if( (Y==0)&&(Z!=0) ) Y=Z*0.1;
587   else Y*=0.1;
588
589   if( (Z==0)&&(X!=0) ) Z=X*0.1;
590   else if( (Z==0)&&(Y!=0) ) Z=Y*0.1;
591   else Z*=0.1;
592
593
594   if( (X==0)&&(Y==0)&&(Z==0) ) {
595
596
597   }
598   XMin-=X; XMax+=X;
599   YMin-=Y; YMax+=Y;
600   ZMin-=Z; ZMax+=Z;
601
602   //Marking of points included in the common
603   const Standard_Integer FinTP1 = TPoints1.NbItems();
604 //  for(Standard_Integer i=0; i<FinTP1; i++) {
605   Standard_Integer i ;
606   for( i=0; i<FinTP1; i++) {
607     IntPolyh_Point & Pt1 = TPoints1[i];
608     Standard_Integer r;
609     if(Pt1.X()<XMin) { 
610       r=1; 
611     }   
612     else { 
613       if(Pt1.X()>XMax) {
614         r=2;
615       } else {
616         r=0; 
617       } 
618     }
619     if(Pt1.Y()<YMin) { 
620       r|=4; 
621     }  
622     else { 
623       if(Pt1.Y()>YMax) {
624         r|=8; 
625       }
626     }
627     if(Pt1.Z()<ZMin) {
628       r|=16; 
629     } else {
630       if(Pt1.Z()>ZMax) {
631         r|=32;
632       }
633     }
634     Pt1.SetPartOfCommon(r);
635   }
636
637   const Standard_Integer FinTP2 = TPoints2.NbItems();
638   for(Standard_Integer ii=0; ii<FinTP2; ii++) {
639     IntPolyh_Point & Pt2 = TPoints2[ii];
640     Standard_Integer rr;
641     if(Pt2.X()<XMin) { 
642       rr=1;
643     }   
644     else { 
645       if(Pt2.X()>XMax) {
646         rr=2;
647       } else {
648         rr=0;
649       } 
650     }
651     if(Pt2.Y()<YMin) {
652       rr|=4; 
653     }  
654     else { 
655       if(Pt2.Y()>YMax) {
656         rr|=8; 
657       }
658     }
659     if(Pt2.Z()<ZMin) {
660       rr|=16;
661     }
662     else {
663       if(Pt2.Z()>ZMax) {
664         rr|=32;
665       }
666     }
667     Pt2.SetPartOfCommon(rr);
668   }
669 }
670 //=======================================================================
671 //function : FillArrayOfEdges
672 //purpose  : Compute edges from the array of points
673 //           FILL THE ARRAY OF EDGES
674 //=======================================================================
675 void IntPolyh_MaillageAffinage::FillArrayOfEdges
676   (const Standard_Integer SurfID)
677 {
678
679   IntPolyh_ArrayOfEdges &TEdges=(SurfID==1)? TEdges1:TEdges2;
680   IntPolyh_ArrayOfTriangles &TTriangles=(SurfID==1)? TTriangles1:TTriangles2;
681   Standard_Integer NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
682   Standard_Integer NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
683
684   //NbEdges = 3 + 3*(NbSamplesV-2) + 3*(NbSamplesU-2) + 
685   //        + 3*(NbSamplesU-2)*(NbSamplesV-2) + (NbSamplesV-1) + (NbSamplesU-1);
686   //NbSamplesU and NbSamples cannot be less than 2, so
687   Standard_Integer NbEdges = 3*NbSamplesU*NbSamplesV - 2*(NbSamplesU+NbSamplesV) + 1;
688   TEdges.Init(NbEdges);
689
690   Standard_Integer CpteurTabEdges=0;
691
692   //maillage u0 v0
693   TEdges[CpteurTabEdges].SetFirstPoint(0);                // U V
694   TEdges[CpteurTabEdges].SetSecondPoint(1);             // U V+1
695   TEdges[CpteurTabEdges].SetSecondTriangle(0);
696   TTriangles[0].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
697   CpteurTabEdges++;
698
699   TEdges[CpteurTabEdges].SetFirstPoint(0);                // U V
700   TEdges[CpteurTabEdges].SetSecondPoint(NbSamplesV);    // U+1 V
701   TEdges[CpteurTabEdges].SetFirstTriangle(1);
702   TTriangles[1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
703   CpteurTabEdges++;
704   
705   TEdges[CpteurTabEdges].SetFirstPoint(0);                // U V
706   TEdges[CpteurTabEdges].SetSecondPoint(NbSamplesV+1);  // U+1 V+1
707   TEdges[CpteurTabEdges].SetFirstTriangle(0);
708   TTriangles[0].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
709   TEdges[CpteurTabEdges].SetSecondTriangle(1);
710   TTriangles[1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
711   CpteurTabEdges++;
712   
713   //maillage surU=u0
714   Standard_Integer PntInit=1;
715   Standard_Integer BoucleMeshV;
716   for(BoucleMeshV=1; BoucleMeshV<NbSamplesV-1;BoucleMeshV++){
717     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
718     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1);             // U V+1
719     //    TEdges[CpteurTabEdges].SetFirstTriangle(-1);
720     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2);
721     TTriangles[BoucleMeshV*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
722     CpteurTabEdges++;
723     
724     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
725     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV+1);    // U+1 V+1
726     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*2);
727     TTriangles[BoucleMeshV*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
728     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2+1);
729     TTriangles[BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
730     CpteurTabEdges++;
731     
732     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
733     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV);  // U+1 V
734     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*2+1);
735     TTriangles[BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
736     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2-2);
737     TTriangles[BoucleMeshV*2-2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
738     CpteurTabEdges++;
739     PntInit++;
740   }  
741         
742   //maillage sur V=v0
743   PntInit=NbSamplesV;
744   for(BoucleMeshV=1; BoucleMeshV<NbSamplesU-1;BoucleMeshV++){      
745     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);    // U V
746     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1); // U V+1
747     TEdges[CpteurTabEdges].SetFirstTriangle((BoucleMeshV-1)*(NbSamplesV-1)*2+1);
748     TTriangles[(BoucleMeshV-1)*(NbSamplesV-1)*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
749     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*(NbSamplesV-1)*2);
750     TTriangles[BoucleMeshV*(NbSamplesV-1)*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
751     CpteurTabEdges++;
752     
753     TEdges[CpteurTabEdges].SetFirstPoint(PntInit); // U V
754     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV+1);     // U+1 V+1
755     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*(NbSamplesV-1)*2);
756     TTriangles[BoucleMeshV*(NbSamplesV-1)*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
757     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*(NbSamplesV-1)*2+1);
758     TTriangles[BoucleMeshV*(NbSamplesV-1)*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
759     CpteurTabEdges++;
760     
761     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);  // U V
762     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV);    // U+1 V
763     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*(NbSamplesV-1)*2+1);
764     TTriangles[BoucleMeshV*(NbSamplesV-1)*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
765     CpteurTabEdges++;
766     PntInit+=NbSamplesV;
767   }
768       
769   PntInit=NbSamplesV+1;
770   //To provide recursion I associate a point with three edges  
771   for(Standard_Integer BoucleMeshU=1; BoucleMeshU<NbSamplesU-1; BoucleMeshU++){
772     for(BoucleMeshV=1; BoucleMeshV<NbSamplesV-1;BoucleMeshV++){
773       TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
774       TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1);             // U V+1
775       TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesV-1)*2*(BoucleMeshU-1)+BoucleMeshV*2+1);
776       TTriangles[(NbSamplesV-1)*2*(BoucleMeshU-1)+BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
777       TEdges[CpteurTabEdges].SetSecondTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2);
778       TTriangles[(NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
779       CpteurTabEdges++;
780       
781       TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
782       TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV+1);    // U+1 V+1 
783       TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2);
784       TTriangles[(NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
785       TEdges[CpteurTabEdges].SetSecondTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2+1);
786       TTriangles[(NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
787       CpteurTabEdges++;
788       
789       TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
790       TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV);  // U+1 V
791       TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2+1);
792       TTriangles[(NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
793       TEdges[CpteurTabEdges].SetSecondTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2-2);
794       TTriangles[(NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2-2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
795       CpteurTabEdges++;            
796       PntInit++;//Pass to the next point
797     }
798     PntInit++;//Pass the last point of the column
799     PntInit++;//Pass the first point of the next column
800   }
801         
802   //close mesh on U=u1
803   PntInit=(NbSamplesU-1)*NbSamplesV; //point U=u1 V=0
804   for(BoucleMeshV=0; BoucleMeshV<NbSamplesV-1; BoucleMeshV++){
805     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);           //U=u1 V
806     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1);        //U=u1 V+1
807     TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesU-2)*(NbSamplesV-1)*2+BoucleMeshV*2+1);
808     TTriangles[(NbSamplesU-2)*(NbSamplesV-1)*2+BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
809     CpteurTabEdges++;
810     PntInit++;
811   }
812   
813   //close mesh on V=v1
814   for(BoucleMeshV=0; BoucleMeshV<NbSamplesU-1;BoucleMeshV++){      
815     TEdges[CpteurTabEdges].SetFirstPoint(NbSamplesV-1+BoucleMeshV*NbSamplesV);       // U V=v1
816     TEdges[CpteurTabEdges].SetSecondPoint(NbSamplesV-1+(BoucleMeshV+1)*NbSamplesV);  //U+1 V=v1
817     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2*(NbSamplesV-1)+(NbSamplesV-2)*2);
818     TTriangles[BoucleMeshV*2*(NbSamplesV-1)+(NbSamplesV-2)*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
819     CpteurTabEdges++;
820   }
821   TEdges.SetNbItems(CpteurTabEdges);
822 }
823
824 //=======================================================================
825 //function : FillArrayOfTriangles
826 //purpose  : Compute triangles from the array of points, and --
827 //           mark the triangles  that use marked points by the
828 //           CommonBox function.
829 //           FILL THE ARRAY OF TRIANGLES
830 //=======================================================================
831 void IntPolyh_MaillageAffinage::FillArrayOfTriangles
832   (const Standard_Integer SurfID)
833 {
834   Standard_Integer CpteurTabTriangles=0;
835   Standard_Integer PntInit=0;
836
837   IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
838   IntPolyh_ArrayOfTriangles &TTriangles=(SurfID==1)? TTriangles1:TTriangles2;
839   Standard_Integer NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
840   Standard_Integer NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
841
842   TTriangles.Init(2*(NbSamplesU-1)*(NbSamplesV-1));
843   //To provide recursion, I associate a point with two triangles  
844   for(Standard_Integer BoucleMeshU=0; BoucleMeshU<NbSamplesU-1; BoucleMeshU++){
845     for(Standard_Integer BoucleMeshV=0; BoucleMeshV<NbSamplesV-1;BoucleMeshV++){
846       
847       // FIRST TRIANGLE
848       TTriangles[CpteurTabTriangles].SetFirstPoint(PntInit);               // U V
849       TTriangles[CpteurTabTriangles].SetSecondPoint(PntInit+1);            // U V+1
850       TTriangles[CpteurTabTriangles].SetThirdPoint(PntInit+NbSamplesV+1); // U+1 V+1
851
852       // IF ITS EDGE CONTACTS WITH THE COMMON BOX IP REMAINS = A 1
853       if( ( (TPoints[PntInit].PartOfCommon()) & (TPoints[PntInit+1].PartOfCommon()) )
854          &&( (TPoints[PntInit+1].PartOfCommon()) & (TPoints[PntInit+NbSamplesV+1].PartOfCommon()))
855          &&( (TPoints[PntInit+NbSamplesV+1].PartOfCommon()) & (TPoints[PntInit].PartOfCommon())) ) 
856         //IF NOT IP=0
857         TTriangles[CpteurTabTriangles].SetIntersectionPossible(Standard_False);
858
859       CpteurTabTriangles++;
860
861       //SECOND TRIANGLE
862       TTriangles[CpteurTabTriangles].SetFirstPoint(PntInit);               // U V
863       TTriangles[CpteurTabTriangles].SetSecondPoint(PntInit+NbSamplesV+1); // U+1 V+1
864       TTriangles[CpteurTabTriangles].SetThirdPoint(PntInit+NbSamplesV);    // U+1 V 
865
866
867       if( ( (TPoints[PntInit].PartOfCommon()) & (TPoints[PntInit+NbSamplesV+1].PartOfCommon()) )
868          &&( (TPoints[PntInit+NbSamplesV+1].PartOfCommon()) & (TPoints[PntInit+NbSamplesV].PartOfCommon()))
869          &&( (TPoints[PntInit+NbSamplesV].PartOfCommon()) & (TPoints[PntInit].PartOfCommon())) ) 
870         TTriangles[CpteurTabTriangles].SetIntersectionPossible(Standard_False);
871
872
873       CpteurTabTriangles++;
874
875       PntInit++;//Pass to the next point
876     }
877     PntInit++;//Pass the last point of the column
878   }
879   TTriangles.SetNbItems(CpteurTabTriangles);
880   const Standard_Integer FinTT = TTriangles.NbItems();
881   if (FinTT==0) {
882   }
883 }
884 //=======================================================================
885 //function : CommonPartRefinement
886 //purpose  : Refine systematicaly all marked triangles of both surfaces
887 //           REFINING OF THE COMMON
888 //=======================================================================
889 void IntPolyh_MaillageAffinage::CommonPartRefinement() 
890 {
891   Standard_Integer FinInit1 = TTriangles1.NbItems();
892   for(Standard_Integer i=0; i<FinInit1; i++) {
893     if(TTriangles1[i].IsIntersectionPossible())
894       TTriangles1[i].MiddleRefinement(i,MaSurface1,TPoints1,TTriangles1,TEdges1);
895   }
896
897   Standard_Integer FinInit2=TTriangles2.NbItems();
898   for(Standard_Integer ii=0; ii<FinInit2; ii++) {
899     if(TTriangles2[ii].IsIntersectionPossible())
900       TTriangles2[ii].MiddleRefinement(ii,MaSurface2,TPoints2,TTriangles2,TEdges2); 
901   }
902
903 }
904 //=======================================================================
905 //function : LocalSurfaceRefinement
906 //purpose  : Refine systematicaly all marked triangles of ONE surface
907 //=======================================================================
908 void IntPolyh_MaillageAffinage::LocalSurfaceRefinement(const Standard_Integer SurfID) {
909 //refine locally, but systematically the chosen surface
910   if (SurfID==1) {
911     const Standard_Integer FinInit1 = TTriangles1.NbItems();
912     for(Standard_Integer i=0; i<FinInit1; i++) {
913       if(TTriangles1[i].IsIntersectionPossible())
914         TTriangles1[i].MiddleRefinement(i,MaSurface1,TPoints1,TTriangles1,TEdges1);
915     }
916   }
917   //
918   if (SurfID==2) {
919     const Standard_Integer FinInit2 = TTriangles2.NbItems();
920     for(Standard_Integer ii=0; ii<FinInit2; ii++) {
921       if(TTriangles2[ii].IsIntersectionPossible()) 
922         TTriangles2[ii].MiddleRefinement(ii,MaSurface2,TPoints2,TTriangles2,TEdges2); 
923     }
924   }
925 }
926 //=======================================================================
927 //function : ComputeDeflections
928 //purpose  : Compute deflection  for   all  triangles  of  one
929 //           surface,and sort min and max of deflections
930 //           REFINING PART
931 //           Calculation of the deflection of all triangles
932 //           --> deflection max
933 //           --> deflection min
934 //=======================================================================
935 void IntPolyh_MaillageAffinage::ComputeDeflections
936   (const Standard_Integer SurfID)
937 {
938   Handle(Adaptor3d_HSurface) aSurface=(SurfID==1)? MaSurface1:MaSurface2;
939   IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
940   IntPolyh_ArrayOfTriangles &TTriangles=(SurfID==1)? TTriangles1:TTriangles2;
941   Standard_Real &FlecheMin=(SurfID==1)? FlecheMin1:FlecheMin2;
942   Standard_Real &FlecheMoy=(SurfID==1)? FlecheMoy1:FlecheMoy2;
943   Standard_Real &FlecheMax=(SurfID==1)? FlecheMax1:FlecheMax2;
944
945   FlecheMax=-RealLast();
946   FlecheMin=RealLast();
947   FlecheMoy=0.0;
948   const Standard_Integer FinTT = TTriangles.NbItems();
949   
950   for(Standard_Integer i = 0; i < FinTT; i++) {
951     IntPolyh_Triangle& aTriangle = TTriangles[i];
952     Standard_Real Fleche = aTriangle.ComputeDeflection(aSurface, TPoints);
953     if (Fleche > FlecheMax)
954       FlecheMax = Fleche;
955     if (Fleche < FlecheMin)
956       FlecheMin = Fleche;
957   }
958 }
959
960 //=======================================================================
961 //function : TrianglesDeflectionsRefinement
962 //purpose  : Refinement of the triangles depending on the deflection
963 //=======================================================================
964 static
965   void TrianglesDeflectionsRefinement(const Handle(Adaptor3d_HSurface)& theS1,
966                                       IntPolyh_ArrayOfTriangles& theTriangles1,
967                                       IntPolyh_ArrayOfEdges& theEdges1,
968                                       IntPolyh_ArrayOfPoints& thePoints1,
969                                       const Standard_Real theFlecheCritique1,
970                                       const Handle(Adaptor3d_HSurface)& theS2,
971                                       IntPolyh_ArrayOfTriangles& theTriangles2,
972                                       IntPolyh_ArrayOfEdges& theEdges2,
973                                       IntPolyh_ArrayOfPoints& thePoints2,
974                                       const Standard_Real theFlecheCritique2)
975 {
976   // Find the intersecting triangles
977   IntPolyh_IndexedDataMapOfIntegerArrayOfInteger aDMILI;
978   GetInterferingTriangles(theTriangles1, thePoints1, theTriangles2, thePoints2, aDMILI);
979   //
980   // Interfering triangles of second surface
981   TColStd_MapOfInteger aMIntS2;
982   // Since the number of the triangles may change during refinement,
983   // save the number of triangles before refinement
984   Standard_Integer FinTT1 = theTriangles1.NbItems();
985   Standard_Integer FinTT2 = theTriangles2.NbItems();
986   //
987   // Analyze interfering triangles
988   for (Standard_Integer i_S1 = 0; i_S1 < FinTT1; i_S1++) {
989     IntPolyh_Triangle& aTriangle1 = theTriangles1[i_S1];
990     if (!aTriangle1.IsIntersectionPossible()) {
991       continue;
992     }
993     //
994     const IntPolyh_ArrayOfInteger *pLI = aDMILI.Seek(i_S1);
995     if (!pLI || !pLI->Length()) {
996       // Mark non-interfering triangles of S1 to avoid their repeated usage
997       aTriangle1.SetIntersectionPossible(Standard_False);
998       continue;
999     }
1000     //
1001     if (aTriangle1.Deflection() > theFlecheCritique1) {
1002       aTriangle1.MiddleRefinement(i_S1, theS1, thePoints1, theTriangles1, theEdges1);
1003     }
1004     //
1005     IntPolyh_ArrayOfInteger::Iterator Iter(*pLI);
1006     for (; Iter.More(); Iter.Next()) {
1007       Standard_Integer i_S2 = Iter.Value();
1008       if (aMIntS2.Add(i_S2)) {
1009         IntPolyh_Triangle & aTriangle2 = theTriangles2[i_S2];
1010         if (aTriangle2.Deflection() > theFlecheCritique2) {
1011           // Refinement of the larger triangles
1012           aTriangle2.MiddleRefinement(i_S2, theS2, thePoints2, theTriangles2, theEdges2);
1013         }
1014       }
1015     }
1016   }
1017   //
1018   // Mark non-interfering triangles of S2 to avoid their repeated usage
1019   if (aMIntS2.Extent() < FinTT2) {
1020     for (Standard_Integer i_S2 = 0; i_S2 < FinTT2; i_S2++) {
1021       if (!aMIntS2.Contains(i_S2)) {
1022         theTriangles2[i_S2].SetIntersectionPossible(Standard_False);
1023       }
1024     }
1025   }
1026 }
1027 //=======================================================================
1028 //function : LargeTrianglesDeflectionsRefinement
1029 //purpose  : Refinement of the large triangles in case one surface is
1030 //           much smaller then the other.
1031 //=======================================================================
1032 static
1033   void LargeTrianglesDeflectionsRefinement(const Handle(Adaptor3d_HSurface)& theS,
1034                                            IntPolyh_ArrayOfTriangles& theTriangles,
1035                                            IntPolyh_ArrayOfEdges& theEdges,
1036                                            IntPolyh_ArrayOfPoints& thePoints,
1037                                            const Bnd_Box& theOppositeBox)
1038 {
1039   // Find all triangles of the bigger surface with bounding boxes
1040   // overlapping the bounding box the other surface
1041   TColStd_ListOfInteger aLIT;
1042   Standard_Integer i, aNbT = theTriangles.NbItems();
1043   for (i = 0; i < aNbT; ++i) {
1044     IntPolyh_Triangle& aTriangle = theTriangles[i];
1045     if (!aTriangle.IsIntersectionPossible() || aTriangle.IsDegenerated()) {
1046       continue;
1047     }
1048     //
1049     const Bnd_Box& aBox = aTriangle.BoundingBox(thePoints);
1050     if (aBox.IsVoid()) {
1051       continue;
1052     }
1053     //
1054     if (aBox.IsOut(theOppositeBox)) {
1055       aTriangle.SetIntersectionPossible(Standard_False);
1056       continue;
1057     }
1058     //
1059     aLIT.Append(i);
1060   }
1061   //
1062   if (aLIT.IsEmpty()) {
1063     return;
1064   }
1065   //
1066   // One surface is very small relatively to the other.
1067   // The criterion of refining for large surface depends on the size of
1068   // the bounding box of the other - since the criterion should be minimized,
1069   // the smallest side of the bounding box is taken
1070   Standard_Real x0, y0, z0, x1, y1, z1;
1071   theOppositeBox.Get(x0, y0, z0, x1, y1, z1);
1072   Standard_Real dx = Abs(x1 - x0);
1073   Standard_Real dy = Abs(y1 - y0);
1074   Standard_Real diag = Abs(z1 - z0);
1075   Standard_Real dd = (dx > dy) ? dy : dx;
1076   if (diag > dd) diag = dd;
1077
1078   // calculation of the criterion of refining
1079   Standard_Real CritereAffinage = 0.0;
1080   Standard_Real DiagPonderation = 0.5;
1081   CritereAffinage = diag*DiagPonderation;
1082   //
1083   // The deflection of the greater is compared to the size of the smaller
1084   TColStd_ListIteratorOfListOfInteger Iter(aLIT);
1085   for (; Iter.More(); Iter.Next()) {
1086     i = Iter.Value();
1087     IntPolyh_Triangle & aTriangle = theTriangles[i];
1088     if (aTriangle.Deflection() > CritereAffinage) {
1089       aTriangle.MultipleMiddleRefinement(CritereAffinage, theOppositeBox, i,
1090                                          theS, thePoints, theTriangles, theEdges);
1091     }
1092     else {
1093       aTriangle.MiddleRefinement(i, theS, thePoints, theTriangles, theEdges);
1094     }
1095   }
1096 }
1097 //=======================================================================
1098 //function : TrianglesDeflectionsRefinementBSB
1099 //purpose  : Refine both surfaces using UBTree as rejection.
1100 //           The criterion used to refine a triangle are:
1101 //           - The deflection;
1102 //           - The  size of the - bounding boxes (one surface may be
1103 //           very small compared to the other).
1104 //=======================================================================
1105 void IntPolyh_MaillageAffinage::TrianglesDeflectionsRefinementBSB()
1106 {
1107   // To estimate a surface in general it can be interesting
1108   // to calculate all deflections
1109   ComputeDeflections(1);
1110   // Check deflection at output
1111   Standard_Real FlecheCritique1;
1112   if (FlecheMin1 > FlecheMax1) {
1113     return;
1114   }
1115   else {
1116     // fleche min + (flechemax-flechemin) * 80/100
1117     FlecheCritique1 = FlecheMin1*0.2 + FlecheMax1*0.8;
1118   }
1119
1120   // Compute deflections for the second surface
1121   ComputeDeflections(2);
1122
1123   //-- Check arrows at output
1124   Standard_Real FlecheCritique2;
1125   if (FlecheMin2 > FlecheMax2) {
1126     return;
1127   }
1128   else {
1129     //fleche min + (flechemax-flechemin) * 80/100
1130     FlecheCritique2 = FlecheMin2*0.2 + FlecheMax2*0.8;
1131   }
1132
1133   // The greatest of two bounding boxes created in FillArrayOfPoints is found.
1134   // Then this value is weighted depending on the discretization 
1135   // (NbSamplesU and NbSamplesV)
1136   Standard_Real diag1, diag2;
1137   Standard_Real x0, y0, z0, x1, y1, z1;
1138
1139   MyBox1.Get(x0, y0, z0, x1, y1, z1);
1140   x0 -= x1; y0 -= y1; z0 -= z1;
1141   diag1 = x0*x0 + y0*y0 + z0*z0;
1142   const Standard_Real NbSamplesUV1 = Standard_Real(NbSamplesU1) * Standard_Real(NbSamplesV1);
1143   diag1 /= NbSamplesUV1;
1144
1145   MyBox2.Get(x0, y0, z0, x1, y1, z1);
1146   x0 -= x1; y0 -= y1; z0 -= z1;
1147   diag2 = x0*x0 + y0*y0 + z0*z0;
1148   const Standard_Real NbSamplesUV2 = Standard_Real(NbSamplesU2) * Standard_Real(NbSamplesV2);
1149   diag2 /= NbSamplesUV2;
1150
1151   // The surface with the greatest bounding box is "discretized"
1152   if (diag1 < diag2) {
1153     // second is discretized
1154     if (FlecheCritique2 < diag1) {
1155       // The corresponding sizes are not too disproportional
1156       TrianglesDeflectionsRefinement(MaSurface1, TTriangles1, TEdges1, TPoints1, FlecheCritique1,
1157                                      MaSurface2, TTriangles2, TEdges2, TPoints2, FlecheCritique2);
1158     }
1159     else {
1160       // second surface is much larger then the first
1161       LargeTrianglesDeflectionsRefinement(MaSurface2, TTriangles2, TEdges2, TPoints2, MyBox1);
1162     }
1163   }
1164   else {
1165     // first is discretized
1166     if (FlecheCritique1 < diag2) {
1167       // The corresponding sizes are not too disproportional
1168       TrianglesDeflectionsRefinement(MaSurface2, TTriangles2, TEdges2, TPoints2, FlecheCritique2,
1169                                      MaSurface1, TTriangles1, TEdges1, TPoints1, FlecheCritique1);
1170     }
1171     else {
1172       // first surface is much larger then the second
1173       LargeTrianglesDeflectionsRefinement(MaSurface1, TTriangles1, TEdges1, TPoints1, MyBox2);
1174     }
1175   }
1176 }
1177 //=======================================================================
1178 //function : maxSR
1179 //purpose  : This function is used for the function project6
1180 //=======================================================================
1181 inline Standard_Real maxSR(const Standard_Real a,
1182                            const Standard_Real b,
1183                            const Standard_Real c)
1184 {
1185   Standard_Real t = a;
1186   if (b > t) t = b;
1187   if (c > t) t = c;
1188   return t;
1189 }
1190 //=======================================================================
1191 //function : minSR
1192 //purpose  : This function is used for the function project6
1193 //=======================================================================
1194 inline Standard_Real minSR(const Standard_Real a,
1195                            const Standard_Real b,
1196                            const Standard_Real c)
1197 {
1198   Standard_Real t = a;
1199   if (b < t) t = b;
1200   if (c < t) t = c;
1201   return t;
1202 }
1203 //=======================================================================
1204 //function : project6
1205 //purpose  : This function is used for the function TriContact
1206 //=======================================================================
1207 Standard_Integer project6(const IntPolyh_Point &ax, 
1208                           const IntPolyh_Point &p1,
1209                           const IntPolyh_Point &p2, 
1210                           const IntPolyh_Point &p3, 
1211                           const IntPolyh_Point &q1,
1212                           const IntPolyh_Point &q2, 
1213                           const IntPolyh_Point &q3)
1214 {
1215   Standard_Real P1 = ax.Dot(p1);
1216   Standard_Real P2 = ax.Dot(p2);
1217   Standard_Real P3 = ax.Dot(p3);
1218   Standard_Real Q1 = ax.Dot(q1);
1219   Standard_Real Q2 = ax.Dot(q2);
1220   Standard_Real Q3 = ax.Dot(q3);
1221   
1222   Standard_Real mx1 = maxSR(P1, P2, P3);
1223   Standard_Real mn1 = minSR(P1, P2, P3);
1224   Standard_Real mx2 = maxSR(Q1, Q2, Q3);
1225   Standard_Real mn2 = minSR(Q1, Q2, Q3);
1226   
1227   if (mn1 > mx2) return 0;
1228   if (mn2 > mx1) return 0;
1229   return 1;
1230 }
1231 //=======================================================================
1232 //function : TriContact
1233 //purpose  : This fonction     Check if two triangles   are  in
1234 //           contact or no,  return 1 if yes, return 0
1235 //           if no.
1236 //=======================================================================
1237 Standard_Integer IntPolyh_MaillageAffinage::TriContact
1238   (const IntPolyh_Point &P1,
1239    const IntPolyh_Point &P2,
1240    const IntPolyh_Point &P3,
1241    const IntPolyh_Point &Q1,
1242    const IntPolyh_Point &Q2,
1243    const IntPolyh_Point &Q3,
1244    Standard_Real &Angle) const
1245 {
1246   /**
1247      The first triangle is (p1,p2,p3).  The other is (q1,q2,q3).
1248      The edges are (e1,e2,e3) and (f1,f2,f3).
1249      The normals are n1 and m1
1250      Outwards are (g1,g2,g3) and (h1,h2,h3).*/
1251   
1252   if(maxSR(P1.X(),P2.X(),P3.X())<minSR(Q1.X(),Q2.X(),Q3.X())) return(0);
1253   if(maxSR(P1.Y(),P2.Y(),P3.Y())<minSR(Q1.Y(),Q2.Y(),Q3.Y())) return(0);
1254   if(maxSR(P1.Z(),P2.Z(),P3.Z())<minSR(Q1.Z(),Q2.Z(),Q3.Z())) return(0);
1255
1256   if(minSR(P1.X(),P2.X(),P3.X())>maxSR(Q1.X(),Q2.X(),Q3.X())) return(0);
1257   if(minSR(P1.Y(),P2.Y(),P3.Y())>maxSR(Q1.Y(),Q2.Y(),Q3.Y())) return(0);
1258   if(minSR(P1.Z(),P2.Z(),P3.Z())>maxSR(Q1.Z(),Q2.Z(),Q3.Z())) return(0);
1259     
1260   IntPolyh_Point p1, p2, p3;
1261   IntPolyh_Point q1, q2, q3;
1262   IntPolyh_Point e1, e2, e3;
1263   IntPolyh_Point f1, f2, f3;
1264   IntPolyh_Point g1, g2, g3;
1265   IntPolyh_Point h1, h2, h3;
1266   IntPolyh_Point n1, m1;
1267   IntPolyh_Point z;
1268
1269   IntPolyh_Point ef11, ef12, ef13;
1270   IntPolyh_Point ef21, ef22, ef23;
1271   IntPolyh_Point ef31, ef32, ef33;
1272
1273   z.SetX(0.0);  z.SetY(0.0);  z.SetZ(0.0);
1274
1275
1276   p1.SetX(P1.X() - P1.X());  p1.SetY(P1.Y() - P1.Y());  p1.SetZ(P1.Z() - P1.Z());
1277   p2.SetX(P2.X() - P1.X());  p2.SetY(P2.Y() - P1.Y());  p2.SetZ(P2.Z() - P1.Z());
1278   p3.SetX(P3.X() - P1.X());  p3.SetY(P3.Y() - P1.Y());  p3.SetZ(P3.Z() - P1.Z());
1279   
1280   q1.SetX(Q1.X() - P1.X());  q1.SetY(Q1.Y() - P1.Y());  q1.SetZ(Q1.Z() - P1.Z());
1281   q2.SetX(Q2.X() - P1.X());  q2.SetY(Q2.Y() - P1.Y());  q2.SetZ(Q2.Z() - P1.Z());
1282   q3.SetX(Q3.X() - P1.X());  q3.SetY(Q3.Y() - P1.Y());  q3.SetZ(Q3.Z() - P1.Z());
1283   
1284   e1.SetX(p2.X() - p1.X());  e1.SetY(p2.Y() - p1.Y());  e1.SetZ(p2.Z() - p1.Z());
1285   e2.SetX(p3.X() - p2.X());  e2.SetY(p3.Y() - p2.Y());  e2.SetZ(p3.Z() - p2.Z());
1286   e3.SetX(p1.X() - p3.X());  e3.SetY(p1.Y() - p3.Y());  e3.SetZ(p1.Z() - p3.Z());
1287
1288   f1.SetX(q2.X() - q1.X());  f1.SetY(q2.Y() - q1.Y());  f1.SetZ(q2.Z() - q1.Z());
1289   f2.SetX(q3.X() - q2.X());  f2.SetY(q3.Y() - q2.Y());  f2.SetZ(q3.Z() - q2.Z());
1290   f3.SetX(q1.X() - q3.X());  f3.SetY(q1.Y() - q3.Y());  f3.SetZ(q1.Z() - q3.Z());
1291   
1292   n1.Cross(e1, e2); //normal to the first triangle
1293   m1.Cross(f1, f2); //normal to the second triangle
1294
1295   g1.Cross(e1, n1); 
1296   g2.Cross(e2, n1);
1297   g3.Cross(e3, n1);
1298   h1.Cross(f1, m1);
1299   h2.Cross(f2, m1);
1300   h3.Cross(f3, m1);
1301
1302   ef11.Cross(e1, f1);
1303   ef12.Cross(e1, f2);
1304   ef13.Cross(e1, f3);
1305   ef21.Cross(e2, f1);
1306   ef22.Cross(e2, f2);
1307   ef23.Cross(e2, f3);
1308   ef31.Cross(e3, f1);
1309   ef32.Cross(e3, f2);
1310   ef33.Cross(e3, f3);
1311   
1312   // Now the testing is done
1313
1314   if (!project6(n1, p1, p2, p3, q1, q2, q3)) return 0; //T2 is not higher or lower than T1
1315   if (!project6(m1, p1, p2, p3, q1, q2, q3)) return 0; //T1 is not higher of lower than T2
1316   
1317   if (!project6(ef11, p1, p2, p3, q1, q2, q3)) return 0; 
1318   if (!project6(ef12, p1, p2, p3, q1, q2, q3)) return 0;
1319   if (!project6(ef13, p1, p2, p3, q1, q2, q3)) return 0;
1320   if (!project6(ef21, p1, p2, p3, q1, q2, q3)) return 0;
1321   if (!project6(ef22, p1, p2, p3, q1, q2, q3)) return 0;
1322   if (!project6(ef23, p1, p2, p3, q1, q2, q3)) return 0;
1323   if (!project6(ef31, p1, p2, p3, q1, q2, q3)) return 0;
1324   if (!project6(ef32, p1, p2, p3, q1, q2, q3)) return 0;
1325   if (!project6(ef33, p1, p2, p3, q1, q2, q3)) return 0;
1326
1327   if (!project6(g1, p1, p2, p3, q1, q2, q3)) return 0; //T2 is outside of T1 in the plane of T1
1328   if (!project6(g2, p1, p2, p3, q1, q2, q3)) return 0; //T2 is outside of T1 in the plane of T1
1329   if (!project6(g3, p1, p2, p3, q1, q2, q3)) return 0; //T2 is outside of T1 in the plane of T1
1330   if (!project6(h1, p1, p2, p3, q1, q2, q3)) return 0; //T1 is outside of T2 in the plane of T2
1331   if (!project6(h2, p1, p2, p3, q1, q2, q3)) return 0; //T1 is outside of T2 in the plane of T2
1332   if (!project6(h3, p1, p2, p3, q1, q2, q3)) return 0; //T1 is outside of T2 in the plane of T2
1333
1334   //Calculation of cosinus angle between two normals
1335   Standard_Real SqModn1=-1.0;
1336   Standard_Real SqModm1=-1.0;
1337   SqModn1=n1.SquareModulus();
1338   if (SqModn1>SquareMyConfusionPrecision){ 
1339     SqModm1=m1.SquareModulus();
1340   }
1341   if (SqModm1>SquareMyConfusionPrecision) {
1342     Angle=(n1.Dot(m1))/(sqrt(SqModn1)*sqrt(SqModm1));
1343   }
1344   return 1;
1345 }
1346 //=======================================================================
1347 //function : TestNbPoints
1348 //purpose  : This function is used by StartingPointsResearch() to control 
1349 //           the  number of points found keep the result in conformity (1 or 2 points) 
1350 //           void TestNbPoints(const Standard_Integer TriSurfID,
1351 //=======================================================================
1352 void TestNbPoints(const Standard_Integer ,
1353                   Standard_Integer &NbPoints,
1354                   Standard_Integer &NbPointsTotal,
1355                   const IntPolyh_StartPoint &Pt1,
1356                   const IntPolyh_StartPoint &Pt2,
1357                   IntPolyh_StartPoint &SP1,
1358                   IntPolyh_StartPoint &SP2)
1359 {
1360   // already checked in TriangleEdgeContact
1361   //  if( (NbPoints==2)&&(Pt1.CheckSameSP(Pt2)) ) NbPoints=1;
1362
1363   if(NbPoints>2) {
1364
1365   }
1366   else {
1367     if ( (NbPoints==1)&&(NbPointsTotal==0) ) { 
1368       SP1=Pt1;
1369       NbPointsTotal=1;
1370     }
1371     else if ( (NbPoints==1)&&(NbPointsTotal==1) ) { 
1372       if(Pt1.CheckSameSP(SP1)!=1) {
1373         SP2=Pt1;
1374         NbPointsTotal=2;
1375       }
1376     }
1377     else if( (NbPoints==1)&&(NbPointsTotal==2) ) {
1378       if ( (SP1.CheckSameSP(Pt1))||(SP2.CheckSameSP(Pt1)) ) 
1379         NbPointsTotal=2;
1380       else NbPointsTotal=3;
1381     }
1382     else if( (NbPoints==2)&&(NbPointsTotal==0) ) {
1383       SP1=Pt1;
1384       SP2=Pt2;
1385       NbPointsTotal=2;
1386     }
1387     else if( (NbPoints==2)&&(NbPointsTotal==1) ) {//there is also Pt1 != Pt2 
1388       if(SP1.CheckSameSP(Pt1)) {
1389         SP2=Pt2;
1390         NbPointsTotal=2;
1391       }
1392       else if (SP1.CheckSameSP(Pt2)) {
1393         SP2=Pt1;
1394         NbPointsTotal=2;
1395       }
1396       else NbPointsTotal=3;///case SP1!=Pt1 && SP1!=Pt2!
1397     }
1398     else if( (NbPoints==2)&&(NbPointsTotal==2) ) {//there is also SP1!=SP2
1399       if( (SP1.CheckSameSP(Pt1))||(SP1.CheckSameSP(Pt2)) ) {
1400         if( (SP2.CheckSameSP(Pt1))||(SP2.CheckSameSP(Pt2)) )
1401           NbPointsTotal=2;
1402         else NbPointsTotal=3;
1403       }
1404       else NbPointsTotal=3;
1405     }
1406   }
1407 }
1408 //=======================================================================
1409 //function : StartingPointsResearch
1410 //purpose  : From  two  triangles compute intersection  points.
1411 //           If I found   more  than two intersection  points
1412 //           it  means that those triangle are coplanar
1413 //=======================================================================
1414 Standard_Integer IntPolyh_MaillageAffinage::StartingPointsResearch
1415   (const Standard_Integer T1,
1416    const Standard_Integer T2,
1417    IntPolyh_StartPoint &SP1, 
1418    IntPolyh_StartPoint &SP2) const 
1419 {
1420   const IntPolyh_Triangle &Tri1=TTriangles1[T1];
1421   const IntPolyh_Triangle &Tri2=TTriangles2[T2];
1422
1423   const IntPolyh_Point &P1=TPoints1[Tri1.FirstPoint()];
1424   const IntPolyh_Point &P2=TPoints1[Tri1.SecondPoint()];
1425   const IntPolyh_Point &P3=TPoints1[Tri1.ThirdPoint()];
1426   const IntPolyh_Point &Q1=TPoints2[Tri2.FirstPoint()];
1427   const IntPolyh_Point &Q2=TPoints2[Tri2.SecondPoint()];
1428   const IntPolyh_Point &Q3=TPoints2[Tri2.ThirdPoint()];
1429
1430
1431
1432   /* The first triangle is (p1,p2,p3).  The other is (q1,q2,q3).
1433      The sides are (e1,e2,e3) and (f1,f2,f3).
1434      The normals are n1 and m1*/
1435
1436   const IntPolyh_Point  e1=P2-P1;
1437   const IntPolyh_Point  e2=P3-P2;
1438   const IntPolyh_Point  e3=P1-P3;
1439
1440   const IntPolyh_Point  f1=Q2-Q1;
1441   const IntPolyh_Point  f2=Q3-Q2;
1442   const IntPolyh_Point  f3=Q1-Q3;
1443   
1444
1445   IntPolyh_Point nn1,mm1;
1446   nn1.Cross(e1, e2); //normal to the first triangle
1447   mm1.Cross(f1, f2); //normal to the second triangle
1448
1449   Standard_Real nn1modulus, mm1modulus;
1450   nn1modulus=sqrt(nn1.SquareModulus());
1451   mm1modulus=sqrt(mm1.SquareModulus());
1452
1453   //-------------------------------------------------
1454   ///calculation of intersections points between triangles
1455   //-------------------------------------------------
1456   Standard_Integer NbPoints=0;
1457   Standard_Integer NbPointsTotal=0;
1458
1459
1460     ///check T1 normal
1461     if(Abs(nn1modulus)<MyConfusionPrecision) {//10.0e-20){
1462
1463     }
1464     else {
1465       const IntPolyh_Point n1=nn1.Divide(nn1modulus);
1466       ///T2 edges with T1
1467       if(NbPointsTotal<3) {
1468         IntPolyh_StartPoint Pt1,Pt2;
1469         NbPoints=TriangleEdgeContact(1,1,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
1470         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1471       }
1472       
1473       if(NbPointsTotal<3) {
1474         IntPolyh_StartPoint Pt1,Pt2;
1475         NbPoints=TriangleEdgeContact(1,2,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
1476         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1477       }
1478       
1479       if(NbPointsTotal<3) {
1480         IntPolyh_StartPoint Pt1,Pt2;
1481         NbPoints=TriangleEdgeContact(1,3,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
1482         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1483       }
1484     }
1485
1486     ///check T2 normal
1487     if(Abs(mm1modulus)<MyConfusionPrecision) {//10.0e-20){
1488
1489     }
1490     else {
1491       const IntPolyh_Point m1=mm1.Divide(mm1modulus);
1492       ///T1 edges with T2
1493       if(NbPointsTotal<3) {
1494         IntPolyh_StartPoint Pt1,Pt2;
1495         NbPoints=TriangleEdgeContact(2,1,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
1496         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1497       }
1498       
1499       if(NbPointsTotal<3) { 
1500         IntPolyh_StartPoint Pt1,Pt2;
1501         NbPoints=TriangleEdgeContact(2,2,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
1502         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1503       }
1504       
1505       if(NbPointsTotal<3) {
1506         IntPolyh_StartPoint Pt1,Pt2;
1507         NbPoints=TriangleEdgeContact(2,3,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
1508         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1509       }
1510     }
1511   //  no need already checked in  TestNbPoints
1512   //  if( (NbPointsTotal==2)&&(SP1.CheckSameSP(SP2)) ) {
1513   //  NbPointsTotal=1;
1514   //SP1.SetCoupleValue(T1,T2);
1515   // }
1516   //  else 
1517   if(NbPointsTotal==2) {
1518     SP1.SetCoupleValue(T1,T2);
1519     SP2.SetCoupleValue(T1,T2);
1520   }
1521   else if(NbPointsTotal==1)
1522     SP1.SetCoupleValue(T1,T2);
1523   else if(NbPointsTotal==3)
1524     SP1.SetCoupleValue(T1,T2);
1525
1526   return (NbPointsTotal);
1527 }
1528 //=======================================================================
1529 //function : NextStartingPointsResearch
1530 //purpose  : from  two triangles  and an intersection   point I
1531 //           seach the other point (if it exist).
1532 //           This function is used by StartPointChain
1533 //=======================================================================
1534 Standard_Integer IntPolyh_MaillageAffinage::NextStartingPointsResearch
1535   (const Standard_Integer T1,
1536    const Standard_Integer T2,
1537    const IntPolyh_StartPoint &SPInit,
1538    IntPolyh_StartPoint &SPNext) const 
1539 {
1540   Standard_Integer NbPointsTotal=0;
1541   Standard_Integer EdgeInit1=SPInit.E1();
1542   Standard_Integer EdgeInit2=SPInit.E2();
1543   if( (T1<0)||(T2<0) ) NbPointsTotal=0;
1544   else {
1545     
1546     const IntPolyh_Triangle &Tri1=TTriangles1[T1];
1547     const IntPolyh_Triangle &Tri2=TTriangles2[T2];
1548
1549     const IntPolyh_Point &P1=TPoints1[Tri1.FirstPoint()];
1550     const IntPolyh_Point &P2=TPoints1[Tri1.SecondPoint()];
1551     const IntPolyh_Point &P3=TPoints1[Tri1.ThirdPoint()];
1552     const IntPolyh_Point &Q1=TPoints2[Tri2.FirstPoint()];
1553     const IntPolyh_Point &Q2=TPoints2[Tri2.SecondPoint()];
1554     const IntPolyh_Point &Q3=TPoints2[Tri2.ThirdPoint()];
1555     
1556   /* The first triangle is (p1,p2,p3).  The other is (q1,q2,q3).
1557      The edges are (e1,e2,e3) and (f1,f2,f3).
1558      The normals are n1 and m1*/
1559
1560     const IntPolyh_Point  e1=P2-P1;
1561     const IntPolyh_Point  e2=P3-P2;
1562     const IntPolyh_Point  e3=P1-P3;
1563     
1564     const IntPolyh_Point  f1=Q2-Q1;
1565     const IntPolyh_Point  f2=Q3-Q2;
1566     const IntPolyh_Point  f3=Q1-Q3;
1567     
1568     IntPolyh_Point nn1,mm1;
1569     nn1.Cross(e1, e2); //normal to the first triangle
1570     mm1.Cross(f1, f2); //normal to the second triangle
1571
1572     Standard_Real nn1modulus, mm1modulus;
1573     nn1modulus=sqrt(nn1.SquareModulus());
1574     mm1modulus=sqrt(mm1.SquareModulus());
1575     
1576     //-------------------------------------------------
1577     ///calculation of intersections points between triangles
1578     //-------------------------------------------------
1579
1580     Standard_Integer NbPoints=0;
1581     IntPolyh_StartPoint SP1,SP2;
1582
1583     ///check T1 normal
1584     if(Abs(nn1modulus)<MyConfusionPrecision) {//10.0e-20){
1585
1586     }
1587     else {
1588       const IntPolyh_Point n1=nn1.Divide(nn1modulus);
1589       ///T2 edges with T1
1590       if( (NbPointsTotal<3)&&(EdgeInit2!=Tri2.FirstEdge()) ) {
1591         IntPolyh_StartPoint Pt1,Pt2;
1592         NbPoints=TriangleEdgeContact(1,1,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
1593         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1594       }
1595       
1596       if( (NbPointsTotal<3)&&(EdgeInit2!=Tri2.SecondEdge()) ) {
1597         IntPolyh_StartPoint Pt1,Pt2;
1598         NbPoints=TriangleEdgeContact(1,2,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
1599         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1600       }
1601       
1602       if( (NbPointsTotal<3)&&(EdgeInit2!=Tri2.ThirdEdge()) ) {
1603         IntPolyh_StartPoint Pt1,Pt2;
1604         NbPoints=TriangleEdgeContact(1,3,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
1605         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1606       }
1607     }
1608     ///check T2 normal
1609     if(Abs(mm1modulus)<MyConfusionPrecision) {//10.0e-20){
1610
1611     }
1612     else {
1613       const IntPolyh_Point m1=mm1.Divide(mm1modulus);
1614       ///T1 edges with T2
1615       if( (NbPointsTotal<3)&&(EdgeInit1!=Tri1.FirstEdge()) ) {
1616         IntPolyh_StartPoint Pt1,Pt2;
1617         NbPoints=TriangleEdgeContact(2,1,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
1618         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1619       }
1620       
1621       if( (NbPointsTotal<3)&&(EdgeInit1!=Tri1.SecondEdge()) ) {
1622         IntPolyh_StartPoint Pt1,Pt2;
1623         NbPoints=TriangleEdgeContact(2,2,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
1624         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1625       }
1626       
1627       if( (NbPointsTotal<3)&&(EdgeInit1!=Tri1.ThirdEdge()) ) {
1628         IntPolyh_StartPoint Pt1,Pt2;
1629         NbPoints=TriangleEdgeContact(2,3,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
1630         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1631       }
1632     }
1633       
1634     if (NbPointsTotal==1) {
1635       if(SP1.CheckSameSP(SPInit))
1636         NbPointsTotal=0;
1637       else {
1638         SPNext=SP1;
1639       }
1640     }
1641     else if( (NbPointsTotal==2)&&(SP1.CheckSameSP(SPInit)) ) {
1642       NbPointsTotal=1;//SP1 et SPInit sont identiques
1643       SPNext=SP2;
1644     }
1645     else if( (NbPointsTotal==2)&&(SP2.CheckSameSP(SPInit)) ) {
1646       NbPointsTotal=1;//SP2 et SPInit sont identiques
1647       SPNext=SP1;
1648     }
1649
1650     else if(NbPointsTotal>1) {
1651
1652     }
1653   }
1654   SPNext.SetCoupleValue(T1,T2);
1655   return (NbPointsTotal);
1656 }
1657 //=======================================================================
1658 //function : CalculPtsInterTriEdgeCoplanaires
1659 //purpose  : 
1660 //=======================================================================
1661 void CalculPtsInterTriEdgeCoplanaires(const Standard_Integer TriSurfID,
1662                                       const IntPolyh_Point &NormaleTri,
1663                                       const IntPolyh_Triangle &Tri1,
1664                                       const IntPolyh_Triangle &Tri2,
1665                                       const IntPolyh_Point &PE1,
1666                                       const IntPolyh_Point &PE2,
1667                                       const IntPolyh_Point &Edge,
1668                                       const Standard_Integer EdgeIndex,
1669                                       const IntPolyh_Point &PT1,
1670                                       const IntPolyh_Point &PT2,
1671                                       const IntPolyh_Point &Cote,
1672                                       const Standard_Integer CoteIndex,
1673                                       IntPolyh_StartPoint &SP1,
1674                                       IntPolyh_StartPoint &SP2,
1675                                       Standard_Integer &NbPoints)
1676 {
1677   Standard_Real aDE, aDC;
1678   //
1679   gp_Vec aVE(Edge.X(), Edge.Y(), Edge.Z());
1680   gp_Vec aVC(Cote.X(), Cote.Y(), Cote.Z());
1681   //
1682   aDE = aVE.SquareMagnitude();
1683   aDC = aVC.SquareMagnitude();
1684   //
1685   if (aDE > SquareMyConfusionPrecision) {
1686     aVE.Divide(aDE);
1687   }
1688   //
1689   if (aDC > SquareMyConfusionPrecision) {
1690     aVC.Divide(aDC);
1691   }
1692   //
1693   if (!aVE.IsParallel(aVC, MyConfusionPrecision)) {
1694     ///Edge and side are not parallel
1695     IntPolyh_Point Per;
1696     Per.Cross(NormaleTri,Cote);
1697     Standard_Real p1p = Per.Dot(PE1);
1698     Standard_Real p2p = Per.Dot(PE2);
1699     Standard_Real p0p = Per.Dot(PT1);
1700     ///The edge are PT1 are projected on the perpendicular of the side in the plane of the triangle
1701     if ( ( ((p1p>=p0p)&&(p0p>=p2p) )||( (p1p<=p0p)&&(p0p<=p2p) )) && (Abs(p1p-p2p) > MyConfusionPrecision)) {
1702       Standard_Real lambda=(p1p-p0p)/(p1p-p2p);
1703       if (lambda<-MyConfusionPrecision) {
1704
1705       }
1706       IntPolyh_Point PIE;
1707       if (Abs(lambda)<MyConfusionPrecision)//lambda=0
1708         PIE=PE1;
1709       else if (Abs(lambda)>1.0-MyConfusionPrecision)//lambda=1
1710         PIE=PE2;
1711       else
1712         PIE=PE1+Edge*lambda;
1713
1714       Standard_Real alpha=RealLast();
1715       if(Cote.X()!=0) alpha=(PIE.X()-PT1.X())/Cote.X();
1716       else if (Cote.Y()!=0) alpha=(PIE.Y()-PT1.Y())/Cote.Y();
1717       else if (Cote.Z()!=0) alpha=(PIE.Z()-PT1.Z())/Cote.Z();
1718       else {
1719
1720       }
1721       
1722       if (alpha<-MyConfusionPrecision) {
1723
1724       }
1725       else {
1726         if (NbPoints==0) {
1727           SP1.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
1728           if (TriSurfID==1) {
1729             if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
1730               SP1.SetUV1(PT1.U(),PT1.V());
1731               SP1.SetUV1(PIE.U(),PIE.V());
1732               SP1.SetEdge1(-1);
1733             }
1734             if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
1735               SP1.SetUV1(PT2.U(),PT2.V());
1736               SP1.SetUV1(PIE.U(),PIE.V());
1737               SP1.SetEdge1(-1);
1738             }
1739             else {
1740               SP1.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
1741               SP1.SetUV2(PIE.U(),PIE.V());
1742               SP1.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
1743               if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha);
1744               else SP1.SetLambda1(1.0-alpha);
1745             }
1746             NbPoints++;
1747           }
1748           else if (TriSurfID==2) {
1749             if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
1750               SP1.SetUV1(PT1.U(),PT1.V());
1751               SP1.SetUV1(PIE.U(),PIE.V());
1752               SP1.SetEdge2(-1);
1753             }
1754             if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
1755               SP1.SetUV1(PT2.U(),PT2.V());
1756               SP1.SetUV1(PIE.U(),PIE.V());
1757               SP1.SetEdge2(-1);
1758             }
1759             else {
1760               SP1.SetUV1(PIE.U(),PIE.V());
1761               SP1.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
1762               SP1.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
1763               if (Tri2.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda2(alpha);
1764               else SP1.SetLambda2(1.0-alpha);
1765             }
1766             NbPoints++;
1767           }
1768           else {
1769
1770           }
1771         }
1772
1773         else if (NbPoints==1) {
1774           SP2.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
1775           if (TriSurfID==1) {
1776             if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
1777               SP2.SetUV1(PT1.U(),PT1.V());
1778               SP2.SetUV1(PIE.U(),PIE.V());
1779               SP2.SetEdge1(-1);
1780             }
1781             if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
1782               SP2.SetUV1(PT2.U(),PT2.V());
1783               SP2.SetUV1(PIE.U(),PIE.V());
1784               SP2.SetEdge1(-1);
1785             }
1786             else {
1787               SP2.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
1788               SP2.SetUV2(PIE.U(),PIE.V());
1789               SP2.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
1790               if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha);
1791               else SP2.SetLambda1(1.0-alpha);
1792             }
1793             NbPoints++;
1794           }
1795           else if (TriSurfID==2) {
1796             if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
1797               SP2.SetUV1(PT1.U(),PT1.V());
1798               SP2.SetUV1(PIE.U(),PIE.V());
1799               SP2.SetEdge2(-1);
1800             }
1801             if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
1802               SP2.SetUV1(PT2.U(),PT2.V());
1803               SP2.SetUV1(PIE.U(),PIE.V());
1804               SP2.SetEdge2(-1);
1805             }
1806             else {
1807               SP2.SetUV1(PIE.U(),PIE.V());
1808               SP2.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
1809               SP2.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
1810               if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda2(alpha);
1811               else SP2.SetLambda2(1.0-alpha);
1812             }
1813             NbPoints++;
1814           }
1815           else {
1816
1817           }
1818         }
1819
1820         else if( (NbPoints>2)||(NbPoints<0) ) {
1821
1822           }
1823       }
1824     }
1825   }
1826   else {    
1827     //Side and Edge are parallel, with previous 
1828     //rejections they are at the same side
1829     //The points are projected on that side
1830     Standard_Real pe1p= Cote.Dot(PE1);
1831     Standard_Real pe2p= Cote.Dot(PE2);
1832     Standard_Real pt1p= Cote.Dot(PT1);
1833     Standard_Real pt2p= Cote.Dot(PT2);
1834     Standard_Real lambda1=0., lambda2=0., alpha1=0., alpha2=0.;
1835     IntPolyh_Point PEP1,PTP1,PEP2,PTP2;
1836
1837     if (pe1p>pe2p) {
1838       if ( (pt1p<pe1p)&&(pe1p<=pt2p) ) {
1839         lambda1=0.0;
1840         PEP1=PE1;
1841         alpha1=((pe1p-pt1p)/(pt2p-pt1p));
1842         PTP1=PT1+Cote*alpha1;
1843         NbPoints=1;
1844         if (pt1p<=pe2p) {
1845           lambda2=1.0;
1846           PEP2=PE2;
1847           alpha2=((pe2p-pt1p)/(pt2p-pt1p));
1848           PTP2=PT1+Cote*alpha2;
1849           NbPoints=2;
1850         }
1851         else {
1852           lambda2=((pt1p-pe1p)/(pe2p-pe1p));
1853           PEP2=PE1+Edge*lambda2;
1854           alpha2=0.0;
1855           PTP2=PT1;
1856           NbPoints=2;
1857         }
1858       }
1859       else if( (pt2p<pe1p)&&(pe1p<=pt1p) ) {
1860         lambda1=0.0;
1861         PEP1=PE1;
1862         alpha1=((pt1p-pe1p)/(pt1p-pt2p));
1863         PTP1=PT1+Cote*alpha1;
1864         NbPoints=1;
1865         if (pt2p<=pe2p) {
1866           lambda2=1.0;
1867           PEP2=PE2;
1868           alpha2=((pe2p-pt1p)/(pt2p-pt1p));
1869           PTP2=PT1+Cote*alpha2;
1870           NbPoints=2;
1871         }
1872         else {
1873           lambda2=((pt2p-pe1p)/(pe2p-pe1p));
1874           PEP2=PE1+Edge*lambda2;
1875           alpha2=1.0;
1876           PTP2=PT2;
1877           NbPoints=2;
1878         }
1879       }
1880     }
1881     
1882     if (pe1p<pe2p) {
1883       if ( (pt1p<pe2p)&&(pe2p<=pt2p) ) {
1884         lambda1=1.0;
1885         PEP1=PE2;
1886         alpha1=((pe2p-pt1p)/(pt2p-pt1p));
1887         PTP1=PT1+Cote*alpha1;
1888         NbPoints=1;
1889         if (pt1p<=pe1p) {
1890           lambda2=0.0;
1891           PEP2=PE1;
1892           alpha2=((pe1p-pt1p)/(pt2p-pt1p));
1893           PTP2=PT1+Cote*alpha2;
1894           NbPoints=2;
1895         }
1896         else {
1897           lambda2=((pt1p-pe1p)/(pe2p-pe1p));
1898           PEP2=PE2+Edge*lambda2;
1899           alpha2=0.0;
1900           PTP2=PT1;
1901           NbPoints=2;
1902         }
1903       }
1904       else if( (pt2p<pe2p)&&(pe2p<=pt1p) ) {
1905         lambda1=1.0;
1906         PEP1=PE2;
1907         alpha1=((pt1p-pe2p)/(pt1p-pt2p));
1908         PTP1=PT1+Cote*alpha1;
1909         NbPoints=1;
1910         if (pt2p<=pe1p) {
1911           lambda2=0.0;
1912           PEP2=PE1;
1913           alpha2=((pe1p-pt1p)/(pt2p-pt1p));
1914           PTP2=PT1+Cote*alpha2;
1915           NbPoints=2;
1916         }
1917         else {
1918           lambda2=((pt2p-pe1p)/(pe2p-pe1p));
1919           PEP2=PE1+Edge*lambda2;
1920           alpha2=1.0;
1921           PTP2=PT2;
1922           NbPoints=2;
1923         }
1924       }
1925     }
1926   
1927     if (NbPoints!=0) {
1928       SP1.SetXYZ(PEP1.X(),PEP1.Y(),PEP1.Z());
1929       if (TriSurfID==1) {///cote appartient a Tri1
1930         SP1.SetUV1(PTP1.U(),PTP1.V());          
1931         SP1.SetUV2(PEP1.U(),PEP1.V());
1932         SP1.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
1933
1934         if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha1);
1935         else SP1.SetLambda1(1.0-alpha1);
1936         
1937         if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP1.SetLambda2(lambda1);
1938         else SP1.SetLambda2(1.0-lambda1);
1939       }
1940       if (TriSurfID==2) {///cote appartient a Tri2
1941         SP1.SetUV1(PEP1.U(),PTP1.V());          
1942         SP1.SetUV2(PTP1.U(),PEP1.V());
1943         SP1.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
1944
1945         if(Tri2.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha1);
1946         else SP1.SetLambda1(1.0-alpha1);
1947         
1948         if(Tri1.GetEdgeOrientation(EdgeIndex)>0) SP1.SetLambda2(lambda1);
1949         else SP1.SetLambda2(1.0-lambda1);
1950       }
1951       
1952       //It is checked if PEP1!=PEP2
1953       if ( (NbPoints==2)&&(Abs(PEP1.U()-PEP2.U())<MyConfusionPrecision)
1954          &&(Abs(PEP1.V()-PEP2.V())<MyConfusionPrecision) ) NbPoints=1;
1955       if (NbPoints==2) {
1956         SP2.SetXYZ(PEP2.X(),PEP2.Y(),PEP2.Z());
1957         if (TriSurfID==1) {
1958           SP2.SetUV1(PTP2.U(),PTP2.V());          
1959           SP2.SetUV2(PEP2.U(),PEP2.V());
1960           SP2.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
1961
1962           if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha1);
1963           else SP2.SetLambda1(1.0-alpha1);
1964           
1965           if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP2.SetLambda2(lambda1);
1966           else SP2.SetLambda2(1.0-lambda1);
1967         }
1968         if (TriSurfID==2) {
1969           SP2.SetUV1(PEP2.U(),PTP2.V());          
1970           SP2.SetUV2(PTP2.U(),PEP2.V());
1971           SP2.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
1972
1973           if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha1);
1974           else SP2.SetLambda1(1.0-alpha1);
1975           
1976           if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP2.SetLambda2(lambda1);
1977           else SP2.SetLambda2(1.0-lambda1);
1978         } 
1979       }
1980     }
1981   }
1982   //Filter if the point is placed on top, the edge is set  to -1
1983   if (NbPoints>0) {
1984     if(Abs(SP1.Lambda1())<MyConfusionPrecision)
1985       SP1.SetEdge1(-1);
1986     if(Abs(SP1.Lambda1()-1)<MyConfusionPrecision)
1987       SP1.SetEdge1(-1);
1988     if(Abs(SP1.Lambda2())<MyConfusionPrecision)
1989       SP1.SetEdge2(-1);
1990     if(Abs(SP1.Lambda2()-1)<MyConfusionPrecision)
1991       SP1.SetEdge2(-1);
1992   }
1993   if (NbPoints==2) {
1994     if(Abs(SP2.Lambda1())<MyConfusionPrecision)
1995       SP2.SetEdge1(-1);
1996     if(Abs(SP2.Lambda1()-1)<MyConfusionPrecision)
1997       SP2.SetEdge1(-1);
1998     if(Abs(SP2.Lambda2())<MyConfusionPrecision)
1999       SP2.SetEdge2(-1);
2000     if(Abs(SP2.Lambda2()-1)<MyConfusionPrecision)
2001       SP2.SetEdge2(-1);
2002   }
2003 }
2004
2005 //=======================================================================
2006 //function : TriangleEdgeContact
2007 //purpose  : 
2008 //=======================================================================
2009 Standard_Integer IntPolyh_MaillageAffinage::TriangleEdgeContact
2010   (const Standard_Integer TriSurfID,
2011    const Standard_Integer EdgeIndex,
2012    const IntPolyh_Triangle &Tri1,
2013    const IntPolyh_Triangle &Tri2,
2014    const IntPolyh_Point &PT1,
2015    const IntPolyh_Point &PT2,
2016    const IntPolyh_Point &PT3,
2017    const IntPolyh_Point &Cote12,
2018    const IntPolyh_Point &Cote23,
2019    const IntPolyh_Point &Cote31,
2020    const IntPolyh_Point &PE1,const IntPolyh_Point &PE2,
2021    const IntPolyh_Point &Edge,
2022    const IntPolyh_Point &NormaleT,
2023    IntPolyh_StartPoint &SP1,IntPolyh_StartPoint &SP2) const 
2024 {
2025
2026   Standard_Real lambda =0., alpha =0., beta =0.;
2027
2028   //One of edges on which the point is located is known
2029   if (TriSurfID==1) {
2030     SP1.SetEdge2(Tri2.GetEdgeNumber(EdgeIndex));
2031     SP2.SetEdge2(Tri2.GetEdgeNumber(EdgeIndex));
2032   }
2033   else if (TriSurfID==2) {
2034     SP1.SetEdge1(Tri1.GetEdgeNumber(EdgeIndex));
2035     SP2.SetEdge1(Tri1.GetEdgeNumber(EdgeIndex));
2036   }
2037   else {
2038
2039   }
2040
2041   /**The edge is projected on the normal in the triangle if yes 
2042   --> a free intersection (point I)--> Start point is found */
2043   Standard_Integer NbPoints=0;
2044   if(NormaleT.SquareModulus()==0) {
2045
2046   }
2047   else if( (Cote12.SquareModulus()==0)
2048        ||(Cote23.SquareModulus()==0)
2049        ||(Cote31.SquareModulus()==0) ) {
2050
2051   }
2052   else if(Edge.SquareModulus()==0) {
2053
2054   }
2055   else {
2056     Standard_Real pe1 = NormaleT.Dot(PE1);
2057     Standard_Real pe2 = NormaleT.Dot(PE2);
2058     Standard_Real pt1 = NormaleT.Dot(PT1);  
2059     
2060     // PE1I = lambda.Edge
2061    
2062     if( (Abs(pe1-pt1)<MyConfusionPrecision)&&(Abs(pe2-pt1)<MyConfusionPrecision)) {
2063       //edge and triangle are coplanar (two contact points at maximum)
2064
2065
2066       //the tops of the triangle are projected on the perpendicular to the edge 
2067       IntPolyh_Point PerpEdge;
2068       PerpEdge.Cross(NormaleT,Edge);
2069       Standard_Real pp1 = PerpEdge.Dot(PT1);
2070       Standard_Real pp2 = PerpEdge.Dot(PT2);
2071       Standard_Real pp3 = PerpEdge.Dot(PT3);
2072       Standard_Real ppe1 = PerpEdge.Dot(PE1);
2073       
2074
2075       if ( (Abs(pp1-pp2)<MyConfusionPrecision)&&(Abs(pp1-pp3)<MyConfusionPrecision) ) {
2076
2077       }
2078       else {
2079         if ( ( (pp1>=ppe1)&&(pp2<=ppe1)&&(pp3<=ppe1) ) || ( (pp1<=ppe1)&&(pp2>=ppe1)&&(pp3>=ppe1) ) ){
2080           //there are two sides (common top PT1) that can cut the edge
2081           
2082           //first side
2083           CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2084                                            PT1,PT2,Cote12,1,SP1,SP2,NbPoints);
2085           
2086           if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2087               &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2088           
2089           //second side
2090           if (NbPoints<2) CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2091                                                            PT3,PT1,Cote31,3,SP1,SP2,NbPoints);
2092         }
2093         
2094         if ( (NbPoints>1)&&(Abs(SP1.U1()-SP2.U1())<MyConfusionPrecision)
2095             &&(Abs(SP1.V2()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2096         if (NbPoints>=2) return(NbPoints);
2097         
2098         else if ( ( ( (pp2>=ppe1)&&(pp1<=ppe1)&&(pp3<=ppe1) ) || ( (pp2<=ppe1)&&(pp1>=ppe1)&&(pp3>=ppe1) ) )
2099                  && (NbPoints<2) ) {
2100           //there are two sides (common top PT2) that can cut the edge
2101           
2102           //first side
2103           CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2104                                            PT1,PT2,Cote12,1,SP1,SP2,NbPoints);
2105           
2106           if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2107               &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2108           
2109           //second side
2110           if(NbPoints<2) CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2111                                                           PT2,PT3,Cote23,2,SP1,SP2,NbPoints);
2112         }
2113         if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2114             &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2115         if (NbPoints>=2) return(NbPoints);
2116         
2117         else if ( (( (pp3>=ppe1)&&(pp1<=ppe1)&&(pp2<=ppe1) ) || ( (pp3<=ppe1)&&(pp1>=ppe1)&&(pp2>=ppe1) ))
2118                 && (NbPoints<2) ) {
2119           //there are two sides (common top PT3) that can cut the edge
2120           
2121           //first side
2122           CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2123                                            PT3,PT1,Cote31,3,SP1,SP2,NbPoints);
2124           
2125           if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2126               &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2127           
2128           //second side
2129           if (NbPoints<2) CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2130                                                            PT2,PT3,Cote23,2,SP1,SP2,NbPoints);
2131         }
2132         if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2133             &&(Abs(SP2.V1()-SP1.V1())<MyConfusionPrecision) ) NbPoints=1;
2134         if (NbPoints>=2) return(NbPoints);
2135       }
2136     }
2137
2138     //------------------------------------------------------
2139     // NON COPLANAR edge and triangle (a contact point)
2140     //------------------------------------------------------
2141     else if(   ( (pe1>=pt1)&&(pt1>=pe2) ) || ( (pe1<=pt1)&&(pt1<=pe2) )  ) { //
2142       lambda=(pe1-pt1)/(pe1-pe2);
2143       IntPolyh_Point PI;
2144       if (lambda<-MyConfusionPrecision) {
2145
2146       }
2147       else if (Abs(lambda)<MyConfusionPrecision) {//lambda==0
2148         PI=PE1;
2149         if(TriSurfID==1) SP1.SetEdge2(-1);
2150         else SP1.SetEdge1(-1);
2151       }
2152       else if (Abs(lambda-1.0)<MyConfusionPrecision) {//lambda==1
2153         PI=PE2;
2154         if(TriSurfID==1) SP1.SetEdge2(-1);
2155         else SP1.SetEdge1(-1);
2156       }
2157       else {
2158         PI=PE1+Edge*lambda;
2159         if(TriSurfID==1) {
2160           if(Tri2.GetEdgeOrientation(EdgeIndex)>0)
2161             SP1.SetLambda2(lambda);
2162           else SP1.SetLambda2(1.0-lambda);
2163         }
2164         if(TriSurfID==2) {
2165           if(Tri1.GetEdgeOrientation(EdgeIndex)>0)
2166             SP1.SetLambda1(lambda);
2167           else SP1.SetLambda1(1.0-lambda);
2168         }
2169       }
2170         
2171       Standard_Real Cote23X=Cote23.X();
2172       Standard_Real D1=0.0;
2173       Standard_Real D3,D4;
2174
2175       //Combination Eq1 Eq2
2176       if(Abs(Cote23X)>MyConfusionPrecision) {
2177         D1=Cote12.Y()-Cote12.X()*Cote23.Y()/Cote23X;
2178       }
2179       if(Abs(D1)>MyConfusionPrecision) {
2180         alpha = ( PI.Y()-PT1.Y()-(PI.X()-PT1.X())*Cote23.Y()/Cote23X )/D1;
2181
2182         ///It is checked if 1.0>=alpha>=0.0
2183         if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
2184         else beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23X;
2185       }
2186       //Combination Eq1 and Eq2 with Cote23.X()==0
2187       else if ( (Abs(Cote12.X())>MyConfusionPrecision)
2188               &&(Abs(Cote23X)<MyConfusionPrecision) ) { //There is Cote23.X()==0
2189         alpha = (PI.X()-PT1.X())/Cote12.X();
2190         
2191         if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
2192         
2193         else if (Abs(Cote23.Y())>MyConfusionPrecision) beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
2194         else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
2195         else  {
2196
2197         }
2198       }
2199       //Combination Eq1 and Eq3
2200       else if ( (Abs(Cote23.X())>MyConfusionPrecision)
2201               &&(Abs( D3= (Cote12.Z()-Cote12.X()*Cote23.Z()/Cote23.X()) ) > MyConfusionPrecision) ) {
2202         
2203         alpha = (PI.Z()-PT1.Z()-(PI.X()-PT1.X())*Cote23.Z()/Cote23.X())/D3;
2204         
2205         if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
2206         else beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23.X();
2207       }
2208       //Combination Eq2 and Eq3
2209       else if ( (Abs(Cote23.Y())>MyConfusionPrecision)
2210               &&(Abs( D4= (Cote12.Z()-Cote12.Y()*Cote23.Z()/Cote23.Y()) ) > MyConfusionPrecision) ) {
2211         
2212         alpha = (PI.Z()-PT1.Z()-(PI.Y()-PT1.Y())*Cote23.Z()/Cote23.Y())/D4;
2213         
2214         if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
2215         else beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
2216       }
2217       //Combination Eq2 and Eq3 with Cote23.Y()==0
2218       else if ( (Abs(Cote12.Y())>MyConfusionPrecision)
2219              && (Abs(Cote23.Y())<MyConfusionPrecision) ) {
2220         alpha = (PI.Y()-PT1.Y())/Cote12.Y();
2221         
2222         if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
2223         
2224         else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
2225         
2226         else {
2227           printf("\nCote PT2PT3 nul1\n");
2228           PT2.Dump(2004);
2229           PT3.Dump(3004);
2230         }
2231       }
2232       //Combination Eq1 and Eq3 with Cote23.Z()==0
2233       else if ( (Abs(Cote12.Z())>MyConfusionPrecision)
2234              && (Abs(Cote23.Z())<MyConfusionPrecision) ) {
2235         alpha = (PI.Z()-PT1.Z())/Cote12.Z();
2236         
2237         if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
2238         
2239         else if (Abs(Cote23.X())>MyConfusionPrecision) beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23.X();
2240         
2241         else {
2242
2243         }
2244       }
2245       
2246       else { //Particular case not processed ?
2247
2248         alpha=RealLast();
2249         beta=RealLast();
2250       }
2251       
2252       if( (beta<-MyConfusionPrecision)||(beta>(alpha+MyConfusionPrecision)) ) return(0);
2253       else {
2254         SP1.SetXYZ(PI.X(),PI.Y(),PI.Z());
2255
2256         if (TriSurfID==1) {
2257           SP1.SetUV2(PI.U(),PI.V());
2258           SP1.SetUV1(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
2259           NbPoints++;
2260           if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
2261             SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
2262             SP1.SetUV1(PT1.U(),PT1.V());
2263             SP1.SetEdge1(-1);
2264           }
2265           else if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
2266             SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
2267             SP1.SetUV1(PT2.U(),PT2.V());
2268             SP1.SetEdge1(-1);
2269           }
2270           else if ( (Abs(beta-1)<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
2271             SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
2272             SP1.SetUV1(PT3.U(),PT3.V());
2273             SP1.SetEdge1(-1);
2274           }
2275           else if (beta<MyConfusionPrecision) {//beta==0 
2276             SP1.SetEdge1(Tri1.GetEdgeNumber(1));
2277             if(Tri1.GetEdgeOrientation(1)>0)
2278               SP1.SetLambda1(alpha);
2279             else SP1.SetLambda1(1.0-alpha);
2280           }
2281           else if (Abs(beta-alpha)<MyConfusionPrecision) {//beta==alpha  
2282             SP1.SetEdge1(Tri1.GetEdgeNumber(3));
2283             if(Tri1.GetEdgeOrientation(3)>0)
2284               SP1.SetLambda1(1.0-alpha);
2285             else SP1.SetLambda1(alpha);
2286           }
2287           else if (Abs(alpha-1)<MyConfusionPrecision) {//alpha==1
2288             SP1.SetEdge1(Tri1.GetEdgeNumber(2));
2289             if(Tri1.GetEdgeOrientation(2)>0)
2290               SP1.SetLambda1(beta);
2291             else SP1.SetLambda1(1.0-beta);
2292           }
2293         }
2294         else if(TriSurfID==2) {
2295           SP1.SetUV1(PI.U(),PI.V());
2296           SP1.SetUV2(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
2297           NbPoints++;
2298           if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
2299             SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
2300             SP1.SetUV2(PT1.U(),PT1.V());
2301             SP1.SetEdge2(-1);
2302           }
2303           else if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
2304             SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
2305             SP1.SetUV2(PT2.U(),PT2.V());
2306             SP1.SetEdge2(-1);
2307           }
2308           else if ( (Abs(beta-1)<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
2309             SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
2310             SP1.SetUV2(PT3.U(),PT3.V());
2311             SP1.SetEdge2(-1);
2312           }
2313           else if (beta<MyConfusionPrecision) { //beta==0
2314             SP1.SetEdge2(Tri2.GetEdgeNumber(1));
2315             if(Tri2.GetEdgeOrientation(1)>0)
2316               SP1.SetLambda2(alpha);
2317             else SP1.SetLambda2(1.0-alpha);
2318           }
2319           else if (Abs(beta-alpha)<MyConfusionPrecision) {//beta==alpha
2320             SP1.SetEdge2(Tri2.GetEdgeNumber(3));
2321             if(Tri2.GetEdgeOrientation(3)>0)
2322               SP1.SetLambda2(1.0-alpha);
2323             else SP1.SetLambda2(alpha);
2324           }
2325           else if (Abs(alpha-1)<MyConfusionPrecision) {//alpha==1
2326             SP1.SetEdge2(Tri2.GetEdgeNumber(2));        
2327             if(Tri2.GetEdgeOrientation(2)>0)
2328               SP1.SetLambda2(alpha);
2329             else SP1.SetLambda2(1.0-alpha);
2330           }
2331         }
2332         else {
2333
2334         }
2335       }
2336     }
2337     else return 0;
2338   }
2339   return (NbPoints);
2340 }
2341 //=======================================================================
2342 //function : TriangleCompare
2343 //purpose  : Analyze  each couple of  triangles from the two --
2344 //           array  of triangles,  to   see  if they are  in
2345 //           contact,  and  compute the  incidence.  Then  put
2346 //           couples  in contact  in  the  array  of  couples
2347 //=======================================================================
2348 Standard_Integer IntPolyh_MaillageAffinage::TriangleCompare ()
2349 {
2350   // Find couples with interfering bounding boxes
2351   IntPolyh_IndexedDataMapOfIntegerArrayOfInteger aDMILI;
2352   GetInterferingTriangles(TTriangles1, TPoints1,
2353                           TTriangles2, TPoints2,
2354                           aDMILI);
2355   if (aDMILI.IsEmpty()) {
2356     return 0;
2357   }
2358   //
2359   Standard_Real CoupleAngle = -2.0;
2360   //
2361   // Intersection of the triangles
2362   Standard_Integer i, aNb = aDMILI.Extent();
2363   for (i = 1; i <= aNb; ++i) {
2364     const Standard_Integer i_S1 = aDMILI.FindKey(i);
2365     IntPolyh_Triangle &Triangle1 =  TTriangles1[i_S1];
2366     const IntPolyh_Point& P1 = TPoints1[Triangle1.FirstPoint()];
2367     const IntPolyh_Point& P2 = TPoints1[Triangle1.SecondPoint()];
2368     const IntPolyh_Point& P3 = TPoints1[Triangle1.ThirdPoint()];
2369     //
2370     const IntPolyh_ArrayOfInteger& aLI2 = aDMILI(i);
2371     IntPolyh_ArrayOfInteger::Iterator aItLI(aLI2);
2372     for (; aItLI.More(); aItLI.Next()) {
2373       const Standard_Integer i_S2 = aItLI.Value();
2374       IntPolyh_Triangle &Triangle2 =  TTriangles2[i_S2];
2375       const IntPolyh_Point& Q1 = TPoints2[Triangle2.FirstPoint()];
2376       const IntPolyh_Point& Q2 = TPoints2[Triangle2.SecondPoint()];
2377       const IntPolyh_Point& Q3 = TPoints2[Triangle2.ThirdPoint()];
2378       //
2379       if (TriContact(P1, P2, P3, Q1, Q2, Q3, CoupleAngle)) {
2380         IntPolyh_Couple aCouple(i_S1, i_S2, CoupleAngle);
2381         TTrianglesContacts.Append(aCouple);
2382         //
2383         Triangle1.SetIntersection(Standard_True);
2384         Triangle2.SetIntersection(Standard_True);
2385       }
2386     }
2387   }
2388   return TTrianglesContacts.Extent();
2389 }
2390
2391 //=======================================================================
2392 //function : CheckCoupleAndGetAngle
2393 //purpose  : 
2394 //=======================================================================
2395 Standard_Boolean CheckCoupleAndGetAngle(const Standard_Integer T1, 
2396                                         const Standard_Integer T2,
2397                                         Standard_Real& Angle, 
2398                                         IntPolyh_ListOfCouples &TTrianglesContacts) 
2399 {
2400   IntPolyh_ListIteratorOfListOfCouples aIt(TTrianglesContacts);
2401   for (; aIt.More(); aIt.Next()) {
2402     IntPolyh_Couple& TestCouple = aIt.ChangeValue();
2403     if (!TestCouple.IsAnalyzed()) {
2404       if (TestCouple.FirstValue() == T1 && TestCouple.SecondValue() == T2) {
2405         TestCouple.SetAnalyzed(Standard_True);
2406         Angle = TestCouple.Angle();
2407         return Standard_True;
2408       }
2409     }
2410   }
2411   return Standard_False;
2412 }
2413 //=======================================================================
2414 //function : CheckCoupleAndGetAngle2
2415 //purpose  : 
2416 //=======================================================================
2417 Standard_Boolean CheckCoupleAndGetAngle2(const Standard_Integer T1,
2418                                          const Standard_Integer T2,
2419                                          const Standard_Integer T11,
2420                                          const Standard_Integer T22,
2421                                          IntPolyh_ListIteratorOfListOfCouples& theItCT11,
2422                                          IntPolyh_ListIteratorOfListOfCouples& theItCT22,
2423                                          Standard_Real & Angle,
2424                                          IntPolyh_ListOfCouples &TTrianglesContacts) 
2425 {
2426   ///couple T1 T2 is found in the list
2427   ///T11 and T22 are two other triangles implied  in the contact edge edge
2428   /// CT11 couple( T1,T22) and CT22 couple (T2,T11)
2429   /// these couples will be marked if there is a start point
2430   Standard_Boolean Test1 , Test2, Test3;
2431   Test1 = Test2 = Test3 = Standard_False;
2432   //
2433   IntPolyh_ListIteratorOfListOfCouples aIt(TTrianglesContacts);
2434   for (; aIt.More(); aIt.Next()) {
2435     IntPolyh_Couple& TestCouple = aIt.ChangeValue();
2436     if (TestCouple.IsAnalyzed()) {
2437       continue;
2438     }
2439     //
2440     if (TestCouple.FirstValue() == T1) {
2441       if (TestCouple.SecondValue() == T2) {
2442         Test1 = Standard_True;
2443         TestCouple.SetAnalyzed(Standard_True);
2444         Angle = TestCouple.Angle();
2445       }
2446       else if (TestCouple.SecondValue() == T22) {
2447         Test2 = Standard_True;
2448         theItCT11 = aIt;
2449         Angle = TestCouple.Angle();
2450       }
2451     }
2452     else if (TestCouple.FirstValue() == T11) {
2453       if (TestCouple.SecondValue() == T2) {
2454         Test3 = Standard_True;
2455         theItCT22 = aIt;
2456         Angle = TestCouple.Angle();
2457       }
2458     }
2459     //
2460     if (Test1 && Test2 && Test3) {
2461       break;
2462     }
2463   }
2464   return Test1;
2465 }
2466 //=======================================================================
2467 //function : CheckNextStartPoint
2468 //purpose  : it is checked if the point is not a top
2469 //           then it is stored in one or several valid arrays with 
2470 // the proper list number
2471 //=======================================================================
2472 Standard_Integer CheckNextStartPoint(IntPolyh_SectionLine & SectionLine,
2473                                      IntPolyh_ArrayOfTangentZones & TTangentZones,
2474                                      IntPolyh_StartPoint & SP,
2475                                      const Standard_Boolean Prepend)//=Standard_False) 
2476 {
2477   Standard_Integer Test=1;
2478   if( (SP.E1()==-1)||(SP.E2()==-1) ) {
2479     //The tops of triangle are analyzed
2480     //It is checked if they are not in the array TTangentZones
2481     Standard_Integer FinTTZ=TTangentZones.NbItems();
2482     for(Standard_Integer uiui=0; uiui<FinTTZ; uiui++) {
2483       IntPolyh_StartPoint TestSP=TTangentZones[uiui];
2484       if ( (Abs(SP.U1()-TestSP.U1())<MyConfusionPrecision)
2485           &&(Abs(SP.V1()-TestSP.V1())<MyConfusionPrecision) ) {
2486         if ( (Abs(SP.U2()-TestSP.U2())<MyConfusionPrecision)
2487             &&(Abs(SP.V2()-TestSP.V2())<MyConfusionPrecision) ) {
2488           Test=0;//SP is already in the list of  tops
2489           uiui=FinTTZ;
2490         }
2491       }
2492     }
2493     if (Test) {//the top does not belong to the list of TangentZones
2494       SP.SetChainList(-1);
2495       TTangentZones[FinTTZ]=SP;
2496       TTangentZones.IncrementNbItems();
2497       Test=0;//the examined point is a top
2498     }
2499   }
2500   else if (Test) {
2501     if (Prepend)
2502       SectionLine.Prepend(SP);
2503     else {
2504       SectionLine[SectionLine.NbStartPoints()]=SP;
2505       SectionLine.IncrementNbStartPoints();
2506     }
2507
2508   }
2509   //if the point is not a top Test=1
2510   //The chain is continued
2511   return(Test); 
2512 }
2513 //=======================================================================
2514 //function : StartPointsChain
2515 //purpose  : Loop on the array of couples. Compute StartPoints.
2516 //           Try to chain  the StartPoints into SectionLines or
2517 //           put  the  point  in  the    ArrayOfTangentZones if
2518 //           chaining it, is not possible.
2519 //=======================================================================
2520 Standard_Integer IntPolyh_MaillageAffinage::StartPointsChain
2521   (IntPolyh_ArrayOfSectionLines& TSectionLines,
2522    IntPolyh_ArrayOfTangentZones& TTangentZones) 
2523 {
2524   //Loop on the array of couples filled in the function COMPARE()
2525   IntPolyh_ListIteratorOfListOfCouples aIt(TTrianglesContacts);
2526   for (; aIt.More(); aIt.Next()) {
2527     IntPolyh_Couple& aCouple = aIt.ChangeValue();
2528     // Check if the couple of triangles has not been already examined.
2529     if(!aCouple.IsAnalyzed()) {
2530
2531       Standard_Integer SectionLineIndex=TSectionLines.NbItems();
2532       // fill last section line if still empty (eap)
2533       if (SectionLineIndex > 0
2534           &&
2535           TSectionLines[SectionLineIndex-1].NbStartPoints() == 0)
2536         SectionLineIndex -= 1;
2537       else
2538         TSectionLines.IncrementNbItems();
2539
2540       IntPolyh_SectionLine &  MySectionLine=TSectionLines[SectionLineIndex];
2541       if (MySectionLine.GetN() == 0) // eap
2542         MySectionLine.Init(10000);//Initialisation of array of StartPoint
2543
2544       Standard_Integer NbPoints=-1;
2545       Standard_Integer T1I, T2I;
2546       T1I = aCouple.FirstValue();
2547       T2I = aCouple.SecondValue();
2548       
2549       // Start points for the current couple are found
2550       IntPolyh_StartPoint SP1, SP2;
2551       NbPoints=StartingPointsResearch(T1I,T2I,SP1, SP2);//first calculation
2552       aCouple.SetAnalyzed(Standard_True);//the couple is marked
2553
2554       if(NbPoints==1) {// particular case top/triangle or edge/edge
2555         //the start point is input in the array
2556         SP1.SetChainList(SectionLineIndex);
2557         SP1.SetAngle(aCouple.Angle());
2558         //it is checked if the point is not atop of the triangle
2559         if(CheckNextStartPoint(MySectionLine,TTangentZones,SP1)) {
2560           IntPolyh_StartPoint SPNext1;
2561           Standard_Integer TestSP1=0;
2562           
2563           //chain of a side
2564           IntPolyh_StartPoint SP11;//=SP1;
2565           if(SP1.E1()>=0) { //&&(SP1.E2()!=-1) already tested if the point is not a top
2566             Standard_Integer NextTriangle1=0;
2567             if (TEdges1[SP1.E1()].FirstTriangle()!=T1I) NextTriangle1=TEdges1[SP1.E1()].FirstTriangle();
2568             else NextTriangle1=TEdges1[SP1.E1()].SecondTriangle();
2569             
2570             Standard_Real Angle=-2.0;
2571             if (CheckCoupleAndGetAngle(NextTriangle1,T2I,Angle,TTrianglesContacts)) {
2572               //it is checked if the couple exists and is marked
2573               Standard_Integer NbPoints11=0;
2574               NbPoints11=NextStartingPointsResearch(NextTriangle1,T2I,SP1,SP11);
2575               if (NbPoints11==1) {
2576                 SP11.SetChainList(SectionLineIndex);
2577                 SP11.SetAngle(Angle);
2578                 
2579                 if(CheckNextStartPoint(MySectionLine,TTangentZones,SP11)) {            
2580                   Standard_Integer EndChainList=1;
2581                   while (EndChainList!=0) {
2582                     TestSP1=GetNextChainStartPoint(SP11,SPNext1,MySectionLine,TTangentZones);
2583                     if(TestSP1==1) {
2584                       SPNext1.SetChainList(SectionLineIndex);
2585                       if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext1))
2586                         SP11=SPNext1;
2587                       else EndChainList=0;
2588                     }
2589                     else EndChainList=0; //There is no next point
2590                   }
2591                 }
2592                
2593               }
2594               else {
2595                 if(NbPoints11>1) {//The point is input in the array TTangentZones
2596                   TTangentZones[TTangentZones.NbItems()]=SP11;//default list number = -1
2597                   TTangentZones.IncrementNbItems();
2598                 }
2599                 else {
2600
2601                 }
2602               }
2603             }
2604           }
2605           else if (SP1.E2()<0){
2606
2607           }
2608           //chain of the other side
2609           IntPolyh_StartPoint SP12;//=SP1;
2610           if (SP1.E2()>=0) { //&&(SP1.E1()!=-1) already tested
2611             Standard_Integer NextTriangle2;
2612             if (TEdges2[SP1.E2()].FirstTriangle()!=T2I) NextTriangle2=TEdges2[SP1.E2()].FirstTriangle();
2613             else NextTriangle2=TEdges2[SP1.E2()].SecondTriangle();
2614             
2615             Standard_Real Angle=-2.0;
2616             if(CheckCoupleAndGetAngle(T1I,NextTriangle2,Angle,TTrianglesContacts)) {
2617               Standard_Integer NbPoints12=0;
2618               NbPoints12=NextStartingPointsResearch(T1I,NextTriangle2,SP1, SP12);
2619               if (NbPoints12==1) {
2620                 
2621                 SP12.SetChainList(SectionLineIndex);
2622                 SP12.SetAngle(Angle);
2623                 Standard_Boolean Prepend = Standard_True; // eap
2624                 
2625                 if(CheckNextStartPoint(MySectionLine,TTangentZones,SP12, Prepend)) {
2626                   Standard_Integer EndChainList=1;
2627                   while (EndChainList!=0) {
2628                     TestSP1=GetNextChainStartPoint(SP12,SPNext1,
2629                                                    MySectionLine,TTangentZones,
2630                                                    Prepend); // eap
2631                     if(TestSP1==1) {
2632                       SPNext1.SetChainList(SectionLineIndex);
2633                       if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext1,Prepend))
2634                         SP12=SPNext1;
2635                       else EndChainList=0;
2636                     }
2637                     else EndChainList=0; //there is no next point
2638                   }
2639                 }
2640                 
2641                 else {
2642                   if(NbPoints12>1) {//The points are input in the array TTangentZones
2643                     TTangentZones[TTangentZones.NbItems()]=SP12;//default list number = -1
2644                     TTangentZones.IncrementNbItems();
2645                   }
2646                   else {
2647
2648                   }
2649                 }
2650               }
2651             }
2652           }
2653           else if(SP1.E1()<0){
2654
2655           }
2656         }
2657       }
2658       else if(NbPoints==2) {
2659         //the start points are input in the array 
2660         IntPolyh_StartPoint SPNext2;
2661         Standard_Integer TestSP2=0;
2662         Standard_Integer EndChainList=1;
2663
2664         SP1.SetChainList(SectionLineIndex);
2665         SP1.SetAngle(aCouple.Angle());
2666         if(CheckNextStartPoint(MySectionLine,TTangentZones,SP1)) {
2667
2668           //chain of a side
2669           while (EndChainList!=0) {
2670             TestSP2=GetNextChainStartPoint(SP1,SPNext2,MySectionLine,TTangentZones);
2671             if(TestSP2==1) {
2672               SPNext2.SetChainList(SectionLineIndex);
2673               if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext2))
2674                 SP1=SPNext2;
2675               else EndChainList=0;
2676             }
2677             else EndChainList=0; //there is no next point
2678           }
2679         }
2680         
2681         SP2.SetChainList(SectionLineIndex);
2682         SP2.SetAngle(aCouple.Angle());
2683         Standard_Boolean Prepend = Standard_True; // eap
2684
2685         if(CheckNextStartPoint(MySectionLine,TTangentZones,SP2,Prepend)) {
2686           
2687           //chain of the other side
2688           EndChainList=1;
2689           while (EndChainList!=0) {
2690             TestSP2=GetNextChainStartPoint(SP2,SPNext2,
2691                                            MySectionLine,TTangentZones,
2692                                            Prepend); // eap
2693             if(TestSP2==1) {
2694               SPNext2.SetChainList(SectionLineIndex);
2695               if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext2,Prepend))
2696                 SP2=SPNext2;
2697               else EndChainList=0;
2698             }
2699             else EndChainList=0; //there is no next point
2700           }
2701         }
2702       }
2703
2704       else if( (NbPoints>2)&&(NbPoints<7) ) {
2705         //More than two start points 
2706         //the start point is input in the table
2707         SP1.SetChainList(SectionLineIndex);
2708         CheckNextStartPoint(MySectionLine,TTangentZones,SP1);
2709       }
2710       
2711       else {
2712
2713       }
2714     }
2715   }
2716
2717   return(1);
2718 }
2719 //=======================================================================
2720 //function : GetNextChainStartPoint
2721 //purpose  : Mainly  used  by StartPointsChain(), this function
2722 //           try to compute the next StartPoint.
2723 //           GetNextChainStartPoint is used only if it is known that there are 2 contact points
2724 //=======================================================================
2725 Standard_Integer IntPolyh_MaillageAffinage::GetNextChainStartPoint
2726   (const IntPolyh_StartPoint & SP,
2727    IntPolyh_StartPoint & SPNext,
2728    IntPolyh_SectionLine & MySectionLine,
2729    IntPolyh_ArrayOfTangentZones & TTangentZones,
2730    const Standard_Boolean Prepend) 
2731 {
2732   Standard_Integer NbPoints=0;
2733   if( (SP.E1()>=0)&&(SP.E2()==-2) ) {
2734     //case if the point is on edge of T1
2735     Standard_Integer NextTriangle1;
2736     if (TEdges1[SP.E1()].FirstTriangle()!=SP.T1()) NextTriangle1=TEdges1[SP.E1()].FirstTriangle();
2737     else 
2738       NextTriangle1=TEdges1[SP.E1()].SecondTriangle();
2739     //If is checked if two triangles intersect
2740     Standard_Real Angle= -2.0;
2741     if (CheckCoupleAndGetAngle(NextTriangle1,SP.T2(),Angle,TTrianglesContacts)) {
2742       NbPoints=NextStartingPointsResearch(NextTriangle1,SP.T2(),SP,SPNext);
2743       if( NbPoints!=1 ) {
2744         if (NbPoints>1)
2745           CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend);
2746         else {
2747
2748         NbPoints=0;
2749         }
2750       }
2751       else 
2752         SPNext.SetAngle(Angle);
2753     }
2754     else NbPoints=0;//this couple does not intersect
2755   }
2756   else if( (SP.E1()==-2)&&(SP.E2()>=0) ) {
2757     //case if the point is on edge of T2
2758     Standard_Integer NextTriangle2;
2759     if (TEdges2[SP.E2()].FirstTriangle()!=SP.T2()) NextTriangle2=TEdges2[SP.E2()].FirstTriangle();
2760     else 
2761       NextTriangle2=TEdges2[SP.E2()].SecondTriangle();
2762     Standard_Real Angle= -2.0;
2763     if (CheckCoupleAndGetAngle(SP.T1(),NextTriangle2,Angle,TTrianglesContacts)) {
2764       NbPoints=NextStartingPointsResearch(SP.T1(),NextTriangle2,SP,SPNext);
2765       if( NbPoints!=1 ) {
2766         if (NbPoints>1)
2767           CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend);
2768         else {
2769
2770         NbPoints=0;
2771         }
2772       }
2773       else
2774         SPNext.SetAngle(Angle);
2775     }
2776     else NbPoints=0;
2777   }
2778   else if( (SP.E1()==-2)&&(SP.E2()==-2) ) { 
2779     ///no edge is touched or cut
2780
2781     NbPoints=0;
2782   }
2783   else if( (SP.E1()>=0)&&(SP.E2()>=0) ) {
2784     ///the point is located on two edges
2785       Standard_Integer NextTriangle1;
2786       if (TEdges1[SP.E1()].FirstTriangle()!=SP.T1()) NextTriangle1=TEdges1[SP.E1()].FirstTriangle();
2787       else 
2788         NextTriangle1=TEdges1[SP.E1()].SecondTriangle();
2789       Standard_Integer NextTriangle2;
2790       if (TEdges2[SP.E2()].FirstTriangle()!=SP.T2()) NextTriangle2=TEdges2[SP.E2()].FirstTriangle();
2791       else 
2792         NextTriangle2=TEdges2[SP.E2()].SecondTriangle();
2793       Standard_Real Angle= -2.0;
2794
2795       IntPolyh_ListIteratorOfListOfCouples aItCT11, aItCT22;
2796       if (CheckCoupleAndGetAngle2(NextTriangle1,NextTriangle2,
2797                                   SP.T1(),SP.T2(), aItCT11, aItCT22,
2798                                   Angle,TTrianglesContacts)) {
2799         NbPoints=NextStartingPointsResearch(NextTriangle1,NextTriangle2,SP,SPNext);
2800         if( NbPoints!=1 ) {
2801           if (NbPoints>1) {
2802             ///The new point is checked
2803             if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend)>0) {
2804             }
2805             else {
2806
2807             }
2808           }
2809           NbPoints=0;
2810         }
2811         else {//NbPoints==1
2812           SPNext.SetAngle(Angle);
2813           if (aItCT11.More()) aItCT11.ChangeValue().SetAnalyzed(Standard_True);
2814           if (aItCT22.More()) aItCT22.ChangeValue().SetAnalyzed(Standard_True);
2815         }
2816       }
2817       else NbPoints=0;
2818     }
2819   else if( (SP.E1()==-1)||(SP.E2()==-1) ) {
2820     ///the points are tops of triangle
2821     ///the point is atored in an intermediary array
2822   }
2823   return(NbPoints);
2824 }
2825 //=======================================================================
2826 //function : GetArrayOfPoints
2827 //purpose  : 
2828 //=======================================================================
2829 const IntPolyh_ArrayOfPoints& IntPolyh_MaillageAffinage::GetArrayOfPoints
2830   (const Standard_Integer SurfID)const
2831 {
2832  if (SurfID==1)
2833  return(TPoints1);
2834  return(TPoints2);
2835 }
2836 //=======================================================================
2837 //function : GetArrayOfEdges
2838 //purpose  : 
2839 //=======================================================================
2840 const IntPolyh_ArrayOfEdges& IntPolyh_MaillageAffinage::GetArrayOfEdges
2841   (const Standard_Integer SurfID)const
2842 {
2843  if (SurfID==1)
2844  return(TEdges1);
2845  return(TEdges2);
2846 }
2847 //=======================================================================
2848 //function : GetArrayOfTriangles
2849 //purpose  : 
2850 //=======================================================================
2851 const IntPolyh_ArrayOfTriangles& 
2852   IntPolyh_MaillageAffinage::GetArrayOfTriangles
2853   (const Standard_Integer SurfID)const{
2854   if (SurfID==1)
2855   return(TTriangles1);
2856   return(TTriangles2);
2857 }
2858
2859 //=======================================================================
2860 //function : GetBox
2861 //purpose  : 
2862 //=======================================================================
2863 Bnd_Box IntPolyh_MaillageAffinage::GetBox(const Standard_Integer SurfID) const
2864 {
2865   if (SurfID==1)
2866   return(MyBox1);
2867   return(MyBox2);
2868 }
2869 //=======================================================================
2870 //function : GetArrayOfCouples
2871 //purpose  : 
2872 //=======================================================================
2873 IntPolyh_ListOfCouples &IntPolyh_MaillageAffinage::GetCouples()
2874 {
2875   return TTrianglesContacts;
2876 }
2877 //=======================================================================
2878 //function : SetEnlargeZone
2879 //purpose  : 
2880 //=======================================================================
2881 void IntPolyh_MaillageAffinage::SetEnlargeZone(Standard_Boolean& EnlargeZone)
2882 {
2883   myEnlargeZone = EnlargeZone;
2884 }
2885 //=======================================================================
2886 //function : GetEnlargeZone
2887 //purpose  : 
2888 //=======================================================================
2889 Standard_Boolean IntPolyh_MaillageAffinage::GetEnlargeZone() const
2890 {
2891   return myEnlargeZone;
2892 }
2893
2894 //=======================================================================
2895 //function : GetMinDeflection
2896 //purpose  : 
2897 //=======================================================================
2898 Standard_Real IntPolyh_MaillageAffinage::GetMinDeflection(const Standard_Integer SurfID) const
2899 {
2900   return (SurfID==1)? FlecheMin1:FlecheMin2;
2901 }
2902
2903 //=======================================================================
2904 //function : GetMaxDeflection
2905 //purpose  : 
2906 //=======================================================================
2907 Standard_Real IntPolyh_MaillageAffinage::GetMaxDeflection(const Standard_Integer SurfID) const
2908 {
2909   return (SurfID==1)? FlecheMax1:FlecheMax2;
2910 }
2911
2912 //=======================================================================
2913 //function : DegeneratedIndex
2914 //purpose  : 
2915 //=======================================================================
2916 void DegeneratedIndex(const TColStd_Array1OfReal& aXpars,
2917                       const Standard_Integer aNbX,
2918                       const Handle(Adaptor3d_HSurface)& aS,
2919                       const Standard_Integer aIsoDirection,
2920                       Standard_Integer& aI1,
2921                       Standard_Integer& aI2)
2922 {
2923   Standard_Integer i;
2924   Standard_Boolean bDegX1, bDegX2;
2925   Standard_Real aDegX1, aDegX2, aTol2, aX;
2926   //
2927   aI1=0;
2928   aI2=0;
2929   aTol2=MyTolerance*MyTolerance;
2930   //
2931   if (aIsoDirection==1){ // V=const
2932     bDegX1=IsDegenerated(aS, 1, aTol2, aDegX1);
2933     bDegX2=IsDegenerated(aS, 2, aTol2, aDegX2);
2934   }
2935   else if (aIsoDirection==2){ // U=const
2936     bDegX1=IsDegenerated(aS, 3, aTol2, aDegX1);
2937     bDegX2=IsDegenerated(aS, 4, aTol2, aDegX2);
2938   }
2939   else {
2940     return;
2941   }
2942   //
2943   if (!(bDegX1 || bDegX2)) {
2944     return;
2945   }
2946   //  
2947   for(i=1; i<=aNbX; ++i) {
2948     aX=aXpars(i);
2949     if (bDegX1) {
2950       if (fabs(aX-aDegX1) < MyTolerance) {
2951         aI1=i;
2952       }
2953     }
2954     if (bDegX2) {
2955       if (fabs(aX-aDegX2) < MyTolerance) {
2956         aI2=i;
2957       }
2958     }
2959   }
2960 }
2961 //=======================================================================
2962 //function : IsDegenerated
2963 //purpose  : 
2964 //=======================================================================
2965 Standard_Boolean IsDegenerated(const Handle(Adaptor3d_HSurface)& aS,
2966                                const Standard_Integer aIndex,
2967                                const Standard_Real aTol2,
2968                                Standard_Real& aDegX)
2969 {
2970   Standard_Boolean bRet;
2971   Standard_Integer i, aNbP;
2972   Standard_Real aU, dU, aU1, aU2, aV, dV, aV1, aV2, aD2;
2973   gp_Pnt aP1, aP2;
2974   //
2975   bRet=Standard_False;
2976   aNbP=3;
2977   aDegX=99;
2978   //
2979   aU1=aS->FirstUParameter();
2980   aU2=aS->LastUParameter();
2981   aV1=aS->FirstVParameter();
2982   aV2=aS->LastVParameter();
2983   //
2984   if (aIndex<3) { // V=const
2985     aV=aV1;
2986     if (aIndex==2) {
2987       aV=aV2;
2988     }
2989     dU=(aU2-aU1)/(aNbP-1);
2990     aU=aU1;
2991     aP1=aS->Value(aU, aV);
2992     for (i=1; i<aNbP; ++i) {
2993       aU=i*dU;
2994       if (i==aNbP-1){
2995         aU=aU2;
2996       }
2997       aP2=aS->Value(aU, aV);
2998       aD2=aP1.SquareDistance(aP2);
2999       if (aD2>aTol2) {
3000         return bRet;
3001       }
3002       aP1=aP2;
3003     }
3004     aDegX=aV;
3005     bRet=!bRet;
3006   }
3007   else {// U=const
3008     aU=aU1;
3009     if (aIndex==4) {
3010       aU=aU2;
3011     }
3012     dV=(aV2-aV1)/(aNbP-1);
3013     aV=aV1;
3014     aP1=aS->Value(aU, aV);
3015     for (i=1; i<aNbP; ++i) {
3016       aV=i*dV;
3017       if (i==aNbP-1){
3018         aV=aV2;
3019       }
3020       aP2=aS->Value(aU, aV);
3021       aD2=aP1.SquareDistance(aP2);
3022       if (aD2>aTol2) {
3023         return bRet;
3024       }
3025       aP1=aP2;
3026     }
3027     bRet=!bRet;
3028     aDegX=aU;
3029   }
3030   //
3031   return bRet;
3032 }
3033 //=======================================================================
3034 //function : EnlargeZone
3035 //purpose  : 
3036 //=======================================================================
3037 void EnlargeZone(const Handle(Adaptor3d_HSurface)& MaSurface,
3038                  Standard_Real &u0, 
3039                  Standard_Real &u1, 
3040                  Standard_Real &v0, 
3041                  Standard_Real &v1) 
3042 {
3043   if(MaSurface->GetType() == GeomAbs_BSplineSurface ||
3044      MaSurface->GetType() == GeomAbs_BezierSurface) {
3045     if((!MaSurface->IsUClosed() && !MaSurface->IsUPeriodic()) &&
3046        (Abs(u0) < 1.e+100 && Abs(u1) < 1.e+100) ) {
3047       Standard_Real delta_u = 0.01*Abs(u1 - u0);
3048       u0 -= delta_u;
3049       u1 += delta_u;
3050     }
3051     if((!MaSurface->IsVClosed() && !MaSurface->IsVPeriodic()) &&
3052        (Abs(v0) < 1.e+100 && Abs(v1) < 1.e+100) ) {
3053       Standard_Real delta_v = 0.01*Abs(v1 - v0);
3054       v0 -= delta_v;
3055       v1 += delta_v;
3056     }
3057   }
3058 }