0024510: Remove unused local variables
[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
9 // under the terms of the GNU Lesser General Public 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
21 //  modified by Edward AGAPOV (eap) Thu Feb 14 2002 (occ139)
22 //  - make Section Line parts rightly connected (prepend 2nd part to the 1st)
23 //  - TriangleShape() for debugging purpose
24
25 //  Modified by skv - Thu Sep 25 17:42:42 2003 OCC567
26 //  modified by ofv Thu Apr  8 14:58:13 2004 fip
27
28
29 #include <IntPolyh_MaillageAffinage.ixx>
30
31 #include <Precision.hxx>
32 #include <gp_Pnt.hxx>
33 #include <gp.hxx>
34
35 #include <TColStd_ListIteratorOfListOfInteger.hxx>
36
37 #include <Bnd_BoundSortBox.hxx>
38 #include <Bnd_HArray1OfBox.hxx> 
39
40 #include <IntCurveSurface_ThePolyhedronOfHInter.hxx>
41
42 #include <IntPolyh_ArrayOfCouples.hxx>
43 #include <IntPolyh_Edge.hxx>
44 #include <IntPolyh_Couple.hxx>
45
46 static Standard_Real MyTolerance=10.0e-7;
47 static Standard_Real MyConfusionPrecision=10.0e-12;
48 static Standard_Real SquareMyConfusionPrecision=10.0e-24;
49 //
50 static
51   inline Standard_Real maxSR(const Standard_Real a,
52                              const Standard_Real b,
53                              const Standard_Real c);
54
55 static
56   inline Standard_Real minSR(const Standard_Real a,
57                              const Standard_Real b,
58                              const Standard_Real c);
59 static
60   Standard_Integer project6(const IntPolyh_Point &ax, 
61                             const IntPolyh_Point &p1,
62                             const IntPolyh_Point &p2, 
63                             const IntPolyh_Point &p3, 
64                             const IntPolyh_Point &q1,
65                             const IntPolyh_Point &q2, 
66                             const IntPolyh_Point &q3);
67 static
68   void TestNbPoints(const Standard_Integer ,
69                     Standard_Integer &NbPoints,
70                     Standard_Integer &NbPointsTotal,
71                     const IntPolyh_StartPoint &Pt1,
72                     const IntPolyh_StartPoint &Pt2,
73                     IntPolyh_StartPoint &SP1,
74                     IntPolyh_StartPoint &SP2);
75 static
76   void CalculPtsInterTriEdgeCoplanaires(const Standard_Integer TriSurfID,
77                                         const IntPolyh_Point &NormaleTri,
78                                         const IntPolyh_Point &PE1,
79                                         const IntPolyh_Point &PE2,
80                                         const IntPolyh_Point &Edge,
81                                         const IntPolyh_Point &PT1,
82                                         const IntPolyh_Point &PT2,
83                                         const IntPolyh_Point &Cote,
84                                         const Standard_Integer CoteIndex,
85                                         IntPolyh_StartPoint &SP1,
86                                         IntPolyh_StartPoint &SP2,
87                                         Standard_Integer &NbPoints);
88 static
89   void CalculPtsInterTriEdgeCoplanaires2(const Standard_Integer TriSurfID,
90                                          const IntPolyh_Point &NormaleTri,
91                                          const IntPolyh_Triangle &Tri1,
92                                          const IntPolyh_Triangle &Tri2,
93                                          const IntPolyh_Point &PE1,
94                                          const IntPolyh_Point &PE2,
95                                          const IntPolyh_Point &Edge,
96                                          const Standard_Integer EdgeIndex,
97                                          const IntPolyh_Point &PT1,
98                                          const IntPolyh_Point &PT2,
99                                          const IntPolyh_Point &Cote,
100                                        const Standard_Integer CoteIndex,
101                                          IntPolyh_StartPoint &SP1,
102                                          IntPolyh_StartPoint &SP2,
103                                          Standard_Integer &NbPoints); 
104 static
105   Standard_Integer CheckCoupleAndGetAngle(const Standard_Integer T1, 
106                                           const Standard_Integer T2,
107                                           Standard_Real& Angle, 
108                                           IntPolyh_ArrayOfCouples &TTrianglesContacts);
109 static
110   Standard_Integer CheckCoupleAndGetAngle2(const Standard_Integer T1,
111                                            const Standard_Integer T2,
112                                            const Standard_Integer T11, 
113                                            const Standard_Integer T22,
114                                            Standard_Integer &CT11,
115                                            Standard_Integer &CT22, 
116                                            Standard_Real & Angle,
117                                            IntPolyh_ArrayOfCouples &TTrianglesContacts);
118 static
119   Standard_Integer CheckNextStartPoint(IntPolyh_SectionLine & SectionLine,
120                                        IntPolyh_ArrayOfTangentZones & TTangentZones,
121                                        IntPolyh_StartPoint & SP,
122                                        const Standard_Boolean Prepend=Standard_False); 
123
124 static
125   Standard_Boolean IsDegenerated(const Handle(Adaptor3d_HSurface)& aS,
126                                  const Standard_Integer aIndex,
127                                  const Standard_Real aTol2,
128                                  Standard_Real& aDegX);
129 static
130   void DegeneratedIndex(const TColStd_Array1OfReal& Xpars,
131                         const Standard_Integer aNbX,
132                         const Handle(Adaptor3d_HSurface)& aS,
133                         const Standard_Integer aIsoDirection,
134                         Standard_Integer& aI1,
135                         Standard_Integer& aI2);
136 static
137   void EnlargeZone(const Handle(Adaptor3d_HSurface)& MaSurface,
138                    Standard_Real &u0, 
139                    Standard_Real &u1, 
140                    Standard_Real &v0, 
141                    Standard_Real &v1);
142  
143 //=======================================================================
144 //function : IntPolyh_MaillageAffinage
145 //purpose  : 
146 //=======================================================================
147 IntPolyh_MaillageAffinage::IntPolyh_MaillageAffinage
148   (const Handle(Adaptor3d_HSurface)& Surface1,
149    const Handle(Adaptor3d_HSurface)& Surface2,
150    const Standard_Integer )
151 :
152   MaSurface1(Surface1), 
153   MaSurface2(Surface2), 
154   NbSamplesU1(10), 
155   NbSamplesU2(10), 
156   NbSamplesV1(10), 
157   NbSamplesV2(10),
158   FlecheMax1(0.0), 
159   FlecheMax2(0.0), 
160   FlecheMin1(0.0), 
161   FlecheMin2(0.0),
162   FlecheMoy1(0.0), 
163   FlecheMoy2(0.0), 
164   myEnlargeZone(Standard_False) 
165
166 }
167 //=======================================================================
168 //function : IntPolyh_MaillageAffinage
169 //purpose  : 
170 //=======================================================================
171 IntPolyh_MaillageAffinage::IntPolyh_MaillageAffinage
172   (const Handle(Adaptor3d_HSurface)& Surface1,
173    const Standard_Integer NbSU1,
174    const Standard_Integer NbSV1,
175    const Handle(Adaptor3d_HSurface)& Surface2,
176    const Standard_Integer NbSU2,
177    const Standard_Integer NbSV2,
178    const Standard_Integer )
179 :
180   MaSurface1(Surface1), 
181   MaSurface2(Surface2), 
182   NbSamplesU1(NbSU1), 
183   NbSamplesU2(NbSU2), 
184   NbSamplesV1(NbSV1), 
185   NbSamplesV2(NbSV2),
186   FlecheMax1(0.0), 
187   FlecheMax2(0.0), 
188   FlecheMin1(0.0), 
189   FlecheMin2(0.0),
190   FlecheMoy1(0.0), 
191   FlecheMoy2(0.0), 
192   myEnlargeZone(Standard_False)
193
194 }
195 //=======================================================================
196 //function : FillArrayOfPnt
197 //purpose  : Compute points on one surface and fill an array of points
198 //=======================================================================
199 void IntPolyh_MaillageAffinage::FillArrayOfPnt
200   (const Standard_Integer SurfID)
201 {
202   Standard_Integer NbSamplesU, NbSamplesV, i, aNbSamplesU1, aNbSamplesV1;
203   Standard_Real u0, u1, v0, v1, aU, aV, dU, dV;
204   //
205   const Handle(Adaptor3d_HSurface&) MaSurface=(SurfID==1)? MaSurface1 : MaSurface2;
206   NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
207   NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
208   //
209   u0 = (MaSurface)->FirstUParameter();  
210   u1 = (MaSurface)->LastUParameter();
211   v0 = (MaSurface)->FirstVParameter();  
212   v1 = (MaSurface)->LastVParameter();  
213
214   if(myEnlargeZone) { 
215     EnlargeZone(MaSurface, u0, u1, v0, v1);
216   }
217   //
218   TColStd_Array1OfReal aUpars(1, NbSamplesU);
219   TColStd_Array1OfReal aVpars(1, NbSamplesV);
220   //
221   aNbSamplesU1=NbSamplesU-1;
222   aNbSamplesV1=NbSamplesV-1;
223   //
224   dU=(u1-u0)/Standard_Real(aNbSamplesU1);
225   dV=(v1-v0)/Standard_Real(aNbSamplesV1);
226   //
227   for (i=0; i<NbSamplesU; ++i) {
228     aU=u0+i*dU;
229     if (i==aNbSamplesU1) {
230       aU=u1;
231     }
232     aUpars.SetValue(i+1, aU);
233   }
234   //
235   for (i=0; i<NbSamplesV; ++i) {
236     aV=v0+i*dV;
237     if (i==aNbSamplesV1) {
238       aV=v1;
239     }
240     aVpars.SetValue(i+1, aV);
241   }
242   //
243   FillArrayOfPnt(SurfID, aUpars, aVpars);
244 }
245 //=======================================================================
246 //function : FillArrayOfPnt
247 //purpose  : Compute points on one surface and fill an array of points
248 //           FILL AN ARRAY OF POINTS
249 //=======================================================================
250 void IntPolyh_MaillageAffinage::FillArrayOfPnt
251   (const Standard_Integer SurfID,
252    const Standard_Boolean isShiftFwd)
253 {
254   Standard_Integer NbSamplesU, NbSamplesV, i, aNbSamplesU1, aNbSamplesV1; 
255   Standard_Real u0, u1, v0, v1, aU, aV, dU, dV;
256   const Handle(Adaptor3d_HSurface)& MaSurface=(SurfID==1)? MaSurface1 : MaSurface2;
257   NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
258   NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
259
260   u0 = (MaSurface)->FirstUParameter();  
261   u1 = (MaSurface)->LastUParameter();
262   v0 = (MaSurface)->FirstVParameter();  
263   v1 = (MaSurface)->LastVParameter();  
264
265   if(myEnlargeZone) {
266     EnlargeZone(MaSurface, u0, u1, v0, v1);
267   }
268   //
269   TColStd_Array1OfReal aUpars(1, NbSamplesU);
270   TColStd_Array1OfReal aVpars(1, NbSamplesV);
271   //
272   aNbSamplesU1=NbSamplesU-1;
273   aNbSamplesV1=NbSamplesV-1;
274   //
275   dU=(u1-u0)/Standard_Real(aNbSamplesU1);
276   dV=(v1-v0)/Standard_Real(aNbSamplesV1);
277   //
278   for (i=0; i<NbSamplesU; ++i) {
279     aU=u0+i*dU;
280     if (i==aNbSamplesU1) {
281       aU=u1;
282     }
283     aUpars.SetValue(i+1, aU);
284   }
285   //
286   for (i=0; i<NbSamplesV; ++i) {
287     aV=v0+i*dV;
288     if (i==aNbSamplesV1) {
289       aV=v1;
290     }
291     aVpars.SetValue(i+1, aV);
292   }
293   //
294   FillArrayOfPnt(SurfID, isShiftFwd, aUpars, aVpars);  
295 }
296 //=======================================================================
297 //function : FillArrayOfPnt
298 //purpose  : Compute points on one surface and fill an array of points
299 //=======================================================================
300 void IntPolyh_MaillageAffinage::FillArrayOfPnt
301   (const Standard_Integer SurfID, 
302    const TColStd_Array1OfReal& Upars,
303    const TColStd_Array1OfReal& Vpars)
304 {
305   Standard_Boolean bDegI, bDeg;
306   Standard_Integer aNbU, aNbV, iCnt, i, j;
307   Standard_Integer aID1, aID2, aJD1, aJD2;
308   Standard_Real aTol, aU, aV, aX, aY, aZ;
309   gp_Pnt aP;
310   //
311   aNbU=(SurfID==1)? NbSamplesU1 : NbSamplesU2;
312   aNbV=(SurfID==1)? NbSamplesV1 : NbSamplesV2;
313   Bnd_Box& aBox = (SurfID==1) ? MyBox1 : MyBox2;
314   Handle(Adaptor3d_HSurface)& aS=(SurfID==1)? MaSurface1:MaSurface2;
315   IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
316   //
317   aJD1=0;
318   aJD2=0;
319   aID1=0;
320   aID2=0;
321   DegeneratedIndex(Vpars, aNbV, aS, 1, aJD1, aJD2);
322   if (!(aJD1 || aJD2)) {
323     DegeneratedIndex(Upars, aNbU, aS, 2, aID1, aID2);
324   }
325   //
326   TPoints.Init(aNbU*aNbV);
327   iCnt=0;
328   for(i=1; i<=aNbU; ++i){
329     bDegI=(aID1==i || aID2==i);
330     aU=Upars(i);
331     for(j=1; j<=aNbV; ++j){
332       aV=Vpars(j);
333       aP=aS->Value(aU, aV);
334       aP.Coord(aX, aY, aZ);
335       IntPolyh_Point& aIP=TPoints[iCnt];
336       aIP.Set(aX, aY, aZ, aU, aV);
337       //
338       bDeg=bDegI || (aJD1==j || aJD2==j);
339       if (bDeg) {
340         aIP.SetDegenerated(bDeg);
341       }
342       ++iCnt;
343       aBox.Add(aP);
344     }
345   }
346   //
347   TPoints.SetNbItems(iCnt);
348   //
349   IntCurveSurface_ThePolyhedronOfHInter polyhedron(aS, Upars, Vpars);
350   //
351   aTol=polyhedron.DeflectionOverEstimation();
352   aTol*=1.2;
353
354   Standard_Real a1,a2,a3,b1,b2,b3;
355   //
356   aBox.Get(a1,a2,a3,b1,b2,b3);
357   aBox.Update(a1-aTol,a2-aTol,a3-aTol,b1+aTol,b2+aTol,b3+aTol);
358   aBox.Enlarge(MyTolerance);
359 }
360
361 //=======================================================================
362 //function : FillArrayOfPnt
363 //purpose  : Compute points on one surface and fill an array of points
364 //           REMPLISSAGE DU TABLEAU DE POINTS
365 //=======================================================================
366 void IntPolyh_MaillageAffinage::FillArrayOfPnt
367   (const Standard_Integer SurfID,
368    const Standard_Boolean isShiftFwd,
369    const TColStd_Array1OfReal& Upars,
370    const TColStd_Array1OfReal& Vpars)
371 {
372   Standard_Boolean bDegI, bDeg;
373   Standard_Integer aNbU, aNbV, iCnt, i, j;
374   Standard_Integer aID1, aID2, aJD1, aJD2;
375   Standard_Real Tol, resol, aU, aV, aMag;
376   Standard_Real aX, aY, aZ;
377   gp_Pnt aP;
378   gp_Vec aDU, aDV, aNorm;
379   //
380   aNbU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
381   aNbV=(SurfID==1)? NbSamplesV1:NbSamplesV2; 
382   Bnd_Box& aBox = (SurfID==1) ? MyBox1 : MyBox2;
383   Handle(Adaptor3d_HSurface) aS=(SurfID==1)? MaSurface1:MaSurface2;
384   IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
385   //
386   resol = gp::Resolution();
387   //
388   IntCurveSurface_ThePolyhedronOfHInter polyhedron(aS, Upars, Vpars);
389   Tol=polyhedron.DeflectionOverEstimation();
390   aJD1=0;
391   aJD2=0;
392   aID1=0;
393   aID2=0;
394   DegeneratedIndex(Vpars, aNbV, aS, 1, aJD1, aJD2);
395   if (!(aJD1 || aJD2)) {
396     DegeneratedIndex(Upars, aNbU, aS, 2, aID1, aID2);
397   }
398   //
399   TPoints.Init(aNbU*aNbV);
400   iCnt=0;
401   for(i=1; i<=aNbU; ++i){
402     bDegI=(aID1==i || aID2==i);
403     aU = Upars(i);
404     for(j=1; j<=aNbV; ++j){
405       aV = Vpars(j);
406       aS->D1(aU, aV, aP, aDU, aDV);
407       
408       aNorm = aDU.Crossed(aDV);
409       aMag = aNorm.Magnitude();
410       if (aMag > resol) {
411         aNorm /= aMag;
412         aNorm.Multiply(Tol*1.5);
413         //
414         if (isShiftFwd) {
415           aP.Translate(aNorm);
416         }
417         else{
418           aP.Translate(aNorm.Reversed());
419         }
420       }
421       
422       IntPolyh_Point& aIP=TPoints[iCnt];
423       aP.Coord(aX, aY, aZ);
424       aIP.Set(aX, aY, aZ, aU, aV);
425       //
426       bDeg=bDegI || (aJD1==j || aJD2==j);
427       if (bDeg) {
428         aIP.SetDegenerated(bDeg);
429       }
430       ++iCnt;
431       aBox.Add(aP);
432     }
433   }
434   //
435   TPoints.SetNbItems(iCnt);
436   //
437   Tol*=1.2;
438   //
439   Standard_Real a1,a2,a3,b1,b2,b3;
440   //
441   aBox.Get(a1,a2,a3,b1,b2,b3);
442   aBox.Update(a1-Tol,a2-Tol,a3-Tol,b1+Tol,b2+Tol,b3+Tol);
443   aBox.Enlarge(MyTolerance);
444 }
445 //=======================================================================
446 //function : CommonBox
447 //purpose  : Compute the common box  witch is the intersection
448 //           of the two bounding boxes,  and mark the points of
449 //           the two surfaces that are inside.
450 //           REJECTION BOUNDING BOXES
451 //           DETERMINATION OF THE COMMON BOX
452 //=======================================================================
453 void IntPolyh_MaillageAffinage::CommonBox (const Bnd_Box &,
454                                            const Bnd_Box &,
455                                            Standard_Real &XMin,
456                                            Standard_Real &YMin,
457                                            Standard_Real &ZMin,
458                                            Standard_Real &XMax,
459                                            Standard_Real &YMax,
460                                            Standard_Real &ZMax) 
461 {
462   Standard_Real x10,y10,z10,x11,y11,z11;
463   Standard_Real x20,y20,z20,x21,y21,z21;
464
465   MyBox1.Get(x10,y10,z10,x11,y11,z11);
466   MyBox2.Get(x20,y20,z20,x21,y21,z21);
467   XMin = 0.;
468   YMin = 0.;
469   ZMin = 0.;
470   XMax = 0.;
471   YMax = 0.;
472   ZMax = 0.;
473
474   if((x10>x21)||(x20>x11)||(y10>y21)||(y20>y11)||(z10>z21)||(z20>z11)) {
475   }
476   else {
477     if(x11<=x21) XMax=x11; else { if(x21<=x11) XMax=x21;}
478     if(x20<=x10) XMin=x10; else { if(x10<=x20) XMin=x20;}
479     if(y11<=y21) YMax=y11; else { if(y21<=y11) YMax=y21;}
480     if(y20<=y10) YMin=y10; else { if(y10<=y20) YMin=y20;}
481     if(z11<=z21) ZMax=z11; else { if(z21<=z11) ZMax=z21;}
482     if(z20<=z10) ZMin=z10; else { if(z10<=z20) ZMin=z20;}
483
484     if(((XMin==XMax)&&(!(YMin==YMax)&&!(ZMin==ZMax)))
485         ||((YMin==YMax)&&(!(XMin==XMax)&&!(ZMin==ZMax)))//ou exclusif ??
486         ||((ZMin==ZMax)&&(!(XMin==XMax)&&!(YMin==YMax)))) {
487     }
488   }
489   //
490   Standard_Real X,Y,Z;
491   X=XMax-XMin; 
492   Y=YMax-YMin; 
493   Z=ZMax-ZMin; 
494   //extension of the box
495   if( (X==0)&&(Y!=0) ) X=Y*0.1;
496   else if( (X==0)&&(Z!=0) ) X=Z*0.1;
497   else X*=0.1;
498
499   if( (Y==0)&&(X!=0) ) Y=X*0.1;
500   else if( (Y==0)&&(Z!=0) ) Y=Z*0.1;
501   else Y*=0.1;
502
503   if( (Z==0)&&(X!=0) ) Z=X*0.1;
504   else if( (Z==0)&&(Y!=0) ) Z=Y*0.1;
505   else Z*=0.1;
506
507
508   if( (X==0)&&(Y==0)&&(Z==0) ) {
509
510
511   }
512   XMin-=X; XMax+=X;
513   YMin-=Y; YMax+=Y;
514   ZMin-=Z; ZMax+=Z;
515
516   //Marking of points included in the common
517   const Standard_Integer FinTP1 = TPoints1.NbItems();
518 //  for(Standard_Integer i=0; i<FinTP1; i++) {
519   Standard_Integer i ;
520   for( i=0; i<FinTP1; i++) {
521     IntPolyh_Point & Pt1 = TPoints1[i];
522     Standard_Integer r;
523     if(Pt1.X()<XMin) { 
524       r=1; 
525     }   
526     else { 
527       if(Pt1.X()>XMax) {
528         r=2;
529       } else {
530         r=0; 
531       } 
532     }
533     if(Pt1.Y()<YMin) { 
534       r|=4; 
535     }  
536     else { 
537       if(Pt1.Y()>YMax) {
538         r|=8; 
539       }
540     }
541     if(Pt1.Z()<ZMin) {
542       r|=16; 
543     } else {
544       if(Pt1.Z()>ZMax) {
545         r|=32;
546       }
547     }
548     Pt1.SetPartOfCommon(r);
549   }
550
551   const Standard_Integer FinTP2 = TPoints2.NbItems();
552   for(Standard_Integer ii=0; ii<FinTP2; ii++) {
553     IntPolyh_Point & Pt2 = TPoints2[ii];
554     Standard_Integer rr;
555     if(Pt2.X()<XMin) { 
556       rr=1;
557     }   
558     else { 
559       if(Pt2.X()>XMax) {
560         rr=2;
561       } else {
562         rr=0;
563       } 
564     }
565     if(Pt2.Y()<YMin) {
566       rr|=4; 
567     }  
568     else { 
569       if(Pt2.Y()>YMax) {
570         rr|=8; 
571       }
572     }
573     if(Pt2.Z()<ZMin) {
574       rr|=16;
575     }
576     else {
577       if(Pt2.Z()>ZMax) {
578         rr|=32;
579       }
580     }
581     Pt2.SetPartOfCommon(rr);
582   }
583 }
584 //=======================================================================
585 //function : FillArrayOfEdges
586 //purpose  : Compute edges from the array of points
587 //           FILL THE ARRAY OF EDGES
588 //=======================================================================
589 void IntPolyh_MaillageAffinage::FillArrayOfEdges
590   (const Standard_Integer SurfID)
591 {
592
593   IntPolyh_ArrayOfEdges &TEdges=(SurfID==1)? TEdges1:TEdges2;
594   Standard_Integer NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
595   Standard_Integer NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
596
597   //NbEdges = 3 + 3*(NbSamplesV-2) + 3*(NbSamplesU-2) + 
598   //        + 3*(NbSamplesU-2)*(NbSamplesV-2) + (NbSamplesV-1) + (NbSamplesU-1);
599   //NbSamplesU and NbSamples cannot be less than 2, so
600   Standard_Integer NbEdges = 3*NbSamplesU*NbSamplesV - 2*(NbSamplesU+NbSamplesV) + 1;
601   TEdges.Init(NbEdges);
602
603   Standard_Integer CpteurTabEdges=0;
604
605   //maillage u0 v0
606   TEdges[CpteurTabEdges].SetFirstPoint(0);                // U V
607   TEdges[CpteurTabEdges].SetSecondPoint(1);             // U V+1
608   //  TEdges[CpteurTabEdges].SetFirstTriangle(-1);
609   TEdges[CpteurTabEdges].SetSecondTriangle(0);
610   CpteurTabEdges++;
611   
612   TEdges[CpteurTabEdges].SetFirstPoint(0);                // U V
613   TEdges[CpteurTabEdges].SetSecondPoint(NbSamplesV);    // U+1 V
614   TEdges[CpteurTabEdges].SetFirstTriangle(0);
615   TEdges[CpteurTabEdges].SetSecondTriangle(1);
616   CpteurTabEdges++;
617   
618   TEdges[CpteurTabEdges].SetFirstPoint(0);                // U V
619   TEdges[CpteurTabEdges].SetSecondPoint(NbSamplesV+1);  // U+1 V+1
620   TEdges[CpteurTabEdges].SetFirstTriangle(1);
621   //  TEdges[CpteurTabEdges].SetSecondTriangle(-1);
622   CpteurTabEdges++;
623   
624   //maillage surU=u0
625   Standard_Integer PntInit=1;
626   Standard_Integer BoucleMeshV;
627   for(BoucleMeshV=1; BoucleMeshV<NbSamplesV-1;BoucleMeshV++){
628     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
629     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1);             // U V+1
630     //    TEdges[CpteurTabEdges].SetFirstTriangle(-1);
631     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2);
632     CpteurTabEdges++;
633     
634     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
635     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV+1);    // U+1 V+1
636     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*2);
637     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2+1);
638     CpteurTabEdges++;
639     
640     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
641     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV);  // U+1 V
642     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*2+1);
643     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2-2);
644     CpteurTabEdges++;
645     PntInit++;
646   }  
647         
648   //maillage sur V=v0
649   PntInit=NbSamplesV;
650   for(BoucleMeshV=1; BoucleMeshV<NbSamplesU-1;BoucleMeshV++){      
651     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);    // U V
652     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1); // U V+1
653     TEdges[CpteurTabEdges].SetFirstTriangle((BoucleMeshV-1)*(NbSamplesV-1)*2+1);
654     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*(NbSamplesV-1)*2);
655     CpteurTabEdges++;
656     
657     TEdges[CpteurTabEdges].SetFirstPoint(PntInit); // U V
658     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV+1);     // U+1 V+1
659     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*(NbSamplesV-1)*2);
660     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*(NbSamplesV-1)*2+1);
661     CpteurTabEdges++;
662     
663     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);  // U V
664     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV);    // U+1 V
665     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*(NbSamplesV-1)*2+1);
666     //    TEdges[CpteurTabEdges].SetSecondTriangle(-1);
667     CpteurTabEdges++;
668     PntInit+=NbSamplesV;
669   }
670       
671   PntInit=NbSamplesV+1;
672   //To provide recursion I associate a point with three edges  
673   for(Standard_Integer BoucleMeshU=1; BoucleMeshU<NbSamplesU-1; BoucleMeshU++){
674     for(Standard_Integer BoucleMeshV=1; BoucleMeshV<NbSamplesV-1;BoucleMeshV++){
675       TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
676       TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1);             // U V+1
677       TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesV-1)*2*(BoucleMeshU-1)+BoucleMeshV*2+1);
678       TEdges[CpteurTabEdges].SetSecondTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2);
679       CpteurTabEdges++;
680       
681       TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
682       TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV+1);    // U+1 V+1 
683       TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2);
684       TEdges[CpteurTabEdges].SetSecondTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2+1);
685       CpteurTabEdges++;
686       
687       TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
688       TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV);  // U+1 V
689       TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2+1);
690       TEdges[CpteurTabEdges].SetSecondTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2-2);
691       CpteurTabEdges++;     
692       PntInit++;//Pass to the next point
693     }
694     PntInit++;//Pass the last point of the column
695     PntInit++;//Pass the first point of the next column
696   }
697         
698   //close mesh on U=u1
699   PntInit=(NbSamplesU-1)*NbSamplesV; //point U=u1 V=0
700   for(BoucleMeshV=0; BoucleMeshV<NbSamplesV-1; BoucleMeshV++){
701     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);           //U=u1 V
702     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1);        //U=u1 V+1
703     TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesU-2)*(NbSamplesV-1)*2+BoucleMeshV*2+1);
704     //    TEdges[CpteurTabEdges].SetSecondTriangle(-1);
705     CpteurTabEdges++;
706     PntInit++;
707   }
708   
709   //close mesh on V=v1
710   for(BoucleMeshV=0; BoucleMeshV<NbSamplesU-1;BoucleMeshV++){      
711     TEdges[CpteurTabEdges].SetFirstPoint(NbSamplesV-1+BoucleMeshV*NbSamplesV);       // U V=v1
712     TEdges[CpteurTabEdges].SetSecondPoint(NbSamplesV-1+(BoucleMeshV+1)*NbSamplesV);  //U+1 V=v1
713     //    TEdges[CpteurTabEdges].SetFirstTriangle(-1);
714     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2*(NbSamplesV-1)+(NbSamplesV-2)*2);
715     CpteurTabEdges++;
716   }
717   TEdges.SetNbItems(CpteurTabEdges);
718 }
719
720 //=======================================================================
721 //function : FillArrayOfTriangles
722 //purpose  : Compute triangles from the array of points, and --
723 //           mark the triangles  that use marked points by the
724 //           CommonBox function.
725 //           FILL THE ARRAY OF TRIANGLES
726 //=======================================================================
727 void IntPolyh_MaillageAffinage::FillArrayOfTriangles
728   (const Standard_Integer SurfID)
729 {
730   Standard_Integer CpteurTabTriangles=0;
731   Standard_Integer PntInit=0;
732
733   IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
734   IntPolyh_ArrayOfTriangles &TTriangles=(SurfID==1)? TTriangles1:TTriangles2;
735   Standard_Integer NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
736   Standard_Integer NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
737
738   TTriangles.Init(2*(NbSamplesU-1)*(NbSamplesV-1));
739   //To provide recursion, I associate a point with two triangles  
740   for(Standard_Integer BoucleMeshU=0; BoucleMeshU<NbSamplesU-1; BoucleMeshU++){
741     for(Standard_Integer BoucleMeshV=0; BoucleMeshV<NbSamplesV-1;BoucleMeshV++){
742       
743       // FIRST TRIANGLE
744       TTriangles[CpteurTabTriangles].SetFirstPoint(PntInit);               // U V
745       TTriangles[CpteurTabTriangles].SetSecondPoint(PntInit+1);            // U V+1
746       TTriangles[CpteurTabTriangles].SetThirdPoint(PntInit+NbSamplesV+1); // U+1 V+1
747
748       // IF ITS EDGE CONTACTS WITH THE COMMON BOX IP REMAINS = A 1
749       if( ( (TPoints[PntInit].PartOfCommon()) & (TPoints[PntInit+1].PartOfCommon()) )
750          &&( (TPoints[PntInit+1].PartOfCommon()) & (TPoints[PntInit+NbSamplesV+1].PartOfCommon()))
751          &&( (TPoints[PntInit+NbSamplesV+1].PartOfCommon()) & (TPoints[PntInit].PartOfCommon())) ) 
752         //IF NOT IP=0
753         TTriangles[CpteurTabTriangles].SetIndiceIntersectionPossible(0);
754
755       CpteurTabTriangles++;
756
757       //SECOND TRIANGLE
758       TTriangles[CpteurTabTriangles].SetFirstPoint(PntInit);               // U V
759       TTriangles[CpteurTabTriangles].SetSecondPoint(PntInit+NbSamplesV+1); // U+1 V+1
760       TTriangles[CpteurTabTriangles].SetThirdPoint(PntInit+NbSamplesV);    // U+1 V 
761
762
763       if( ( (TPoints[PntInit].PartOfCommon()) & (TPoints[PntInit+NbSamplesV+1].PartOfCommon()) )
764          &&( (TPoints[PntInit+NbSamplesV+1].PartOfCommon()) & (TPoints[PntInit+NbSamplesV].PartOfCommon()))
765          &&( (TPoints[PntInit+NbSamplesV].PartOfCommon()) & (TPoints[PntInit].PartOfCommon())) ) 
766         TTriangles[CpteurTabTriangles].SetIndiceIntersectionPossible(0);
767
768
769       CpteurTabTriangles++;
770
771       PntInit++;//Pass to the next point
772     }
773     PntInit++;//Pass the last point of the column
774   }
775   TTriangles.SetNbItems(CpteurTabTriangles);
776   const Standard_Integer FinTT = TTriangles.NbItems();
777   if (FinTT==0) {
778   }
779 }
780 //=======================================================================
781 //function : LinkEdges2Triangles
782 //purpose  : fill the  edge fields in  Triangle object  for the
783 //           two array of triangles.
784 //=======================================================================
785 void IntPolyh_MaillageAffinage::LinkEdges2Triangles() 
786 {
787   const Standard_Integer FinTT1 = TTriangles1.NbItems();
788   const Standard_Integer FinTT2 = TTriangles2.NbItems();
789
790   for(Standard_Integer uiui1=0; uiui1<FinTT1; uiui1++) {
791     IntPolyh_Triangle & MyTriangle1=TTriangles1[uiui1];
792     if ( (MyTriangle1.FirstEdge()) == -1 ) {
793       MyTriangle1.SetEdgeandOrientation(1,TEdges1);
794       MyTriangle1.SetEdgeandOrientation(2,TEdges1);
795       MyTriangle1.SetEdgeandOrientation(3,TEdges1);
796     }
797   }
798   for(Standard_Integer uiui2=0; uiui2<FinTT2; uiui2++) {
799     IntPolyh_Triangle & MyTriangle2=TTriangles2[uiui2];
800     if ( (MyTriangle2.FirstEdge()) == -1 ) {
801       MyTriangle2.SetEdgeandOrientation(1,TEdges2);
802       MyTriangle2.SetEdgeandOrientation(2,TEdges2);
803       MyTriangle2.SetEdgeandOrientation(3,TEdges2);
804     }
805   }
806 }
807 //=======================================================================
808 //function : CommonPartRefinement
809 //purpose  : Refine systematicaly all marked triangles of both surfaces
810 //           REFINING OF THE COMMON
811 //=======================================================================
812 void IntPolyh_MaillageAffinage::CommonPartRefinement() 
813 {
814   Standard_Integer FinInit1 = TTriangles1.NbItems();
815   for(Standard_Integer i=0; i<FinInit1; i++) {
816     if(TTriangles1[i].IndiceIntersectionPossible()!=0) 
817       TTriangles1[i].MiddleRefinement(i,MaSurface1,TPoints1,TTriangles1,TEdges1);
818   }
819
820   Standard_Integer FinInit2=TTriangles2.NbItems();
821   for(Standard_Integer ii=0; ii<FinInit2; ii++) {
822     if(TTriangles2[ii].IndiceIntersectionPossible()!=0) 
823       TTriangles2[ii].MiddleRefinement(ii,MaSurface2,TPoints2,TTriangles2,TEdges2); 
824   }
825
826 }
827 //=======================================================================
828 //function : LocalSurfaceRefinement
829 //purpose  : Refine systematicaly all marked triangles of ONE surface
830 //=======================================================================
831 void IntPolyh_MaillageAffinage::LocalSurfaceRefinement(const Standard_Integer SurfID) {
832 //refine locally, but systematically the chosen surface
833   if (SurfID==1) {
834     const Standard_Integer FinInit1 = TTriangles1.NbItems();
835     for(Standard_Integer i=0; i<FinInit1; i++) {
836       if(TTriangles1[i].IndiceIntersectionPossible()!=0)
837         TTriangles1[i].MiddleRefinement(i,MaSurface1,TPoints1,TTriangles1,TEdges1);
838     }
839   }
840   //
841   if (SurfID==2) {
842     const Standard_Integer FinInit2 = TTriangles2.NbItems();
843     for(Standard_Integer ii=0; ii<FinInit2; ii++) {
844       if(TTriangles2[ii].IndiceIntersectionPossible()!=0) 
845         TTriangles2[ii].MiddleRefinement(ii,MaSurface2,TPoints2,TTriangles2,TEdges2); 
846     }
847   }
848 }
849 //=======================================================================
850 //function : ComputeDeflections
851 //purpose  : Compute deflection  for   all  triangles  of  one
852 //           surface,and sort min and max of deflections
853 //           REFINING PART
854 //           Calculation of the deflection of all triangles
855 //           --> deflection max
856 //           --> deflection min
857 //=======================================================================
858 void IntPolyh_MaillageAffinage::ComputeDeflections
859   (const Standard_Integer SurfID)
860 {
861   Handle(Adaptor3d_HSurface) MaSurface=(SurfID==1)? MaSurface1:MaSurface2;
862   IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
863   IntPolyh_ArrayOfTriangles &TTriangles=(SurfID==1)? TTriangles1:TTriangles2;
864   Standard_Real &FlecheMin=(SurfID==1)? FlecheMin1:FlecheMin2;
865   Standard_Real &FlecheMoy=(SurfID==1)? FlecheMoy1:FlecheMoy2;
866   Standard_Real &FlecheMax=(SurfID==1)? FlecheMax1:FlecheMax2;
867
868   Standard_Integer CpteurTabFleche=0;
869   FlecheMax=-RealLast();
870   FlecheMin=RealLast();
871   FlecheMoy=0.0;
872   const Standard_Integer FinTT = TTriangles.NbItems();
873   
874   for(CpteurTabFleche=0; CpteurTabFleche<FinTT; CpteurTabFleche++) {
875     IntPolyh_Triangle &Triangle = TTriangles[CpteurTabFleche];
876     if ( Triangle.GetFleche() < 0) { //pas normal
877
878     }
879     else{
880       Triangle.TriangleDeflection(MaSurface, TPoints);
881       Standard_Real Fleche=Triangle.GetFleche();
882
883       if (Fleche > FlecheMax)
884         FlecheMax=Fleche;
885       if (Fleche < FlecheMin)
886         FlecheMin=Fleche;
887     }
888   }
889 }
890 //=======================================================================
891 //function : TrianglesDeflectionsRefinementBSB
892 //purpose  : Refine  both  surfaces using  BoundSortBox  as --
893 //           rejection.  The  criterions  used to refine a  --
894 //           triangle are:  The deflection The  size of the --
895 //           bounding boxes   (one surface may be   very small
896 //           compared to the other)
897 //=======================================================================
898 void IntPolyh_MaillageAffinage::TrianglesDeflectionsRefinementBSB() 
899 {
900   const Standard_Integer FinTT1 = TTriangles1.NbItems();
901   const Standard_Integer FinTT2 = TTriangles2.NbItems();
902   
903   ComputeDeflections(1);
904   // To estimate a surface in general it can be interesting 
905   //to calculate all deflections
906   //-- Check deflection at output
907                                                                  
908   Standard_Real FlecheCritique1;
909   if(FlecheMin1>FlecheMax1) {
910     return;    
911   }
912   else {//fleche min + (flechemax-flechemin) * 80/100
913     FlecheCritique1 = FlecheMin1*0.2+FlecheMax1*0.8;
914   }
915   
916   ComputeDeflections(2);
917   //-- Check arrows at output
918   
919   Standard_Real FlecheCritique2;
920   if(FlecheMin2>FlecheMax2) { 
921
922     return;
923   }
924   else {//fleche min + (flechemax-flechemin) * 80/100
925     FlecheCritique2 = FlecheMin2*0.2+FlecheMax2*0.8;
926   }
927   
928   //Bounding boxes
929   Bnd_BoundSortBox BndBSB;
930   Standard_Real diag1,diag2;
931   Standard_Real x0,y0,z0,x1,y1,z1;
932   
933   //The greatest of two bounding boxes created in FillArrayOfPoints is found.
934   //Then this value is weighted depending on the discretization 
935   //(NbSamplesU and NbSamplesV)
936   MyBox1.Get(x0,y0,z0,x1,y1,z1);
937   x0-=x1; y0-=y1; z0-=z1;
938   diag1=x0*x0+y0*y0+z0*z0;
939   const Standard_Real NbSamplesUV1=Standard_Real(NbSamplesU1) * Standard_Real(NbSamplesV1);
940   diag1/=NbSamplesUV1;
941
942   MyBox2.Get(x0,y0,z0,x1,y1,z1);
943   x0-=x1; y0-=y1; z0-=z1;
944   diag2=x0*x0+y0*y0+z0*z0;
945   const Standard_Real NbSamplesUV2=Standard_Real(NbSamplesU2) * Standard_Real(NbSamplesV2);
946   diag2/=NbSamplesUV2;
947   
948   //-- The surface with the greatest bounding box is "discretized"
949   
950   //Standard_Integer NbInterTentees=0;
951   
952   if(diag1<diag2) { 
953     
954     if(FlecheCritique2<diag1) {//the corresponding sizes are not too disproportional
955
956       Handle(Bnd_HArray1OfBox) HBnd = new  Bnd_HArray1OfBox(1,FinTT2);
957         
958       for(Standard_Integer i=0; i<FinTT2; i++){
959         if (TTriangles2[i].IndiceIntersectionPossible()!=0) {
960           Bnd_Box b;
961           const IntPolyh_Triangle& T=TTriangles2[i];
962           const IntPolyh_Point&    PA=TPoints2[T.FirstPoint()];
963           const IntPolyh_Point&    PB=TPoints2[T.SecondPoint()]; 
964           const IntPolyh_Point&    PC=TPoints2[T.ThirdPoint()]; 
965           gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
966           gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
967           gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
968           b.Add(pntA);//Box b, which contains triangle i of surface 2 is created./
969           b.Add(pntB);
970           b.Add(pntC);
971           b.Enlarge(T.GetFleche()+MyTolerance);
972           HBnd->SetValue(i+1,b);//Box b is added in the array HBnd
973         }
974       }
975       
976       //Inititalization of the boundary, sorting of boxes
977       BndBSB.Initialize(HBnd);//contains boxes of 2
978       
979       Standard_Integer FinTT1Init=FinTT1;
980       for(Standard_Integer i_S1=0; i_S1<FinTT1Init; i_S1++) {
981         if(TTriangles1[i_S1].IndiceIntersectionPossible()!=0) {
982           //-- Loop on the boxes of mesh 1 
983           Bnd_Box b;
984           const IntPolyh_Triangle& T=TTriangles1[i_S1];
985           const IntPolyh_Point&    PA=TPoints1[T.FirstPoint()]; 
986           const IntPolyh_Point&    PB=TPoints1[T.SecondPoint()]; 
987           const IntPolyh_Point&    PC=TPoints1[T.ThirdPoint()]; 
988           gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
989           gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
990           gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
991           b.Add(pntA);
992           b.Add(pntB);
993           b.Add(pntC);
994           b.Enlarge(T.GetFleche());
995           //-- List of boxes of 2, which touch this box (of 1)
996           const TColStd_ListOfInteger& ListeOf2 = BndBSB.Compare(b);
997           
998           if((ListeOf2.IsEmpty())==0) {
999             IntPolyh_Triangle &Triangle1 = TTriangles1[i_S1];
1000             if(Triangle1.GetFleche()>FlecheCritique1)
1001               Triangle1.MiddleRefinement(i_S1,MaSurface1,TPoints1,
1002                                          TTriangles1, TEdges1);
1003             
1004             for (TColStd_ListIteratorOfListOfInteger Iter(ListeOf2); 
1005                  Iter.More(); 
1006                  Iter.Next()) {
1007               Standard_Integer i_S2=Iter.Value()-1;
1008               //if the box of s1 contacts with the boxes of s2 
1009               //the arrow of the triangle is checked
1010               IntPolyh_Triangle & Triangle2 = TTriangles2[i_S2];
1011               if(Triangle2.IndiceIntersectionPossible()!=0)
1012                 if(Triangle2.GetFleche()>FlecheCritique2)
1013                   Triangle2.MiddleRefinement( i_S2, MaSurface2, TPoints2,
1014                                              TTriangles2, TEdges2);
1015             }
1016           }
1017         }
1018       }
1019     }
1020
1021     //--------------------------------------------------------------------
1022     //FlecheCritique2 > diag1
1023     else {
1024       //2 is discretized
1025
1026       Handle(Bnd_HArray1OfBox) HBnd = new  Bnd_HArray1OfBox(1,FinTT2);
1027     
1028       for(Standard_Integer i=0; i<FinTT2; i++){
1029         if (TTriangles2[i].IndiceIntersectionPossible()!=0) {
1030           Bnd_Box b;
1031           const IntPolyh_Triangle& T=TTriangles2[i];
1032           const IntPolyh_Point&    PA=TPoints2[T.FirstPoint()];
1033           const IntPolyh_Point&    PB=TPoints2[T.SecondPoint()]; 
1034           const IntPolyh_Point&    PC=TPoints2[T.ThirdPoint()]; 
1035           gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
1036           gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
1037           gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
1038           b.Add(pntA);//Box b, which contains triangle i of surface 2 is created/
1039           b.Add(pntB);
1040           b.Add(pntC);
1041           b.Enlarge(T.GetFleche()+MyTolerance);
1042           //-- BndBSB.Add(b,i+1);
1043           HBnd->SetValue(i+1,b);//Box b is added in array HBnd
1044         }
1045       }
1046       
1047       //Inititalization of the ouput bounding box
1048       BndBSB.Initialize(HBnd);//contains boxes of 2
1049       
1050
1051       //The bounding box Be1 of surface1 is compared BSB of surface2
1052       const TColStd_ListOfInteger& ListeOf2 = BndBSB.Compare(MyBox1);
1053       
1054       if((ListeOf2.IsEmpty())==0) {
1055         //if the bounding box Be1 of s1 contacts with 
1056         //the boxes of s2 the deflection of triangle of s2 is checked
1057
1058         // Be1 is very small in relation to Be2
1059         //The criterion of refining for surface2 depends on the size of Be1
1060         //As it is known that this criterion should be minimized, 
1061         //the smallest side of the bounding box is taken
1062         Standard_Real x0,x1,y0,y1,z0,z1;
1063         MyBox1.Get(x0,y0,z0,x1,y1,z1);
1064         Standard_Real dx=Abs(x1-x0);
1065         Standard_Real dy=Abs(y1-y0);
1066         Standard_Real diag=Abs(z1-z0);
1067         Standard_Real dd=-1.0;
1068         if (dx>dy)
1069           dd=dy;
1070         else
1071           dd=dx;
1072         if (diag>dd) diag=dd;
1073         
1074         //if Be1 contacts with the boxes of s2, the deflection 
1075         //of the triangles of s2 is checked (greater)
1076         //in relation to the size of Be1 (smaller)
1077         for (TColStd_ListIteratorOfListOfInteger Iter(ListeOf2); 
1078              Iter.More(); 
1079              Iter.Next()) {
1080           Standard_Integer i_S2=Iter.Value()-1;
1081           
1082           IntPolyh_Triangle & Triangle2=TTriangles2[i_S2];
1083           if(Triangle2.IndiceIntersectionPossible()) {
1084             
1085             //calculation of the criterion of refining
1086             //The deflection of the greater is compared to the size of the smaller
1087             Standard_Real CritereAffinage=0.0;
1088             Standard_Real DiagPonderation=0.5;
1089             CritereAffinage = diag*DiagPonderation;
1090             if(Triangle2.GetFleche()>CritereAffinage)
1091               Triangle2.MultipleMiddleRefinement2(CritereAffinage, MyBox1, i_S2,
1092                                                   MaSurface2, TPoints2,
1093                                                   TTriangles2,TEdges2);
1094             
1095             else Triangle2.MiddleRefinement(i_S2,MaSurface2,TPoints2,
1096                                             TTriangles2, TEdges2);
1097           }
1098         }
1099       }
1100     }
1101   }
1102   
1103   
1104   else {     //-- The greater is discretised
1105
1106     if(FlecheCritique1<diag2) {//the respective sizes are not to much disproportional
1107     
1108       Handle(Bnd_HArray1OfBox) HBnd = new  Bnd_HArray1OfBox(1,FinTT1);
1109       
1110       for(Standard_Integer i=0; i<FinTT1; i++){
1111         if(TTriangles1[i].IndiceIntersectionPossible()!=0) {
1112           Bnd_Box b;
1113           const IntPolyh_Triangle& T=TTriangles1[i];
1114           const IntPolyh_Point&    PA=TPoints1[T.FirstPoint()];
1115           const IntPolyh_Point&    PB=TPoints1[T.SecondPoint()]; 
1116           const IntPolyh_Point&    PC=TPoints1[T.ThirdPoint()]; 
1117           gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
1118           gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
1119           gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
1120           b.Add(pntA);//Box b, which contains triangle i of surface 2 is created.
1121           b.Add(pntB);
1122           b.Add(pntC);
1123           b.Enlarge(T.GetFleche()+MyTolerance);
1124           HBnd->SetValue(i+1,b);//Boite b is added in the array HBnd
1125         }
1126       }
1127       BndBSB.Initialize(HBnd);
1128       
1129       Standard_Integer FinTT2init=FinTT2;
1130       for(Standard_Integer i_S2=0; i_S2<FinTT2init; i_S2++) {
1131         if (TTriangles2[i_S2].IndiceIntersectionPossible()!=0) {      
1132           //-- Loop on the boxes of mesh 2 
1133           Bnd_Box b;
1134           const IntPolyh_Triangle& T=TTriangles2[i_S2];
1135           const IntPolyh_Point&    PA=TPoints2[T.FirstPoint()]; 
1136           const IntPolyh_Point&    PB=TPoints2[T.SecondPoint()]; 
1137           const IntPolyh_Point&    PC=TPoints2[T.ThirdPoint()]; 
1138           gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
1139           gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
1140           gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
1141           b.Add(pntA);
1142           b.Add(pntB);
1143           b.Add(pntC);
1144           b.Enlarge(T.GetFleche()+MyTolerance);
1145           //-- List of boxes of 1 touching this box (of 2)
1146           const TColStd_ListOfInteger& ListeOf1 = BndBSB.Compare(b);
1147           IntPolyh_Triangle & Triangle2=TTriangles2[i_S2];
1148           if((ListeOf1.IsEmpty())==0) {
1149             
1150             if(Triangle2.GetFleche()>FlecheCritique2)
1151               Triangle2.MiddleRefinement(i_S2,MaSurface2,TPoints2,
1152                                          TTriangles2, TEdges2);
1153             
1154             for (TColStd_ListIteratorOfListOfInteger Iter(ListeOf1); 
1155                  Iter.More(); 
1156                  Iter.Next()) {
1157               Standard_Integer i_S1=Iter.Value()-1;
1158               IntPolyh_Triangle & Triangle1=TTriangles1[i_S1];
1159               if (Triangle1.IndiceIntersectionPossible())
1160                 if(Triangle1.GetFleche()>FlecheCritique1)
1161                   Triangle1.MiddleRefinement(i_S1,MaSurface1,TPoints1,
1162                                              TTriangles1, TEdges1); 
1163             }
1164           }
1165         }
1166       }
1167     }
1168     //-----------------------------------------------------------------------------
1169     else {// FlecheCritique1>diag2
1170       // 1 is discretized
1171
1172       Handle(Bnd_HArray1OfBox) HBnd = new  Bnd_HArray1OfBox(1,FinTT1);
1173     
1174       for(Standard_Integer i=0; i<FinTT1; i++){
1175         if (TTriangles1[i].IndiceIntersectionPossible()!=0) {
1176           Bnd_Box b;
1177           const IntPolyh_Triangle& T=TTriangles1[i];
1178           const IntPolyh_Point&    PA=TPoints1[T.FirstPoint()];
1179           const IntPolyh_Point&    PB=TPoints1[T.SecondPoint()]; 
1180           const IntPolyh_Point&    PC=TPoints1[T.ThirdPoint()]; 
1181           gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
1182           gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
1183           gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
1184           b.Add(pntA);//Box b, which contains triangle i of surface 1 is created./
1185           b.Add(pntB);
1186           b.Add(pntC);
1187           b.Enlarge(T.GetFleche()+MyTolerance);
1188           HBnd->SetValue(i+1,b);//Box b is added in the array HBnd
1189         }
1190       }
1191       
1192       //Inititalisation of the boundary output box
1193       BndBSB.Initialize(HBnd);//contains boxes of 1
1194       
1195       //Bounding box Be2 of surface2 is compared to BSB of surface1
1196       const TColStd_ListOfInteger& ListeOf1 = BndBSB.Compare(MyBox2);
1197       
1198       if((ListeOf1.IsEmpty())==0) {
1199         //if the bounding box Be2 of s2 contacts 
1200         //with boxes of s1 the deflection of the triangle of s1 is checked
1201         
1202         // Be2 is very small compared to Be1
1203         //The criterion of refining for surface1 depends on the size of Be2
1204         //As this criterion should be minimized, 
1205         //the smallest side of the bounding box is taken
1206         Standard_Real x0,x1,y0,y1,z0,z1;
1207         MyBox2.Get(x0,y0,z0,x1,y1,z1);
1208         Standard_Real dx=Abs(x1-x0);
1209         Standard_Real dy=Abs(y1-y0);
1210         Standard_Real diag=Abs(z1-z0);
1211         Standard_Real dd=-1.0;
1212         if (dx>dy)
1213           dd=dy;
1214         else
1215           dd=dx;
1216         if (diag>dd) diag=dd;
1217         
1218         //if Be2 contacts with boxes of s1, the deflection of 
1219         //triangles of s1 (greater) is checked
1220         //comparatively to the size of Be2 (smaller).
1221         for (TColStd_ListIteratorOfListOfInteger Iter(ListeOf1); 
1222              Iter.More(); 
1223              Iter.Next()) {
1224           Standard_Integer i_S1=Iter.Value()-1;
1225
1226           IntPolyh_Triangle & Triangle1=TTriangles1[i_S1];
1227           if(Triangle1.IndiceIntersectionPossible()) {
1228
1229             //calculation of the criterion of refining
1230             //The deflection of the greater is compared 
1231             //with the size of the smaller.
1232             Standard_Real CritereAffinage=0.0;
1233             Standard_Real DiagPonderation=0.5;
1234             CritereAffinage = diag*DiagPonderation;;
1235             if(Triangle1.GetFleche()>CritereAffinage)
1236               Triangle1.MultipleMiddleRefinement2(CritereAffinage,MyBox2, i_S1,
1237                                                   MaSurface1, TPoints1,
1238                                                   TTriangles1, TEdges1); 
1239             
1240             else Triangle1.MiddleRefinement(i_S1,MaSurface1,TPoints1,
1241                                             TTriangles1, TEdges1);
1242             
1243           }
1244         }
1245       }
1246     }
1247   }
1248 }
1249 //=======================================================================
1250 //function : maxSR
1251 //purpose  : This function is used for the function project6
1252 //=======================================================================
1253 inline Standard_Real maxSR(const Standard_Real a,
1254                            const Standard_Real b,
1255                            const Standard_Real c)
1256 {
1257   Standard_Real t = a;
1258   if (b > t) t = b;
1259   if (c > t) t = c;
1260   return t;
1261 }
1262 //=======================================================================
1263 //function : minSR
1264 //purpose  : This function is used for the function project6
1265 //=======================================================================
1266 inline Standard_Real minSR(const Standard_Real a,
1267                            const Standard_Real b,
1268                            const Standard_Real c)
1269 {
1270   Standard_Real t = a;
1271   if (b < t) t = b;
1272   if (c < t) t = c;
1273   return t;
1274 }
1275 //=======================================================================
1276 //function : project6
1277 //purpose  : This function is used for the function TriContact
1278 //=======================================================================
1279 Standard_Integer project6(const IntPolyh_Point &ax, 
1280                           const IntPolyh_Point &p1,
1281                           const IntPolyh_Point &p2, 
1282                           const IntPolyh_Point &p3, 
1283                           const IntPolyh_Point &q1,
1284                           const IntPolyh_Point &q2, 
1285                           const IntPolyh_Point &q3)
1286 {
1287   Standard_Real P1 = ax.Dot(p1);
1288   Standard_Real P2 = ax.Dot(p2);
1289   Standard_Real P3 = ax.Dot(p3);
1290   Standard_Real Q1 = ax.Dot(q1);
1291   Standard_Real Q2 = ax.Dot(q2);
1292   Standard_Real Q3 = ax.Dot(q3);
1293   
1294   Standard_Real mx1 = maxSR(P1, P2, P3);
1295   Standard_Real mn1 = minSR(P1, P2, P3);
1296   Standard_Real mx2 = maxSR(Q1, Q2, Q3);
1297   Standard_Real mn2 = minSR(Q1, Q2, Q3);
1298   
1299   if (mn1 > mx2) return 0;
1300   if (mn2 > mx1) return 0;
1301   return 1;
1302 }
1303 //=======================================================================
1304 //function : TriContact
1305 //purpose  : This fonction     Check if two triangles   are  in
1306 //           contact or no,  return 1 if yes, return 0
1307 //           if no.
1308 //=======================================================================
1309 Standard_Integer IntPolyh_MaillageAffinage::TriContact
1310   (const IntPolyh_Point &P1,
1311    const IntPolyh_Point &P2,
1312    const IntPolyh_Point &P3,
1313    const IntPolyh_Point &Q1,
1314    const IntPolyh_Point &Q2,
1315    const IntPolyh_Point &Q3,
1316    Standard_Real &Angle) const
1317 {
1318   /**
1319      The first triangle is (p1,p2,p3).  The other is (q1,q2,q3).
1320      The edges are (e1,e2,e3) and (f1,f2,f3).
1321      The normals are n1 and m1
1322      Outwards are (g1,g2,g3) and (h1,h2,h3).*/
1323
1324   IntPolyh_Point p1, p2, p3;
1325   IntPolyh_Point q1, q2, q3;
1326   IntPolyh_Point e1, e2, e3;
1327   IntPolyh_Point f1, f2, f3;
1328   IntPolyh_Point g1, g2, g3;
1329   IntPolyh_Point h1, h2, h3;
1330   IntPolyh_Point n1, m1;
1331   IntPolyh_Point z;
1332
1333   IntPolyh_Point ef11, ef12, ef13;
1334   IntPolyh_Point ef21, ef22, ef23;
1335   IntPolyh_Point ef31, ef32, ef33;
1336   
1337   z.SetX(0.0);  z.SetY(0.0);  z.SetZ(0.0);
1338   
1339   if(maxSR(P1.X(),P2.X(),P3.X())<minSR(Q1.X(),Q2.X(),Q3.X())) return(0);
1340   if(maxSR(P1.Y(),P2.Y(),P3.Y())<minSR(Q1.Y(),Q2.Y(),Q3.Y())) return(0);
1341   if(maxSR(P1.Z(),P2.Z(),P3.Z())<minSR(Q1.Z(),Q2.Z(),Q3.Z())) return(0);
1342
1343   if(minSR(P1.X(),P2.X(),P3.X())>maxSR(Q1.X(),Q2.X(),Q3.X())) return(0);
1344   if(minSR(P1.Y(),P2.Y(),P3.Y())>maxSR(Q1.Y(),Q2.Y(),Q3.Y())) return(0);
1345   if(minSR(P1.Z(),P2.Z(),P3.Z())>maxSR(Q1.Z(),Q2.Z(),Q3.Z())) return(0);
1346
1347   p1.SetX(P1.X() - P1.X());  p1.SetY(P1.Y() - P1.Y());  p1.SetZ(P1.Z() - P1.Z());
1348   p2.SetX(P2.X() - P1.X());  p2.SetY(P2.Y() - P1.Y());  p2.SetZ(P2.Z() - P1.Z());
1349   p3.SetX(P3.X() - P1.X());  p3.SetY(P3.Y() - P1.Y());  p3.SetZ(P3.Z() - P1.Z());
1350   
1351   q1.SetX(Q1.X() - P1.X());  q1.SetY(Q1.Y() - P1.Y());  q1.SetZ(Q1.Z() - P1.Z());
1352   q2.SetX(Q2.X() - P1.X());  q2.SetY(Q2.Y() - P1.Y());  q2.SetZ(Q2.Z() - P1.Z());
1353   q3.SetX(Q3.X() - P1.X());  q3.SetY(Q3.Y() - P1.Y());  q3.SetZ(Q3.Z() - P1.Z());
1354   
1355   e1.SetX(p2.X() - p1.X());  e1.SetY(p2.Y() - p1.Y());  e1.SetZ(p2.Z() - p1.Z());
1356   e2.SetX(p3.X() - p2.X());  e2.SetY(p3.Y() - p2.Y());  e2.SetZ(p3.Z() - p2.Z());
1357   e3.SetX(p1.X() - p3.X());  e3.SetY(p1.Y() - p3.Y());  e3.SetZ(p1.Z() - p3.Z());
1358
1359   f1.SetX(q2.X() - q1.X());  f1.SetY(q2.Y() - q1.Y());  f1.SetZ(q2.Z() - q1.Z());
1360   f2.SetX(q3.X() - q2.X());  f2.SetY(q3.Y() - q2.Y());  f2.SetZ(q3.Z() - q2.Z());
1361   f3.SetX(q1.X() - q3.X());  f3.SetY(q1.Y() - q3.Y());  f3.SetZ(q1.Z() - q3.Z());
1362   
1363   n1.Cross(e1, e2); //normal to the first triangle
1364   m1.Cross(f1, f2); //normal to the second triangle
1365
1366   g1.Cross(e1, n1); 
1367   g2.Cross(e2, n1);
1368   g3.Cross(e3, n1);
1369   h1.Cross(f1, m1);
1370   h2.Cross(f2, m1);
1371   h3.Cross(f3, m1);
1372
1373   ef11.Cross(e1, f1);
1374   ef12.Cross(e1, f2);
1375   ef13.Cross(e1, f3);
1376   ef21.Cross(e2, f1);
1377   ef22.Cross(e2, f2);
1378   ef23.Cross(e2, f3);
1379   ef31.Cross(e3, f1);
1380   ef32.Cross(e3, f2);
1381   ef33.Cross(e3, f3);
1382   
1383   // Now the testing is done
1384
1385   if (!project6(n1, p1, p2, p3, q1, q2, q3)) return 0; //T2 is not higher or lower than T1
1386   if (!project6(m1, p1, p2, p3, q1, q2, q3)) return 0; //T1 is not higher of lower than T2
1387   
1388   if (!project6(ef11, p1, p2, p3, q1, q2, q3)) return 0; 
1389   if (!project6(ef12, p1, p2, p3, q1, q2, q3)) return 0;
1390   if (!project6(ef13, p1, p2, p3, q1, q2, q3)) return 0;
1391   if (!project6(ef21, p1, p2, p3, q1, q2, q3)) return 0;
1392   if (!project6(ef22, p1, p2, p3, q1, q2, q3)) return 0;
1393   if (!project6(ef23, p1, p2, p3, q1, q2, q3)) return 0;
1394   if (!project6(ef31, p1, p2, p3, q1, q2, q3)) return 0;
1395   if (!project6(ef32, p1, p2, p3, q1, q2, q3)) return 0;
1396   if (!project6(ef33, p1, p2, p3, q1, q2, q3)) return 0;
1397
1398   if (!project6(g1, p1, p2, p3, q1, q2, q3)) return 0; //T2 is outside of T1 in the plane of T1
1399   if (!project6(g2, p1, p2, p3, q1, q2, q3)) return 0; //T2 is outside of T1 in the plane of T1
1400   if (!project6(g3, p1, p2, p3, q1, q2, q3)) return 0; //T2 is outside of T1 in the plane of T1
1401   if (!project6(h1, p1, p2, p3, q1, q2, q3)) return 0; //T1 is outside of T2 in the plane of T2
1402   if (!project6(h2, p1, p2, p3, q1, q2, q3)) return 0; //T1 is outside of T2 in the plane of T2
1403   if (!project6(h3, p1, p2, p3, q1, q2, q3)) return 0; //T1 is outside of T2 in the plane of T2
1404
1405   //Calculation of cosinus angle between two normals
1406   Standard_Real SqModn1=-1.0;
1407   Standard_Real SqModm1=-1.0;
1408   SqModn1=n1.SquareModulus();
1409   if (SqModn1>SquareMyConfusionPrecision){ 
1410     SqModm1=m1.SquareModulus();
1411   }
1412   if (SqModm1>SquareMyConfusionPrecision) {
1413     Angle=(n1.Dot(m1))/(sqrt(SqModn1)*sqrt(SqModm1));
1414   }
1415   return 1;
1416 }
1417 //=======================================================================
1418 //function : TestNbPoints
1419 //purpose  : This function is used by StartingPointsResearch() to control 
1420 //           the  number of points found keep the result in conformity (1 or 2 points) 
1421 //           void TestNbPoints(const Standard_Integer TriSurfID,
1422 //=======================================================================
1423 void TestNbPoints(const Standard_Integer ,
1424                   Standard_Integer &NbPoints,
1425                   Standard_Integer &NbPointsTotal,
1426                   const IntPolyh_StartPoint &Pt1,
1427                   const IntPolyh_StartPoint &Pt2,
1428                   IntPolyh_StartPoint &SP1,
1429                   IntPolyh_StartPoint &SP2)
1430 {
1431   // already checked in TriangleEdgeContact2
1432   //  if( (NbPoints==2)&&(Pt1.CheckSameSP(Pt2)) ) NbPoints=1;
1433
1434   if(NbPoints>2) {
1435
1436   }
1437   else {
1438     if ( (NbPoints==1)&&(NbPointsTotal==0) ) { 
1439       SP1=Pt1;
1440       NbPointsTotal=1;
1441     }
1442     else if ( (NbPoints==1)&&(NbPointsTotal==1) ) { 
1443       if(Pt1.CheckSameSP(SP1)!=1) {
1444         SP2=Pt1;
1445         NbPointsTotal=2;
1446       }
1447     }
1448     else if( (NbPoints==1)&&(NbPointsTotal==2) ) {
1449       if ( (SP1.CheckSameSP(Pt1))||(SP2.CheckSameSP(Pt1)) ) 
1450         NbPointsTotal=2;
1451       else NbPointsTotal=3;
1452     }
1453     else if( (NbPoints==2)&&(NbPointsTotal==0) ) {
1454       SP1=Pt1;
1455       SP2=Pt2;
1456       NbPointsTotal=2;
1457     }
1458     else if( (NbPoints==2)&&(NbPointsTotal==1) ) {//there is also Pt1 != Pt2 
1459       if(SP1.CheckSameSP(Pt1)) {
1460         SP2=Pt2;
1461         NbPointsTotal=2;
1462       }
1463       else if (SP1.CheckSameSP(Pt2)) {
1464         SP2=Pt1;
1465         NbPointsTotal=2;
1466       }
1467       else NbPointsTotal=3;///case SP1!=Pt1 && SP1!=Pt2!
1468     }
1469     else if( (NbPoints==2)&&(NbPointsTotal==2) ) {//there is also SP1!=SP2
1470       if( (SP1.CheckSameSP(Pt1))||(SP1.CheckSameSP(Pt2)) ) {
1471         if( (SP2.CheckSameSP(Pt1))||(SP2.CheckSameSP(Pt2)) )
1472           NbPointsTotal=2;
1473         else NbPointsTotal=3;
1474       }
1475       else NbPointsTotal=3;
1476     }
1477   }
1478 }
1479 //=======================================================================
1480 //function : StartingPointsResearch
1481 //purpose  : 
1482 //=======================================================================
1483 Standard_Integer IntPolyh_MaillageAffinage::StartingPointsResearch
1484   (const Standard_Integer T1,
1485    const Standard_Integer T2,
1486    IntPolyh_StartPoint &SP1, 
1487    IntPolyh_StartPoint &SP2) const 
1488 {
1489   const IntPolyh_Point &P1=TPoints1[TTriangles1[T1].FirstPoint()];
1490   const IntPolyh_Point &P2=TPoints1[TTriangles1[T1].SecondPoint()];
1491   const IntPolyh_Point &P3=TPoints1[TTriangles1[T1].ThirdPoint()];
1492   const IntPolyh_Point &Q1=TPoints2[TTriangles2[T2].FirstPoint()];
1493   const IntPolyh_Point &Q2=TPoints2[TTriangles2[T2].SecondPoint()];
1494   const IntPolyh_Point &Q3=TPoints2[TTriangles2[T2].ThirdPoint()];
1495
1496
1497   /* The first triangle is (p1,p2,p3).  The other is (q1,q2,q3).
1498      The sides are (e1,e2,e3) and (f1,f2,f3).
1499      The normals are n1 and m1*/
1500
1501   const IntPolyh_Point  e1=P2-P1;
1502   const IntPolyh_Point  e2=P3-P2;
1503   const IntPolyh_Point  e3=P1-P3;
1504
1505   const IntPolyh_Point  f1=Q2-Q1;
1506   const IntPolyh_Point  f2=Q3-Q2;
1507   const IntPolyh_Point  f3=Q1-Q3;
1508   
1509
1510   IntPolyh_Point nn1,mm1;
1511   nn1.Cross(e1, e2); //normal of the first triangle
1512   mm1.Cross(f1, f2); //normal of the  second triangle
1513
1514   Standard_Real nn1modulus, mm1modulus;
1515   nn1modulus=sqrt(nn1.SquareModulus());
1516   mm1modulus=sqrt(mm1.SquareModulus());
1517
1518   //-------------------------------------------------------
1519   ///calculation of intersection points between two triangles
1520   //-------------------------------------------------------
1521   Standard_Integer NbPoints=0;
1522   Standard_Integer NbPointsTotal=0;
1523   IntPolyh_StartPoint Pt1,Pt2;
1524
1525
1526     ///check T1 normal
1527     if(Abs(nn1modulus)<MyConfusionPrecision){//10.0e-20) {
1528
1529     }
1530     else {
1531       const IntPolyh_Point n1=nn1.Divide(nn1modulus);
1532       ///T2 edges with T1
1533       if(NbPointsTotal<2) {
1534         NbPoints=TriangleEdgeContact(1,1,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
1535         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1536       }
1537       
1538       if(NbPointsTotal<2) {
1539         NbPoints=TriangleEdgeContact(1,2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
1540         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1541       }
1542       
1543       if(NbPointsTotal<2) {
1544         NbPoints=TriangleEdgeContact(1,3,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
1545         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1546       }
1547     }
1548
1549     ///check T2 normal
1550     if(Abs(mm1modulus)<MyConfusionPrecision) {//10.0e-20){
1551
1552     }
1553     else {
1554       const IntPolyh_Point m1=mm1.Divide(mm1modulus);
1555       ///T1 edges with T2  
1556       if(NbPointsTotal<2) {
1557         NbPoints=TriangleEdgeContact(2,1,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
1558         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1559       }
1560       
1561       if(NbPointsTotal<2) { 
1562         NbPoints=TriangleEdgeContact(2,2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
1563         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1564       }
1565       
1566       if(NbPointsTotal<2) {
1567         NbPoints=TriangleEdgeContact(2,3,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
1568         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1569       }
1570     }
1571
1572   /*  if( (NbPointsTotal >1)&&( Abs(SP1.U1()-SP2.U1())<MyConfusionPrecision)
1573       &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) )*/
1574   if( (NbPoints)&&(SP1.CheckSameSP(SP2)) )
1575     NbPointsTotal=1;
1576   
1577   SP1.SetCoupleValue(T1,T2);
1578   SP2.SetCoupleValue(T1,T2);
1579   return (NbPointsTotal);
1580 }
1581 //=======================================================================
1582 //function : StartingPointsResearch2
1583 //purpose  : From  two  triangles compute intersection  points.
1584 //           If I found   more  than two intersection  points
1585 //           it  means that those triangle are coplanar
1586 //=======================================================================
1587 Standard_Integer IntPolyh_MaillageAffinage::StartingPointsResearch2
1588   (const Standard_Integer T1,
1589    const Standard_Integer T2,
1590    IntPolyh_StartPoint &SP1, 
1591    IntPolyh_StartPoint &SP2) const 
1592 {
1593   const IntPolyh_Triangle &Tri1=TTriangles1[T1];
1594   const IntPolyh_Triangle &Tri2=TTriangles2[T2];
1595
1596   const IntPolyh_Point &P1=TPoints1[Tri1.FirstPoint()];
1597   const IntPolyh_Point &P2=TPoints1[Tri1.SecondPoint()];
1598   const IntPolyh_Point &P3=TPoints1[Tri1.ThirdPoint()];
1599   const IntPolyh_Point &Q1=TPoints2[Tri2.FirstPoint()];
1600   const IntPolyh_Point &Q2=TPoints2[Tri2.SecondPoint()];
1601   const IntPolyh_Point &Q3=TPoints2[Tri2.ThirdPoint()];
1602
1603
1604
1605   /* The first triangle is (p1,p2,p3).  The other is (q1,q2,q3).
1606      The sides are (e1,e2,e3) and (f1,f2,f3).
1607      The normals are n1 and m1*/
1608
1609   const IntPolyh_Point  e1=P2-P1;
1610   const IntPolyh_Point  e2=P3-P2;
1611   const IntPolyh_Point  e3=P1-P3;
1612
1613   const IntPolyh_Point  f1=Q2-Q1;
1614   const IntPolyh_Point  f2=Q3-Q2;
1615   const IntPolyh_Point  f3=Q1-Q3;
1616   
1617
1618   IntPolyh_Point nn1,mm1;
1619   nn1.Cross(e1, e2); //normal to the first triangle
1620   mm1.Cross(f1, f2); //normal to the second triangle
1621
1622   Standard_Real nn1modulus, mm1modulus;
1623   nn1modulus=sqrt(nn1.SquareModulus());
1624   mm1modulus=sqrt(mm1.SquareModulus());
1625
1626   //-------------------------------------------------
1627   ///calculation of intersections points between triangles
1628   //-------------------------------------------------
1629   Standard_Integer NbPoints=0;
1630   Standard_Integer NbPointsTotal=0;
1631
1632
1633     ///check T1 normal
1634     if(Abs(nn1modulus)<MyConfusionPrecision) {//10.0e-20){
1635
1636     }
1637     else {
1638       const IntPolyh_Point n1=nn1.Divide(nn1modulus);
1639       ///T2 edges with T1
1640       if(NbPointsTotal<3) {
1641         IntPolyh_StartPoint Pt1,Pt2;
1642         NbPoints=TriangleEdgeContact2(1,1,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
1643         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1644       }
1645       
1646       if(NbPointsTotal<3) {
1647         IntPolyh_StartPoint Pt1,Pt2;
1648         NbPoints=TriangleEdgeContact2(1,2,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
1649         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1650       }
1651       
1652       if(NbPointsTotal<3) {
1653         IntPolyh_StartPoint Pt1,Pt2;
1654         NbPoints=TriangleEdgeContact2(1,3,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
1655         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1656       }
1657     }
1658
1659     ///check T2 normal
1660     if(Abs(mm1modulus)<MyConfusionPrecision) {//10.0e-20){
1661
1662     }
1663     else {
1664       const IntPolyh_Point m1=mm1.Divide(mm1modulus);
1665       ///T1 edges with T2
1666       if(NbPointsTotal<3) {
1667         IntPolyh_StartPoint Pt1,Pt2;
1668         NbPoints=TriangleEdgeContact2(2,1,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
1669         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1670       }
1671       
1672       if(NbPointsTotal<3) { 
1673         IntPolyh_StartPoint Pt1,Pt2;
1674         NbPoints=TriangleEdgeContact2(2,2,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
1675         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1676       }
1677       
1678       if(NbPointsTotal<3) {
1679         IntPolyh_StartPoint Pt1,Pt2;
1680         NbPoints=TriangleEdgeContact2(2,3,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
1681         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1682       }
1683     }
1684   //  no need already checked in  TestNbPoints
1685   //  if( (NbPointsTotal==2)&&(SP1.CheckSameSP(SP2)) ) {
1686   //  NbPointsTotal=1;
1687   //SP1.SetCoupleValue(T1,T2);
1688   // }
1689   //  else 
1690   if(NbPointsTotal==2) {
1691     SP1.SetCoupleValue(T1,T2);
1692     SP2.SetCoupleValue(T1,T2);
1693   }
1694   else if(NbPointsTotal==1)
1695     SP1.SetCoupleValue(T1,T2);
1696   else if(NbPointsTotal==3)
1697     SP1.SetCoupleValue(T1,T2);
1698
1699   return (NbPointsTotal);
1700 }
1701 //=======================================================================
1702 //function : NextStartingPointsResearch
1703 //purpose  : 
1704 //=======================================================================
1705 Standard_Integer IntPolyh_MaillageAffinage::NextStartingPointsResearch
1706   (const Standard_Integer T1,
1707    const Standard_Integer T2,
1708    const IntPolyh_StartPoint &SPInit,
1709    IntPolyh_StartPoint &SPNext) const 
1710 {
1711   Standard_Integer NbPointsTotal=0;
1712   if( (T1<0)||(T2<0) ) NbPointsTotal=0;
1713   else {
1714     const IntPolyh_Point &P1=TPoints1[TTriangles1[T1].FirstPoint()];
1715     const IntPolyh_Point &P2=TPoints1[TTriangles1[T1].SecondPoint()];
1716     const IntPolyh_Point &P3=TPoints1[TTriangles1[T1].ThirdPoint()];
1717     const IntPolyh_Point &Q1=TPoints2[TTriangles2[T2].FirstPoint()];
1718     const IntPolyh_Point &Q2=TPoints2[TTriangles2[T2].SecondPoint()];
1719     const IntPolyh_Point &Q3=TPoints2[TTriangles2[T2].ThirdPoint()];
1720     
1721   /* The first triangle is (p1,p2,p3).  The other is (q1,q2,q3).
1722      The sides are (e1,e2,e3) and (f1,f2,f3).
1723      The normals are n1 and m1*/
1724
1725     const IntPolyh_Point  e1=P2-P1;
1726     const IntPolyh_Point  e2=P3-P2;
1727     const IntPolyh_Point  e3=P1-P3;
1728     
1729     const IntPolyh_Point  f1=Q2-Q1;
1730     const IntPolyh_Point  f2=Q3-Q2;
1731     const IntPolyh_Point  f3=Q1-Q3;
1732     
1733     IntPolyh_Point nn1,mm1;
1734     nn1.Cross(e1, e2); //normal to the first triangle
1735     mm1.Cross(f1, f2); //normal to the second triangle
1736
1737     Standard_Real nn1modulus, mm1modulus;
1738     nn1modulus=sqrt(nn1.SquareModulus());
1739     mm1modulus=sqrt(mm1.SquareModulus());
1740      
1741     //-------------------------------------------------
1742     ///calculation of intersections points between triangles
1743     //-------------------------------------------------
1744     Standard_Integer NbPoints=0;
1745     IntPolyh_StartPoint SP1,SP2;
1746             
1747     ///check T1 normal
1748     if(Abs(nn1modulus)<MyConfusionPrecision) {//10.0e-20){
1749
1750     }
1751     else {
1752       const IntPolyh_Point n1=nn1.Divide(nn1modulus);
1753       ///T2 edges with T1
1754       if(NbPointsTotal<3) {
1755         IntPolyh_StartPoint Pt1,Pt2;
1756         NbPoints=TriangleEdgeContact(1,1,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
1757         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1758       }
1759       
1760       if(NbPointsTotal<3) {
1761         IntPolyh_StartPoint Pt1,Pt2;
1762         NbPoints=TriangleEdgeContact(1,2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
1763         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1764       }
1765       
1766       if(NbPointsTotal<3) {
1767         IntPolyh_StartPoint Pt1,Pt2;
1768         NbPoints=TriangleEdgeContact(1,3,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
1769         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1770       }
1771     }
1772
1773     ///check T2 normal
1774     if(Abs(mm1modulus)<MyConfusionPrecision) {//10.0e-20){
1775
1776     }
1777     else {
1778       const IntPolyh_Point m1=mm1.Divide(mm1modulus);
1779       ///T1 edges with T2
1780       if(NbPointsTotal<3) {
1781         IntPolyh_StartPoint Pt1,Pt2;
1782         NbPoints=TriangleEdgeContact(2,1,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
1783         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1784       }
1785       
1786       if(NbPointsTotal<3) {
1787         IntPolyh_StartPoint Pt1,Pt2;
1788         NbPoints=TriangleEdgeContact(2,2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
1789         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1790       }
1791       
1792       if(NbPointsTotal<3) {
1793         IntPolyh_StartPoint Pt1,Pt2;
1794         NbPoints=TriangleEdgeContact(2,3,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
1795         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1796       }
1797     }
1798
1799     if (NbPointsTotal==1) {
1800       /*      if( (Abs(SP1.U1()-SPInit.U1())<MyConfusionPrecision)
1801               &&(Abs(SP1.V1()-SPInit.V1())<MyConfusionPrecision) ) */
1802       if(SP1.CheckSameSP(SP2))
1803         NbPointsTotal=0;
1804       else {
1805
1806         NbPointsTotal=0;
1807       }
1808     }
1809
1810     //    if ( (NbPointsTotal==2)&&( Abs(SP1.U1()-SPInit.U1())<MyConfusionPrecision)
1811     //&&( Abs(SP1.V1()-SPInit.V1())<MyConfusionPrecision) ) {
1812     if( (NbPointsTotal==2)&&(SP1.CheckSameSP(SPInit)) ) {
1813       NbPointsTotal=1;//SP1 et SPInit sont identiques
1814       SPNext=SP2;
1815     }
1816     //    if( (NbPointsTotal==2)&&( Abs(SP2.U1()-SPInit.U1())<MyConfusionPrecision)
1817     //&&( Abs(SP2.V1()-SPInit.V1())<MyConfusionPrecision) ) {
1818     if( (NbPointsTotal==2)&&(SP2.CheckSameSP(SPInit)) ) {
1819       NbPointsTotal=1;//SP2 et SPInit sont identiques
1820       SPNext=SP1;
1821     }
1822     if(NbPointsTotal>1) {
1823
1824     }
1825   }
1826   SPNext.SetCoupleValue(T1,T2);
1827   return (NbPointsTotal);
1828 }
1829 //=======================================================================
1830 //function : NextStartingPointsResearch2
1831 //purpose  : from  two triangles  and an intersection   point I
1832 //           seach the other point (if it exist).
1833 //           This function is used by StartPointChain
1834 //=======================================================================
1835 Standard_Integer IntPolyh_MaillageAffinage::NextStartingPointsResearch2
1836   (const Standard_Integer T1,
1837    const Standard_Integer T2,
1838    const IntPolyh_StartPoint &SPInit,
1839    IntPolyh_StartPoint &SPNext) const 
1840 {
1841   Standard_Integer NbPointsTotal=0;
1842   Standard_Integer EdgeInit1=SPInit.E1();
1843   Standard_Integer EdgeInit2=SPInit.E2();
1844   if( (T1<0)||(T2<0) ) NbPointsTotal=0;
1845   else {
1846     
1847     const IntPolyh_Triangle &Tri1=TTriangles1[T1];
1848     const IntPolyh_Triangle &Tri2=TTriangles2[T2];
1849
1850     const IntPolyh_Point &P1=TPoints1[Tri1.FirstPoint()];
1851     const IntPolyh_Point &P2=TPoints1[Tri1.SecondPoint()];
1852     const IntPolyh_Point &P3=TPoints1[Tri1.ThirdPoint()];
1853     const IntPolyh_Point &Q1=TPoints2[Tri2.FirstPoint()];
1854     const IntPolyh_Point &Q2=TPoints2[Tri2.SecondPoint()];
1855     const IntPolyh_Point &Q3=TPoints2[Tri2.ThirdPoint()];
1856     
1857   /* The first triangle is (p1,p2,p3).  The other is (q1,q2,q3).
1858      The edges are (e1,e2,e3) and (f1,f2,f3).
1859      The normals are n1 and m1*/
1860
1861     const IntPolyh_Point  e1=P2-P1;
1862     const IntPolyh_Point  e2=P3-P2;
1863     const IntPolyh_Point  e3=P1-P3;
1864     
1865     const IntPolyh_Point  f1=Q2-Q1;
1866     const IntPolyh_Point  f2=Q3-Q2;
1867     const IntPolyh_Point  f3=Q1-Q3;
1868     
1869     IntPolyh_Point nn1,mm1;
1870     nn1.Cross(e1, e2); //normal to the first triangle
1871     mm1.Cross(f1, f2); //normal to the second triangle
1872
1873     Standard_Real nn1modulus, mm1modulus;
1874     nn1modulus=sqrt(nn1.SquareModulus());
1875     mm1modulus=sqrt(mm1.SquareModulus());
1876     
1877     //-------------------------------------------------
1878     ///calculation of intersections points between triangles
1879     //-------------------------------------------------
1880
1881     Standard_Integer NbPoints=0;
1882     IntPolyh_StartPoint SP1,SP2;
1883
1884     ///check T1 normal
1885     if(Abs(nn1modulus)<MyConfusionPrecision) {//10.0e-20){
1886
1887     }
1888     else {
1889       const IntPolyh_Point n1=nn1.Divide(nn1modulus);
1890       ///T2 edges with T1
1891       if( (NbPointsTotal<3)&&(EdgeInit2!=Tri2.FirstEdge()) ) {
1892         IntPolyh_StartPoint Pt1,Pt2;
1893         NbPoints=TriangleEdgeContact2(1,1,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
1894         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1895       }
1896       
1897       if( (NbPointsTotal<3)&&(EdgeInit2!=Tri2.SecondEdge()) ) {
1898         IntPolyh_StartPoint Pt1,Pt2;
1899         NbPoints=TriangleEdgeContact2(1,2,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
1900         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1901       }
1902       
1903       if( (NbPointsTotal<3)&&(EdgeInit2!=Tri2.ThirdEdge()) ) {
1904         IntPolyh_StartPoint Pt1,Pt2;
1905         NbPoints=TriangleEdgeContact2(1,3,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
1906         TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1907       }
1908     }
1909     ///check T2 normal
1910     if(Abs(mm1modulus)<MyConfusionPrecision) {//10.0e-20){
1911
1912     }
1913     else {
1914       const IntPolyh_Point m1=mm1.Divide(mm1modulus);
1915       ///T1 edges with T2
1916       if( (NbPointsTotal<3)&&(EdgeInit1!=Tri1.FirstEdge()) ) {
1917         IntPolyh_StartPoint Pt1,Pt2;
1918         NbPoints=TriangleEdgeContact2(2,1,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
1919         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1920       }
1921       
1922       if( (NbPointsTotal<3)&&(EdgeInit1!=Tri1.SecondEdge()) ) {
1923         IntPolyh_StartPoint Pt1,Pt2;
1924         NbPoints=TriangleEdgeContact2(2,2,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
1925         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1926       }
1927       
1928       if( (NbPointsTotal<3)&&(EdgeInit1!=Tri1.ThirdEdge()) ) {
1929         IntPolyh_StartPoint Pt1,Pt2;
1930         NbPoints=TriangleEdgeContact2(2,3,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
1931         TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
1932       }
1933     }
1934       
1935     if (NbPointsTotal==1) {
1936       if(SP1.CheckSameSP(SPInit))
1937         NbPointsTotal=0;
1938       else {
1939         SPNext=SP1;
1940       }
1941     }
1942     else if( (NbPointsTotal==2)&&(SP1.CheckSameSP(SPInit)) ) {
1943       NbPointsTotal=1;//SP1 et SPInit sont identiques
1944       SPNext=SP2;
1945     }
1946     else if( (NbPointsTotal==2)&&(SP2.CheckSameSP(SPInit)) ) {
1947       NbPointsTotal=1;//SP2 et SPInit sont identiques
1948       SPNext=SP1;
1949     }
1950
1951     else if(NbPointsTotal>1) {
1952
1953     }
1954   }
1955   SPNext.SetCoupleValue(T1,T2);
1956   return (NbPointsTotal);
1957 }
1958 //=======================================================================
1959 //function : CalculPtsInterTriEdgeCoplanaires
1960 //purpose  : 
1961 //=======================================================================
1962 void CalculPtsInterTriEdgeCoplanaires(const Standard_Integer TriSurfID,
1963                                       const IntPolyh_Point &NormaleTri,
1964                                       const IntPolyh_Point &PE1,
1965                                       const IntPolyh_Point &PE2,
1966                                       const IntPolyh_Point &Edge,
1967                                       const IntPolyh_Point &PT1,
1968                                       const IntPolyh_Point &PT2,
1969                                       const IntPolyh_Point &Cote,
1970                                       const Standard_Integer CoteIndex,
1971                                       IntPolyh_StartPoint &SP1,
1972                                       IntPolyh_StartPoint &SP2,
1973                                       Standard_Integer &NbPoints) 
1974 {
1975   IntPolyh_Point TestParalleles;
1976   TestParalleles.Cross(Edge,Cote);
1977   if(sqrt(TestParalleles.SquareModulus())<MyConfusionPrecision) {
1978     IntPolyh_Point Per;
1979     Per.Cross(NormaleTri,Cote);
1980     Standard_Real p1p = Per.Dot(PE1);
1981     Standard_Real p2p = Per.Dot(PE2);
1982     Standard_Real p0p = Per.Dot(PT1);    
1983     if ( ( (p1p>=p0p)&&(p2p<=p0p) )||( (p1p<=p0p)&&(p2p>=p0p) ) ) {
1984       Standard_Real lambda=(p1p-p0p)/(p1p-p2p);
1985       if (lambda<-MyConfusionPrecision) {
1986
1987       }
1988       IntPolyh_Point PIE=PE1+Edge*lambda;       
1989       Standard_Real alpha=RealLast();
1990       if(Cote.X()!=0) alpha=(PIE.X()-PT1.X())/Cote.X();
1991       else if (Cote.Y()!=0) alpha=(PIE.Y()-PT1.Y())/Cote.Y();
1992       else if (Cote.Z()!=0) alpha=(PIE.Z()-PT1.Z())/Cote.Z();
1993       else {
1994
1995       }
1996       if (alpha<-MyConfusionPrecision) {
1997
1998       }
1999       else {
2000         if (NbPoints==0) {
2001           SP1.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
2002           if (TriSurfID==1) {
2003             SP1.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
2004             SP1.SetUV2(PIE.U(),PIE.V());
2005             SP1.SetEdge1(CoteIndex);
2006             NbPoints++;
2007           }
2008           else if (TriSurfID==2) {
2009             SP1.SetUV1(PIE.U(),PIE.V());
2010             SP1.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
2011             SP1.SetEdge2(CoteIndex);
2012             NbPoints++;
2013           }
2014           else {
2015
2016           }
2017         }
2018
2019         else if (NbPoints==1) {
2020           SP2.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
2021           if (TriSurfID==1) {
2022             SP2.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
2023             SP2.SetUV2(PIE.U(),PIE.V());
2024             SP2.SetEdge1(CoteIndex);
2025             NbPoints++;
2026           }
2027           else if (TriSurfID==2) {
2028             SP2.SetUV1(PIE.U(),PIE.V());
2029             SP2.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
2030             SP2.SetEdge2(CoteIndex);
2031             NbPoints++;
2032           }
2033           else {
2034
2035           }
2036         }
2037
2038         else if( (NbPoints>2)||(NbPoints<0) ) {
2039
2040         }
2041       }
2042     }
2043   }
2044   else {    //Cote et Edge paralleles, avec les rejections precedentes ils sont sur la meme droite
2045     //On projette les points sur cette droite
2046     Standard_Real pe1p= Cote.Dot(PE1);
2047     Standard_Real pe2p= Cote.Dot(PE2);
2048     Standard_Real pt1p= Cote.Dot(PT1);
2049     Standard_Real pt2p= Cote.Dot(PT2);
2050     
2051     IntPolyh_Point PEP1,PTP1,PEP2,PTP2;
2052
2053     //PEP1 et PEP2 sont les points de contact entre le triangle et l'edge dans le repere UV de l'edge
2054     //PTP1 et PTP2 sont les correspondants respectifs a PEP1 et PEP2 dans le repere UV du triangle
2055
2056
2057     if (pe1p>pe2p) {
2058       if ( (pt1p<pe1p)&&(pe1p<=pt2p) ) {
2059         PEP1=PE1;
2060         PTP1=PT1+Cote*((pe1p-pt1p)/(pt2p-pt1p));
2061         NbPoints=1;
2062         if (pt1p<=pe2p) {
2063           PEP2=PE2;
2064           PTP2=PT1+Cote*((pe2p-pt1p)/(pt2p-pt1p));
2065           NbPoints=2;
2066         }
2067         else {
2068           PEP2=PE1+Edge*((pt1p-pe1p)/(pe2p-pe1p));
2069           PTP2=PT1;
2070           NbPoints=2;
2071         }
2072       }
2073       else if( (pt2p<pe1p)&&(pe1p<=pt1p) ) {
2074         PEP1=PE1;
2075         PTP1=PT1+Cote*((pt1p-pe1p)/(pt1p-pt2p));
2076         NbPoints=1;
2077         if (pt2p<=pe2p) {
2078           PEP2=PE2;
2079           PTP2=PT1+Cote*((pe2p-pt1p)/(pt2p-pt1p));
2080           NbPoints=2;
2081         }
2082         else {
2083           PEP2=PE1+Edge*((pt2p-pe1p)/(pe2p-pe1p));
2084           PTP2=PT2;
2085           NbPoints=2;
2086         }
2087       }
2088     }
2089     
2090     if (pe1p<pe2p) {
2091       if ( (pt1p<pe2p)&&(pe2p<=pt2p) ) {
2092         PEP1=PE2;
2093         PTP1=PT1+Cote*((pe2p-pt1p)/(pt2p-pt1p));
2094         NbPoints=1;
2095         if (pt1p<=pe1p) {
2096           PEP2=PE1;
2097           PTP2=PT1+Cote*((pe1p-pt1p)/(pt2p-pt1p));
2098           NbPoints=2;
2099         }
2100         else {
2101           PEP2=PE2+Edge*((pt1p-pe1p)/(pe2p-pe1p));
2102           PTP2=PT1;
2103           NbPoints=2;
2104         }
2105       }
2106       else if( (pt2p<pe2p)&&(pe2p<=pt1p) ) {
2107         PEP1=PE2;
2108         PTP1=PT1+Cote*((pt1p-pe2p)/(pt1p-pt2p));
2109         NbPoints=1;
2110         if (pt2p<=pe1p) {
2111           PEP2=PE1;
2112           PTP2=PT1+Cote*((pe1p-pt1p)/(pt2p-pt1p));
2113           NbPoints=2;
2114         }
2115         else {
2116           PEP2=PE1+Edge*((pt2p-pe1p)/(pe2p-pe1p));
2117           PTP2=PT2;
2118           NbPoints=2;
2119         }
2120       }
2121     }
2122   
2123     if (NbPoints!=0) {
2124       if (Abs(PEP1.U()-PEP2.U())<MyConfusionPrecision
2125           &&(Abs(PEP1.V()-PEP2.V())<MyConfusionPrecision) ) NbPoints=1;
2126       
2127       SP1.SetXYZ(PEP1.X(),PEP1.Y(),PEP1.Z());
2128       if (TriSurfID==1) {
2129         SP1.SetUV1(PTP1.U(),PTP1.V());    
2130         SP1.SetUV2(PEP1.U(),PEP1.V());
2131         SP1.SetEdge1(CoteIndex);
2132       }
2133       if (TriSurfID==2) {
2134         SP1.SetUV1(PEP1.U(),PTP1.V());    
2135         SP1.SetUV2(PTP1.U(),PEP1.V());
2136         SP1.SetEdge2(CoteIndex);
2137       }
2138       
2139       if (NbPoints==2) {
2140         SP2.SetXYZ(PEP2.X(),PEP2.Y(),PEP2.Z());
2141         if (TriSurfID==1) {
2142           SP2.SetUV1(PTP2.U(),PTP2.V());          
2143           SP2.SetUV2(PEP2.U(),PEP2.V());
2144           SP2.SetEdge1(CoteIndex);
2145         }
2146         if (TriSurfID==2) {
2147           SP2.SetUV1(PEP2.U(),PTP2.V());          
2148           SP2.SetUV2(PTP2.U(),PEP2.V());
2149           SP2.SetEdge2(CoteIndex);
2150         } 
2151       }
2152     }
2153   }
2154 }
2155 //=======================================================================
2156 //function : CalculPtsInterTriEdgeCoplanaires2
2157 //purpose  : 
2158 //=======================================================================
2159 void CalculPtsInterTriEdgeCoplanaires2(const Standard_Integer TriSurfID,
2160                                        const IntPolyh_Point &NormaleTri,
2161                                        const IntPolyh_Triangle &Tri1,
2162                                        const IntPolyh_Triangle &Tri2,
2163                                        const IntPolyh_Point &PE1,
2164                                        const IntPolyh_Point &PE2,
2165                                        const IntPolyh_Point &Edge,
2166                                        const Standard_Integer EdgeIndex,
2167                                        const IntPolyh_Point &PT1,
2168                                        const IntPolyh_Point &PT2,
2169                                        const IntPolyh_Point &Cote,
2170                                        const Standard_Integer CoteIndex,
2171                                        IntPolyh_StartPoint &SP1,
2172                                        IntPolyh_StartPoint &SP2,
2173                                        Standard_Integer &NbPoints)
2174  {
2175   IntPolyh_Point TestParalleles;
2176   TestParalleles.Cross(Edge,Cote);
2177
2178   if(sqrt(TestParalleles.SquareModulus())>MyConfusionPrecision) {
2179     ///Edge and side are not parallel
2180     IntPolyh_Point Per;
2181     Per.Cross(NormaleTri,Cote);
2182     Standard_Real p1p = Per.Dot(PE1);
2183     Standard_Real p2p = Per.Dot(PE2);
2184     Standard_Real p0p = Per.Dot(PT1);
2185     ///The edge are PT1 are projected on the perpendicular of the side in the plane of the triangle
2186     if ( ( (p1p>=p0p)&&(p0p>=p2p) )||( (p1p<=p0p)&&(p0p<=p2p) ) ) {
2187       Standard_Real lambda=(p1p-p0p)/(p1p-p2p);
2188       if (lambda<-MyConfusionPrecision) {
2189
2190       }
2191       IntPolyh_Point PIE;
2192       if (Abs(lambda)<MyConfusionPrecision)//lambda=0
2193         PIE=PE1;
2194       else if (Abs(lambda)>1.0-MyConfusionPrecision)//lambda=1
2195         PIE=PE2;
2196       else
2197         PIE=PE1+Edge*lambda;
2198
2199       Standard_Real alpha=RealLast();
2200       if(Cote.X()!=0) alpha=(PIE.X()-PT1.X())/Cote.X();
2201       else if (Cote.Y()!=0) alpha=(PIE.Y()-PT1.Y())/Cote.Y();
2202       else if (Cote.Z()!=0) alpha=(PIE.Z()-PT1.Z())/Cote.Z();
2203       else {
2204
2205       }
2206       
2207       if (alpha<-MyConfusionPrecision) {
2208
2209       }
2210       else {
2211         if (NbPoints==0) {
2212           SP1.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
2213           if (TriSurfID==1) {
2214             if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
2215               SP1.SetUV1(PT1.U(),PT1.V());
2216               SP1.SetUV1(PIE.U(),PIE.V());
2217               SP1.SetEdge1(-1);
2218             }
2219             if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
2220               SP1.SetUV1(PT2.U(),PT2.V());
2221               SP1.SetUV1(PIE.U(),PIE.V());
2222               SP1.SetEdge1(-1);
2223             }
2224             else {
2225               SP1.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
2226               SP1.SetUV2(PIE.U(),PIE.V());
2227               SP1.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
2228               if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha);
2229               else SP1.SetLambda1(1.0-alpha);
2230             }
2231             NbPoints++;
2232           }
2233           else if (TriSurfID==2) {
2234             if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
2235               SP1.SetUV1(PT1.U(),PT1.V());
2236               SP1.SetUV1(PIE.U(),PIE.V());
2237               SP1.SetEdge2(-1);
2238             }
2239             if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
2240               SP1.SetUV1(PT2.U(),PT2.V());
2241               SP1.SetUV1(PIE.U(),PIE.V());
2242               SP1.SetEdge2(-1);
2243             }
2244             else {
2245               SP1.SetUV1(PIE.U(),PIE.V());
2246               SP1.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
2247               SP1.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
2248               if (Tri2.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda2(alpha);
2249               else SP1.SetLambda2(1.0-alpha);
2250             }
2251             NbPoints++;
2252           }
2253           else {
2254
2255           }
2256         }
2257
2258         else if (NbPoints==1) {
2259           SP2.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
2260           if (TriSurfID==1) {
2261             if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
2262               SP2.SetUV1(PT1.U(),PT1.V());
2263               SP2.SetUV1(PIE.U(),PIE.V());
2264               SP2.SetEdge1(-1);
2265             }
2266             if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
2267               SP2.SetUV1(PT2.U(),PT2.V());
2268               SP2.SetUV1(PIE.U(),PIE.V());
2269               SP2.SetEdge1(-1);
2270             }
2271             else {
2272               SP2.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
2273               SP2.SetUV2(PIE.U(),PIE.V());
2274               SP2.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
2275               if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha);
2276               else SP2.SetLambda1(1.0-alpha);
2277             }
2278             NbPoints++;
2279           }
2280           else if (TriSurfID==2) {
2281             if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
2282               SP2.SetUV1(PT1.U(),PT1.V());
2283               SP2.SetUV1(PIE.U(),PIE.V());
2284               SP2.SetEdge2(-1);
2285             }
2286             if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
2287               SP2.SetUV1(PT2.U(),PT2.V());
2288               SP2.SetUV1(PIE.U(),PIE.V());
2289               SP2.SetEdge2(-1);
2290             }
2291             else {
2292               SP2.SetUV1(PIE.U(),PIE.V());
2293               SP2.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
2294               SP2.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
2295               if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda2(alpha);
2296               else SP2.SetLambda2(1.0-alpha);
2297             }
2298             NbPoints++;
2299           }
2300           else {
2301
2302           }
2303         }
2304
2305         else if( (NbPoints>2)||(NbPoints<0) ) {
2306
2307           }
2308       }
2309     }
2310   }
2311   else {    
2312     //Side and Edge are parallel, with previous 
2313     //rejections they are at the same side
2314     //The points are projected on that side
2315     Standard_Real pe1p= Cote.Dot(PE1);
2316     Standard_Real pe2p= Cote.Dot(PE2);
2317     Standard_Real pt1p= Cote.Dot(PT1);
2318     Standard_Real pt2p= Cote.Dot(PT2);
2319     Standard_Real lambda1=0., lambda2=0., alpha1=0., alpha2=0.;
2320     IntPolyh_Point PEP1,PTP1,PEP2,PTP2;
2321
2322     if (pe1p>pe2p) {
2323       if ( (pt1p<pe1p)&&(pe1p<=pt2p) ) {
2324         lambda1=0.0;
2325         PEP1=PE1;
2326         alpha1=((pe1p-pt1p)/(pt2p-pt1p));
2327         PTP1=PT1+Cote*alpha1;
2328         NbPoints=1;
2329         if (pt1p<=pe2p) {
2330           lambda2=1.0;
2331           PEP2=PE2;
2332           alpha2=((pe2p-pt1p)/(pt2p-pt1p));
2333           PTP2=PT1+Cote*alpha2;
2334           NbPoints=2;
2335         }
2336         else {
2337           lambda2=((pt1p-pe1p)/(pe2p-pe1p));
2338           PEP2=PE1+Edge*lambda2;
2339           alpha2=0.0;
2340           PTP2=PT1;
2341           NbPoints=2;
2342         }
2343       }
2344       else if( (pt2p<pe1p)&&(pe1p<=pt1p) ) {
2345         lambda1=0.0;
2346         PEP1=PE1;
2347         alpha1=((pt1p-pe1p)/(pt1p-pt2p));
2348         PTP1=PT1+Cote*alpha1;
2349         NbPoints=1;
2350         if (pt2p<=pe2p) {
2351           lambda2=1.0;
2352           PEP2=PE2;
2353           alpha2=((pe2p-pt1p)/(pt2p-pt1p));
2354           PTP2=PT1+Cote*alpha2;
2355           NbPoints=2;
2356         }
2357         else {
2358           lambda2=((pt2p-pe1p)/(pe2p-pe1p));
2359           PEP2=PE1+Edge*lambda2;
2360           alpha2=1.0;
2361           PTP2=PT2;
2362           NbPoints=2;
2363         }
2364       }
2365     }
2366     
2367     if (pe1p<pe2p) {
2368       if ( (pt1p<pe2p)&&(pe2p<=pt2p) ) {
2369         lambda1=1.0;
2370         PEP1=PE2;
2371         alpha1=((pe2p-pt1p)/(pt2p-pt1p));
2372         PTP1=PT1+Cote*alpha1;
2373         NbPoints=1;
2374         if (pt1p<=pe1p) {
2375           lambda2=0.0;
2376           PEP2=PE1;
2377           alpha2=((pe1p-pt1p)/(pt2p-pt1p));
2378           PTP2=PT1+Cote*alpha2;
2379           NbPoints=2;
2380         }
2381         else {
2382           lambda2=((pt1p-pe1p)/(pe2p-pe1p));
2383           PEP2=PE2+Edge*lambda2;
2384           alpha2=0.0;
2385           PTP2=PT1;
2386           NbPoints=2;
2387         }
2388       }
2389       else if( (pt2p<pe2p)&&(pe2p<=pt1p) ) {
2390         lambda1=1.0;
2391         PEP1=PE2;
2392         alpha1=((pt1p-pe2p)/(pt1p-pt2p));
2393         PTP1=PT1+Cote*alpha1;
2394         NbPoints=1;
2395         if (pt2p<=pe1p) {
2396           lambda2=0.0;
2397           PEP2=PE1;
2398           alpha2=((pe1p-pt1p)/(pt2p-pt1p));
2399           PTP2=PT1+Cote*alpha2;
2400           NbPoints=2;
2401         }
2402         else {
2403           lambda2=((pt2p-pe1p)/(pe2p-pe1p));
2404           PEP2=PE1+Edge*lambda2;
2405           alpha2=1.0;
2406           PTP2=PT2;
2407           NbPoints=2;
2408         }
2409       }
2410     }
2411   
2412     if (NbPoints!=0) {
2413       SP1.SetXYZ(PEP1.X(),PEP1.Y(),PEP1.Z());
2414       if (TriSurfID==1) {///cote appartient a Tri1
2415         SP1.SetUV1(PTP1.U(),PTP1.V());    
2416         SP1.SetUV2(PEP1.U(),PEP1.V());
2417         SP1.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
2418
2419         if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha1);
2420         else SP1.SetLambda1(1.0-alpha1);
2421         
2422         if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP1.SetLambda2(lambda1);
2423         else SP1.SetLambda2(1.0-lambda1);
2424       }
2425       if (TriSurfID==2) {///cote appartient a Tri2
2426         SP1.SetUV1(PEP1.U(),PTP1.V());    
2427         SP1.SetUV2(PTP1.U(),PEP1.V());
2428         SP1.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
2429
2430         if(Tri2.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha1);
2431         else SP1.SetLambda1(1.0-alpha1);
2432         
2433         if(Tri1.GetEdgeOrientation(EdgeIndex)>0) SP1.SetLambda2(lambda1);
2434         else SP1.SetLambda2(1.0-lambda1);
2435       }
2436       
2437       //It is checked if PEP1!=PEP2
2438       if ( (NbPoints==2)&&(Abs(PEP1.U()-PEP2.U())<MyConfusionPrecision)
2439          &&(Abs(PEP1.V()-PEP2.V())<MyConfusionPrecision) ) NbPoints=1;
2440       if (NbPoints==2) {
2441         SP2.SetXYZ(PEP2.X(),PEP2.Y(),PEP2.Z());
2442         if (TriSurfID==1) {
2443           SP2.SetUV1(PTP2.U(),PTP2.V());          
2444           SP2.SetUV2(PEP2.U(),PEP2.V());
2445           SP2.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
2446
2447           if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha1);
2448           else SP2.SetLambda1(1.0-alpha1);
2449           
2450           if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP2.SetLambda2(lambda1);
2451           else SP2.SetLambda2(1.0-lambda1);
2452         }
2453         if (TriSurfID==2) {
2454           SP2.SetUV1(PEP2.U(),PTP2.V());          
2455           SP2.SetUV2(PTP2.U(),PEP2.V());
2456           SP2.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
2457
2458           if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha1);
2459           else SP2.SetLambda1(1.0-alpha1);
2460           
2461           if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP2.SetLambda2(lambda1);
2462           else SP2.SetLambda2(1.0-lambda1);
2463         } 
2464       }
2465     }
2466   }
2467   //Filter if the point is placed on top, the edge is set  to -1
2468   if (NbPoints>0) {
2469     if(Abs(SP1.Lambda1())<MyConfusionPrecision)
2470       SP1.SetEdge1(-1);
2471     if(Abs(SP1.Lambda1()-1)<MyConfusionPrecision)
2472       SP1.SetEdge1(-1);
2473     if(Abs(SP1.Lambda2())<MyConfusionPrecision)
2474       SP1.SetEdge2(-1);
2475     if(Abs(SP1.Lambda2()-1)<MyConfusionPrecision)
2476       SP1.SetEdge2(-1);
2477   }
2478   if (NbPoints==2) {
2479     if(Abs(SP2.Lambda1())<MyConfusionPrecision)
2480       SP2.SetEdge1(-1);
2481     if(Abs(SP2.Lambda1()-1)<MyConfusionPrecision)
2482       SP2.SetEdge1(-1);
2483     if(Abs(SP2.Lambda2())<MyConfusionPrecision)
2484       SP2.SetEdge2(-1);
2485     if(Abs(SP2.Lambda2()-1)<MyConfusionPrecision)
2486       SP2.SetEdge2(-1);
2487   }
2488 }
2489 //=======================================================================
2490 //function : TriangleEdgeContact
2491 //purpose  : 
2492 //=======================================================================
2493 Standard_Integer IntPolyh_MaillageAffinage::TriangleEdgeContact
2494   (const Standard_Integer TriSurfID,
2495    const Standard_Integer EdgeIndex,
2496    const IntPolyh_Point &PT1,
2497    const IntPolyh_Point &PT2,
2498    const IntPolyh_Point &PT3,
2499    const IntPolyh_Point &Cote12,
2500    const IntPolyh_Point &Cote23,
2501    const IntPolyh_Point &Cote31,
2502    const IntPolyh_Point &PE1,const IntPolyh_Point &PE2,
2503    const IntPolyh_Point &Edge,
2504    const IntPolyh_Point &NormaleT,
2505    IntPolyh_StartPoint &SP1,IntPolyh_StartPoint &SP2) const 
2506 {
2507
2508   Standard_Real lambda =0.;
2509   Standard_Real alpha =0.;
2510   Standard_Real beta =0.;
2511
2512     
2513   //The edge, on which the point is located, is known.
2514   if (TriSurfID==1) {
2515     SP1.SetEdge2(EdgeIndex);
2516     SP2.SetEdge2(EdgeIndex);
2517   }
2518   else if (TriSurfID==2) {
2519     SP1.SetEdge1(EdgeIndex);
2520     SP2.SetEdge1(EdgeIndex);
2521   }
2522   else {
2523
2524   }
2525
2526 /**The edge is projected on the normal of the triangle if yes 
2527   --> free intersection (point I)--> start point is found*/
2528   Standard_Integer NbPoints=0;
2529
2530   if(NormaleT.SquareModulus()==0) {
2531
2532   }
2533   else if( (Cote12.SquareModulus()==0)
2534        ||(Cote23.SquareModulus()==0)
2535        ||(Cote31.SquareModulus()==0) ) {
2536
2537   }
2538   else if(Edge.SquareModulus()==0) {
2539
2540   }
2541   else {
2542     Standard_Real pe1 = NormaleT.Dot(PE1);
2543     Standard_Real pe2 = NormaleT.Dot(PE2);
2544     Standard_Real pt1 = NormaleT.Dot(PT1);  
2545     
2546     // PE1I = lambda.Edge
2547     
2548     if( (Abs(pe1-pe2)<MyConfusionPrecision)&&(Abs(pe1-pt1)<MyConfusionPrecision) ) {
2549       //edge and triangle are coplanar (two contact points maximum)
2550
2551       //The tops of the triangle are projected on the perpendicular of the edge 
2552       
2553       IntPolyh_Point PerpEdge;
2554       PerpEdge.Cross(NormaleT,Edge);
2555       Standard_Real pp1 = PerpEdge.Dot(PT1);
2556       Standard_Real pp2 = PerpEdge.Dot(PT2);
2557       Standard_Real pp3 = PerpEdge.Dot(PT3);
2558       Standard_Real ppe1 = PerpEdge.Dot(PE1);
2559                
2560       if ( ( (pp1>ppe1)&&(pp2<=ppe1)&&(pp3<=ppe1) ) || ( (pp1<ppe1)&&(pp2>=ppe1)&&(pp3>=ppe1) ) ){
2561         //there are two sides (commun top PT1) that can cut the edge
2562         
2563         //first side
2564         CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
2565                                          PE1,PE2,Edge,PT1,PT2,
2566                                          Cote12,1,SP1,SP2,NbPoints);
2567
2568         if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2569             &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2570
2571         //second side
2572         if (NbPoints<2) {
2573           CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
2574                                            PE1,PE2,Edge,PT3,PT1,
2575                                            Cote31,3,SP1,SP2,NbPoints);
2576         }
2577       }
2578
2579       if ( (NbPoints>1)&&(Abs(SP1.U1()-SP2.U1())<MyConfusionPrecision)
2580           &&(Abs(SP1.V2()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2581       if (NbPoints>=2) return(NbPoints);
2582               
2583       else if ( ( ( (pp2>ppe1)&&(pp1<=ppe1)&&(pp3<=ppe1) ) || ( (pp2<ppe1)&&(pp1>=ppe1)&&(pp3>=ppe1) ) )
2584                && (NbPoints<2) ) {
2585         //there are two sides (common top PT2) that can cut the edge
2586         
2587         //first side
2588         CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
2589                                          PE1,PE2,Edge,PT1,PT2,
2590                                          Cote12,1,SP1,SP2,NbPoints);
2591
2592         if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2593             &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2594
2595         //second side
2596         if(NbPoints<2) CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
2597                                                         PE1,PE2,Edge,PT2,PT3,
2598                                                         Cote23,2,SP1,SP2,NbPoints);
2599       }
2600       if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2601           &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2602       if (NbPoints>=2) return(NbPoints);
2603                      //= remove
2604       else if ( (( (pp3>ppe1)&&(pp1<=ppe1)&&(pp2<=ppe1) ) || ( (pp3<ppe1)&&(pp1>=ppe1)&&(pp2>=ppe1) ))
2605                && (NbPoints<2) ) {
2606         //there are two sides (common top PT3) that can cut the edge
2607         
2608         //first side
2609         CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
2610                                          PE1,PE2,Edge,PT3,PT1,Cote31,
2611                                          3,SP1,SP2,NbPoints);
2612         
2613         if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2614             &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2615
2616         //second side
2617         if (NbPoints<2) CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
2618                                                          PE1,PE2,Edge,PT2,PT3,
2619                                                          Cote23,2,SP1,SP2,NbPoints);
2620       }
2621       if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2622           &&(Abs(SP2.V1()-SP1.V1())<MyConfusionPrecision) ) NbPoints=1;
2623       if (NbPoints>=2) return(NbPoints);
2624     }
2625     
2626
2627     //------------------------------------------------------
2628     // edge and triangle NON COPLANAR (a contact point)
2629     //------------------------------------------------------
2630     else if(   ( (pe1>=pt1)&&(pe2<=pt1) ) || ( (pe1<=pt1)&&(pe2>=pt1) )  ) {
2631       lambda=(pe1-pt1)/(pe1-pe2);
2632       IntPolyh_Point PI;
2633       if (lambda<-MyConfusionPrecision) {
2634
2635   }
2636       else if (Abs(lambda)<MyConfusionPrecision) {//lambda==0
2637         PI=PE1;
2638         if(TriSurfID==1) SP1.SetEdge2(0);
2639         else SP1.SetEdge1(0);
2640       }
2641       else if (Abs(lambda-1.0)<MyConfusionPrecision) {//lambda==1
2642         PI=PE2;
2643         if(TriSurfID==1) SP1.SetEdge2(0);
2644         else SP1.SetEdge1(0);
2645       }
2646       else {
2647         PI=PE1+Edge*lambda;
2648         if(TriSurfID==1) SP1.SetEdge2(EdgeIndex);
2649         else SP1.SetEdge1(EdgeIndex);
2650       }
2651      
2652       if(Abs(Cote23.X())>MyConfusionPrecision) {
2653         Standard_Real D=(Cote12.Y()-Cote12.X()*Cote23.Y()/Cote23.X());
2654         if(D!=0) alpha = (PI.Y()-PT1.Y()-(PI.X()-PT1.X())*Cote23.Y()/Cote23.X())/D;
2655         else { 
2656
2657         }
2658         if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
2659         else beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23.X();
2660       }
2661       
2662       else if (Abs(Cote12.X())>MyConfusionPrecision) { //On a Cote23.X()==0
2663         alpha = (PI.X()-PT1.X())/Cote12.X();
2664         
2665         if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
2666         
2667         else if (Abs(Cote23.Y())>MyConfusionPrecision) beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
2668         else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
2669         else  {
2670
2671         }
2672       }
2673       
2674       else if (Abs(Cote23.Y())>MyConfusionPrecision) {
2675         //On a Cote23.X()==0 et Cote12.X()==0 ==> equation can't be used
2676         Standard_Real D=(Cote12.Z()-Cote12.Y()*Cote23.Z()/Cote23.Y());
2677         
2678         if(D!=0) alpha = (PI.Z()-PT1.Z()-(PI.Y()-PT1.Y())*Cote23.Z()/Cote23.Y())/D;
2679         else{ 
2680
2681         }
2682   
2683         if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
2684         else beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
2685       }
2686       
2687       else if (Abs(Cote12.Y())>MyConfusionPrecision) {
2688         //On a Cote23.X()==0, Cote12.X()==0 et Cote23.Y()==0
2689         alpha = (PI.Y()-PT1.Y())/Cote12.Y();
2690         
2691         if ((Abs(alpha)<MyConfusionPrecision)||(Abs(alpha-1.0)<MyConfusionPrecision)) return(0);
2692         
2693         else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
2694         
2695         else {
2696
2697         }
2698       }    
2699       
2700       else { //two equations of three can't be used
2701
2702         alpha=RealLast();
2703         beta=RealLast();
2704       }
2705       
2706       if( (beta<-MyConfusionPrecision)||(beta>(alpha+MyConfusionPrecision)) ) return(0);
2707       else {
2708         SP1.SetXYZ(PI.X(),PI.Y(),PI.Z());
2709
2710         if (TriSurfID==1) {
2711           SP1.SetUV2(PI.U(),PI.V());
2712           SP1.SetUV1(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
2713           NbPoints++;
2714           if (beta<MyConfusionPrecision) {//beta==0 && alpha
2715             SP1.SetEdge1(1);
2716             SP1.SetLambda1(alpha);
2717           }
2718           if (Abs(beta-alpha)<MyConfusionPrecision) {//beta==alpha  
2719             SP1.SetEdge1(3);
2720             SP1.SetLambda1(1.0-alpha);
2721           }
2722           if (Abs(alpha-1)<MyConfusionPrecision)//alpha==1
2723             SP1.SetEdge1(2);
2724           if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
2725             SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
2726             SP1.SetUV1(PT1.U(),PT1.V());
2727             SP1.SetEdge1(0);
2728           }
2729           if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
2730             SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
2731             SP1.SetUV1(PT2.U(),PT2.V());
2732             SP1.SetEdge1(0);
2733           }
2734           if ( (Abs(beta-1)<MyConfusionPrecision)||(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
2735             SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
2736             SP1.SetUV1(PT3.U(),PT3.V());
2737             SP1.SetEdge1(0);
2738           }
2739         }
2740         else if(TriSurfID==2) {
2741           SP1.SetUV1(PI.U(),PI.V());
2742           SP1.SetUV2(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
2743           NbPoints++;
2744           if (beta<MyConfusionPrecision) { //beta==0
2745             SP1.SetEdge2(1);
2746           }
2747           if (Abs(beta-alpha)<MyConfusionPrecision)//beta==alpha
2748             SP1.SetEdge2(3);
2749           if (Abs(alpha-1)<MyConfusionPrecision)//alpha==1
2750             SP1.SetEdge2(2);    
2751           if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
2752             SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
2753             SP1.SetUV2(PT1.U(),PT1.V());
2754             SP1.SetEdge2(0);
2755           }
2756           if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
2757             SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
2758             SP1.SetUV2(PT2.U(),PT2.V());
2759             SP1.SetEdge2(0);
2760           }
2761           if ( (Abs(beta-1)<MyConfusionPrecision)||(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
2762             SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
2763             SP1.SetUV2(PT3.U(),PT3.V());
2764             SP1.SetEdge2(0);
2765           }
2766         }
2767         else{ 
2768
2769         }
2770       }
2771     }
2772     else return 0;
2773   }
2774   return (NbPoints);
2775 }
2776
2777 //=======================================================================
2778 //function : TriangleEdgeContact2
2779 //purpose  : 
2780 //=======================================================================
2781 Standard_Integer IntPolyh_MaillageAffinage::TriangleEdgeContact2
2782   (const Standard_Integer TriSurfID,
2783    const Standard_Integer EdgeIndex,
2784    const IntPolyh_Triangle &Tri1,
2785    const IntPolyh_Triangle &Tri2,
2786    const IntPolyh_Point &PT1,
2787    const IntPolyh_Point &PT2,
2788    const IntPolyh_Point &PT3,
2789    const IntPolyh_Point &Cote12,
2790    const IntPolyh_Point &Cote23,
2791    const IntPolyh_Point &Cote31,
2792    const IntPolyh_Point &PE1,const IntPolyh_Point &PE2,
2793    const IntPolyh_Point &Edge,
2794    const IntPolyh_Point &NormaleT,
2795    IntPolyh_StartPoint &SP1,IntPolyh_StartPoint &SP2) const 
2796 {
2797
2798   Standard_Real lambda =0., alpha =0., beta =0.;
2799
2800   //One of edges on which the point is located is known
2801   if (TriSurfID==1) {
2802     SP1.SetEdge2(Tri2.GetEdgeNumber(EdgeIndex));
2803     SP2.SetEdge2(Tri2.GetEdgeNumber(EdgeIndex));
2804   }
2805   else if (TriSurfID==2) {
2806     SP1.SetEdge1(Tri1.GetEdgeNumber(EdgeIndex));
2807     SP2.SetEdge1(Tri1.GetEdgeNumber(EdgeIndex));
2808   }
2809   else {
2810
2811   }
2812
2813   /**The edge is projected on the normal in the triangle if yes 
2814   --> a free intersection (point I)--> Start point is found */
2815   Standard_Integer NbPoints=0;
2816   if(NormaleT.SquareModulus()==0) {
2817
2818   }
2819   else if( (Cote12.SquareModulus()==0)
2820        ||(Cote23.SquareModulus()==0)
2821        ||(Cote31.SquareModulus()==0) ) {
2822
2823   }
2824   else if(Edge.SquareModulus()==0) {
2825
2826   }
2827   else {
2828     Standard_Real pe1 = NormaleT.Dot(PE1);
2829     Standard_Real pe2 = NormaleT.Dot(PE2);
2830     Standard_Real pt1 = NormaleT.Dot(PT1);  
2831     
2832     // PE1I = lambda.Edge
2833    
2834     if( (Abs(pe1-pt1)<MyConfusionPrecision)&&(Abs(pe2-pt1)<MyConfusionPrecision)) {
2835       //edge and triangle are coplanar (two contact points at maximum)
2836
2837
2838       //the tops of the triangle are projected on the perpendicular to the edge 
2839       IntPolyh_Point PerpEdge;
2840       PerpEdge.Cross(NormaleT,Edge);
2841       Standard_Real pp1 = PerpEdge.Dot(PT1);
2842       Standard_Real pp2 = PerpEdge.Dot(PT2);
2843       Standard_Real pp3 = PerpEdge.Dot(PT3);
2844       Standard_Real ppe1 = PerpEdge.Dot(PE1);
2845       
2846
2847       if ( (Abs(pp1-pp2)<MyConfusionPrecision)&&(Abs(pp1-pp3)<MyConfusionPrecision) ) {
2848
2849       }
2850       else {
2851         if ( ( (pp1>=ppe1)&&(pp2<=ppe1)&&(pp3<=ppe1) ) || ( (pp1<=ppe1)&&(pp2>=ppe1)&&(pp3>=ppe1) ) ){
2852           //there are two sides (common top PT1) that can cut the edge
2853           
2854           //first side
2855           CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2856                                             PT1,PT2,Cote12,1,SP1,SP2,NbPoints);
2857           
2858           if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2859               &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2860           
2861           //second side
2862           if (NbPoints<2) CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2863                                                             PT3,PT1,Cote31,3,SP1,SP2,NbPoints);
2864         }
2865         
2866         if ( (NbPoints>1)&&(Abs(SP1.U1()-SP2.U1())<MyConfusionPrecision)
2867             &&(Abs(SP1.V2()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2868         if (NbPoints>=2) return(NbPoints);
2869         
2870         else if ( ( ( (pp2>=ppe1)&&(pp1<=ppe1)&&(pp3<=ppe1) ) || ( (pp2<=ppe1)&&(pp1>=ppe1)&&(pp3>=ppe1) ) )
2871                  && (NbPoints<2) ) {
2872           //there are two sides (common top PT2) that can cut the edge
2873           
2874           //first side
2875           CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2876                                             PT1,PT2,Cote12,1,SP1,SP2,NbPoints);
2877           
2878           if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2879               &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2880           
2881           //second side
2882           if(NbPoints<2) CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2883                                                            PT2,PT3,Cote23,2,SP1,SP2,NbPoints);
2884         }
2885         if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2886             &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2887         if (NbPoints>=2) return(NbPoints);
2888         
2889         else if ( (( (pp3>=ppe1)&&(pp1<=ppe1)&&(pp2<=ppe1) ) || ( (pp3<=ppe1)&&(pp1>=ppe1)&&(pp2>=ppe1) ))
2890                 && (NbPoints<2) ) {
2891           //there are two sides (common top PT3) that can cut the edge
2892           
2893           //first side
2894           CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2895                                             PT3,PT1,Cote31,3,SP1,SP2,NbPoints);
2896           
2897           if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2898               &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
2899           
2900           //second side
2901           if (NbPoints<2) CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
2902                                                             PT2,PT3,Cote23,2,SP1,SP2,NbPoints);
2903         }
2904         if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
2905             &&(Abs(SP2.V1()-SP1.V1())<MyConfusionPrecision) ) NbPoints=1;
2906         if (NbPoints>=2) return(NbPoints);
2907       }
2908     }
2909
2910     //------------------------------------------------------
2911     // NON COPLANAR edge and triangle (a contact point)
2912     //------------------------------------------------------
2913     else if(   ( (pe1>=pt1)&&(pt1>=pe2) ) || ( (pe1<=pt1)&&(pt1<=pe2) )  ) { //
2914       lambda=(pe1-pt1)/(pe1-pe2);
2915       IntPolyh_Point PI;
2916       if (lambda<-MyConfusionPrecision) {
2917
2918       }
2919       else if (Abs(lambda)<MyConfusionPrecision) {//lambda==0
2920         PI=PE1;
2921         if(TriSurfID==1) SP1.SetEdge2(-1);
2922         else SP1.SetEdge1(-1);
2923       }
2924       else if (Abs(lambda-1.0)<MyConfusionPrecision) {//lambda==1
2925         PI=PE2;
2926         if(TriSurfID==1) SP1.SetEdge2(-1);
2927         else SP1.SetEdge1(-1);
2928       }
2929       else {
2930         PI=PE1+Edge*lambda;
2931         if(TriSurfID==1) {
2932           if(Tri2.GetEdgeOrientation(EdgeIndex)>0)
2933             SP1.SetLambda2(lambda);
2934           else SP1.SetLambda2(1.0-lambda);
2935         }
2936         if(TriSurfID==2) {
2937           if(Tri1.GetEdgeOrientation(EdgeIndex)>0)
2938             SP1.SetLambda1(lambda);
2939           else SP1.SetLambda1(1.0-lambda);
2940         }
2941       }
2942         
2943       Standard_Real Cote23X=Cote23.X();
2944       Standard_Real D1=0.0;
2945       Standard_Real D3,D4;
2946
2947       //Combination Eq1 Eq2
2948       if(Abs(Cote23X)>MyConfusionPrecision) {
2949         D1=Cote12.Y()-Cote12.X()*Cote23.Y()/Cote23X;
2950       }
2951       if(Abs(D1)>MyConfusionPrecision) {
2952         alpha = ( PI.Y()-PT1.Y()-(PI.X()-PT1.X())*Cote23.Y()/Cote23X )/D1;
2953
2954         ///It is checked if 1.0>=alpha>=0.0
2955         if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
2956         else beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23X;
2957       }
2958       //Combination Eq1 and Eq2 with Cote23.X()==0
2959       else if ( (Abs(Cote12.X())>MyConfusionPrecision)
2960               &&(Abs(Cote23X)<MyConfusionPrecision) ) { //There is Cote23.X()==0
2961         alpha = (PI.X()-PT1.X())/Cote12.X();
2962         
2963         if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
2964         
2965         else if (Abs(Cote23.Y())>MyConfusionPrecision) beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
2966         else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
2967         else  {
2968
2969         }
2970       }
2971       //Combination Eq1 and Eq3
2972       else if ( (Abs(Cote23.X())>MyConfusionPrecision)
2973               &&(Abs( D3= (Cote12.Z()-Cote12.X()*Cote23.Z()/Cote23.X()) ) > MyConfusionPrecision) ) {
2974         
2975         alpha = (PI.Z()-PT1.Z()-(PI.X()-PT1.X())*Cote23.Z()/Cote23.X())/D3;
2976         
2977         if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
2978         else beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23.X();
2979       }
2980       //Combination Eq2 and Eq3
2981       else if ( (Abs(Cote23.Y())>MyConfusionPrecision)
2982               &&(Abs( D4= (Cote12.Z()-Cote12.Y()*Cote23.Z()/Cote23.Y()) ) > MyConfusionPrecision) ) {
2983         
2984         alpha = (PI.Z()-PT1.Z()-(PI.Y()-PT1.Y())*Cote23.Z()/Cote23.Y())/D4;
2985         
2986         if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
2987         else beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
2988       }
2989       //Combination Eq2 and Eq3 with Cote23.Y()==0
2990       else if ( (Abs(Cote12.Y())>MyConfusionPrecision)
2991              && (Abs(Cote23.Y())<MyConfusionPrecision) ) {
2992         alpha = (PI.Y()-PT1.Y())/Cote12.Y();
2993         
2994         if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
2995         
2996         else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
2997         
2998         else {
2999           printf("\nCote PT2PT3 nul1\n");
3000           PT2.Dump(2004);
3001           PT3.Dump(3004);
3002         }
3003       }
3004       //Combination Eq1 and Eq3 with Cote23.Z()==0
3005       else if ( (Abs(Cote12.Z())>MyConfusionPrecision)
3006              && (Abs(Cote23.Z())<MyConfusionPrecision) ) {
3007         alpha = (PI.Z()-PT1.Z())/Cote12.Z();
3008         
3009         if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
3010         
3011         else if (Abs(Cote23.X())>MyConfusionPrecision) beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23.X();
3012         
3013         else {
3014
3015         }
3016       }
3017       
3018       else { //Particular case not processed ?
3019
3020         alpha=RealLast();
3021         beta=RealLast();
3022       }
3023       
3024       if( (beta<-MyConfusionPrecision)||(beta>(alpha+MyConfusionPrecision)) ) return(0);
3025       else {
3026         SP1.SetXYZ(PI.X(),PI.Y(),PI.Z());
3027
3028         if (TriSurfID==1) {
3029           SP1.SetUV2(PI.U(),PI.V());
3030           SP1.SetUV1(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
3031           NbPoints++;
3032           if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
3033             SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
3034             SP1.SetUV1(PT1.U(),PT1.V());
3035             SP1.SetEdge1(-1);
3036           }
3037           else if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
3038             SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
3039             SP1.SetUV1(PT2.U(),PT2.V());
3040             SP1.SetEdge1(-1);
3041           }
3042           else if ( (Abs(beta-1)<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
3043             SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
3044             SP1.SetUV1(PT3.U(),PT3.V());
3045             SP1.SetEdge1(-1);
3046           }
3047           else if (beta<MyConfusionPrecision) {//beta==0 
3048             SP1.SetEdge1(Tri1.GetEdgeNumber(1));
3049             if(Tri1.GetEdgeOrientation(1)>0)
3050               SP1.SetLambda1(alpha);
3051             else SP1.SetLambda1(1.0-alpha);
3052           }
3053           else if (Abs(beta-alpha)<MyConfusionPrecision) {//beta==alpha  
3054             SP1.SetEdge1(Tri1.GetEdgeNumber(3));
3055             if(Tri1.GetEdgeOrientation(3)>0)
3056               SP1.SetLambda1(1.0-alpha);
3057             else SP1.SetLambda1(alpha);
3058           }
3059           else if (Abs(alpha-1)<MyConfusionPrecision) {//alpha==1
3060             SP1.SetEdge1(Tri1.GetEdgeNumber(2));
3061             if(Tri1.GetEdgeOrientation(2)>0)
3062               SP1.SetLambda1(beta);
3063             else SP1.SetLambda1(1.0-beta);
3064           }
3065         }
3066         else if(TriSurfID==2) {
3067           SP1.SetUV1(PI.U(),PI.V());
3068           SP1.SetUV2(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
3069           NbPoints++;
3070           if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
3071             SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
3072             SP1.SetUV2(PT1.U(),PT1.V());
3073             SP1.SetEdge2(-1);
3074           }
3075           else if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
3076             SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
3077             SP1.SetUV2(PT2.U(),PT2.V());
3078             SP1.SetEdge2(-1);
3079           }
3080           else if ( (Abs(beta-1)<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
3081             SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
3082             SP1.SetUV2(PT3.U(),PT3.V());
3083             SP1.SetEdge2(-1);
3084           }
3085           else if (beta<MyConfusionPrecision) { //beta==0
3086             SP1.SetEdge2(Tri2.GetEdgeNumber(1));
3087             if(Tri2.GetEdgeOrientation(1)>0)
3088               SP1.SetLambda2(alpha);
3089             else SP1.SetLambda2(1.0-alpha);
3090           }
3091           else if (Abs(beta-alpha)<MyConfusionPrecision) {//beta==alpha
3092             SP1.SetEdge2(Tri2.GetEdgeNumber(3));
3093             if(Tri2.GetEdgeOrientation(3)>0)
3094               SP1.SetLambda2(1.0-alpha);
3095             else SP1.SetLambda2(alpha);
3096           }
3097           else if (Abs(alpha-1)<MyConfusionPrecision) {//alpha==1
3098             SP1.SetEdge2(Tri2.GetEdgeNumber(2));        
3099             if(Tri2.GetEdgeOrientation(2)>0)
3100               SP1.SetLambda2(alpha);
3101             else SP1.SetLambda2(1.0-alpha);
3102           }
3103         }
3104         else {
3105
3106         }
3107       }
3108     }
3109     else return 0;
3110   }
3111   return (NbPoints);
3112 }
3113 //=======================================================================
3114 //function : TriangleComparePSP
3115 //purpose  : The   same as   TriangleCompare, plus compute the
3116 //           StartPoints without chaining them.
3117 //=======================================================================
3118 Standard_Integer IntPolyh_MaillageAffinage::TriangleComparePSP () 
3119 {
3120   Standard_Integer CpteurTab=0;
3121   Standard_Integer CpteurTabSP=0;
3122   Standard_Real CoupleAngle=-2.0;
3123   const Standard_Integer FinTT1 = TTriangles1.NbItems();
3124   const Standard_Integer FinTT2 = TTriangles2.NbItems();
3125
3126   for(Standard_Integer i_S1=0; i_S1<FinTT1; i_S1++) {
3127     IntPolyh_Triangle &Triangle1 =  TTriangles1[i_S1];
3128     if ((Triangle1.IndiceIntersectionPossible() == 0) ||
3129         (Triangle1.GetFleche() < 0.))
3130       continue;
3131     for(Standard_Integer i_S2=0; i_S2<FinTT2; i_S2++){
3132       IntPolyh_Triangle &Triangle2 =  TTriangles2[i_S2];
3133       if ((Triangle2.IndiceIntersectionPossible() != 0) && 
3134           (Triangle2.GetFleche() >= 0.)) {
3135         IntPolyh_StartPoint SP1, SP2;
3136         //If a triangle is dead or not in BSB, comparison is not possible
3137         //
3138         Standard_Integer iDeg1, iDeg2, iDeg3, iDeg;
3139         //
3140         const IntPolyh_Point& P1=TPoints1[Triangle1.FirstPoint()];
3141         const IntPolyh_Point& P2=TPoints1[Triangle1.SecondPoint()];
3142         const IntPolyh_Point& P3=TPoints1[Triangle1.ThirdPoint()];
3143         iDeg1=(P1.Degenerated()) ? 1 : 0;
3144         iDeg2=(P2.Degenerated()) ? 1 : 0;
3145         iDeg3=(P3.Degenerated()) ? 1 : 0;
3146         iDeg=iDeg1+iDeg2+iDeg3;
3147         if (iDeg>1) {
3148           continue;
3149         }
3150         //
3151         const IntPolyh_Point& Q1=TPoints2[Triangle2.FirstPoint()];
3152         const IntPolyh_Point& Q2=TPoints2[Triangle2.SecondPoint()];
3153         const IntPolyh_Point& Q3=TPoints2[Triangle2.ThirdPoint()];
3154         iDeg1=(Q1.Degenerated()) ? 1 : 0;
3155         iDeg2=(Q2.Degenerated()) ? 1 : 0;
3156         iDeg3=(Q3.Degenerated()) ? 1 : 0;
3157         iDeg=iDeg1+iDeg2+iDeg3;
3158         if (iDeg>1) {
3159           continue;
3160         }
3161         //
3162         if (TriContact(P1, P2, P3, Q1, Q2, Q3, CoupleAngle)) {
3163           Triangle1.SetIndiceIntersection(1);//The triangle is cut by another
3164           Triangle2.SetIndiceIntersection(1);
3165           
3166           Standard_Integer NbPoints;
3167           NbPoints=StartingPointsResearch(i_S1,i_S2,SP1, SP2);
3168           
3169           if (NbPoints==0) {
3170             
3171           }
3172           
3173           if ( (NbPoints>0)&&(NbPoints<3) ) {
3174             SP1.SetCoupleValue(i_S1,i_S2);
3175             TStartPoints[CpteurTabSP]=SP1;
3176             CpteurTabSP++;
3177             
3178             
3179           }
3180           
3181           if(NbPoints==2) {       
3182             SP2.SetCoupleValue(i_S1,i_S2);
3183             TStartPoints[CpteurTabSP]=SP2;
3184             CpteurTabSP++;
3185             
3186             
3187           }
3188           
3189           if(NbPoints>2) {
3190             
3191           }
3192           CpteurTab++;
3193         }
3194       }
3195     }
3196   }
3197   return(CpteurTabSP);
3198 }
3199 //=======================================================================
3200 //function : TriangleCompare
3201 //purpose  : Analyze  each couple of  triangles from the two --
3202 //           array  of triangles,  to   see  if they are  in
3203 //           contact,  and  compute the  incidence.  Then  put
3204 //           couples  in contact  in  the  array  of  couples
3205 //=======================================================================
3206 Standard_Integer IntPolyh_MaillageAffinage::TriangleCompare ()
3207 {
3208   Standard_Integer CpteurTab=0;
3209
3210   const Standard_Integer FinTT1 = TTriangles1.NbItems();
3211   const Standard_Integer FinTT2 = TTriangles2.NbItems();
3212
3213   Standard_Integer TTClimit = 200;
3214   Standard_Integer NbTTC = FinTT1 * FinTT2 / 10;
3215   if (NbTTC < TTClimit)
3216     NbTTC = TTClimit;
3217   TTrianglesContacts.Init(NbTTC);
3218   // eap
3219   //TTrianglesContacts.Init(FinTT1 * FinTT2 / 10);
3220
3221   Standard_Real CoupleAngle=-2.0;
3222   for(Standard_Integer i_S1=0; i_S1<FinTT1; i_S1++) {
3223     IntPolyh_Triangle &Triangle1 =  TTriangles1[i_S1];
3224     if ((Triangle1.IndiceIntersectionPossible() == 0) || 
3225         (Triangle1.GetFleche() < 0.))
3226       continue;
3227     for(Standard_Integer i_S2=0; i_S2<FinTT2; i_S2++){
3228       IntPolyh_Triangle &Triangle2 =  TTriangles2[i_S2];
3229       if ((Triangle2.IndiceIntersectionPossible() != 0) && 
3230           (Triangle2.GetFleche() >= 0.)) {
3231         //If a triangle is dead or not in BSB, comparison is not possible
3232         Standard_Integer iDeg1, iDeg2, iDeg3, iDeg;
3233         //
3234         const IntPolyh_Point& P1=TPoints1[Triangle1.FirstPoint()];
3235         const IntPolyh_Point& P2=TPoints1[Triangle1.SecondPoint()];
3236         const IntPolyh_Point& P3=TPoints1[Triangle1.ThirdPoint()];
3237         iDeg1=(P1.Degenerated()) ? 1 : 0;
3238         iDeg2=(P2.Degenerated()) ? 1 : 0;
3239         iDeg3=(P3.Degenerated()) ? 1 : 0;
3240         iDeg=iDeg1+iDeg2+iDeg3;
3241         if (iDeg>1) {
3242           continue;
3243         }
3244         //
3245         const IntPolyh_Point& Q1=TPoints2[Triangle2.FirstPoint()];
3246         const IntPolyh_Point& Q2=TPoints2[Triangle2.SecondPoint()];
3247         const IntPolyh_Point& Q3=TPoints2[Triangle2.ThirdPoint()];
3248         iDeg1=(Q1.Degenerated()) ? 1 : 0;
3249         iDeg2=(Q2.Degenerated()) ? 1 : 0;
3250         iDeg3=(Q3.Degenerated()) ? 1 : 0;
3251         iDeg=iDeg1+iDeg2+iDeg3;
3252         if (iDeg>1) {
3253           continue;
3254         }
3255         //
3256         if (TriContact(P1, P2, P3, Q1, Q2, Q3, CoupleAngle)) {
3257           if (CpteurTab >= NbTTC)
3258             {
3259               TTrianglesContacts.SetNbItems(CpteurTab);
3260               return(CpteurTab);
3261             }
3262           TTrianglesContacts[CpteurTab].SetCoupleValue(i_S1, i_S2);
3263           TTrianglesContacts[CpteurTab].SetAngleValue(CoupleAngle);
3264           //test  TTrianglesContacts[CpteurTab].Dump(CpteurTab);
3265           
3266           Triangle1.SetIndiceIntersection(1);//The triangle is cut by another
3267           Triangle2.SetIndiceIntersection(1);
3268           CpteurTab++;
3269         }
3270       }
3271     }
3272   }
3273   TTrianglesContacts.SetNbItems(CpteurTab);
3274
3275   return(CpteurTab);
3276 }
3277
3278 //=======================================================================
3279 //function : StartPointsCalcul
3280 //purpose  : From the array  of couples compute  all the start
3281 //           points and display them on the screen
3282 //=======================================================================
3283 void IntPolyh_MaillageAffinage::StartPointsCalcul() const
3284 {
3285   const Standard_Integer FinTTC = TTrianglesContacts.NbItems();
3286 //   printf("StartPointsCalcul() from IntPolyh_MaillageAffinage.cxx : StartPoints:\n");
3287   for(Standard_Integer ii=0; ii<FinTTC; ii++) {
3288     IntPolyh_StartPoint SP1,SP2;
3289     Standard_Integer T1,T2;
3290     T1=TTrianglesContacts[ii].FirstValue();
3291     T2=TTrianglesContacts[ii].SecondValue();
3292     StartingPointsResearch(T1,T2,SP1,SP2);
3293     if ( (SP1.E1()!=-1)&&(SP1.E2()!=-1) ) SP1.Dump(ii);
3294     if ( (SP2.E1()!=-1)&&(SP2.E2()!=-1) ) SP2.Dump(ii);
3295   }
3296 }
3297 //=======================================================================
3298 //function : CheckCoupleAndGetAngle
3299 //purpose  : 
3300 //=======================================================================
3301 Standard_Integer CheckCoupleAndGetAngle(const Standard_Integer T1, 
3302                                         const Standard_Integer T2,
3303                                         Standard_Real& Angle, 
3304                                         IntPolyh_ArrayOfCouples &TTrianglesContacts) 
3305 {
3306   Standard_Integer Test=0;
3307   const Standard_Integer FinTTC = TTrianglesContacts.NbItems();
3308   for (Standard_Integer oioi=0; oioi<FinTTC; oioi++) {
3309     IntPolyh_Couple TestCouple = TTrianglesContacts[oioi];
3310     if ( (TestCouple.FirstValue()==T1)&&(TestCouple.AnalyseFlagValue()!=1) ) {
3311       if (TestCouple.SecondValue()==T2) {
3312         Test=oioi;
3313         TTrianglesContacts[oioi].SetAnalyseFlag(1);
3314         Angle=TTrianglesContacts[oioi].AngleValue();
3315         oioi=FinTTC;
3316       }
3317     }
3318   }
3319   return(Test);
3320 }
3321 //=======================================================================
3322 //function : CheckCoupleAndGetAngle2
3323 //purpose  : 
3324 //=======================================================================
3325 Standard_Integer CheckCoupleAndGetAngle2(const Standard_Integer T1,
3326                                          const Standard_Integer T2,
3327                                          const Standard_Integer T11, 
3328                                          const Standard_Integer T22,
3329                                          Standard_Integer &CT11,
3330                                          Standard_Integer &CT22, 
3331                                          Standard_Real & Angle,
3332                                          IntPolyh_ArrayOfCouples &TTrianglesContacts) 
3333 {
3334   ///couple T1 T2 is found in the list
3335   ///T11 and T22 are two other triangles implied  in the contact edge edge
3336   /// CT11 couple( T1,T22) and CT22 couple (T2,T11)
3337   /// these couples will be marked if there is a start point
3338   Standard_Integer Test1=0;
3339   Standard_Integer Test2=0;
3340   Standard_Integer Test3=0;
3341   const Standard_Integer FinTTC = TTrianglesContacts.NbItems();
3342   for (Standard_Integer oioi=0; oioi<FinTTC; oioi++) {
3343     IntPolyh_Couple TestCouple = TTrianglesContacts[oioi];
3344     if( (Test1==0)||(Test2==0)||(Test3==0) ) {
3345       if ( (TestCouple.FirstValue()==T1)&&(TestCouple.AnalyseFlagValue()!=1) ) {
3346         if (TestCouple.SecondValue()==T2) {
3347           Test1=1;
3348           TTrianglesContacts[oioi].SetAnalyseFlag(1);
3349           Angle=TTrianglesContacts[oioi].AngleValue();
3350         }
3351         else if (TestCouple.SecondValue()==T22) {
3352           Test2=1;
3353           CT11=oioi;
3354           Angle=TTrianglesContacts[oioi].AngleValue();
3355         }
3356       }
3357       else if( (TestCouple.FirstValue()==T11)&&(TestCouple.AnalyseFlagValue()!=1) ) {
3358         if (TestCouple.SecondValue()==T2) {
3359           Test3=1;
3360           CT22=oioi;
3361           Angle=TTrianglesContacts[oioi].AngleValue();
3362         }
3363       }
3364     }
3365     else
3366       oioi=FinTTC;
3367   }
3368   return(Test1);
3369 }
3370 //=======================================================================
3371 //function : CheckNextStartPoint
3372 //purpose  : it is checked if the point is not a top
3373 //           then it is stored in one or several valid arrays with 
3374 // the proper list number
3375 //=======================================================================
3376 Standard_Integer CheckNextStartPoint(IntPolyh_SectionLine & SectionLine,
3377                                      IntPolyh_ArrayOfTangentZones & TTangentZones,
3378                                      IntPolyh_StartPoint & SP,
3379                                      const Standard_Boolean Prepend)//=Standard_False) 
3380 {
3381   Standard_Integer Test=1;
3382   if( (SP.E1()==-1)||(SP.E2()==-1) ) {
3383     //The tops of triangle are analyzed
3384     //It is checked if they are not in the array TTangentZones
3385     Standard_Integer FinTTZ=TTangentZones.NbItems();
3386     for(Standard_Integer uiui=0; uiui<FinTTZ; uiui++) {
3387       IntPolyh_StartPoint TestSP=TTangentZones[uiui];
3388       if ( (Abs(SP.U1()-TestSP.U1())<MyConfusionPrecision)
3389           &&(Abs(SP.V1()-TestSP.V1())<MyConfusionPrecision) ) {
3390         if ( (Abs(SP.U2()-TestSP.U2())<MyConfusionPrecision)
3391             &&(Abs(SP.V2()-TestSP.V2())<MyConfusionPrecision) ) {
3392           Test=0;//SP is already in the list of  tops
3393           uiui=FinTTZ;
3394         }
3395       }
3396     }
3397     if (Test) {//the top does not belong to the list of TangentZones
3398       SP.SetChainList(-1);
3399       TTangentZones[FinTTZ]=SP;
3400       TTangentZones.IncrementNbItems();
3401       Test=0;//the examined point is a top
3402     }
3403   }
3404   else if (Test) {
3405     if (Prepend)
3406       SectionLine.Prepend(SP);
3407     else {
3408       SectionLine[SectionLine.NbStartPoints()]=SP;
3409       SectionLine.IncrementNbStartPoints();
3410     }
3411
3412   }
3413   //if the point is not a top Test=1
3414   //The chain is continued
3415   return(Test); 
3416 }
3417 //=======================================================================
3418 //function : StartPointsChain
3419 //purpose  : Loop on the array of couples. Compute StartPoints.
3420 //           Try to chain  the StartPoints into SectionLines or
3421 //           put  the  point  in  the    ArrayOfTangentZones if
3422 //           chaining it, is not possible.
3423 //=======================================================================
3424 Standard_Integer IntPolyh_MaillageAffinage::StartPointsChain
3425   (IntPolyh_ArrayOfSectionLines& TSectionLines,
3426    IntPolyh_ArrayOfTangentZones& TTangentZones) 
3427 {
3428 //Loop on the array of couples filled in the function COMPARE()
3429   const Standard_Integer FinTTC = TTrianglesContacts.NbItems();
3430
3431 //Array of tops of triangles
3432   for(Standard_Integer IndexA=0; IndexA<FinTTC; IndexA++) {
3433     //First couple of triangles.
3434     //It is checked if the couple of triangles has not been already examined.
3435     if(TTrianglesContacts[IndexA].AnalyseFlagValue()!=1) {
3436
3437       Standard_Integer SectionLineIndex=TSectionLines.NbItems();
3438       // fill last section line if still empty (eap)
3439       if (SectionLineIndex > 0
3440           &&
3441           TSectionLines[SectionLineIndex-1].NbStartPoints() == 0)
3442         SectionLineIndex -= 1;
3443       else
3444         TSectionLines.IncrementNbItems();
3445
3446       IntPolyh_SectionLine &  MySectionLine=TSectionLines[SectionLineIndex];
3447       if (MySectionLine.GetN() == 0) // eap
3448         MySectionLine.Init(10000);//Initialisation of array of StartPoint
3449
3450       Standard_Integer NbPoints=-1;
3451       Standard_Integer T1I, T2I;
3452       T1I = TTrianglesContacts[IndexA].FirstValue();
3453       T2I = TTrianglesContacts[IndexA].SecondValue();
3454       
3455       // Start points for the current couple are found
3456       IntPolyh_StartPoint SP1, SP2;
3457       NbPoints=StartingPointsResearch2(T1I,T2I,SP1, SP2);//first calculation
3458       TTrianglesContacts[IndexA].SetAnalyseFlag(1);//the couple is marked
3459
3460       if(NbPoints==1) {// particular case top/triangle or edge/edge
3461         //the start point is input in the array
3462         SP1.SetChainList(SectionLineIndex);
3463         SP1.SetAngle(TTrianglesContacts[IndexA].AngleValue());
3464         //it is checked if the point is not atop of the triangle
3465         if(CheckNextStartPoint(MySectionLine,TTangentZones,SP1)) {
3466           IntPolyh_StartPoint SPNext1;
3467           Standard_Integer TestSP1=0;
3468           
3469           //chain of a side
3470           IntPolyh_StartPoint SP11;//=SP1;
3471           if(SP1.E1()>=0) { //&&(SP1.E2()!=-1) already tested if the point is not a top
3472             Standard_Integer NextTriangle1=0;
3473             if (TEdges1[SP1.E1()].FirstTriangle()!=T1I) NextTriangle1=TEdges1[SP1.E1()].FirstTriangle();
3474             else NextTriangle1=TEdges1[SP1.E1()].SecondTriangle();
3475             
3476             Standard_Real Angle=-2.0;
3477             if (CheckCoupleAndGetAngle(NextTriangle1,T2I,Angle,TTrianglesContacts)) {
3478               //it is checked if the couple exists and is marked
3479               Standard_Integer NbPoints11=0;
3480               NbPoints11=NextStartingPointsResearch2(NextTriangle1,T2I,SP1,SP11);
3481               if (NbPoints11==1) {
3482                 SP11.SetChainList(SectionLineIndex);
3483                 SP11.SetAngle(Angle);
3484                 
3485                 if(CheckNextStartPoint(MySectionLine,TTangentZones,SP11)) {         
3486                   Standard_Integer EndChainList=1;
3487                   while (EndChainList!=0) {
3488                     TestSP1=GetNextChainStartPoint(SP11,SPNext1,MySectionLine,TTangentZones);
3489                     if(TestSP1==1) {
3490                       SPNext1.SetChainList(SectionLineIndex);
3491                       if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext1))
3492                         SP11=SPNext1;
3493                       else EndChainList=0;
3494                     }
3495                     else EndChainList=0; //There is no next point
3496                   }
3497                 }
3498                
3499               }
3500               else {
3501                 if(NbPoints11>1) {//The point is input in the array TTangentZones
3502                   TTangentZones[TTangentZones.NbItems()]=SP11;//default list number = -1
3503                   TTangentZones.IncrementNbItems();
3504                 }
3505                 else {
3506
3507                 }
3508               }
3509             }
3510           }
3511           else if (SP1.E2()<0){
3512
3513           }
3514           //chain of the other side
3515           IntPolyh_StartPoint SP12;//=SP1;
3516           if (SP1.E2()>=0) { //&&(SP1.E1()!=-1) already tested
3517             Standard_Integer NextTriangle2;
3518             if (TEdges2[SP1.E2()].FirstTriangle()!=T2I) NextTriangle2=TEdges2[SP1.E2()].FirstTriangle();
3519             else NextTriangle2=TEdges2[SP1.E2()].SecondTriangle();
3520             
3521             Standard_Real Angle=-2.0;
3522             if(CheckCoupleAndGetAngle(T1I,NextTriangle2,Angle,TTrianglesContacts)) {
3523               Standard_Integer NbPoints12=0;
3524               NbPoints12=NextStartingPointsResearch2(T1I,NextTriangle2,SP1, SP12);
3525               if (NbPoints12==1) {
3526                 
3527                 SP12.SetChainList(SectionLineIndex);
3528                 SP12.SetAngle(Angle);
3529                 Standard_Boolean Prepend = Standard_True; // eap
3530                 
3531                 if(CheckNextStartPoint(MySectionLine,TTangentZones,SP12, Prepend)) {
3532                   Standard_Integer EndChainList=1;
3533                   while (EndChainList!=0) {
3534                     TestSP1=GetNextChainStartPoint(SP12,SPNext1,
3535                                                    MySectionLine,TTangentZones,
3536                                                    Prepend); // eap
3537                     if(TestSP1==1) {
3538                       SPNext1.SetChainList(SectionLineIndex);
3539                       if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext1,Prepend))
3540                         SP12=SPNext1;
3541                       else EndChainList=0;
3542                     }
3543                     else EndChainList=0; //there is no next point
3544                   }
3545                 }
3546                 
3547                 else {
3548                   if(NbPoints12>1) {//The points are input in the array TTangentZones
3549                     TTangentZones[TTangentZones.NbItems()]=SP12;//default list number = -1
3550                     TTangentZones.IncrementNbItems();
3551                   }
3552                   else {
3553
3554                   }
3555                 }
3556               }
3557             }
3558           }
3559           else if(SP1.E1()<0){
3560
3561           }
3562         }
3563       }
3564       else if(NbPoints==2) {
3565         //the start points are input in the array 
3566         IntPolyh_StartPoint SPNext2;
3567         Standard_Integer TestSP2=0;
3568         Standard_Integer EndChainList=1;
3569
3570         SP1.SetChainList(SectionLineIndex);
3571         SP1.SetAngle(TTrianglesContacts[IndexA].AngleValue());
3572         if(CheckNextStartPoint(MySectionLine,TTangentZones,SP1)) {
3573
3574           //chain of a side
3575           while (EndChainList!=0) {
3576             TestSP2=GetNextChainStartPoint(SP1,SPNext2,MySectionLine,TTangentZones);
3577             if(TestSP2==1) {
3578               SPNext2.SetChainList(SectionLineIndex);
3579               if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext2))
3580                 SP1=SPNext2;
3581               else EndChainList=0;
3582             }
3583             else EndChainList=0; //there is no next point
3584           }
3585         }
3586         
3587         SP2.SetChainList(SectionLineIndex);
3588         SP2.SetAngle(TTrianglesContacts[IndexA].AngleValue());
3589         Standard_Boolean Prepend = Standard_True; // eap
3590
3591         if(CheckNextStartPoint(MySectionLine,TTangentZones,SP2,Prepend)) {
3592           
3593           //chain of the other side
3594           EndChainList=1;
3595           while (EndChainList!=0) {
3596             TestSP2=GetNextChainStartPoint(SP2,SPNext2,
3597                                            MySectionLine,TTangentZones,
3598                                            Prepend); // eap
3599             if(TestSP2==1) {
3600               SPNext2.SetChainList(SectionLineIndex);
3601               if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext2,Prepend))
3602                 SP2=SPNext2;
3603               else EndChainList=0;
3604             }
3605             else EndChainList=0; //there is no next point
3606           }
3607         }
3608       }
3609
3610       else if( (NbPoints>2)&&(NbPoints<7) ) {
3611         //More than two start points 
3612         //the start point is input in the table
3613         SP1.SetChainList(SectionLineIndex);
3614         CheckNextStartPoint(MySectionLine,TTangentZones,SP1);
3615       }
3616       
3617       else {
3618
3619       }
3620     }
3621   }
3622
3623   return(1);
3624 }
3625 //=======================================================================
3626 //function : GetNextChainStartPoint
3627 //purpose  : Mainly  used  by StartPointsChain(), this function
3628 //           try to compute the next StartPoint.
3629 //           GetNextChainStartPoint is used only if it is known that there are 2 contact points
3630 //=======================================================================
3631 Standard_Integer IntPolyh_MaillageAffinage::GetNextChainStartPoint
3632   (const IntPolyh_StartPoint & SP,
3633    IntPolyh_StartPoint & SPNext,
3634    IntPolyh_SectionLine & MySectionLine,
3635    IntPolyh_ArrayOfTangentZones & TTangentZones,
3636    const Standard_Boolean Prepend) 
3637 {
3638   Standard_Integer NbPoints=0;
3639   if( (SP.E1()>=0)&&(SP.E2()==-2) ) {
3640     //case if the point is on edge of T1
3641     Standard_Integer NextTriangle1;
3642     if (TEdges1[SP.E1()].FirstTriangle()!=SP.T1()) NextTriangle1=TEdges1[SP.E1()].FirstTriangle();
3643     else 
3644       NextTriangle1=TEdges1[SP.E1()].SecondTriangle();
3645     //If is checked if two triangles intersect
3646     Standard_Real Angle= -2.0;
3647     if (CheckCoupleAndGetAngle(NextTriangle1,SP.T2(),Angle,TTrianglesContacts)) {
3648       NbPoints=NextStartingPointsResearch2(NextTriangle1,SP.T2(),SP,SPNext);
3649       if( NbPoints!=1 ) {
3650         if (NbPoints>1)
3651           CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend);
3652         else {
3653
3654         NbPoints=0;
3655         }
3656       }
3657       else 
3658         SPNext.SetAngle(Angle);
3659     }
3660     else NbPoints=0;//this couple does not intersect
3661   }
3662   else if( (SP.E1()==-2)&&(SP.E2()>=0) ) {
3663     //case if the point is on edge of T2
3664     Standard_Integer NextTriangle2;
3665     if (TEdges2[SP.E2()].FirstTriangle()!=SP.T2()) NextTriangle2=TEdges2[SP.E2()].FirstTriangle();
3666     else 
3667       NextTriangle2=TEdges2[SP.E2()].SecondTriangle();
3668     Standard_Real Angle= -2.0;
3669     if (CheckCoupleAndGetAngle(SP.T1(),NextTriangle2,Angle,TTrianglesContacts)) {
3670       NbPoints=NextStartingPointsResearch2(SP.T1(),NextTriangle2,SP,SPNext);
3671       if( NbPoints!=1 ) {
3672         if (NbPoints>1)
3673           CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend);
3674         else {
3675
3676         NbPoints=0;
3677         }
3678       }
3679       else
3680         SPNext.SetAngle(Angle);
3681     }
3682     else NbPoints=0;
3683   }
3684   else if( (SP.E1()==-2)&&(SP.E2()==-2) ) { 
3685     ///no edge is touched or cut
3686
3687     NbPoints=0;
3688   }
3689   else if( (SP.E1()>=0)&&(SP.E2()>=0) ) {
3690     ///the point is located on two edges
3691       Standard_Integer NextTriangle1;
3692       Standard_Integer CpleT11=-1;
3693       Standard_Integer CpleT22=-1;
3694       if (TEdges1[SP.E1()].FirstTriangle()!=SP.T1()) NextTriangle1=TEdges1[SP.E1()].FirstTriangle();
3695       else 
3696         NextTriangle1=TEdges1[SP.E1()].SecondTriangle();
3697       Standard_Integer NextTriangle2;
3698       if (TEdges2[SP.E2()].FirstTriangle()!=SP.T2()) NextTriangle2=TEdges2[SP.E2()].FirstTriangle();
3699       else 
3700         NextTriangle2=TEdges2[SP.E2()].SecondTriangle();
3701       Standard_Real Angle= -2.0;
3702       if (CheckCoupleAndGetAngle2(NextTriangle1,NextTriangle2,
3703                                   SP.T1(),SP.T2(),CpleT11,CpleT22,
3704                                   Angle,TTrianglesContacts)) {
3705         NbPoints=NextStartingPointsResearch2(NextTriangle1,NextTriangle2,SP,SPNext);
3706         if( NbPoints!=1 ) {
3707           if (NbPoints>1) {
3708             ///The new point is checked
3709             if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend)>0) {
3710             }
3711             else {
3712
3713             }
3714           }
3715           NbPoints=0;
3716         }
3717         else {//NbPoints==1
3718           SPNext.SetAngle(Angle);     
3719           //The couples (Ti,Tj) (Ti',Tj') are marked
3720           if (CpleT11>=0) TTrianglesContacts[CpleT11].SetAnalyseFlag(1);
3721           else {
3722
3723           }
3724           if (CpleT22>=0) TTrianglesContacts[CpleT22].SetAnalyseFlag(1);
3725           else {
3726
3727           }
3728         }
3729       }
3730       else NbPoints=0;
3731     }
3732   else if( (SP.E1()==-1)||(SP.E2()==-1) ) {
3733     ///the points are tops of triangle
3734     ///the point is atored in an intermediary array
3735   }
3736   return(NbPoints);
3737 }
3738 //=======================================================================
3739 //function : GetArrayOfPoints
3740 //purpose  : 
3741 //=======================================================================
3742 const IntPolyh_ArrayOfPoints& IntPolyh_MaillageAffinage::GetArrayOfPoints
3743   (const Standard_Integer SurfID)const
3744 {
3745  if (SurfID==1)
3746  return(TPoints1);
3747  return(TPoints2);
3748 }
3749 //=======================================================================
3750 //function : GetArrayOfEdges
3751 //purpose  : 
3752 //=======================================================================
3753 const IntPolyh_ArrayOfEdges& IntPolyh_MaillageAffinage::GetArrayOfEdges
3754   (const Standard_Integer SurfID)const
3755 {
3756  if (SurfID==1)
3757  return(TEdges1);
3758  return(TEdges2);
3759 }
3760 //=======================================================================
3761 //function : GetArrayOfTriangles
3762 //purpose  : 
3763 //=======================================================================
3764 const IntPolyh_ArrayOfTriangles& 
3765   IntPolyh_MaillageAffinage::GetArrayOfTriangles
3766   (const Standard_Integer SurfID)const{
3767   if (SurfID==1)
3768   return(TTriangles1);
3769   return(TTriangles2);
3770 }
3771
3772 //=======================================================================
3773 //function : GetBox
3774 //purpose  : 
3775 //=======================================================================
3776 Bnd_Box IntPolyh_MaillageAffinage::GetBox(const Standard_Integer SurfID) const
3777 {
3778   if (SurfID==1)
3779   return(MyBox1);
3780   return(MyBox2);
3781 }
3782
3783 //=======================================================================
3784 //function : GetBoxDraw
3785 //purpose  : 
3786 //=======================================================================
3787 void IntPolyh_MaillageAffinage::GetBoxDraw(const Standard_Integer SurfID)const
3788 {
3789 Standard_Real x0,y0,z0,x1,y1,z1;
3790   if (SurfID==1) {
3791     MyBox1.Get(x0,y0,z0,x1,y1,z1);
3792   }
3793   else {
3794     MyBox2.Get(x0,y0,z0,x1,y1,z1);
3795   }
3796 }
3797 //=======================================================================
3798 //function : GetArrayOfCouples
3799 //purpose  : 
3800 //=======================================================================
3801 IntPolyh_ArrayOfCouples &IntPolyh_MaillageAffinage::GetArrayOfCouples()
3802 {
3803   return TTrianglesContacts;
3804 }
3805 //=======================================================================
3806 //function : SetEnlargeZone
3807 //purpose  : 
3808 //=======================================================================
3809 void IntPolyh_MaillageAffinage::SetEnlargeZone(Standard_Boolean& EnlargeZone)
3810 {
3811   myEnlargeZone = EnlargeZone;
3812 }
3813 //=======================================================================
3814 //function : GetEnlargeZone
3815 //purpose  : 
3816 //=======================================================================
3817 Standard_Boolean IntPolyh_MaillageAffinage::GetEnlargeZone() const
3818 {
3819   return myEnlargeZone;
3820 }
3821 //=======================================================================
3822 //function : DegeneratedIndex
3823 //purpose  : 
3824 //=======================================================================
3825 void DegeneratedIndex(const TColStd_Array1OfReal& aXpars,
3826                       const Standard_Integer aNbX,
3827                       const Handle(Adaptor3d_HSurface)& aS,
3828                       const Standard_Integer aIsoDirection,
3829                       Standard_Integer& aI1,
3830                       Standard_Integer& aI2)
3831 {
3832   Standard_Integer i;
3833   Standard_Boolean bDegX1, bDegX2;
3834   Standard_Real aDegX1, aDegX2, aTol2, aX;
3835   //
3836   aI1=0;
3837   aI2=0;
3838   aTol2=MyTolerance*MyTolerance;
3839   //
3840   if (aIsoDirection==1){ // V=const
3841     bDegX1=IsDegenerated(aS, 1, aTol2, aDegX1);
3842     bDegX2=IsDegenerated(aS, 2, aTol2, aDegX2);
3843   }
3844   else if (aIsoDirection==2){ // U=const
3845     bDegX1=IsDegenerated(aS, 3, aTol2, aDegX1);
3846     bDegX2=IsDegenerated(aS, 4, aTol2, aDegX2);
3847   }
3848   else {
3849     return;
3850   }
3851   //
3852   if (!(bDegX1 || bDegX2)) {
3853     return;
3854   }
3855   //  
3856   for(i=1; i<=aNbX; ++i) {
3857     aX=aXpars(i);
3858     if (bDegX1) {
3859       if (fabs(aX-aDegX1) < MyTolerance) {
3860         aI1=i;
3861       }
3862     }
3863     if (bDegX2) {
3864       if (fabs(aX-aDegX2) < MyTolerance) {
3865         aI2=i;
3866       }
3867     }