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