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