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