src/IntPolyh/IntPolyh_MaillageAffinage.cxx

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