1 #include <mystdlib.h>
2
3 #include <linalg.hpp>
4 #include <csg.hpp>
5
6 namespace netgen
7 {
8
Face(int pi1,int pi2,int pi3,const ARRAY<Point<3>> & points,int ainputnr)9 Polyhedra::Face::Face (int pi1, int pi2, int pi3,
10 const ARRAY<Point<3> > & points,
11 int ainputnr)
12 {
13 inputnr = ainputnr;
14
15 pnums[0] = pi1;
16 pnums[1] = pi2;
17 pnums[2] = pi3;
18
19
20 bbox.Set (points[pi1]);
21 bbox.Add (points[pi2]);
22 bbox.Add (points[pi3]);
23
24 v1 = points[pi2] - points[pi1];
25 v2 = points[pi3] - points[pi1];
26
27 n = Cross (v1, v2);
28
29 nn = n;
30 nn.Normalize();
31 // PseudoInverse (v1, v2, w1, w2);
32
33 Mat<2,3> mat;
34 Mat<3,2> inv;
35 for (int i = 0; i < 3; i++)
36 {
37 mat(0,i) = v1(i);
38 mat(1,i) = v2(i);
39 }
40 CalcInverse (mat, inv);
41 for (int i = 0; i < 3; i++)
42 {
43 w1(i) = inv(i,0);
44 w2(i) = inv(i,1);
45 }
46 }
47
48
Polyhedra()49 Polyhedra :: Polyhedra ()
50 {
51 surfaceactive.SetSize(0);
52 surfaceids.SetSize(0);
53 eps_base1 = 1e-8;
54 }
55
~Polyhedra()56 Polyhedra :: ~Polyhedra ()
57 {
58 ;
59 }
60
CreateDefault()61 Primitive * Polyhedra :: CreateDefault ()
62 {
63 return new Polyhedra();
64 }
65
BoxInSolid(const BoxSphere<3> & box) const66 INSOLID_TYPE Polyhedra :: BoxInSolid (const BoxSphere<3> & box) const
67 {
68 /*
69 for (i = 1; i <= faces.Size(); i++)
70 if (FaceBoxIntersection (i, box))
71 return DOES_INTERSECT;
72 */
73 for (int i = 0; i < faces.Size(); i++)
74 {
75 if (!faces[i].bbox.Intersect (box))
76 continue;
77 //(*testout) << "face " << i << endl;
78
79 const Point<3> & p1 = points[faces[i].pnums[0]];
80 const Point<3> & p2 = points[faces[i].pnums[1]];
81 const Point<3> & p3 = points[faces[i].pnums[2]];
82
83 if (fabs (faces[i].nn * (p1 - box.Center())) > box.Diam()/2)
84 continue;
85
86 //(*testout) << "still in loop" << endl;
87
88 double dist2 = MinDistTP2 (p1, p2, p3, box.Center());
89 //(*testout) << "p1 " << p1 << " p2 " << p2 << " p3 " << p3 << endl
90 // << " box.Center " << box.Center() << " box.Diam() " << box.Diam() << endl
91 // << " dist2 " << dist2 << " sqr(box.Diam()/2) " << sqr(box.Diam()/2) << endl;
92 if (dist2 < sqr (box.Diam()/2))
93 {
94 //(*testout) << "DOES_INTERSECT" << endl;
95 return DOES_INTERSECT;
96 }
97 };
98
99 return PointInSolid (box.Center(), 1e-3 * box.Diam());
100 }
101
102
PointInSolid(const Point<3> & p,double eps) const103 INSOLID_TYPE Polyhedra :: PointInSolid (const Point<3> & p,
104 double eps) const
105 {
106 //(*testout) << "PointInSolid p " << p << " eps " << eps << endl;
107 //(*testout) << "bbox " << poly_bbox << endl;
108
109 if((p(0) > poly_bbox.PMax()(0) + eps) || (p(0) < poly_bbox.PMin()(0) - eps) ||
110 (p(1) > poly_bbox.PMax()(1) + eps) || (p(1) < poly_bbox.PMin()(1) - eps) ||
111 (p(2) > poly_bbox.PMax()(2) + eps) || (p(2) < poly_bbox.PMin()(2) - eps))
112 {
113 //(*testout) << "returning IS_OUTSIDE" << endl;
114 return IS_OUTSIDE;
115 }
116
117 Vec<3> n, v1, v2;
118
119 // random (?) numbers:
120 n(0) = -0.424621;
121 n(1) = 0.15432;
122 n(2) = 0.89212238;
123
124 int cnt = 0;
125
126 for (int i = 0; i < faces.Size(); i++)
127 {
128 const Point<3> & p1 = points[faces[i].pnums[0]];
129
130 Vec<3> v0 = p - p1;
131
132 double lam3 = faces[i].nn * v0;
133
134 if(fabs(lam3) < eps)
135 {
136 double lam1 = (faces[i].w1 * v0);
137 double lam2 = (faces[i].w2 * v0);
138 if (lam1 >= -eps_base1 && lam2 >= -eps_base1 && lam1+lam2 <= 1+eps_base1)
139 {
140 //(*testout) << "returning DOES_INTERSECT" << endl;
141 return DOES_INTERSECT;
142 }
143 }
144 else
145 {
146 lam3 = -(faces[i].n * v0) / (faces[i].n * n);
147
148 if (lam3 < 0) continue;
149
150 Vec<3> rs = v0 + lam3 * n;
151
152 double lam1 = (faces[i].w1 * rs);
153 double lam2 = (faces[i].w2 * rs);
154 if (lam1 >= 0 && lam2 >= 0 && lam1+lam2 <= 1)
155 cnt++;
156 }
157
158 }
159
160 //(*testout) << " cnt = " << cnt%2 << endl;
161 return (cnt % 2) ? IS_INSIDE : IS_OUTSIDE;
162 }
163
164
165
166
GetTangentialSurfaceIndices(const Point<3> & p,ARRAY<int> & surfind,double eps) const167 void Polyhedra :: GetTangentialSurfaceIndices (const Point<3> & p,
168 ARRAY<int> & surfind, double eps) const
169 {
170 for (int i = 0; i < faces.Size(); i++)
171 {
172 const Point<3> & p1 = points[faces[i].pnums[0]];
173
174 Vec<3> v0 = p - p1;
175 double lam3 = -(faces[i].nn * v0); // n->nn
176
177 if (fabs (lam3) > eps) continue;
178
179 double lam1 = (faces[i].w1 * v0);
180 double lam2 = (faces[i].w2 * v0);
181
182 if (lam1 >= -eps_base1 && lam2 >= -eps_base1 && lam1+lam2 <= 1+eps_base1)
183 if (!surfind.Contains (GetSurfaceId(i)))
184 surfind.Append (GetSurfaceId(i));
185 }
186
187 }
188
189
190
VecInSolid(const Point<3> & p,const Vec<3> & v,double eps) const191 INSOLID_TYPE Polyhedra :: VecInSolid (const Point<3> & p,
192 const Vec<3> & v,
193 double eps) const
194 {
195 ARRAY<int> point_on_faces;
196 INSOLID_TYPE res(DOES_INTERSECT);
197
198 Vec<3> vn = v;
199 vn.Normalize();
200 for (int i = 0; i < faces.Size(); i++)
201 {
202 const Point<3> & p1 = points[faces[i].pnums[0]];
203
204 Vec<3> v0 = p - p1;
205 double lam3 = -(faces[i].nn * v0); // n->nn
206
207
208 if (fabs (lam3) > eps) continue;
209 //(*testout) << "lam3 <= eps" << endl;
210
211 double lam1 = (faces[i].w1 * v0);
212 double lam2 = (faces[i].w2 * v0);
213
214 if (lam1 >= -eps_base1 && lam2 >= -eps_base1 && lam1+lam2 <= 1+eps_base1)
215 {
216 point_on_faces.Append(i);
217
218 double scal = vn * faces[i].nn; // n->nn
219
220 res = DOES_INTERSECT;
221 if (scal > eps_base1) res = IS_OUTSIDE;
222 if (scal < -eps_base1) res = IS_INSIDE;
223 }
224 }
225
226 //(*testout) << "point_on_faces.Size() " << point_on_faces.Size()
227 // << " res " << res << endl;
228
229 if (point_on_faces.Size() == 0)
230 return PointInSolid (p, 0);
231 if (point_on_faces.Size() == 1)
232 return res;
233
234
235
236
237 double mindist(0);
238 bool first = true;
239
240 for(int i=0; i<point_on_faces.Size(); i++)
241 {
242 for(int j=0; j<3; j++)
243 {
244 double dist = Dist(p,points[faces[point_on_faces[i]].pnums[j]]);
245 if(dist > eps && (first || dist < mindist))
246 {
247 mindist = dist;
248 first = false;
249 }
250 }
251 }
252
253 Point<3> p2 = p + (1e-2*mindist) * vn;
254 res = PointInSolid (p2, eps);
255
256 // (*testout) << "mindist " << mindist << " res " << res << endl;
257
258 return res;
259
260
261 }
262
263
264 /*
265 INSOLID_TYPE Polyhedra :: VecInSolid2 (const Point<3> & p,
266 const Vec<3> & v1,
267 const Vec<3> & v2,
268 double eps) const
269 {
270 INSOLID_TYPE res;
271
272 res = VecInSolid(p,v1,eps);
273 if(res != DOES_INTERSECT)
274 return res;
275
276 int point_on_n_faces = 0;
277
278 Vec<3> v1n = v1;
279 v1n.Normalize();
280 Vec<3> v2n = v2;
281 v2n.Normalize();
282
283
284 for (int i = 0; i < faces.Size(); i++)
285 {
286 const Point<3> & p1 = points[faces[i].pnums[0]];
287
288 Vec<3> v0 = p - p1;
289 double lam3 = -(faces[i].n * v0);
290
291 if (fabs (lam3) > eps) continue;
292
293 double lam1 = (faces[i].w1 * v0);
294 double lam2 = (faces[i].w2 * v0);
295
296 if (lam1 >= -eps && lam2 >= -eps && lam1+lam2 <= 1+eps)
297 {
298 double scal1 = v1n * faces[i].n;
299 if (fabs (scal1) > eps) continue;
300
301
302 point_on_n_faces++;
303
304 double scal2 = v2n * faces[i].n;
305 res = DOES_INTERSECT;
306 if (scal2 > eps) res = IS_OUTSIDE;
307 if (scal2 < -eps) res = IS_INSIDE;
308 }
309 }
310
311 if (point_on_n_faces == 1)
312 return res;
313
314 cerr << "primitive::vecinsolid2 makes nonsense for polyhedra" << endl;
315
316 return Primitive :: VecInSolid2 (p, v1, v2, eps);
317 }
318 */
319
320
321
VecInSolid2(const Point<3> & p,const Vec<3> & v1,const Vec<3> & v2,double eps) const322 INSOLID_TYPE Polyhedra :: VecInSolid2 (const Point<3> & p,
323 const Vec<3> & v1,
324 const Vec<3> & v2,
325 double eps) const
326 {
327 //(*testout) << "VecInSolid2 eps " << eps << endl;
328 INSOLID_TYPE res = VecInSolid(p,v1,eps);
329 //(*testout) << "VecInSolid = " <<res <<endl;
330
331 if(res != DOES_INTERSECT)
332 return res;
333
334 int point_on_n_faces = 0;
335
336 Vec<3> v1n = v1;
337 v1n.Normalize();
338 Vec<3> v2n = v2 - (v2 * v1n) * v1n;
339 v2n.Normalize();
340
341 double cosv2, cosv2max = -1;
342
343
344 for (int i = 0; i < faces.Size(); i++)
345 {
346 const Point<3> & p1 = points[faces[i].pnums[0]];
347
348 Vec<3> v0 = p - p1;
349 if (fabs (faces[i].nn * v0) > eps) continue; // n->nn
350 if (fabs (v1n * faces[i].nn) > eps_base1) continue; // n->nn
351
352 double lam1 = (faces[i].w1 * v0);
353 double lam2 = (faces[i].w2 * v0);
354
355 if (lam1 >= -eps_base1 && lam2 >= -eps_base1 && lam1+lam2 <= 1+eps_base1)
356 {
357 // v1 is in face
358
359 Point<3> fc = Center (points[faces[i].pnums[0]],
360 points[faces[i].pnums[1]],
361 points[faces[i].pnums[2]]);
362
363 Vec<3> vpfc = fc - p;
364 cosv2 = (v2n * vpfc) / vpfc.Length();
365 if (cosv2 > cosv2max)
366 {
367 cosv2max = cosv2;
368 point_on_n_faces++;
369
370 double scal2 = v2n * faces[i].nn; // n->nn
371 res = DOES_INTERSECT;
372 if (scal2 > eps_base1) res = IS_OUTSIDE;
373 if (scal2 < -eps_base1) res = IS_INSIDE;
374
375 }
376 }
377 }
378
379 if (point_on_n_faces >= 1)
380 return res;
381
382 (*testout) << "primitive::vecinsolid2 makes nonsense for polyhedra" << endl;
383 cerr << "primitive::vecinsolid2 makes nonsense for polyhedra" << endl;
384
385 return Primitive :: VecInSolid2 (p, v1, v2, eps);
386 }
387
388
389
GetTangentialVecSurfaceIndices2(const Point<3> & p,const Vec<3> & v1,const Vec<3> & v2,ARRAY<int> & surfind,double eps) const390 void Polyhedra :: GetTangentialVecSurfaceIndices2 (const Point<3> & p, const Vec<3> & v1, const Vec<3> & v2,
391 ARRAY<int> & surfind, double eps) const
392 {
393 Vec<3> v1n = v1;
394 v1n.Normalize();
395 Vec<3> v2n = v2; // - (v2 * v1n) * v1n;
396 v2n.Normalize();
397
398
399 for (int i = 0; i < faces.Size(); i++)
400 {
401 const Point<3> & p1 = points[faces[i].pnums[0]];
402
403 Vec<3> v0 = p - p1;
404 if (fabs (v0 * faces[i].nn) > eps) continue; // n->nn
405 if (fabs (v1n * faces[i].nn) > eps_base1) continue; // n->nn
406 if (fabs (v2n * faces[i].nn) > eps_base1) continue; // n->nn
407
408 double lam01 = (faces[i].w1 * v0);
409 double lam02 = (faces[i].w2 * v0);
410 double lam03 = 1-lam01-lam02;
411 double lam11 = (faces[i].w1 * v1);
412 double lam12 = (faces[i].w2 * v1);
413 double lam13 = -lam11-lam12;
414 double lam21 = (faces[i].w1 * v2);
415 double lam22 = (faces[i].w2 * v2);
416 double lam23 = -lam21-lam22;
417
418 bool ok1 = lam01 > eps_base1 ||
419 lam01 > -eps_base1 && lam11 > eps_base1 ||
420 lam01 > -eps_base1 && lam11 > -eps_base1 && lam21 > eps_base1;
421
422 bool ok2 = lam02 > eps_base1 ||
423 lam02 > -eps_base1 && lam12 > eps_base1 ||
424 lam02 > -eps_base1 && lam12 > -eps_base1 && lam22 > eps_base1;
425
426 bool ok3 = lam03 > eps_base1 ||
427 lam03 > -eps_base1 && lam13 > eps_base1 ||
428 lam03 > -eps_base1 && lam13 > -eps_base1 && lam23 > eps_base1;
429
430 if (ok1 && ok2 && ok3)
431 {
432 if (!surfind.Contains (GetSurfaceId(faces[i].planenr)))
433 surfind.Append (GetSurfaceId(faces[i].planenr));
434 }
435 }
436 }
437
438
439
440
441
442
443
444
445
446
447
448
GetPrimitiveData(const char * & classname,ARRAY<double> & coeffs) const449 void Polyhedra :: GetPrimitiveData (const char *& classname,
450 ARRAY<double> & coeffs) const
451 {
452 classname = "Polyhedra";
453 coeffs.SetSize(0);
454 coeffs.Append (points.Size());
455 coeffs.Append (faces.Size());
456 coeffs.Append (planes.Size());
457
458 /*
459 int i, j;
460 for (i = 1; i <= planes.Size(); i++)
461 {
462 planes.Elem(i)->Print (*testout);
463 }
464 for (i = 1; i <= faces.Size(); i++)
465 {
466 (*testout) << "face " << i << " has plane " << faces.Get(i).planenr << endl;
467 for (j = 1; j <= 3; j++)
468 (*testout) << points.Get(faces.Get(i).pnums[j-1]);
469 (*testout) << endl;
470 }
471 */
472 }
473
SetPrimitiveData(ARRAY<double> & coeffs)474 void Polyhedra :: SetPrimitiveData (ARRAY<double> & coeffs)
475 {
476 ;
477 }
478
Reduce(const BoxSphere<3> & box)479 void Polyhedra :: Reduce (const BoxSphere<3> & box)
480 {
481 for (int i = 0; i < planes.Size(); i++)
482 surfaceactive[i] = 0;
483
484 for (int i = 0; i < faces.Size(); i++)
485 if (FaceBoxIntersection (i, box))
486 surfaceactive[faces[i].planenr] = 1;
487 }
488
UnReduce()489 void Polyhedra :: UnReduce ()
490 {
491 for (int i = 0; i < planes.Size(); i++)
492 surfaceactive[i] = 1;
493 }
494
AddPoint(const Point<3> & p)495 int Polyhedra :: AddPoint (const Point<3> & p)
496 {
497 if(points.Size() == 0)
498 poly_bbox.Set(p);
499 else
500 poly_bbox.Add(p);
501
502 return points.Append (p);
503 }
504
AddFace(int pi1,int pi2,int pi3,int inputnum)505 int Polyhedra :: AddFace (int pi1, int pi2, int pi3, int inputnum)
506 {
507 (*testout) << "polyhedra, add face " << pi1 << ", " << pi2 << ", " << pi3 << endl;
508
509 if(pi1 == pi2 || pi2 == pi3 || pi3 == pi1)
510 {
511 ostringstream msg;
512 msg << "Illegal point numbers for polyhedron face: " << pi1+1 << ", " << pi2+1 << ", " << pi3+1;
513 throw NgException(msg.str());
514 }
515
516 faces.Append (Face (pi1, pi2, pi3, points, inputnum));
517
518 Point<3> p1 = points[pi1];
519 Point<3> p2 = points[pi2];
520 Point<3> p3 = points[pi3];
521
522 Vec<3> v1 = p2 - p1;
523 Vec<3> v2 = p3 - p1;
524
525 Vec<3> n = Cross (v1, v2);
526 n.Normalize();
527
528 Plane pl (p1, n);
529 // int inverse;
530 // int identicto = -1;
531 // for (int i = 0; i < planes.Size(); i++)
532 // if (pl.IsIdentic (*planes[i], inverse, 1e-9*max3(v1.Length(),v2.Length(),Dist(p2,p3))))
533 // {
534 // if (!inverse)
535 // identicto = i;
536 // }
537 // // cout << "is identic = " << identicto << endl;
538 // identicto = -1; // changed April 10, JS
539
540 // if (identicto != -1)
541 // faces.Last().planenr = identicto;
542 // else
543 {
544 planes.Append (new Plane (p1, n));
545 surfaceactive.Append (1);
546 surfaceids.Append (0);
547 faces.Last().planenr = planes.Size()-1;
548 }
549
550 // (*testout) << "is plane nr " << faces.Last().planenr << endl;
551
552 return faces.Size();
553 }
554
555
556
FaceBoxIntersection(int fnr,const BoxSphere<3> & box) const557 int Polyhedra :: FaceBoxIntersection (int fnr, const BoxSphere<3> & box) const
558 {
559 /*
560 (*testout) << "check face box intersection, fnr = " << fnr << endl;
561 (*testout) << "box = " << box << endl;
562 (*testout) << "face-box = " << faces[fnr].bbox << endl;
563 */
564
565 if (!faces[fnr].bbox.Intersect (box))
566 return 0;
567
568 const Point<3> & p1 = points[faces[fnr].pnums[0]];
569 const Point<3> & p2 = points[faces[fnr].pnums[1]];
570 const Point<3> & p3 = points[faces[fnr].pnums[2]];
571
572 double dist2 = MinDistTP2 (p1, p2, p3, box.Center());
573 /*
574 (*testout) << "p1 = " << p1 << endl;
575 (*testout) << "p2 = " << p2 << endl;
576 (*testout) << "p3 = " << p3 << endl;
577
578 (*testout) << "box.Center() = " << box.Center() << endl;
579 (*testout) << "center = " << box.Center() << endl;
580 (*testout) << "dist2 = " << dist2 << endl;
581 (*testout) << "diam = " << box.Diam() << endl;
582 */
583 if (dist2 < sqr (box.Diam()/2))
584 {
585 // (*testout) << "intersect" << endl;
586 return 1;
587 }
588 return 0;
589 }
590
591
GetPolySurfs(ARRAY<ARRAY<int> * > & polysurfs)592 void Polyhedra :: GetPolySurfs(ARRAY < ARRAY<int> * > & polysurfs)
593 {
594 int maxnum = -1;
595
596 for(int i = 0; i<faces.Size(); i++)
597 {
598 if(faces[i].inputnr > maxnum)
599 maxnum = faces[i].inputnr;
600 }
601
602 polysurfs.SetSize(maxnum+1);
603 for(int i=0; i<polysurfs.Size(); i++)
604 polysurfs[i] = new ARRAY<int>;
605
606 for(int i = 0; i<faces.Size(); i++)
607 polysurfs[faces[i].inputnr]->Append(faces[i].planenr);
608 }
609
610
CalcSpecialPoints(ARRAY<Point<3>> & pts) const611 void Polyhedra::CalcSpecialPoints (ARRAY<Point<3> > & pts) const
612 {
613 for (int i = 0; i < points.Size(); i++)
614 pts.Append (points[i]);
615 }
616
617
AnalyzeSpecialPoint(const Point<3> & pt,ARRAY<Point<3>> & specpts) const618 void Polyhedra :: AnalyzeSpecialPoint (const Point<3> & pt,
619 ARRAY<Point<3> > & specpts) const
620 {
621 ;
622 }
623
SpecialPointTangentialVector(const Point<3> & p,int s1,int s2) const624 Vec<3> Polyhedra :: SpecialPointTangentialVector (const Point<3> & p, int s1, int s2) const
625 {
626 const double eps = 1e-10*poly_bbox.Diam();
627
628 for (int fi1 = 0; fi1 < faces.Size(); fi1++)
629 for (int fi2 = 0; fi2 < faces.Size(); fi2++)
630 {
631 int si1 = faces[fi1].planenr;
632 int si2 = faces[fi2].planenr;
633
634 if (surfaceids[si1] != s1 || surfaceids[si2] != s2) continue;
635
636 //(*testout) << "check pair fi1/fi2 " << fi1 << "/" << fi2 << endl;
637
638 Vec<3> n1 = GetSurface(si1) . GetNormalVector (p);
639 Vec<3> n2 = GetSurface(si2) . GetNormalVector (p);
640 Vec<3> t = Cross (n1, n2);
641
642 //(*testout) << "t = " << t << endl;
643
644
645 /*
646 int samepts = 0;
647 for (int j = 0; j < 3; j++)
648 for (int k = 0; k < 3; k++)
649 if (Dist(points[faces[fi1].pnums[j]],
650 points[faces[fi2].pnums[k]]) < eps)
651 samepts++;
652 if (samepts < 2) continue;
653 */
654
655 bool shareedge = false;
656 for(int j = 0; !shareedge && j < 3; j++)
657 {
658 Vec<3> v1 = points[faces[fi1].pnums[(j+1)%3]] - points[faces[fi1].pnums[j]];
659 double smax = v1.Length();
660 v1 *= 1./smax;
661
662 int pospos;
663 if(fabs(v1(0)) > 0.5)
664 pospos = 0;
665 else if(fabs(v1(1)) > 0.5)
666 pospos = 1;
667 else
668 pospos = 2;
669
670 double sp = (p(pospos) - points[faces[fi1].pnums[j]](pospos)) / v1(pospos);
671 if(sp < -eps || sp > smax+eps)
672 continue;
673
674 for (int k = 0; !shareedge && k < 3; k ++)
675 {
676 Vec<3> v2 = points[faces[fi2].pnums[(k+1)%3]] - points[faces[fi2].pnums[k]];
677 v2.Normalize();
678 if(v2 * v1 > 0)
679 v2 -= v1;
680 else
681 v2 += v1;
682
683 //(*testout) << "v2.Length2() " << v2.Length2() << endl;
684
685 if(v2.Length2() > 1e-18)
686 continue;
687
688 double sa,sb;
689
690 sa = (points[faces[fi2].pnums[k]](pospos) - points[faces[fi1].pnums[j]](pospos)) / v1(pospos);
691 sb = (points[faces[fi2].pnums[(k+1)%3]](pospos) - points[faces[fi1].pnums[j]](pospos)) / v1(pospos);
692
693
694 if(Dist(points[faces[fi1].pnums[j]] + sa*v1, points[faces[fi2].pnums[k]]) > eps)
695 continue;
696
697 if(sa > sb)
698 {
699 double aux = sa; sa = sb; sb = aux;
700 }
701
702 //testout->precision(20);
703 //(*testout) << "sa " << sa << " sb " << sb << " smax " << smax << " sp " << sp << " v1 " << v1 << endl;
704 //testout->precision(8);
705
706
707 shareedge = (sa < -eps && sb > eps) ||
708 (sa < smax-eps && sb > smax+eps) ||
709 (sa > -eps && sb < smax+eps);
710
711 if(!shareedge)
712 continue;
713
714 sa = max2(sa,0.);
715 sb = min2(sb,smax);
716
717 if(sp < sa+eps)
718 shareedge = (t * v1 > 0);
719 else if (sp > sb-eps)
720 shareedge = (t * v1 < 0);
721
722 }
723 }
724 if (!shareedge) continue;
725
726 t.Normalize();
727
728
729 return t;
730 }
731
732 return Vec<3> (0,0,0);
733 }
734
735
736 }
737
738
739