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