src/IntPolyh/IntPolyh_MaillageAffinage.cxx

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