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