Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4PolyPhiFace.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// Implementation of G4PolyPhiFace, the face that bounds a polycone or
27// polyhedra at its phi opening.
28//
29// Author: David C. Williams ([email protected])
30// --------------------------------------------------------------------
31
32#include "G4PolyPhiFace.hh"
33#include "G4ClippablePolygon.hh"
34#include "G4ReduciblePolygon.hh"
35#include "G4AffineTransform.hh"
36#include "G4SolidExtentList.hh"
38
39#include "Randomize.hh"
40#include "G4TwoVector.hh"
41
42// Constructor
43//
44// Points r,z should be supplied in clockwise order in r,z. For example:
45//
46// [1]---------[2] ^ R
47// | | |
48// | | +--> z
49// [0]---------[3]
50//
52 G4double phi,
53 G4double deltaPhi,
54 G4double phiOther )
55{
57
58 numEdges = rz->NumVertices();
59
60 rMin = rz->Amin();
61 rMax = rz->Amax();
62 zMin = rz->Bmin();
63 zMax = rz->Bmax();
64
65 //
66 // Is this the "starting" phi edge of the two?
67 //
68 G4bool start = (phiOther > phi);
69
70 //
71 // Build radial vector
72 //
73 radial = G4ThreeVector( std::cos(phi), std::sin(phi), 0.0 );
74
75 //
76 // Build normal
77 //
78 G4double zSign = start ? 1 : -1;
79 normal = G4ThreeVector( zSign*radial.y(), -zSign*radial.x(), 0 );
80
81 //
82 // Is allBehind?
83 //
84 allBehind = (zSign*(std::cos(phiOther)*radial.y() - std::sin(phiOther)*radial.x()) < 0);
85
86 //
87 // Adjacent edges
88 //
89 G4double midPhi = phi + (start ? +0.5 : -0.5)*deltaPhi;
90 G4double cosMid = std::cos(midPhi),
91 sinMid = std::sin(midPhi);
92 //
93 // Allocate corners
94 //
96 //
97 // Fill them
98 //
100
103
104 iterRZ.Begin();
105 do // Loop checking, 13.08.2015, G.Cosmo
106 {
107 corn->r = iterRZ.GetA();
108 corn->z = iterRZ.GetB();
109 corn->x = corn->r*radial.x();
110 corn->y = corn->r*radial.y();
111
112 // Add pointer on prev corner
113 //
114 if( corn == corners )
115 { corn->prev = corners+numEdges-1;}
116 else
117 { corn->prev = helper; }
118
119 // Add pointer on next corner
120 //
121 if( corn < corners+numEdges-1 )
122 { corn->next = corn+1;}
123 else
124 { corn->next = corners; }
125
126 helper = corn;
127 } while( ++corn, iterRZ.Next() );
128
129 //
130 // Allocate edges
131 //
133
134 //
135 // Fill them
136 //
137 G4double rFact = std::cos(0.5*deltaPhi);
138 G4double rFactNormalize = 1.0/std::sqrt(1.0+rFact*rFact);
139
141 * here = corners;
142 G4PolyPhiFaceEdge* edge = edges;
143 do // Loop checking, 13.08.2015, G.Cosmo
144 {
145 G4ThreeVector sideNorm;
146
147 edge->v0 = prev;
148 edge->v1 = here;
149
150 G4double dr = here->r - prev->r,
151 dz = here->z - prev->z;
152
153 edge->length = std::sqrt( dr*dr + dz*dz );
154
155 edge->tr = dr/edge->length;
156 edge->tz = dz/edge->length;
157
158 if ((here->r < DBL_MIN) && (prev->r < DBL_MIN))
159 {
160 //
161 // Sigh! Always exceptions!
162 // This edge runs at r==0, so its adjoing surface is not a
163 // PolyconeSide or PolyhedraSide, but the opposite PolyPhiFace.
164 //
165 G4double zSignOther = start ? -1 : 1;
166 sideNorm = G4ThreeVector( zSignOther*std::sin(phiOther),
167 -zSignOther*std::cos(phiOther), 0 );
168 }
169 else
170 {
171 sideNorm = G4ThreeVector( edge->tz*cosMid,
172 edge->tz*sinMid,
173 -edge->tr*rFact );
174 sideNorm *= rFactNormalize;
175 }
176 sideNorm += normal;
177
178 edge->norm3D = sideNorm.unit();
179 } while( edge++, prev=here, ++here < corners+numEdges );
180
181 //
182 // Go back and fill in corner "normals"
183 //
184 G4PolyPhiFaceEdge* prevEdge = edges+numEdges-1;
185 edge = edges;
186 do // Loop checking, 13.08.2015, G.Cosmo
187 {
188 //
189 // Calculate vertex 2D normals (on the phi surface)
190 //
191 G4double rPart = prevEdge->tr + edge->tr;
192 G4double zPart = prevEdge->tz + edge->tz;
193 G4double norm = std::sqrt( rPart*rPart + zPart*zPart );
194 G4double rNorm = +zPart/norm;
195 G4double zNorm = -rPart/norm;
196
197 edge->v0->rNorm = rNorm;
198 edge->v0->zNorm = zNorm;
199
200 //
201 // Calculate the 3D normals.
202 //
203 // Find the vector perpendicular to the z axis
204 // that defines the plane that contains the vertex normal
205 //
206 G4ThreeVector xyVector;
207
208 if (edge->v0->r < DBL_MIN)
209 {
210 //
211 // This is a vertex at r==0, which is a special
212 // case. The normal we will construct lays in the
213 // plane at the center of the phi opening.
214 //
215 // We also know that rNorm < 0
216 //
217 G4double zSignOther = start ? -1 : 1;
218 G4ThreeVector normalOther( zSignOther*std::sin(phiOther),
219 -zSignOther*std::cos(phiOther), 0 );
220
221 xyVector = - normal - normalOther;
222 }
223 else
224 {
225 //
226 // This is a vertex at r > 0. The plane
227 // is the average of the normal and the
228 // normal of the adjacent phi face
229 //
230 xyVector = G4ThreeVector( cosMid, sinMid, 0 );
231 if (rNorm < 0)
232 xyVector -= normal;
233 else
234 xyVector += normal;
235 }
236
237 //
238 // Combine it with the r/z direction from the face
239 //
240 edge->v0->norm3D = rNorm*xyVector.unit() + G4ThreeVector( 0, 0, zNorm );
241 } while( prevEdge=edge, ++edge < edges+numEdges );
242
243 //
244 // Build point on surface
245 //
246 G4double rAve = 0.5*(rMax-rMin),
247 zAve = 0.5*(zMax-zMin);
248 surface = G4ThreeVector( rAve*radial.x(), rAve*radial.y(), zAve );
249}
250
251// Diagnose
252//
253// Throw an exception if something is found inconsistent with
254// the solid.
255//
256// For debugging purposes only
257//
259{
261 do // Loop checking, 13.08.2015, G.Cosmo
262 {
263 G4ThreeVector test(corner->x, corner->y, corner->z);
264 test -= 1E-6*corner->norm3D;
265
266 if (owner->Inside(test) != kInside)
267 G4Exception( "G4PolyPhiFace::Diagnose()", "GeomSolids0002",
268 FatalException, "Bad vertex normal found." );
269 } while( ++corner < corners+numEdges );
270}
271
272// Fake default constructor - sets only member data and allocates memory
273// for usage restricted to object persistency.
274//
276 : numEdges(0), rMin(0.), rMax(0.), zMin(0.), zMax(0.), kCarTolerance(0.)
277{
278}
279
280
281//
282// Destructor
283//
285{
286 delete [] edges;
287 delete [] corners;
288}
289
290
291//
292// Copy constructor
293//
295 : G4VCSGface()
296{
297 CopyStuff( source );
298}
299
300
301//
302// Assignment operator
303//
305{
306 if (this == &source) { return *this; }
307
308 delete [] edges;
309 delete [] corners;
310
311 CopyStuff( source );
312
313 return *this;
314}
315
316
317//
318// CopyStuff (protected)
319//
321{
322 //
323 // The simple stuff
324 //
325 numEdges = source.numEdges;
326 normal = source.normal;
327 radial = source.radial;
328 surface = source.surface;
329 rMin = source.rMin;
330 rMax = source.rMax;
331 zMin = source.zMin;
332 zMax = source.zMax;
333 allBehind = source.allBehind;
334 triangles = nullptr;
335
337 fSurfaceArea = source.fSurfaceArea;
338
339 //
340 // Corner dynamic array
341 //
344 *sourceCorn = source.corners;
345 do // Loop checking, 13.08.2015, G.Cosmo
346 {
347 *corn = *sourceCorn;
348 } while( ++sourceCorn, ++corn < corners+numEdges );
349
350 //
351 // Edge dynamic array
352 //
354
356 * here = corners;
357 G4PolyPhiFaceEdge* edge = edges,
358 * sourceEdge = source.edges;
359 do // Loop checking, 13.08.2015, G.Cosmo
360 {
361 *edge = *sourceEdge;
362 edge->v0 = prev;
363 edge->v1 = here;
364 } while( ++sourceEdge, ++edge, prev=here, ++here < corners+numEdges );
365}
366
367// Intersect
368//
370 const G4ThreeVector& v,
371 G4bool outgoing,
372 G4double surfTolerance,
373 G4double& distance,
374 G4double& distFromSurface,
375 G4ThreeVector& aNormal,
376 G4bool& isAllBehind )
377{
378 G4double normSign = outgoing ? +1 : -1;
379
380 //
381 // These don't change
382 //
383 isAllBehind = allBehind;
384 aNormal = normal;
385
386 //
387 // Correct normal? Here we have straight sides, and can safely ignore
388 // intersections where the dot product with the normal is zero.
389 //
390 G4double dotProd = normSign*normal.dot(v);
391
392 if (dotProd <= 0) return false;
393
394 //
395 // Calculate distance to surface. If the side is too far
396 // behind the point, we must reject it.
397 //
398 G4ThreeVector ps = p - surface;
399 distFromSurface = -normSign*ps.dot(normal);
400
401 if (distFromSurface < -surfTolerance) return false;
402
403 //
404 // Calculate precise distance to intersection with the side
405 // (along the trajectory, not normal to the surface)
406 //
407 distance = distFromSurface/dotProd;
408
409 //
410 // Calculate intersection point in r,z
411 //
412 G4ThreeVector ip = p + distance*v;
413
414 G4double r = radial.dot(ip);
415
416 //
417 // And is it inside the r/z extent?
418 //
419 return InsideEdgesExact( r, ip.z(), normSign, p, v );
420}
421
422// Distance
423//
425{
426 G4double normSign = outgoing ? +1 : -1;
427 //
428 // Correct normal?
429 //
430 G4ThreeVector ps = p - surface;
431 G4double distPhi = -normSign*normal.dot(ps);
432
433 if (distPhi < -0.5*kCarTolerance)
434 return kInfinity;
435 else if (distPhi < 0)
436 distPhi = 0.0;
437
438 //
439 // Calculate projected point in r,z
440 //
441 G4double r = radial.dot(p);
442
443 //
444 // Are we inside the face?
445 //
446 G4double distRZ2;
447
448 if (InsideEdges( r, p.z(), &distRZ2, 0 ))
449 {
450 //
451 // Yup, answer is just distPhi
452 //
453 return distPhi;
454 }
455 else
456 {
457 //
458 // Nope. Penalize by distance out
459 //
460 return std::sqrt( distPhi*distPhi + distRZ2 );
461 }
462}
463
464// Inside
465//
467 G4double tolerance,
468 G4double* bestDistance )
469{
470 //
471 // Get distance along phi, which if negative means the point
472 // is nominally inside the shape.
473 //
474 G4ThreeVector ps = p - surface;
475 G4double distPhi = normal.dot(ps);
476
477 //
478 // Calculate projected point in r,z
479 //
480 G4double r = radial.dot(p);
481
482 //
483 // Are we inside the face?
484 //
485 G4double distRZ2;
486 G4PolyPhiFaceVertex* base3Dnorm = nullptr;
487 G4ThreeVector* head3Dnorm = nullptr;
488
489 if (InsideEdges( r, p.z(), &distRZ2, &base3Dnorm, &head3Dnorm ))
490 {
491 //
492 // Looks like we're inside. Distance is distance in phi.
493 //
494 *bestDistance = std::fabs(distPhi);
495
496 //
497 // Use distPhi to decide fate
498 //
499 if (distPhi < -tolerance) return kInside;
500 if (distPhi < tolerance) return kSurface;
501 return kOutside;
502 }
503 else
504 {
505 //
506 // We're outside the extent of the face,
507 // so the distance is penalized by distance from edges in RZ
508 //
509 *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
510
511 //
512 // Use edge normal to decide fate
513 //
514 G4ThreeVector cc( base3Dnorm->r*radial.x(),
515 base3Dnorm->r*radial.y(),
516 base3Dnorm->z );
517 cc = p - cc;
518 G4double normDist = head3Dnorm->dot(cc);
519 if ( distRZ2 > tolerance*tolerance )
520 {
521 //
522 // We're far enough away that kSurface is not possible
523 //
524 return normDist < 0 ? kInside : kOutside;
525 }
526
527 if (normDist < -tolerance) return kInside;
528 if (normDist < tolerance) return kSurface;
529 return kOutside;
530 }
531}
532
533// Normal
534//
535// This virtual member is simple for our planer shape,
536// which has only one normal
537//
539 G4double* bestDistance )
540{
541 //
542 // Get distance along phi, which if negative means the point
543 // is nominally inside the shape.
544 //
545 G4double distPhi = normal.dot(p);
546
547 //
548 // Calculate projected point in r,z
549 //
550 G4double r = radial.dot(p);
551
552 //
553 // Are we inside the face?
554 //
555 G4double distRZ2;
556
557 if (InsideEdges( r, p.z(), &distRZ2, 0 ))
558 {
559 //
560 // Yup, answer is just distPhi
561 //
562 *bestDistance = std::fabs(distPhi);
563 }
564 else
565 {
566 //
567 // Nope. Penalize by distance out
568 //
569 *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
570 }
571
572 return normal;
573}
574
575// Extent
576//
577// This actually isn't needed by polycone or polyhedra...
578//
580{
581 G4double max = -kInfinity;
582
584 do // Loop checking, 13.08.2015, G.Cosmo
585 {
586 G4double here = axis.x()*corner->r*radial.x()
587 + axis.y()*corner->r*radial.y()
588 + axis.z()*corner->z;
589 if (here > max) max = here;
590 } while( ++corner < corners + numEdges );
591
592 return max;
593}
594
595// CalculateExtent
596//
597// See notes in G4VCSGface
598//
600 const G4VoxelLimits& voxelLimit,
601 const G4AffineTransform& transform,
602 G4SolidExtentList& extentList )
603{
604 //
605 // Construct a (sometimes big) clippable polygon,
606 //
607 // Perform the necessary transformations while doing so
608 //
609 G4ClippablePolygon polygon;
610
612 do // Loop checking, 13.08.2015, G.Cosmo
613 {
614 G4ThreeVector point( 0, 0, corner->z );
615 point += radial*corner->r;
616
617 polygon.AddVertexInOrder( transform.TransformPoint( point ) );
618 } while( ++corner < corners + numEdges );
619
620 //
621 // Clip away
622 //
623 if (polygon.PartialClip( voxelLimit, axis ))
624 {
625 //
626 // Add it to the list
627 //
628 polygon.SetNormal( transform.TransformAxis(normal) );
629 extentList.AddSurface( polygon );
630 }
631}
632
633// InsideEdgesExact
634//
635// Decide if the point in r,z is inside the edges of our face,
636// **but** do so consistently with other faces.
637//
638// This routine has functionality similar to InsideEdges, but uses
639// an algorithm to decide if a trajectory falls inside or outside the
640// face that uses only the trajectory p,v values and the three dimensional
641// points representing the edges of the polygon. The objective is to plug up
642// any leaks between touching G4PolyPhiFaces (at r==0) and any other face
643// that uses the same convention.
644//
645// See: "Computational Geometry in C (Second Edition)"
646// http://cs.smith.edu/~orourke/
647//
649 G4double normSign,
650 const G4ThreeVector& p,
651 const G4ThreeVector& v )
652{
653 //
654 // Quick check of extent
655 //
656 if ( (r < rMin-kCarTolerance)
657 || (r > rMax+kCarTolerance) ) return false;
658
659 if ( (z < zMin-kCarTolerance)
660 || (z > zMax+kCarTolerance) ) return false;
661
662 //
663 // Exact check: loop over all vertices
664 //
665 G4double qx = p.x() + v.x(),
666 qy = p.y() + v.y(),
667 qz = p.z() + v.z();
668
669 G4int answer = 0;
671 *prev = corners+numEdges-1;
672
673 G4double cornZ, prevZ;
674
675 prevZ = ExactZOrder( z, qx, qy, qz, v, normSign, prev );
676 do // Loop checking, 13.08.2015, G.Cosmo
677 {
678 //
679 // Get z order of this vertex, and compare to previous vertex
680 //
681 cornZ = ExactZOrder( z, qx, qy, qz, v, normSign, corn );
682
683 if (cornZ < 0)
684 {
685 if (prevZ < 0) continue;
686 }
687 else if (cornZ > 0)
688 {
689 if (prevZ > 0) continue;
690 }
691 else
692 {
693 //
694 // By chance, we overlap exactly (within precision) with
695 // the current vertex. Continue if the same happened previously
696 // (e.g. the previous vertex had the same z value)
697 //
698 if (prevZ == 0) continue;
699
700 //
701 // Otherwise, to decide what to do, we need to know what is
702 // coming up next. Specifically, we need to find the next vertex
703 // with a non-zero z order.
704 //
705 // One might worry about infinite loops, but the above conditional
706 // should prevent it
707 //
708 G4PolyPhiFaceVertex *next = corn;
709 G4double nextZ;
710 do // Loop checking, 13.08.2015, G.Cosmo
711 {
712 ++next;
713 if (next == corners+numEdges) next = corners;
714
715 nextZ = ExactZOrder( z, qx, qy, qz, v, normSign, next );
716 } while( nextZ == 0 );
717
718 //
719 // If we won't be changing direction, go to the next vertex
720 //
721 if (nextZ*prevZ < 0) continue;
722 }
723
724
725 //
726 // We overlap in z with the side of the face that stretches from
727 // vertex "prev" to "corn". On which side (left or right) do
728 // we lay with respect to this segment?
729 //
730 G4ThreeVector qa( qx - prev->x, qy - prev->y, qz - prev->z ),
731 qb( qx - corn->x, qy - corn->y, qz - corn->z );
732
733 G4double aboveOrBelow = normSign*qa.cross(qb).dot(v);
734
735 if (aboveOrBelow > 0)
736 ++answer;
737 else if (aboveOrBelow < 0)
738 --answer;
739 else
740 {
741 //
742 // A precisely zero answer here means we exactly
743 // intersect (within roundoff) the edge of the face.
744 // Return true in this case.
745 //
746 return true;
747 }
748 } while( prevZ = cornZ, prev=corn, ++corn < corners+numEdges );
749
750 return answer!=0;
751}
752
753// InsideEdges (don't care about distance)
754//
755// Decide if the point in r,z is inside the edges of our face
756//
757// This routine can be made a zillion times quicker by implementing
758// better code, for example:
759//
760// int pnpoly(int npol, float *xp, float *yp, float x, float y)
761// {
762// int i, j, c = 0;
763// for (i = 0, j = npol-1; i < npol; j = i++) {
764// if ((((yp[i]<=y) && (y<yp[j])) ||
765// ((yp[j]<=y) && (y<yp[i]))) &&
766// (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
767//
768// c = !c;
769// }
770// return c;
771// }
772//
773// See "Point in Polyon Strategies", Eric Haines [Graphic Gems IV] pp. 24-46
774//
775// The algorithm below is rather unique, but is based on code needed to
776// calculate the distance to the shape. Left here for testing...
777//
779{
780 //
781 // Quick check of extent
782 //
783 if ( r < rMin || r > rMax ) return false;
784 if ( z < zMin || z > zMax ) return false;
785
786 //
787 // More thorough check
788 //
789 G4double notUsed;
790
791 return InsideEdges( r, z, &notUsed, 0 );
792}
793
794// InsideEdges (care about distance)
795//
796// Decide if the point in r,z is inside the edges of our face
797//
799 G4double* bestDist2,
800 G4PolyPhiFaceVertex** base3Dnorm,
801 G4ThreeVector** head3Dnorm )
802{
803 G4double bestDistance2 = kInfinity;
804 G4bool answer = false;
805
806 G4PolyPhiFaceEdge* edge = edges;
807 do // Loop checking, 13.08.2015, G.Cosmo
808 {
809 G4PolyPhiFaceVertex* testMe;
810 //
811 // Get distance perpendicular to the edge
812 //
813 G4double dr = (r-edge->v0->r), dz = (z-edge->v0->z);
814
815 G4double distOut = dr*edge->tz - dz*edge->tr;
816 G4double distance2 = distOut*distOut;
817 if (distance2 > bestDistance2) continue; // No hope!
818
819 //
820 // Check to see if normal intersects edge within the edge's boundary
821 //
822 G4double q = dr*edge->tr + dz*edge->tz;
823
824 //
825 // If it doesn't, penalize distance2 appropriately
826 //
827 if (q < 0)
828 {
829 distance2 += q*q;
830 testMe = edge->v0;
831 }
832 else if (q > edge->length)
833 {
834 G4double s2 = q-edge->length;
835 distance2 += s2*s2;
836 testMe = edge->v1;
837 }
838 else
839 {
840 testMe = nullptr;
841 }
842
843 //
844 // Closest edge so far?
845 //
846 if (distance2 < bestDistance2)
847 {
848 bestDistance2 = distance2;
849 if (testMe)
850 {
851 G4double distNorm = dr*testMe->rNorm + dz*testMe->zNorm;
852 answer = (distNorm <= 0);
853 if (base3Dnorm)
854 {
855 *base3Dnorm = testMe;
856 *head3Dnorm = &testMe->norm3D;
857 }
858 }
859 else
860 {
861 answer = (distOut <= 0);
862 if (base3Dnorm)
863 {
864 *base3Dnorm = edge->v0;
865 *head3Dnorm = &edge->norm3D;
866 }
867 }
868 }
869 } while( ++edge < edges + numEdges );
870
871 *bestDist2 = bestDistance2;
872 return answer;
873}
874
875// Calculation of Surface Area of a Triangle
876// In the same time Random Point in Triangle is given
877//
879 G4ThreeVector p2,
880 G4ThreeVector p3,
881 G4ThreeVector* p4 )
882{
883 G4ThreeVector v, w;
884
885 v = p3 - p1;
886 w = p1 - p2;
887 G4double lambda1 = G4UniformRand();
888 G4double lambda2 = lambda1*G4UniformRand();
889
890 *p4=p2 + lambda1*w + lambda2*v;
891 return 0.5*(v.cross(w)).mag();
892}
893
894// Compute surface area
895//
897{
898 if ( fSurfaceArea==0. ) { Triangulate(); }
899 return fSurfaceArea;
900}
901
902// Return random point on face
903//
905{
906 Triangulate();
907 return surface_point;
908}
909
910//
911// Auxiliary Functions used for Finding the PointOnFace using Triangulation
912//
913
914// Calculation of 2*Area of Triangle with Sign
915//
917 G4TwoVector b,
918 G4TwoVector c )
919{
920 return ((b.x()-a.x())*(c.y()-a.y())-
921 (c.x()-a.x())*(b.y()-a.y()));
922}
923
924// Boolean function for sign of Surface
925//
927 G4TwoVector b,
928 G4TwoVector c )
929{
930 return Area2(a,b,c)>0;
931}
932
933// Boolean function for sign of Surface
934//
936 G4TwoVector b,
937 G4TwoVector c )
938{
939 return Area2(a,b,c)>=0;
940}
941
942// Boolean function for sign of Surface
943//
945 G4TwoVector b,
946 G4TwoVector c )
947{
948 return Area2(a,b,c)==0;
949}
950
951// Boolean function for finding "Proper" Intersection
952// That means Intersection of two lines segments (a,b) and (c,d)
953//
955 G4TwoVector b,
957{
958 if( Collinear(a,b,c) || Collinear(a,b,d)||
959 Collinear(c,d,a) || Collinear(c,d,b) ) { return false; }
960
961 G4bool Positive;
962 Positive = !(Left(a,b,c))^!(Left(a,b,d));
963 return Positive && (!Left(c,d,a)^!Left(c,d,b));
964}
965
966// Boolean function for determining if Point c is between a and b
967// For the tree points(a,b,c) on the same line
968//
970{
971 if( !Collinear(a,b,c) ) { return false; }
972
973 if(a.x()!=b.x())
974 {
975 return ((a.x()<=c.x())&&(c.x()<=b.x()))||
976 ((a.x()>=c.x())&&(c.x()>=b.x()));
977 }
978 else
979 {
980 return ((a.y()<=c.y())&&(c.y()<=b.y()))||
981 ((a.y()>=c.y())&&(c.y()>=b.y()));
982 }
983}
984
985// Boolean function for finding Intersection "Proper" or not
986// Between two line segments (a,b) and (c,d)
987//
989 G4TwoVector b,
991{
992 if( IntersectProp(a,b,c,d) )
993 { return true; }
994 else if( Between(a,b,c)||
995 Between(a,b,d)||
996 Between(c,d,a)||
997 Between(c,d,b) )
998 { return true; }
999 else
1000 { return false; }
1001}
1002
1003// Boolean Diagonalie help to determine
1004// if diagonal s of segment (a,b) is convex or reflex
1005//
1008{
1010 G4PolyPhiFaceVertex *corner_next=triangles;
1011
1012 // For each Edge (corner,corner_next)
1013 do // Loop checking, 13.08.2015, G.Cosmo
1014 {
1015 corner_next=corner->next;
1016
1017 // Skip edges incident to a of b
1018 //
1019 if( (corner!=a)&&(corner_next!=a)
1020 &&(corner!=b)&&(corner_next!=b) )
1021 {
1022 G4TwoVector rz1,rz2,rz3,rz4;
1023 rz1 = G4TwoVector(a->r,a->z);
1024 rz2 = G4TwoVector(b->r,b->z);
1025 rz3 = G4TwoVector(corner->r,corner->z);
1026 rz4 = G4TwoVector(corner_next->r,corner_next->z);
1027 if( Intersect(rz1,rz2,rz3,rz4) ) { return false; }
1028 }
1029 corner=corner->next;
1030
1031 } while( corner != triangles );
1032
1033 return true;
1034}
1035
1036// Boolean function that determine if b is Inside Cone (a0,a,a1)
1037// being a the center of the Cone
1038//
1040{
1041 // a0,a and a1 are consecutive vertices
1042 //
1044 a1=a->next;
1045 a0=a->prev;
1046
1047 G4TwoVector arz,arz0,arz1,brz;
1048 arz=G4TwoVector(a->r,a->z);arz0=G4TwoVector(a0->r,a0->z);
1049 arz1=G4TwoVector(a1->r,a1->z);brz=G4TwoVector(b->r,b->z);
1050
1051
1052 if(LeftOn(arz,arz1,arz0)) // If a is convex vertex
1053 {
1054 return Left(arz,brz,arz0)&&Left(brz,arz,arz1);
1055 }
1056 else // Else a is reflex
1057 {
1058 return !( LeftOn(arz,brz,arz1)&&LeftOn(brz,arz,arz0));
1059 }
1060}
1061
1062// Boolean function finding if Diagonal is possible
1063// inside Polycone or PolyHedra
1064//
1066{
1067 return InCone(a,b) && InCone(b,a) && Diagonalie(a,b);
1068}
1069
1070// Initialisation for Triangulisation by ear tips
1071// For details see "Computational Geometry in C" by Joseph O'Rourke
1072//
1074{
1076 G4PolyPhiFaceVertex* c_prev, *c_next;
1077
1078 do // Loop checking, 13.08.2015, G.Cosmo
1079 {
1080 // We need to determine three consecutive vertices
1081 //
1082 c_next=corner->next;
1083 c_prev=corner->prev;
1084
1085 // Calculation of ears
1086 //
1087 corner->ear=Diagonal(c_prev,c_next);
1088 corner=corner->next;
1089
1090 } while( corner!=triangles );
1091}
1092
1093// Triangulisation by ear tips for Polycone or Polyhedra
1094// For details see "Computational Geometry in C" by Joseph O'Rourke
1095//
1097{
1098 // The copy of Polycone is made and this copy is reordered in order to
1099 // have a list of triangles. This list is used for GetPointOnFace().
1100
1102 triangles = tri_help;
1104
1105 std::vector<G4double> areas;
1106 std::vector<G4ThreeVector> points;
1107 G4double area=0.;
1108 G4PolyPhiFaceVertex *v0,*v1,*v2,*v3,*v4;
1109 v2=triangles;
1110
1111 // Make copy for prev/next for triang=corners
1112 //
1113 G4PolyPhiFaceVertex* helper = corners;
1114 G4PolyPhiFaceVertex* helper2 = corners;
1115 do // Loop checking, 13.08.2015, G.Cosmo
1116 {
1117 triang->r = helper->r;
1118 triang->z = helper->z;
1119 triang->x = helper->x;
1120 triang->y= helper->y;
1121
1122 // add pointer on prev corner
1123 //
1124 if( helper==corners )
1125 { triang->prev=triangles+numEdges-1; }
1126 else
1127 { triang->prev=helper2; }
1128
1129 // add pointer on next corner
1130 //
1131 if( helper<corners+numEdges-1 )
1132 { triang->next=triang+1; }
1133 else
1134 { triang->next=triangles; }
1135 helper2=triang;
1136 helper=helper->next;
1137 triang=triang->next;
1138
1139 } while( helper!=corners );
1140
1141 EarInit();
1142
1143 G4int n=numEdges;
1144 G4int i=0;
1145 G4ThreeVector p1,p2,p3,p4;
1146 const G4int max_n_loops=numEdges*10000; // protection against infinite loop
1147
1148 // Each step of outer loop removes one ear
1149 //
1150 while(n>3) // Loop checking, 13.08.2015, G.Cosmo
1151 { // Inner loop searches for one ear
1152 v2=triangles;
1153 do // Loop checking, 13.08.2015, G.Cosmo
1154 {
1155 if(v2->ear) // Ear found. Fill variables
1156 {
1157 // (v1,v3) is diagonal
1158 //
1159 v3=v2->next; v4=v3->next;
1160 v1=v2->prev; v0=v1->prev;
1161
1162 // Calculate areas and points
1163
1164 p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
1165 p2=G4ThreeVector((v1)->x,(v1)->y,(v1)->z);
1166 p3=G4ThreeVector((v3)->x,(v3)->y,(v3)->z);
1167
1168 G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
1169 points.push_back(p4);
1170 areas.push_back(result1);
1171 area=area+result1;
1172
1173 // Update earity of diagonal endpoints
1174 //
1175 v1->ear=Diagonal(v0,v3);
1176 v3->ear=Diagonal(v1,v4);
1177
1178 // Cut off the ear v2
1179 // Has to be done for a copy and not for real PolyPhiFace
1180 //
1181 v1->next=v3;
1182 v3->prev=v1;
1183 triangles=v3; // In case the head was v2
1184 --n;
1185
1186 break; // out of inner loop
1187 } // end if ear found
1188
1189 v2=v2->next;
1190
1191 } while( v2!=triangles );
1192
1193 ++i;
1194 if(i>=max_n_loops)
1195 {
1196 G4Exception( "G4PolyPhiFace::Triangulation()",
1197 "GeomSolids0003", FatalException,
1198 "Maximum number of steps is reached for triangulation!" );
1199 }
1200 } // end outer while loop
1201
1202 if(v2->next)
1203 {
1204 // add last triangle
1205 //
1206 v2=v2->next;
1207 p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
1208 p2=G4ThreeVector((v2->next)->x,(v2->next)->y,(v2->next)->z);
1209 p3=G4ThreeVector((v2->prev)->x,(v2->prev)->y,(v2->prev)->z);
1210 G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
1211 points.push_back(p4);
1212 areas.push_back(result1);
1213 area=area+result1;
1214 }
1215
1216 // Surface Area is stored
1217 //
1218 fSurfaceArea = area;
1219
1220 // Second Step: choose randomly one surface
1221 //
1222 G4double chose = area*G4UniformRand();
1223
1224 // Third Step: Get a point on choosen surface
1225 //
1226 G4double Achose1, Achose2;
1227 Achose1=0; Achose2=0.;
1228 i=0;
1229 do // Loop checking, 13.08.2015, G.Cosmo
1230 {
1231 Achose2+=areas[i];
1232 if(chose>=Achose1 && chose<Achose2)
1233 {
1234 G4ThreeVector point;
1235 point=points[i] ;
1236 surface_point=point;
1237 break;
1238 }
1239 i++; Achose1=Achose2;
1240 } while( i<numEdges-2 );
1241
1242 delete [] tri_help;
1243 tri_help = nullptr;
1244}
const G4double kCarTolerance
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
const G4double a0
CLHEP::Hep3Vector G4ThreeVector
CLHEP::Hep2Vector G4TwoVector
Definition: G4TwoVector.hh:36
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
#define G4UniformRand()
Definition: Randomize.hh:52
double x() const
double y() const
double z() const
Hep3Vector unit() const
double x() const
double y() const
Hep3Vector cross(const Hep3Vector &) const
double dot(const Hep3Vector &) const
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
virtual G4bool PartialClip(const G4VoxelLimits &voxelLimit, const EAxis IgnoreMe)
virtual void AddVertexInOrder(const G4ThreeVector vertex)
void SetNormal(const G4ThreeVector &newNormal)
G4double GetSurfaceTolerance() const
static G4GeometryTolerance * GetInstance()
G4double SurfaceArea()
G4bool Left(G4TwoVector a, G4TwoVector b, G4TwoVector c)
G4bool Collinear(G4TwoVector a, G4TwoVector b, G4TwoVector c)
void Diagnose(G4VSolid *solid)
G4PolyPhiFace & operator=(const G4PolyPhiFace &source)
virtual ~G4PolyPhiFace()
G4bool Diagonal(G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b)
G4ThreeVector surface_point
G4bool InsideEdges(G4double r, G4double z)
G4bool Diagonalie(G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b)
G4bool IntersectProp(G4TwoVector a, G4TwoVector b, G4TwoVector c, G4TwoVector d)
G4ThreeVector surface
G4bool LeftOn(G4TwoVector a, G4TwoVector b, G4TwoVector c)
G4bool Intersect(const G4ThreeVector &p, const G4ThreeVector &v, G4bool outgoing, G4double surfTolerance, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal, G4bool &allBehind)
G4ThreeVector Normal(const G4ThreeVector &p, G4double *bestDistance)
G4double fSurfaceArea
void CopyStuff(const G4PolyPhiFace &source)
void CalculateExtent(const EAxis axis, const G4VoxelLimits &voxelLimit, const G4AffineTransform &tranform, G4SolidExtentList &extentList)
G4double Area2(G4TwoVector a, G4TwoVector b, G4TwoVector c)
G4double ExactZOrder(G4double z, G4double qx, G4double qy, G4double qz, const G4ThreeVector &v, G4double normSign, const G4PolyPhiFaceVertex *vert) const
G4double Distance(const G4ThreeVector &p, G4bool outgoing)
G4double Extent(const G4ThreeVector axis)
G4PolyPhiFaceEdge * edges
G4ThreeVector radial
G4bool InCone(G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b)
G4bool Between(G4TwoVector a, G4TwoVector b, G4TwoVector c)
G4PolyPhiFaceVertex * corners
G4PolyPhiFaceVertex * triangles
G4ThreeVector normal
EInside Inside(const G4ThreeVector &p, G4double tolerance, G4double *bestDistance)
G4double kCarTolerance
G4ThreeVector GetPointOnFace()
G4double SurfaceTriangle(G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3, G4ThreeVector *p4)
G4bool InsideEdgesExact(G4double r, G4double z, G4double normSign, const G4ThreeVector &p, const G4ThreeVector &v)
G4PolyPhiFace(const G4ReduciblePolygon *rz, G4double phi, G4double deltaPhi, G4double phiOther)
G4double Amin() const
G4int NumVertices() const
G4double Bmin() const
G4double Bmax() const
G4double Amax() const
void AddSurface(const G4ClippablePolygon &surface)
virtual EInside Inside(const G4ThreeVector &p) const =0
EAxis
Definition: geomdefs.hh:54
EInside
Definition: geomdefs.hh:67
@ kInside
Definition: geomdefs.hh:70
@ kOutside
Definition: geomdefs.hh:68
@ kSurface
Definition: geomdefs.hh:69
G4PolyPhiFaceVertex * v1
G4ThreeVector norm3D
G4PolyPhiFaceVertex * v0
G4ThreeVector norm3D
G4PolyPhiFaceVertex * next
G4PolyPhiFaceVertex * prev
#define DBL_MIN
Definition: templates.hh:54