| Geant4 Cross Reference |
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 (davidw@scipp.ucsc.edu)
30 // --------------------------------------------------------------------
31
32 #include "G4PolyPhiFace.hh"
33 #include "G4ClippablePolygon.hh"
34 #include "G4ReduciblePolygon.hh"
35 #include "G4AffineTransform.hh"
36 #include "G4SolidExtentList.hh"
37 #include "G4GeometryTolerance.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 //
51 G4PolyPhiFace::G4PolyPhiFace( const G4ReduciblePolygon* rz,
52 G4double phi,
53 G4double deltaPhi,
54 G4double phiOther )
55 {
56 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
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 //
95 const std::size_t maxEdges = numEdges>0 ? numEdges : 1;
96 corners = new G4PolyPhiFaceVertex[maxEdges];
97 //
98 // Fill them
99 //
100 G4ReduciblePolygonIterator iterRZ(rz);
101
102 G4PolyPhiFaceVertex* corn = corners;
103 G4PolyPhiFaceVertex* helper = corners;
104
105 iterRZ.Begin();
106 do // Loop checking, 13.08.2015, G.Cosmo
107 {
108 corn->r = iterRZ.GetA();
109 corn->z = iterRZ.GetB();
110 corn->x = corn->r*radial.x();
111 corn->y = corn->r*radial.y();
112
113 // Add pointer on prev corner
114 //
115 if( corn == corners )
116 { corn->prev = corners+maxEdges-1;}
117 else
118 { corn->prev = helper; }
119
120 // Add pointer on next corner
121 //
122 if( corn < corners+maxEdges-1 )
123 { corn->next = corn+1;}
124 else
125 { corn->next = corners; }
126
127 helper = corn;
128 } while( ++corn, iterRZ.Next() );
129
130 //
131 // Allocate edges
132 //
133 edges = new G4PolyPhiFaceEdge[maxEdges];
134
135 //
136 // Fill them
137 //
138 G4double rFact = std::cos(0.5*deltaPhi);
139 G4double rFactNormalize = 1.0/std::sqrt(1.0+rFact*rFact);
140
141 G4PolyPhiFaceVertex* prev = corners+maxEdges-1,
142 * here = corners;
143 G4PolyPhiFaceEdge* edge = edges;
144 do // Loop checking, 13.08.2015, G.Cosmo
145 {
146 G4ThreeVector sideNorm;
147
148 edge->v0 = prev;
149 edge->v1 = here;
150
151 G4double dr = here->r - prev->r,
152 dz = here->z - prev->z;
153
154 edge->length = std::sqrt( dr*dr + dz*dz );
155
156 edge->tr = dr/edge->length;
157 edge->tz = dz/edge->length;
158
159 if ((here->r < DBL_MIN) && (prev->r < DBL_MIN))
160 {
161 //
162 // Sigh! Always exceptions!
163 // This edge runs at r==0, so its adjoing surface is not a
164 // PolyconeSide or PolyhedraSide, but the opposite PolyPhiFace.
165 //
166 G4double zSignOther = start ? -1 : 1;
167 sideNorm = G4ThreeVector( zSignOther*std::sin(phiOther),
168 -zSignOther*std::cos(phiOther), 0 );
169 }
170 else
171 {
172 sideNorm = G4ThreeVector( edge->tz*cosMid,
173 edge->tz*sinMid,
174 -edge->tr*rFact );
175 sideNorm *= rFactNormalize;
176 }
177 sideNorm += normal;
178
179 edge->norm3D = sideNorm.unit();
180 } while( edge++, prev=here, ++here < corners+maxEdges );
181
182 //
183 // Go back and fill in corner "normals"
184 //
185 G4PolyPhiFaceEdge* prevEdge = edges+maxEdges-1;
186 edge = edges;
187 do // Loop checking, 13.08.2015, G.Cosmo
188 {
189 //
190 // Calculate vertex 2D normals (on the phi surface)
191 //
192 G4double rPart = prevEdge->tr + edge->tr;
193 G4double zPart = prevEdge->tz + edge->tz;
194 G4double norm = std::sqrt( rPart*rPart + zPart*zPart );
195 G4double rNorm = +zPart/norm;
196 G4double zNorm = -rPart/norm;
197
198 edge->v0->rNorm = rNorm;
199 edge->v0->zNorm = zNorm;
200
201 //
202 // Calculate the 3D normals.
203 //
204 // Find the vector perpendicular to the z axis
205 // that defines the plane that contains the vertex normal
206 //
207 G4ThreeVector xyVector;
208
209 if (edge->v0->r < DBL_MIN)
210 {
211 //
212 // This is a vertex at r==0, which is a special
213 // case. The normal we will construct lays in the
214 // plane at the center of the phi opening.
215 //
216 // We also know that rNorm < 0
217 //
218 G4double zSignOther = start ? -1 : 1;
219 G4ThreeVector normalOther( zSignOther*std::sin(phiOther),
220 -zSignOther*std::cos(phiOther), 0 );
221
222 xyVector = - normal - normalOther;
223 }
224 else
225 {
226 //
227 // This is a vertex at r > 0. The plane
228 // is the average of the normal and the
229 // normal of the adjacent phi face
230 //
231 xyVector = G4ThreeVector( cosMid, sinMid, 0 );
232 if (rNorm < 0)
233 xyVector -= normal;
234 else
235 xyVector += normal;
236 }
237
238 //
239 // Combine it with the r/z direction from the face
240 //
241 edge->v0->norm3D = rNorm*xyVector.unit() + G4ThreeVector( 0, 0, zNorm );
242 } while( prevEdge=edge, ++edge < edges+maxEdges );
243
244 //
245 // Build point on surface
246 //
247 G4double rAve = 0.5*(rMax-rMin),
248 zAve = 0.5*(zMax-zMin);
249 surface = G4ThreeVector( rAve*radial.x(), rAve*radial.y(), zAve );
250 }
251
252 // Diagnose
253 //
254 // Throw an exception if something is found inconsistent with
255 // the solid.
256 //
257 // For debugging purposes only
258 //
259 void G4PolyPhiFace::Diagnose( G4VSolid* owner )
260 {
261 G4PolyPhiFaceVertex* corner = corners;
262 do // Loop checking, 13.08.2015, G.Cosmo
263 {
264 G4ThreeVector test(corner->x, corner->y, corner->z);
265 test -= 1E-6*corner->norm3D;
266
267 if (owner->Inside(test) != kInside)
268 G4Exception( "G4PolyPhiFace::Diagnose()", "GeomSolids0002",
269 FatalException, "Bad vertex normal found." );
270 } while( ++corner < corners+numEdges );
271 }
272
273 // Fake default constructor - sets only member data and allocates memory
274 // for usage restricted to object persistency.
275 //
276 G4PolyPhiFace::G4PolyPhiFace( __void__&)
277 : numEdges(0), rMin(0.), rMax(0.), zMin(0.), zMax(0.), kCarTolerance(0.)
278 {
279 }
280
281
282 //
283 // Destructor
284 //
285 G4PolyPhiFace::~G4PolyPhiFace()
286 {
287 delete [] edges;
288 delete [] corners;
289 }
290
291
292 //
293 // Copy constructor
294 //
295 G4PolyPhiFace::G4PolyPhiFace( const G4PolyPhiFace& source )
296 {
297 CopyStuff( source );
298 }
299
300
301 //
302 // Assignment operator
303 //
304 G4PolyPhiFace& G4PolyPhiFace::operator=( const G4PolyPhiFace& source )
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 //
320 void G4PolyPhiFace::CopyStuff( const G4PolyPhiFace& source )
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
336 kCarTolerance = source.kCarTolerance;
337 fSurfaceArea = source.fSurfaceArea;
338
339 const std::size_t maxEdges = (numEdges > 0) ? numEdges : 1;
340
341 //
342 // Corner dynamic array
343 //
344 corners = new G4PolyPhiFaceVertex[maxEdges];
345 G4PolyPhiFaceVertex *corn = corners,
346 *sourceCorn = source.corners;
347 do // Loop checking, 13.08.2015, G.Cosmo
348 {
349 *corn = *sourceCorn;
350 } while( ++sourceCorn, ++corn < corners+maxEdges );
351
352 //
353 // Edge dynamic array
354 //
355 edges = new G4PolyPhiFaceEdge[maxEdges];
356
357 G4PolyPhiFaceVertex* prev = corners+maxEdges-1,
358 * here = corners;
359 G4PolyPhiFaceEdge* edge = edges,
360 * sourceEdge = source.edges;
361 do // Loop checking, 13.08.2015, G.Cosmo
362 {
363 *edge = *sourceEdge;
364 edge->v0 = prev;
365 edge->v1 = here;
366 } while( ++sourceEdge, ++edge, prev=here, ++here < corners+maxEdges );
367 }
368
369 // Intersect
370 //
371 G4bool G4PolyPhiFace::Intersect( const G4ThreeVector& p,
372 const G4ThreeVector& v,
373 G4bool outgoing,
374 G4double surfTolerance,
375 G4double& distance,
376 G4double& distFromSurface,
377 G4ThreeVector& aNormal,
378 G4bool& isAllBehind )
379 {
380 G4double normSign = outgoing ? +1 : -1;
381
382 //
383 // These don't change
384 //
385 isAllBehind = allBehind;
386 aNormal = normal;
387
388 //
389 // Correct normal? Here we have straight sides, and can safely ignore
390 // intersections where the dot product with the normal is zero.
391 //
392 G4double dotProd = normSign*normal.dot(v);
393
394 if (dotProd <= 0) return false;
395
396 //
397 // Calculate distance to surface. If the side is too far
398 // behind the point, we must reject it.
399 //
400 G4ThreeVector ps = p - surface;
401 distFromSurface = -normSign*ps.dot(normal);
402
403 if (distFromSurface < -surfTolerance) return false;
404
405 //
406 // Calculate precise distance to intersection with the side
407 // (along the trajectory, not normal to the surface)
408 //
409 distance = distFromSurface/dotProd;
410
411 //
412 // Calculate intersection point in r,z
413 //
414 G4ThreeVector ip = p + distance*v;
415
416 G4double r = radial.dot(ip);
417
418 //
419 // And is it inside the r/z extent?
420 //
421 return InsideEdgesExact( r, ip.z(), normSign, p, v );
422 }
423
424 // Distance
425 //
426 G4double G4PolyPhiFace::Distance( const G4ThreeVector& p, G4bool outgoing )
427 {
428 G4double normSign = outgoing ? +1 : -1;
429 //
430 // Correct normal?
431 //
432 G4ThreeVector ps = p - surface;
433 G4double distPhi = -normSign*normal.dot(ps);
434
435 if (distPhi < -0.5*kCarTolerance)
436 return kInfinity;
437 else if (distPhi < 0)
438 distPhi = 0.0;
439
440 //
441 // Calculate projected point in r,z
442 //
443 G4double r = radial.dot(p);
444
445 //
446 // Are we inside the face?
447 //
448 G4double distRZ2;
449
450 if (InsideEdges( r, p.z(), &distRZ2, nullptr ))
451 {
452 //
453 // Yup, answer is just distPhi
454 //
455 return distPhi;
456 }
457 else
458 {
459 //
460 // Nope. Penalize by distance out
461 //
462 return std::sqrt( distPhi*distPhi + distRZ2 );
463 }
464 }
465
466 // Inside
467 //
468 EInside G4PolyPhiFace::Inside( const G4ThreeVector& p,
469 G4double tolerance,
470 G4double* bestDistance )
471 {
472 //
473 // Get distance along phi, which if negative means the point
474 // is nominally inside the shape.
475 //
476 G4ThreeVector ps = p - surface;
477 G4double distPhi = normal.dot(ps);
478
479 //
480 // Calculate projected point in r,z
481 //
482 G4double r = radial.dot(p);
483
484 //
485 // Are we inside the face?
486 //
487 G4double distRZ2;
488 G4PolyPhiFaceVertex* base3Dnorm = nullptr;
489 G4ThreeVector* head3Dnorm = nullptr;
490
491 if (InsideEdges( r, p.z(), &distRZ2, &base3Dnorm, &head3Dnorm ))
492 {
493 //
494 // Looks like we're inside. Distance is distance in phi.
495 //
496 *bestDistance = std::fabs(distPhi);
497
498 //
499 // Use distPhi to decide fate
500 //
501 if (distPhi < -tolerance) return kInside;
502 if (distPhi < tolerance) return kSurface;
503 return kOutside;
504 }
505 else
506 {
507 //
508 // We're outside the extent of the face,
509 // so the distance is penalized by distance from edges in RZ
510 //
511 *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
512
513 //
514 // Use edge normal to decide fate
515 //
516 G4ThreeVector cc( base3Dnorm->r*radial.x(),
517 base3Dnorm->r*radial.y(),
518 base3Dnorm->z );
519 cc = p - cc;
520 G4double normDist = head3Dnorm->dot(cc);
521 if ( distRZ2 > tolerance*tolerance )
522 {
523 //
524 // We're far enough away that kSurface is not possible
525 //
526 return normDist < 0 ? kInside : kOutside;
527 }
528
529 if (normDist < -tolerance) return kInside;
530 if (normDist < tolerance) return kSurface;
531 return kOutside;
532 }
533 }
534
535 // Normal
536 //
537 // This virtual member is simple for our planer shape,
538 // which has only one normal
539 //
540 G4ThreeVector G4PolyPhiFace::Normal( const G4ThreeVector& p,
541 G4double* bestDistance )
542 {
543 //
544 // Get distance along phi, which if negative means the point
545 // is nominally inside the shape.
546 //
547 G4double distPhi = normal.dot(p);
548
549 //
550 // Calculate projected point in r,z
551 //
552 G4double r = radial.dot(p);
553
554 //
555 // Are we inside the face?
556 //
557 G4double distRZ2;
558
559 if (InsideEdges( r, p.z(), &distRZ2, nullptr ))
560 {
561 //
562 // Yup, answer is just distPhi
563 //
564 *bestDistance = std::fabs(distPhi);
565 }
566 else
567 {
568 //
569 // Nope. Penalize by distance out
570 //
571 *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
572 }
573
574 return normal;
575 }
576
577 // Extent
578 //
579 // This actually isn't needed by polycone or polyhedra...
580 //
581 G4double G4PolyPhiFace::Extent( const G4ThreeVector axis )
582 {
583 G4double max = -kInfinity;
584
585 G4PolyPhiFaceVertex* corner = corners;
586 do // Loop checking, 13.08.2015, G.Cosmo
587 {
588 G4double here = axis.x()*corner->r*radial.x()
589 + axis.y()*corner->r*radial.y()
590 + axis.z()*corner->z;
591 if (here > max) max = here;
592 } while( ++corner < corners + numEdges );
593
594 return max;
595 }
596
597 // CalculateExtent
598 //
599 // See notes in G4VCSGface
600 //
601 void G4PolyPhiFace::CalculateExtent( const EAxis axis,
602 const G4VoxelLimits& voxelLimit,
603 const G4AffineTransform& transform,
604 G4SolidExtentList& extentList )
605 {
606 //
607 // Construct a (sometimes big) clippable polygon,
608 //
609 // Perform the necessary transformations while doing so
610 //
611 G4ClippablePolygon polygon;
612
613 G4PolyPhiFaceVertex* corner = corners;
614 do // Loop checking, 13.08.2015, G.Cosmo
615 {
616 G4ThreeVector point( 0, 0, corner->z );
617 point += radial*corner->r;
618
619 polygon.AddVertexInOrder( transform.TransformPoint( point ) );
620 } while( ++corner < corners + numEdges );
621
622 //
623 // Clip away
624 //
625 if (polygon.PartialClip( voxelLimit, axis ))
626 {
627 //
628 // Add it to the list
629 //
630 polygon.SetNormal( transform.TransformAxis(normal) );
631 extentList.AddSurface( polygon );
632 }
633 }
634
635 // InsideEdgesExact
636 //
637 // Decide if the point in r,z is inside the edges of our face,
638 // **but** do so consistently with other faces.
639 //
640 // This routine has functionality similar to InsideEdges, but uses
641 // an algorithm to decide if a trajectory falls inside or outside the
642 // face that uses only the trajectory p,v values and the three dimensional
643 // points representing the edges of the polygon. The objective is to plug up
644 // any leaks between touching G4PolyPhiFaces (at r==0) and any other face
645 // that uses the same convention.
646 //
647 // See: "Computational Geometry in C (Second Edition)"
648 // http://cs.smith.edu/~orourke/
649 //
650 G4bool G4PolyPhiFace::InsideEdgesExact( G4double r, G4double z,
651 G4double normSign,
652 const G4ThreeVector& p,
653 const G4ThreeVector& v )
654 {
655 //
656 // Quick check of extent
657 //
658 if ( (r < rMin-kCarTolerance)
659 || (r > rMax+kCarTolerance) ) return false;
660
661 if ( (z < zMin-kCarTolerance)
662 || (z > zMax+kCarTolerance) ) return false;
663
664 //
665 // Exact check: loop over all vertices
666 //
667 G4double qx = p.x() + v.x(),
668 qy = p.y() + v.y(),
669 qz = p.z() + v.z();
670
671 G4int answer = 0;
672 G4PolyPhiFaceVertex *corn = corners,
673 *prev = corners+numEdges-1;
674
675 G4double cornZ, prevZ;
676
677 prevZ = ExactZOrder( z, qx, qy, qz, v, normSign, prev );
678 do // Loop checking, 13.08.2015, G.Cosmo
679 {
680 //
681 // Get z order of this vertex, and compare to previous vertex
682 //
683 cornZ = ExactZOrder( z, qx, qy, qz, v, normSign, corn );
684
685 if (cornZ < 0)
686 {
687 if (prevZ < 0) continue;
688 }
689 else if (cornZ > 0)
690 {
691 if (prevZ > 0) continue;
692 }
693 else
694 {
695 //
696 // By chance, we overlap exactly (within precision) with
697 // the current vertex. Continue if the same happened previously
698 // (e.g. the previous vertex had the same z value)
699 //
700 if (prevZ == 0) continue;
701
702 //
703 // Otherwise, to decide what to do, we need to know what is
704 // coming up next. Specifically, we need to find the next vertex
705 // with a non-zero z order.
706 //
707 // One might worry about infinite loops, but the above conditional
708 // should prevent it
709 //
710 G4PolyPhiFaceVertex *next = corn;
711 G4double nextZ;
712 do // Loop checking, 13.08.2015, G.Cosmo
713 {
714 ++next;
715 if (next == corners+numEdges) next = corners;
716
717 nextZ = ExactZOrder( z, qx, qy, qz, v, normSign, next );
718 } while( nextZ == 0 );
719
720 //
721 // If we won't be changing direction, go to the next vertex
722 //
723 if (nextZ*prevZ < 0) continue;
724 }
725
726
727 //
728 // We overlap in z with the side of the face that stretches from
729 // vertex "prev" to "corn". On which side (left or right) do
730 // we lay with respect to this segment?
731 //
732 G4ThreeVector qa( qx - prev->x, qy - prev->y, qz - prev->z ),
733 qb( qx - corn->x, qy - corn->y, qz - corn->z );
734
735 G4double aboveOrBelow = normSign*qa.cross(qb).dot(v);
736
737 if (aboveOrBelow > 0)
738 ++answer;
739 else if (aboveOrBelow < 0)
740 --answer;
741 else
742 {
743 //
744 // A precisely zero answer here means we exactly
745 // intersect (within roundoff) the edge of the face.
746 // Return true in this case.
747 //
748 return true;
749 }
750 } while( prevZ = cornZ, prev=corn, ++corn < corners+numEdges );
751
752 return answer!=0;
753 }
754
755 // InsideEdges (don't care about distance)
756 //
757 // Decide if the point in r,z is inside the edges of our face
758 //
759 // This routine can be made a zillion times quicker by implementing
760 // better code, for example:
761 //
762 // int pnpoly(int npol, float *xp, float *yp, float x, float y)
763 // {
764 // int i, j, c = 0;
765 // for (i = 0, j = npol-1; i < npol; j = i++) {
766 // if ((((yp[i]<=y) && (y<yp[j])) ||
767 // ((yp[j]<=y) && (y<yp[i]))) &&
768 // (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
769 //
770 // c = !c;
771 // }
772 // return c;
773 // }
774 //
775 // See "Point in Polyon Strategies", Eric Haines [Graphic Gems IV] pp. 24-46
776 //
777 // The algorithm below is rather unique, but is based on code needed to
778 // calculate the distance to the shape. Left here for testing...
779 //
780 G4bool G4PolyPhiFace::InsideEdges( G4double r, G4double z )
781 {
782 //
783 // Quick check of extent
784 //
785 if ( r < rMin || r > rMax ) return false;
786 if ( z < zMin || z > zMax ) return false;
787
788 //
789 // More thorough check
790 //
791 G4double notUsed;
792
793 return InsideEdges( r, z, ¬Used, nullptr );
794 }
795
796 // InsideEdges (care about distance)
797 //
798 // Decide if the point in r,z is inside the edges of our face
799 //
800 G4bool G4PolyPhiFace::InsideEdges( G4double r, G4double z,
801 G4double* bestDist2,
802 G4PolyPhiFaceVertex** base3Dnorm,
803 G4ThreeVector** head3Dnorm )
804 {
805 G4double bestDistance2 = kInfinity;
806 G4bool answer = false;
807
808 G4PolyPhiFaceEdge* edge = edges;
809 do // Loop checking, 13.08.2015, G.Cosmo
810 {
811 G4PolyPhiFaceVertex* testMe = nullptr;
812 G4PolyPhiFaceVertex* v0_vertex = edge->v0;
813 G4PolyPhiFaceVertex* v1_vertex = edge->v1;
814 //
815 // Get distance perpendicular to the edge
816 //
817 G4double dr = (r-v0_vertex->r), dz = (z-v0_vertex->z);
818
819 G4double distOut = dr*edge->tz - dz*edge->tr;
820 G4double distance2 = distOut*distOut;
821 if (distance2 > bestDistance2) continue; // No hope!
822
823 //
824 // Check to see if normal intersects edge within the edge's boundary
825 //
826 G4double q = dr*edge->tr + dz*edge->tz;
827
828 //
829 // If it doesn't, penalize distance2 appropriately
830 //
831 if (q < 0)
832 {
833 distance2 += q*q;
834 testMe = v0_vertex;
835 }
836 else if (q > edge->length)
837 {
838 G4double s2 = q-edge->length;
839 distance2 += s2*s2;
840 testMe = v1_vertex;
841 }
842
843 //
844 // Closest edge so far?
845 //
846 if (distance2 < bestDistance2)
847 {
848 bestDistance2 = distance2;
849 if (testMe != nullptr)
850 {
851 G4double distNorm = dr*testMe->rNorm + dz*testMe->zNorm;
852 answer = (distNorm <= 0);
853 if (base3Dnorm != nullptr)
854 {
855 *base3Dnorm = testMe;
856 *head3Dnorm = &testMe->norm3D;
857 }
858 }
859 else
860 {
861 answer = (distOut <= 0);
862 if (base3Dnorm != nullptr)
863 {
864 *base3Dnorm = v0_vertex;
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 //
878 G4double G4PolyPhiFace::SurfaceTriangle( const G4ThreeVector& p1,
879 const G4ThreeVector& p2,
880 const 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 //
896 G4double G4PolyPhiFace::SurfaceArea()
897 {
898 if ( fSurfaceArea==0. ) { Triangulate(); }
899 return fSurfaceArea;
900 }
901
902 // Return random point on face
903 //
904 G4ThreeVector G4PolyPhiFace::GetPointOnFace()
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 //
916 G4double G4PolyPhiFace::Area2( const G4TwoVector& a,
917 const G4TwoVector& b,
918 const 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 //
926 G4bool G4PolyPhiFace::Left( const G4TwoVector& a,
927 const G4TwoVector& b,
928 const G4TwoVector& c )
929 {
930 return Area2(a,b,c)>0;
931 }
932
933 // Boolean function for sign of Surface
934 //
935 G4bool G4PolyPhiFace::LeftOn( const G4TwoVector& a,
936 const G4TwoVector& b,
937 const G4TwoVector& c )
938 {
939 return Area2(a,b,c)>=0;
940 }
941
942 // Boolean function for sign of Surface
943 //
944 G4bool G4PolyPhiFace::Collinear( const G4TwoVector& a,
945 const G4TwoVector& b,
946 const 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 //
954 G4bool G4PolyPhiFace::IntersectProp( const G4TwoVector& a,
955 const G4TwoVector& b,
956 const G4TwoVector& c, const G4TwoVector& d )
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)) != 0);
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 //
969 G4bool G4PolyPhiFace::Between( const G4TwoVector& a, const G4TwoVector& b, const G4TwoVector& c )
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 //
988 G4bool G4PolyPhiFace::Intersect( const G4TwoVector& a,
989 const G4TwoVector& b,
990 const G4TwoVector& c, const G4TwoVector& d )
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 //
1006 G4bool G4PolyPhiFace::Diagonalie( G4PolyPhiFaceVertex* a,
1007 G4PolyPhiFaceVertex* b )
1008 {
1009 G4PolyPhiFaceVertex *corner = triangles;
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 //
1039 G4bool G4PolyPhiFace::InCone( G4PolyPhiFaceVertex* a, G4PolyPhiFaceVertex* b )
1040 {
1041 // a0,a and a1 are consecutive vertices
1042 //
1043 G4PolyPhiFaceVertex *a0,*a1;
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 //
1065 G4bool G4PolyPhiFace::Diagonal( G4PolyPhiFaceVertex* a, G4PolyPhiFaceVertex* b )
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 //
1073 void G4PolyPhiFace::EarInit()
1074 {
1075 G4PolyPhiFaceVertex* corner = triangles;
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 //
1096 void G4PolyPhiFace::Triangulate()
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
1101 const std::size_t maxEdges = (numEdges > 0) ? numEdges : 1;
1102 auto tri_help = new G4PolyPhiFaceVertex[maxEdges];
1103 triangles = tri_help;
1104 G4PolyPhiFaceVertex* triang = triangles;
1105
1106 std::vector<G4double> areas;
1107 std::vector<G4ThreeVector> points;
1108 G4double area=0.;
1109 G4PolyPhiFaceVertex *v0,*v1,*v2,*v3,*v4;
1110 v2=triangles;
1111
1112 // Make copy for prev/next for triang=corners
1113 //
1114 G4PolyPhiFaceVertex* helper = corners;
1115 G4PolyPhiFaceVertex* helper2 = corners;
1116 do // Loop checking, 13.08.2015, G.Cosmo
1117 {
1118 triang->r = helper->r;
1119 triang->z = helper->z;
1120 triang->x = helper->x;
1121 triang->y = helper->y;
1122
1123 // add pointer on prev corner
1124 //
1125 if( helper==corners )
1126 { triang->prev=triangles+maxEdges-1; }
1127 else
1128 { triang->prev=helper2; }
1129
1130 // add pointer on next corner
1131 //
1132 if( helper<corners+maxEdges-1 )
1133 { triang->next=triang+1; }
1134 else
1135 { triang->next=triangles; }
1136 helper2=triang;
1137 helper=helper->next;
1138 triang=triang->next;
1139
1140 } while( helper!=corners );
1141
1142 EarInit();
1143
1144 G4int n=numEdges;
1145 G4int i=0;
1146 G4ThreeVector p1,p2,p3,p4;
1147 const G4int max_n_loops=numEdges*10000; // protection against infinite loop
1148
1149 // Each step of outer loop removes one ear
1150 //
1151 while(n>3) // Loop checking, 13.08.2015, G.Cosmo
1152 { // Inner loop searches for one ear
1153 v2=triangles;
1154 do // Loop checking, 13.08.2015, G.Cosmo
1155 {
1156 if(v2->ear) // Ear found. Fill variables
1157 {
1158 // (v1,v3) is diagonal
1159 //
1160 v3=v2->next; v4=v3->next;
1161 v1=v2->prev; v0=v1->prev;
1162
1163 // Calculate areas and points
1164
1165 p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
1166 p2=G4ThreeVector((v1)->x,(v1)->y,(v1)->z);
1167 p3=G4ThreeVector((v3)->x,(v3)->y,(v3)->z);
1168
1169 G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
1170 points.push_back(p4);
1171 areas.push_back(result1);
1172 area=area+result1;
1173
1174 // Update earity of diagonal endpoints
1175 //
1176 v1->ear=Diagonal(v0,v3);
1177 v3->ear=Diagonal(v1,v4);
1178
1179 // Cut off the ear v2
1180 // Has to be done for a copy and not for real PolyPhiFace
1181 //
1182 v1->next=v3;
1183 v3->prev=v1;
1184 triangles=v3; // In case the head was v2
1185 --n;
1186
1187 break; // out of inner loop
1188 } // end if ear found
1189
1190 v2=v2->next;
1191
1192 } while( v2!=triangles );
1193
1194 ++i;
1195 if(i>=max_n_loops)
1196 {
1197 G4Exception( "G4PolyPhiFace::Triangulation()",
1198 "GeomSolids0003", FatalException,
1199 "Maximum number of steps is reached for triangulation!" );
1200 }
1201 } // end outer while loop
1202
1203 if(v2->next != nullptr)
1204 {
1205 // add last triangle
1206 //
1207 v2=v2->next;
1208 p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
1209 p2=G4ThreeVector((v2->next)->x,(v2->next)->y,(v2->next)->z);
1210 p3=G4ThreeVector((v2->prev)->x,(v2->prev)->y,(v2->prev)->z);
1211 G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
1212 points.push_back(p4);
1213 areas.push_back(result1);
1214 area=area+result1;
1215 }
1216
1217 // Surface Area is stored
1218 //
1219 fSurfaceArea = area;
1220
1221 // Second Step: choose randomly one surface
1222 //
1223 G4double chose = area*G4UniformRand();
1224
1225 // Third Step: Get a point on choosen surface
1226 //
1227 G4double Achose1, Achose2;
1228 Achose1=0; Achose2=0.;
1229 i=0;
1230 do // Loop checking, 13.08.2015, G.Cosmo
1231 {
1232 Achose2+=areas[i];
1233 if(chose>=Achose1 && chose<Achose2)
1234 {
1235 G4ThreeVector point;
1236 point=points[i];
1237 surface_point=point;
1238 break;
1239 }
1240 ++i; Achose1=Achose2;
1241 } while( i<numEdges-2 );
1242
1243 delete [] tri_help;
1244 tri_help = nullptr;
1245 }
1246