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