Geant4  10.02.p02
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 92024 2015-08-13 14:16:00Z 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 // Loop checking, 13.08.2015, G.Cosmo
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 // Loop checking, 13.08.2015, G.Cosmo
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 // Loop checking, 13.08.2015, G.Cosmo
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 // Loop checking, 13.08.2015, G.Cosmo
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 // Loop checking, 13.08.2015, G.Cosmo
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 // Loop checking, 13.08.2015, G.Cosmo
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 // Loop checking, 13.08.2015, G.Cosmo
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 // Loop checking, 13.08.2015, G.Cosmo
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 // Loop checking, 13.08.2015, G.Cosmo
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 // Loop checking, 13.08.2015, G.Cosmo
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  return answer!=0;
786 }
787 
788 
789 //
790 // InsideEdges (don't care aboud distance)
791 //
792 // Decide if the point in r,z is inside the edges of our face
793 //
794 // This routine can be made a zillion times quicker by implementing
795 // better code, for example:
796 //
797 // int pnpoly(int npol, float *xp, float *yp, float x, float y)
798 // {
799 // int i, j, c = 0;
800 // for (i = 0, j = npol-1; i < npol; j = i++) {
801 // if ((((yp[i]<=y) && (y<yp[j])) ||
802 // ((yp[j]<=y) && (y<yp[i]))) &&
803 // (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
804 //
805 // c = !c;
806 // }
807 // return c;
808 // }
809 //
810 // See "Point in Polyon Strategies", Eric Haines [Graphic Gems IV] pp. 24-46
811 //
812 // My algorithm below is rather unique, but is based on code needed to
813 // calculate the distance to the shape. I left it in here because ...
814 // well ... to test it better.
815 //
817 {
818  //
819  // Quick check of extent
820  //
821  if ( r < rMin || r > rMax ) return false;
822  if ( z < zMin || z > zMax ) return false;
823 
824  //
825  // More thorough check
826  //
827  G4double notUsed;
828 
829  return InsideEdges( r, z, &notUsed, 0 );
830 }
831 
832 
833 //
834 // InsideEdges (care about distance)
835 //
836 // Decide if the point in r,z is inside the edges of our face
837 //
839  G4double *bestDist2,
840  G4PolyPhiFaceVertex **base3Dnorm,
841  G4ThreeVector **head3Dnorm )
842 {
843  G4double bestDistance2 = kInfinity;
844  G4bool answer = 0;
845 
846  G4PolyPhiFaceEdge *edge = edges;
847  do // Loop checking, 13.08.2015, G.Cosmo
848  {
849  G4PolyPhiFaceVertex *testMe;
850  //
851  // Get distance perpendicular to the edge
852  //
853  G4double dr = (r-edge->v0->r), dz = (z-edge->v0->z);
854 
855  G4double distOut = dr*edge->tz - dz*edge->tr;
856  G4double distance2 = distOut*distOut;
857  if (distance2 > bestDistance2) continue; // No hope!
858 
859  //
860  // Check to see if normal intersects edge within the edge's boundary
861  //
862  G4double q = dr*edge->tr + dz*edge->tz;
863 
864  //
865  // If it doesn't, penalize distance2 appropriately
866  //
867  if (q < 0)
868  {
869  distance2 += q*q;
870  testMe = edge->v0;
871  }
872  else if (q > edge->length)
873  {
874  G4double s2 = q-edge->length;
875  distance2 += s2*s2;
876  testMe = edge->v1;
877  }
878  else
879  {
880  testMe = 0;
881  }
882 
883  //
884  // Closest edge so far?
885  //
886  if (distance2 < bestDistance2)
887  {
888  bestDistance2 = distance2;
889  if (testMe)
890  {
891  G4double distNorm = dr*testMe->rNorm + dz*testMe->zNorm;
892  answer = (distNorm <= 0);
893  if (base3Dnorm)
894  {
895  *base3Dnorm = testMe;
896  *head3Dnorm = &testMe->norm3D;
897  }
898  }
899  else
900  {
901  answer = (distOut <= 0);
902  if (base3Dnorm)
903  {
904  *base3Dnorm = edge->v0;
905  *head3Dnorm = &edge->norm3D;
906  }
907  }
908  }
909  } while( ++edge < edges + numEdges );
910 
911  *bestDist2 = bestDistance2;
912  return answer;
913 }
914 
915 //
916 // Calculation of Surface Area of a Triangle
917 // In the same time Random Point in Triangle is given
918 //
920  G4ThreeVector p2,
921  G4ThreeVector p3,
922  G4ThreeVector *p4 )
923 {
924  G4ThreeVector v, w;
925 
926  v = p3 - p1;
927  w = p1 - p2;
928  G4double lambda1 = G4UniformRand();
929  G4double lambda2 = lambda1*G4UniformRand();
930 
931  *p4=p2 + lambda1*w + lambda2*v;
932  return 0.5*(v.cross(w)).mag();
933 }
934 
935 //
936 // Compute surface area
937 //
939 {
940  if ( fSurfaceArea==0. ) { Triangulate(); }
941  return fSurfaceArea;
942 }
943 
944 //
945 // Return random point on face
946 //
948 {
949  Triangulate();
950  return surface_point;
951 }
952 
953 //
954 // Auxiliary Functions used for Finding the PointOnFace using Triangulation
955 //
956 
957 //
958 // Calculation of 2*Area of Triangle with Sign
959 //
961  G4TwoVector b,
962  G4TwoVector c )
963 {
964  return ((b.x()-a.x())*(c.y()-a.y())-
965  (c.x()-a.x())*(b.y()-a.y()));
966 }
967 
968 //
969 // Boolean function for sign of Surface
970 //
972  G4TwoVector b,
973  G4TwoVector c )
974 {
975  return Area2(a,b,c)>0;
976 }
977 
978 //
979 // Boolean function for sign of Surface
980 //
982  G4TwoVector b,
983  G4TwoVector c )
984 {
985  return Area2(a,b,c)>=0;
986 }
987 
988 //
989 // Boolean function for sign of Surface
990 //
992  G4TwoVector b,
993  G4TwoVector c )
994 {
995  return Area2(a,b,c)==0;
996 }
997 
998 //
999 // Boolean function for finding "Proper" Intersection
1000 // That means Intersection of two lines segments (a,b) and (c,d)
1001 //
1003  G4TwoVector b,
1004  G4TwoVector c, G4TwoVector d )
1005 {
1006  if( Collinear(a,b,c) || Collinear(a,b,d)||
1007  Collinear(c,d,a) || Collinear(c,d,b) ) { return false; }
1008 
1009  G4bool Positive;
1010  Positive = !(Left(a,b,c))^!(Left(a,b,d));
1011  return Positive && (!Left(c,d,a)^!Left(c,d,b));
1012 }
1013 
1014 //
1015 // Boolean function for determining if Point c is between a and b
1016 // For the tree points(a,b,c) on the same line
1017 //
1019 {
1020  if( !Collinear(a,b,c) ) { return false; }
1021 
1022  if(a.x()!=b.x())
1023  {
1024  return ((a.x()<=c.x())&&(c.x()<=b.x()))||
1025  ((a.x()>=c.x())&&(c.x()>=b.x()));
1026  }
1027  else
1028  {
1029  return ((a.y()<=c.y())&&(c.y()<=b.y()))||
1030  ((a.y()>=c.y())&&(c.y()>=b.y()));
1031  }
1032 }
1033 
1034 //
1035 // Boolean function for finding Intersection "Proper" or not
1036 // Between two line segments (a,b) and (c,d)
1037 //
1039  G4TwoVector b,
1040  G4TwoVector c, G4TwoVector d )
1041 {
1042  if( IntersectProp(a,b,c,d) )
1043  { return true; }
1044  else if( Between(a,b,c)||
1045  Between(a,b,d)||
1046  Between(c,d,a)||
1047  Between(c,d,b) )
1048  { return true; }
1049  else
1050  { return false; }
1051 }
1052 
1053 //
1054 // Boolean Diagonalie help to determine
1055 // if diagonal s of segment (a,b) is convex or reflex
1056 //
1058  G4PolyPhiFaceVertex *b )
1059 {
1060  G4PolyPhiFaceVertex *corner = triangles;
1061  G4PolyPhiFaceVertex *corner_next=triangles;
1062 
1063  // For each Edge (corner,corner_next)
1064  do // Loop checking, 13.08.2015, G.Cosmo
1065  {
1066  corner_next=corner->next;
1067 
1068  // Skip edges incident to a of b
1069  //
1070  if( (corner!=a)&&(corner_next!=a)
1071  &&(corner!=b)&&(corner_next!=b) )
1072  {
1073  G4TwoVector rz1,rz2,rz3,rz4;
1074  rz1 = G4TwoVector(a->r,a->z);
1075  rz2 = G4TwoVector(b->r,b->z);
1076  rz3 = G4TwoVector(corner->r,corner->z);
1077  rz4 = G4TwoVector(corner_next->r,corner_next->z);
1078  if( Intersect(rz1,rz2,rz3,rz4) ) { return false; }
1079  }
1080  corner=corner->next;
1081 
1082  } while( corner != triangles );
1083 
1084  return true;
1085 }
1086 
1087 //
1088 // Boolean function that determine if b is Inside Cone (a0,a,a1)
1089 // being a the center of the Cone
1090 //
1092 {
1093  // a0,a and a1 are consecutive vertices
1094  //
1096  a1=a->next;
1097  a0=a->prev;
1098 
1099  G4TwoVector arz,arz0,arz1,brz;
1100  arz=G4TwoVector(a->r,a->z);arz0=G4TwoVector(a0->r,a0->z);
1101  arz1=G4TwoVector(a1->r,a1->z);brz=G4TwoVector(b->r,b->z);
1102 
1103 
1104  if(LeftOn(arz,arz1,arz0)) // If a is convex vertex
1105  {
1106  return Left(arz,brz,arz0)&&Left(brz,arz,arz1);
1107  }
1108  else // Else a is reflex
1109  {
1110  return !( LeftOn(arz,brz,arz1)&&LeftOn(brz,arz,arz0));
1111  }
1112 }
1113 
1114 //
1115 // Boolean function finding if Diagonal is possible
1116 // inside Polycone or PolyHedra
1117 //
1119 {
1120  return InCone(a,b) && InCone(b,a) && Diagonalie(a,b);
1121 }
1122 
1123 //
1124 // Initialisation for Triangulisation by ear tips
1125 // For details see "Computational Geometry in C" by Joseph O'Rourke
1126 //
1128 {
1129  G4PolyPhiFaceVertex *corner = triangles;
1130  G4PolyPhiFaceVertex *c_prev,*c_next;
1131 
1132  do // Loop checking, 13.08.2015, G.Cosmo
1133  {
1134  // We need to determine three consecutive vertices
1135  //
1136  c_next=corner->next;
1137  c_prev=corner->prev;
1138 
1139  // Calculation of ears
1140  //
1141  corner->ear=Diagonal(c_prev,c_next);
1142  corner=corner->next;
1143 
1144  } while( corner!=triangles );
1145 }
1146 
1147 //
1148 // Triangulisation by ear tips for Polycone or Polyhedra
1149 // For details see "Computational Geometry in C" by Joseph O'Rourke
1150 //
1152 {
1153  // The copy of Polycone is made and this copy is reordered in order to
1154  // have a list of triangles. This list is used for GetPointOnFace().
1155 
1157  triangles = tri_help;
1158  G4PolyPhiFaceVertex *triang = triangles;
1159 
1160  std::vector<G4double> areas;
1161  std::vector<G4ThreeVector> points;
1162  G4double area=0.;
1163  G4PolyPhiFaceVertex *v0,*v1,*v2,*v3,*v4;
1164  v2=triangles;
1165 
1166  // Make copy for prev/next for triang=corners
1167  //
1168  G4PolyPhiFaceVertex *helper = corners;
1169  G4PolyPhiFaceVertex *helper2 = corners;
1170  do // Loop checking, 13.08.2015, G.Cosmo
1171  {
1172  triang->r = helper->r;
1173  triang->z = helper->z;
1174  triang->x = helper->x;
1175  triang->y= helper->y;
1176 
1177  // add pointer on prev corner
1178  //
1179  if( helper==corners )
1180  { triang->prev=triangles+numEdges-1; }
1181  else
1182  { triang->prev=helper2; }
1183 
1184  // add pointer on next corner
1185  //
1186  if( helper<corners+numEdges-1 )
1187  { triang->next=triang+1; }
1188  else
1189  { triang->next=triangles; }
1190  helper2=triang;
1191  helper=helper->next;
1192  triang=triang->next;
1193 
1194  } while( helper!=corners );
1195 
1196  EarInit();
1197 
1198  G4int n=numEdges;
1199  G4int i=0;
1200  G4ThreeVector p1,p2,p3,p4;
1201  const G4int max_n_loops=numEdges*10000; // protection against infinite loop
1202 
1203  // Each step of outer loop removes one ear
1204  //
1205  while(n>3) // Loop checking, 13.08.2015, G.Cosmo
1206  { // Inner loop searches for one ear
1207  v2=triangles;
1208  do // Loop checking, 13.08.2015, G.Cosmo
1209  {
1210  if(v2->ear) // Ear found. Fill variables
1211  {
1212  // (v1,v3) is diagonal
1213  //
1214  v3=v2->next; v4=v3->next;
1215  v1=v2->prev; v0=v1->prev;
1216 
1217  // Calculate areas and points
1218 
1219  p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
1220  p2=G4ThreeVector((v1)->x,(v1)->y,(v1)->z);
1221  p3=G4ThreeVector((v3)->x,(v3)->y,(v3)->z);
1222 
1223  G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
1224  points.push_back(p4);
1225  areas.push_back(result1);
1226  area=area+result1;
1227 
1228  // Update earity of diagonal endpoints
1229  //
1230  v1->ear=Diagonal(v0,v3);
1231  v3->ear=Diagonal(v1,v4);
1232 
1233  // Cut off the ear v2
1234  // Has to be done for a copy and not for real PolyPhiFace
1235  //
1236  v1->next=v3;
1237  v3->prev=v1;
1238  triangles=v3; // In case the head was v2
1239  n--;
1240 
1241  break; // out of inner loop
1242  } // end if ear found
1243 
1244  v2=v2->next;
1245 
1246  } while( v2!=triangles );
1247 
1248  i++;
1249  if(i>=max_n_loops)
1250  {
1251  G4Exception( "G4PolyPhiFace::Triangulation()",
1252  "GeomSolids0003", FatalException,
1253  "Maximum number of steps is reached for triangulation!" );
1254  }
1255  } // end outer while loop
1256 
1257  if(v2->next)
1258  {
1259  // add last triangle
1260  //
1261  v2=v2->next;
1262  p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
1263  p2=G4ThreeVector((v2->next)->x,(v2->next)->y,(v2->next)->z);
1264  p3=G4ThreeVector((v2->prev)->x,(v2->prev)->y,(v2->prev)->z);
1265  G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
1266  points.push_back(p4);
1267  areas.push_back(result1);
1268  area=area+result1;
1269  }
1270 
1271  // Surface Area is stored
1272  //
1273  fSurfaceArea = area;
1274 
1275  // Second Step: choose randomly one surface
1276  //
1277  G4double chose = area*G4UniformRand();
1278 
1279  // Third Step: Get a point on choosen surface
1280  //
1281  G4double Achose1, Achose2;
1282  Achose1=0; Achose2=0.;
1283  i=0;
1284  do // Loop checking, 13.08.2015, G.Cosmo
1285  {
1286  Achose2+=areas[i];
1287  if(chose>=Achose1 && chose<Achose2)
1288  {
1289  G4ThreeVector point;
1290  point=points[i] ;
1291  surface_point=point;
1292  break;
1293  }
1294  i++; Achose1=Achose2;
1295  } while( i<numEdges-2 );
1296 
1297  delete [] tri_help;
1298  tri_help = 0;
1299 }
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 G4float tolerance
static const G4double a1
G4PolyPhiFaceVertex * v0
G4double GetSurfaceTolerance() const
G4ThreeVector norm3D
const G4double w[NPOINTSGL]
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
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:97
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
const G4double x[NPOINTSGL]
#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