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