| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4TriangularFacet implementation
27 //
28 // 31.10.2004, P R Truscott, QinetiQ Ltd, UK -
29 // 01.08.2007, P R Truscott, QinetiQ Ltd, UK
30 // Significant modification t
31 // based on patches/observati
32 // Holmberg.
33 // 12.10.2012, M Gayer, CERN
34 // New implementation reducin
35 // and considerable CPU speed
36 // implementation of G4Tessel
37 // 23.02.2016, E Tcherniaev, CERN
38 // Improved test to detect de
39 // too narrow) triangles.
40 // -------------------------------------------
41
42 #include "G4TriangularFacet.hh"
43
44 #include "Randomize.hh"
45 #include "G4TessellatedGeometryAlgorithms.hh"
46
47 using namespace std;
48
49 //////////////////////////////////////////////
50 //
51 // Definition of triangular facet using absolu
52 // From this for first vector is retained to d
53 // two relative vectors (E0 and E1) define the
54 // the outward surface normal.
55 //
56 G4TriangularFacet::G4TriangularFacet (const G4
57 const G4
58 const G4
59 G4
60 {
61 fVertices = new vector<G4ThreeVector>(3);
62
63 SetVertex(0, vt0);
64 if (vertexType == ABSOLUTE)
65 {
66 SetVertex(1, vt1);
67 SetVertex(2, vt2);
68 fE1 = vt1 - vt0;
69 fE2 = vt2 - vt0;
70 }
71 else
72 {
73 SetVertex(1, vt0 + vt1);
74 SetVertex(2, vt0 + vt2);
75 fE1 = vt1;
76 fE2 = vt2;
77 }
78
79 G4ThreeVector E1xE2 = fE1.cross(fE2);
80 fArea = 0.5 * E1xE2.mag();
81 for (G4int i = 0; i < 3; ++i) fIndices[i] =
82
83 fIsDefined = true;
84 G4double delta = kCarTolerance; // Set toler
85
86 // Check length of edges
87 //
88 G4double leng1 = fE1.mag();
89 G4double leng2 = (fE2-fE1).mag();
90 G4double leng3 = fE2.mag();
91 if (leng1 <= delta || leng2 <= delta || leng
92 {
93 fIsDefined = false;
94 }
95
96 // Check min height of triangle
97 //
98 if (fIsDefined)
99 {
100 if (2.*fArea/std::max(std::max(leng1,leng2
101 {
102 fIsDefined = false;
103 }
104 }
105
106 // Define facet
107 //
108 if (!fIsDefined)
109 {
110 ostringstream message;
111 message << "Facet is too small or too narr
112 << "Triangle area = " << fArea <<
113 << "P0 = " << GetVertex(0) << G4en
114 << "P1 = " << GetVertex(1) << G4en
115 << "P2 = " << GetVertex(2) << G4en
116 << "Side1 length (P0->P1) = " << l
117 << "Side2 length (P1->P2) = " << l
118 << "Side3 length (P2->P0) = " << l
119 G4Exception("G4TriangularFacet::G4Triangul
120 "GeomSolids1001", JustWarning, message);
121 fSurfaceNormal.set(0,0,0);
122 fA = fB = fC = 0.0;
123 fDet = 0.0;
124 fCircumcentre = vt0 + 0.5*fE1 + 0.5*fE2;
125 fArea = fRadius = 0.0;
126 }
127 else
128 {
129 fSurfaceNormal = E1xE2.unit();
130 fA = fE1.mag2();
131 fB = fE1.dot(fE2);
132 fC = fE2.mag2();
133 fDet = std::fabs(fA*fC - fB*fB);
134
135 fCircumcentre =
136 vt0 + (E1xE2.cross(fE1)*fC + fE2.cross(E
137 fRadius = (fCircumcentre - vt0).mag();
138 }
139 }
140
141 //////////////////////////////////////////////
142 //
143 G4TriangularFacet::G4TriangularFacet ()
144 {
145 fVertices = new vector<G4ThreeVector>(3);
146 G4ThreeVector zero(0,0,0);
147 SetVertex(0, zero);
148 SetVertex(1, zero);
149 SetVertex(2, zero);
150 for (G4int i = 0; i < 3; ++i) fIndices[i] =
151 fIsDefined = false;
152 fSurfaceNormal.set(0,0,0);
153 fA = fB = fC = 0;
154 fE1 = zero;
155 fE2 = zero;
156 fDet = 0.0;
157 fArea = fRadius = 0.0;
158 }
159
160 //////////////////////////////////////////////
161 //
162 G4TriangularFacet::~G4TriangularFacet ()
163 {
164 SetVertices(nullptr);
165 }
166
167 //////////////////////////////////////////////
168 //
169 void G4TriangularFacet::CopyFrom (const G4Tria
170 {
171 auto p = (char *) &rhs;
172 copy(p, p + sizeof(*this), (char *)this);
173
174 if (fIndices[0] < 0 && fVertices == nullptr)
175 {
176 fVertices = new vector<G4ThreeVector>(3);
177 for (G4int i = 0; i < 3; ++i) (*fVertices)
178 }
179 }
180
181 //////////////////////////////////////////////
182 //
183 void G4TriangularFacet::MoveFrom (G4Triangular
184 {
185 fSurfaceNormal = std::move(rhs.fSurfaceNorma
186 fArea = rhs.fArea;
187 fCircumcentre = std::move(rhs.fCircumcentre)
188 fRadius = rhs.fRadius;
189 fIndices = rhs.fIndices;
190 fA = rhs.fA; fB = rhs.fB; fC = rhs.fC;
191 fDet = rhs.fDet;
192 fSqrDist = rhs.fSqrDist;
193 fE1 = std::move(rhs.fE1); fE2 = std::move(rh
194 fIsDefined = rhs.fIsDefined;
195 fVertices = rhs.fVertices;
196 rhs.fVertices = nullptr;
197 }
198
199 //////////////////////////////////////////////
200 //
201 G4TriangularFacet::G4TriangularFacet (const G4
202 : G4VFacet(rhs)
203 {
204 CopyFrom(rhs);
205 }
206
207 //////////////////////////////////////////////
208 //
209 G4TriangularFacet::G4TriangularFacet (G4Triang
210 : G4VFacet(rhs)
211 {
212 MoveFrom(rhs);
213 }
214
215 //////////////////////////////////////////////
216 //
217 G4TriangularFacet&
218 G4TriangularFacet::operator=(const G4Triangula
219 {
220 SetVertices(nullptr);
221
222 if (this != &rhs)
223 {
224 delete fVertices;
225 CopyFrom(rhs);
226 }
227
228 return *this;
229 }
230
231 //////////////////////////////////////////////
232 //
233 G4TriangularFacet&
234 G4TriangularFacet::operator=(G4TriangularFacet
235 {
236 SetVertices(nullptr);
237
238 if (this != &rhs)
239 {
240 delete fVertices;
241 MoveFrom(rhs);
242 }
243
244 return *this;
245 }
246
247 //////////////////////////////////////////////
248 //
249 // GetClone
250 //
251 // Simple member function to generate fA dupli
252 //
253 G4VFacet* G4TriangularFacet::GetClone ()
254 {
255 auto fc = new G4TriangularFacet (GetVertex(0
256 GetVertex(2
257 return fc;
258 }
259
260 //////////////////////////////////////////////
261 //
262 // GetFlippedFacet
263 //
264 // Member function to generate an identical fa
265 // pointing at 180 degrees.
266 //
267 G4TriangularFacet* G4TriangularFacet::GetFlipp
268 {
269 auto flipped = new G4TriangularFacet (GetVer
270 GetVer
271 return flipped;
272 }
273
274 //////////////////////////////////////////////
275 //
276 // Distance (G4ThreeVector)
277 //
278 // Determines the vector between p and the clo
279 // This is based on the algorithm published in
280 // Graphics," Philip J Scheider and David H Eb
281 // 2003. at the time of writing, the algorith
282 // technical note "Distance between point and
283 // at http://www.geometrictools.com/Documentat
284 //
285 // The by-product is the square-distance fSqrD
286 // in case needed by the other "Distance" memb
287 //
288 G4ThreeVector G4TriangularFacet::Distance (con
289 {
290 G4ThreeVector D = GetVertex(0) - p;
291 G4double d = fE1.dot(D);
292 G4double e = fE2.dot(D);
293 G4double f = D.mag2();
294 G4double q = fB*e - fC*d;
295 G4double t = fB*d - fA*e;
296 fSqrDist = 0.;
297
298 if (q+t <= fDet)
299 {
300 if (q < 0.0)
301 {
302 if (t < 0.0)
303 {
304 //
305 // We are in region 4.
306 //
307 if (d < 0.0)
308 {
309 t = 0.0;
310 if (-d >= fA) {q = 1.0; fSqrDist = f
311 else {q = -d/fA; fSqrDist =
312 }
313 else
314 {
315 q = 0.0;
316 if (e >= 0.0) {t = 0.0; fSqrDi
317 else if (-e >= fC) {t = 1.0; fSqrD
318 else {t = -e/fC; fSqr
319 }
320 }
321 else
322 {
323 //
324 // We are in region 3.
325 //
326 q = 0.0;
327 if (e >= 0.0) {t = 0.0; fSqrDist
328 else if (-e >= fC) {t = 1.0; fSqrDist
329 else {t = -e/fC; fSqrDis
330 }
331 }
332 else if (t < 0.0)
333 {
334 //
335 // We are in region 5.
336 //
337 t = 0.0;
338 if (d >= 0.0) {q = 0.0; fSqrDist =
339 else if (-d >= fA) {q = 1.0; fSqrDist =
340 else {q = -d/fA; fSqrDist
341 }
342 else
343 {
344 //
345 // We are in region 0.
346 //
347 G4double dist = fSurfaceNormal.dot(D);
348 fSqrDist = dist*dist;
349 return fSurfaceNormal*dist;
350 }
351 }
352 else
353 {
354 if (q < 0.0)
355 {
356 //
357 // We are in region 2.
358 //
359 G4double tmp0 = fB + d;
360 G4double tmp1 = fC + e;
361 if (tmp1 > tmp0)
362 {
363 G4double numer = tmp1 - tmp0;
364 G4double denom = fA - 2.0*fB + fC;
365 if (numer >= denom) {q = 1.0; t = 0.0;
366 else
367 {
368 q = numer/denom;
369 t = 1.0 - q;
370 fSqrDist = q*(fA*q + fB*t +2.0*d) +
371 }
372 }
373 else
374 {
375 q = 0.0;
376 if (tmp1 <= 0.0) {t = 1.0; fSqrDi
377 else if (e >= 0.0) {t = 0.0; fSqrDi
378 else {t = -e/fC; fSqr
379 }
380 }
381 else if (t < 0.0)
382 {
383 //
384 // We are in region 6.
385 //
386 G4double tmp0 = fB + e;
387 G4double tmp1 = fA + d;
388 if (tmp1 > tmp0)
389 {
390 G4double numer = tmp1 - tmp0;
391 G4double denom = fA - 2.0*fB + fC;
392 if (numer >= denom) {t = 1.0; q = 0.0;
393 else
394 {
395 t = numer/denom;
396 q = 1.0 - t;
397 fSqrDist = q*(fA*q + fB*t +2.0*d) +
398 }
399 }
400 else
401 {
402 t = 0.0;
403 if (tmp1 <= 0.0) {q = 1.0; fSqrDi
404 else if (d >= 0.0) {q = 0.0; fSqrDi
405 else {q = -d/fA; fSqr
406 }
407 }
408 else
409 //
410 // We are in region 1.
411 //
412 {
413 G4double numer = fC + e - fB - d;
414 if (numer <= 0.0)
415 {
416 q = 0.0;
417 t = 1.0;
418 fSqrDist = fC + 2.0*e + f;
419 }
420 else
421 {
422 G4double denom = fA - 2.0*fB + fC;
423 if (numer >= denom) {q = 1.0; t = 0.0;
424 else
425 {
426 q = numer/denom;
427 t = 1.0 - q;
428 fSqrDist = q*(fA*q + fB*t + 2.0*d) +
429 }
430 }
431 }
432 }
433 //
434 //
435 // Do fA check for rounding errors in the di
436 // the conventional methods for calculating
437 // near to or at the surface (as required by
438 // We'll therefore also use the magnitude-sq
439 // (Note that I've also tried to get around
440 // existing equations for
441 //
442 // fSqrDist = function(fA,fB,fC,d,q,t)
443 //
444 // and use fA more accurate addition process
445 // breakdown of cummutitivity [where (A+B)+C
446 // doesn't work.
447 // Calculation from u = D + q*fE1 + t*fE2 is
448 // more robust.
449 //
450 if (fSqrDist < 0.0) fSqrDist = 0.;
451 G4ThreeVector u = D + q*fE1 + t*fE2;
452 G4double u2 = u.mag2();
453 //
454 // The following (part of the roundoff error
455 // updates.
456 //
457 if (fSqrDist > u2) fSqrDist = u2;
458
459 return u;
460 }
461
462 //////////////////////////////////////////////
463 //
464 // Distance (G4ThreeVector, G4double)
465 //
466 // Determines the closest distance between poi
467 // use of G4ThreeVector G4TriangularFacet::Dis
468 // square of the distance in variable fSqrDist
469 // the distance is to be greater than minDist,
470 // computation and return fA very large number
471 //
472 G4double G4TriangularFacet::Distance (const G4
473 G4
474 {
475 //
476 // Start with quicky test to determine if th
477 // the triangle is any closer to p than minD
478 // about more accurate test.
479 //
480 G4double dist = kInfinity;
481 if ((p-fCircumcentre).mag()-fRadius < minDis
482 {
483 //
484 // It's possible that the triangle is clos
485 // so do more accurate assessment.
486 //
487 dist = Distance(p).mag();
488 }
489 return dist;
490 }
491
492 //////////////////////////////////////////////
493 //
494 // Distance (G4ThreeVector, G4double, G4bool)
495 //
496 // Determine the distance to point p. kInfini
497 // (1) outgoing is TRUE and the dot product of
498 // and the displacement vector from p to t
499 // (2) outgoing is FALSE and the dot product o
500 // and the displacement vector from p to t
501 // If approximate methods show the distance is
502 // forget about further computation and return
503 //
504 // This method has been heavily modified thank
505 // corrections of Rickard Holmberg.
506 //
507 G4double G4TriangularFacet::Distance (const G4
508 G4
509 const G4
510 {
511 //
512 // Start with quicky test to determine if th
513 // the triangle is any closer to p than minD
514 // about more accurate test.
515 //
516 G4double dist = kInfinity;
517 if ((p-fCircumcentre).mag()-fRadius < minDis
518 {
519 //
520 // It's possible that the triangle is clos
521 // so do more accurate assessment.
522 //
523 G4ThreeVector v = Distance(p);
524 G4double dist1 = sqrt(fSqrDist);
525 G4double dir = v.dot(fSurfaceNormal);
526 G4bool wrongSide = (dir > 0.0 && !outgoing
527 if (dist1 <= kCarTolerance)
528 {
529 //
530 // Point p is very close to triangle. C
531 // in which case return distance of 0.0
532 //
533 if (wrongSide) dist = 0.0;
534 else dist = dist1;
535 }
536 else if (!wrongSide) dist = dist1;
537 }
538 return dist;
539 }
540
541 //////////////////////////////////////////////
542 //
543 // Extent
544 //
545 // Calculates the furthest the triangle extend
546 // defined by the vector axis.
547 //
548 G4double G4TriangularFacet::Extent (const G4Th
549 {
550 G4double ss = GetVertex(0).dot(axis);
551 G4double sp = GetVertex(1).dot(axis);
552 if (sp > ss) ss = sp;
553 sp = GetVertex(2).dot(axis);
554 if (sp > ss) ss = sp;
555 return ss;
556 }
557
558 //////////////////////////////////////////////
559 //
560 // Intersect
561 //
562 // Member function to find the next intersecti
563 // direction of v. If:
564 // (1) "outgoing" is TRUE, only consider the f
565 // the face.
566 // (2) "outgoing" is FALSE, only consider the
567 // the face.
568 // Member functions returns TRUE if there is a
569 // Sets the distance (distance along w), distF
570 // and normal.
571 //
572 // Also considers intersections that happen wi
573 // distances of distFromSurface = 0.5*kCarTole
574 // This is to detect kSurface without doing fA
575 // G4TessellatedSolid::Distance(p,v) calculati
576 //
577 // This member function is thanks the valuable
578 // However, "gotos" are the Work of the Devil
579 // extreme prejudice!!
580 //
581 // IMPORTANT NOTE: These calculations are pre
582 // vector. If G4TessellatedSolid or other cla
583 // with |v| != 1 then there will be errors.
584 //
585 G4bool G4TriangularFacet::Intersect (const G4T
586 const G4T
587 G4b
588 G4d
589 G4d
590 G4T
591 {
592 //
593 // Check whether the direction of the facet
594 // and the need to be outgoing or ingoing.
595 // return false.
596 //
597 G4double w = v.dot(fSurfaceNormal);
598 if ((outgoing && w < -dirTolerance) || (!out
599 {
600 distance = kInfinity;
601 distFromSurface = kInfinity;
602 normal.set(0,0,0);
603 return false;
604 }
605 //
606 // Calculate the orthogonal distance from p
607 // triangle. Then determine if we're on the
608 // surface (at fA distance greater than kCar
609 // "outgoing".
610 //
611 const G4ThreeVector& p0 = GetVertex(0);
612 G4ThreeVector D = p0 - p;
613 distFromSurface = D.dot(fSurfaceNormal);
614 G4bool wrongSide = (outgoing && distFromSurf
615 (!outgoing && distFromSurface > 0.5*kCarT
616
617 if (wrongSide)
618 {
619 distance = kInfinity;
620 distFromSurface = kInfinity;
621 normal.set(0,0,0);
622 return false;
623 }
624
625 wrongSide = (outgoing && distFromSurface < 0
626 || (!outgoing && distFromSurface >
627 if (wrongSide)
628 {
629 //
630 // We're slightly on the wrong side of the
631 // enough using fA precise distance calcul
632 //
633 G4ThreeVector u = Distance(p);
634 if (fSqrDist <= kCarTolerance*kCarToleranc
635 {
636 //
637 // We're very close. Therefore return f
638 // to pretend we intersect.
639 //
640 // distance = -0.5*kCarTolerance
641 distance = 0.0;
642 normal = fSurfaceNormal;
643 return true;
644 }
645 else
646 {
647 //
648 // We're close to the surface containing
649 // far from the triangle, and on the wro
650 // of the surface normal and v. There i
651 //
652 distance = kInfinity;
653 distFromSurface = kInfinity;
654 normal.set(0,0,0);
655 return false;
656 }
657 }
658 if (w < dirTolerance && w > -dirTolerance)
659 {
660 //
661 // The ray is within the plane of the tria
662 // in the plane of the triangle. First try
663 // mu and nu, where mu is fE1/|fE1|. This
664 // the original algorithm due to Rickard H
665 // mathematical justification than the ori
666 // beware Rickard's was less time-consumin
667 //
668 // Note that vprime is not fA unit vector.
669 // since the values of distance along vpri
670 // with the triangle will be used to deter
671 // same time.
672 //
673 G4ThreeVector mu = fE1.unit();
674 G4ThreeVector nu = fSurfaceNormal.cross(mu
675 G4TwoVector pprime(p.dot(mu), p.dot(nu));
676 G4TwoVector vprime(v.dot(mu), v.dot(nu));
677 G4TwoVector P0prime(p0.dot(mu), p0.dot(nu)
678 G4TwoVector E0prime(fE1.mag(), 0.0);
679 G4TwoVector E1prime(fE2.dot(mu), fE2.dot(n
680 G4TwoVector loc[2];
681 if (G4TessellatedGeometryAlgorithms::Inter
682 vprime, P0
683 {
684 //
685 // There is an intersection between the
686 // Now check which part of the line inte
687 // containing the triangle in 3D.
688 //
689 G4double vprimemag = vprime.mag();
690 G4double s0 = (loc[0] - pprime).m
691 G4double s1 = (loc[1] - pprime).m
692 G4double normDist0 = fSurfaceNormal.dot(
693 G4double normDist1 = fSurfaceNormal.dot(
694
695 if ((normDist0 < 0.0 && normDist1 < 0.0)
696 || (normDist0 > 0.0 && normDist1 > 0.0)
697 || (normDist0 == 0.0 && normDist1 == 0.
698 {
699 distance = kInfinity;
700 distFromSurface = kInfinity;
701 normal.set(0,0,0);
702 return false;
703 }
704 else
705 {
706 G4double dnormDist = normDist1 - normD
707 if (fabs(dnormDist) < DBL_EPSILON)
708 {
709 distance = s0;
710 normal = fSurfaceNormal;
711 if (!outgoing) distFromSurface = -di
712 return true;
713 }
714 else
715 {
716 distance = s0 - normDist0*(s1-s0)/dn
717 normal = fSurfaceNormal;
718 if (!outgoing) distFromSurface = -di
719 return true;
720 }
721 }
722 }
723 else
724 {
725 distance = kInfinity;
726 distFromSurface = kInfinity;
727 normal.set(0,0,0);
728 return false;
729 }
730 }
731 //
732 //
733 // Use conventional algorithm to determine t
734 // intersection. This involves determining
735 // line with the plane containing the triang
736 // point is within the triangle.
737 //
738 distance = distFromSurface / w;
739 G4ThreeVector pp = p + v*distance;
740 G4ThreeVector DD = p0 - pp;
741 G4double d = fE1.dot(DD);
742 G4double e = fE2.dot(DD);
743 G4double ss = fB*e - fC*d;
744 G4double t = fB*d - fA*e;
745
746 G4double sTolerance = (fabs(fB)+ fabs(fC) +
747 G4double tTolerance = (fabs(fA)+ fabs(fB) +
748 G4double detTolerance = (fabs(fA)+ fabs(fC)
749
750 //if (ss < 0.0 || t < 0.0 || ss+t > fDet)
751 if (ss < -sTolerance || t < -tTolerance || (
752 {
753 //
754 // The intersection is outside of the tria
755 //
756 distance = distFromSurface = kInfinity;
757 normal.set(0,0,0);
758 return false;
759 }
760 else
761 {
762 //
763 // There is an intersection. Now we only
764 //
765 normal = fSurfaceNormal;
766 if (!outgoing) distFromSurface = -distFrom
767 return true;
768 }
769 }
770
771 //////////////////////////////////////////////
772 //
773 // GetPointOnFace
774 //
775 // Auxiliary method, returns a uniform random
776 //
777 G4ThreeVector G4TriangularFacet::GetPointOnFac
778 {
779 G4double u = G4UniformRand();
780 G4double v = G4UniformRand();
781 if (u+v > 1.) { u = 1. - u; v = 1. - v; }
782 return GetVertex(0) + u*fE1 + v*fE2;
783 }
784
785 //////////////////////////////////////////////
786 //
787 // GetArea
788 //
789 // Auxiliary method for returning the surface
790 //
791 G4double G4TriangularFacet::GetArea() const
792 {
793 return fArea;
794 }
795
796 //////////////////////////////////////////////
797 //
798 G4GeometryType G4TriangularFacet::GetEntityTyp
799 {
800 return "G4TriangularFacet";
801 }
802
803 //////////////////////////////////////////////
804 //
805 G4ThreeVector G4TriangularFacet::GetSurfaceNor
806 {
807 return fSurfaceNormal;
808 }
809
810 //////////////////////////////////////////////
811 //
812 void G4TriangularFacet::SetSurfaceNormal (cons
813 {
814 fSurfaceNormal = normal;
815 }
816