Geant4  10.00.p02
G4GenericPolycone.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: G4GenericPolycone.cc 72131 2013-07-11 07:12:24Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------------
31 // GEANT 4 class source file
32 //
33 //
34 // G4GenericPolycone.cc
35 //
36 // Implementation of a CSG polycone
37 //
38 // --------------------------------------------------------------------
39 
40 #include "G4GenericPolycone.hh"
41 
42 #if !defined(G4GEOM_USE_UGENERICPOLYCONE)
43 
44 #include "G4PolyconeSide.hh"
45 #include "G4PolyPhiFace.hh"
46 
47 #include "Randomize.hh"
48 
49 #include "G4Polyhedron.hh"
50 #include "G4EnclosingCylinder.hh"
51 #include "G4ReduciblePolygon.hh"
52 #include "G4VPVParameterisation.hh"
53 
54 using namespace CLHEP;
55 
56 //
57 // Constructor (generic parameters)
58 //
60  G4double phiStart,
61  G4double phiTotal,
62  G4int numRZ,
63  const G4double r[],
64  const G4double z[] )
65  : G4VCSGfaceted( name )
66 {
67 
68  G4ReduciblePolygon *rz = new G4ReduciblePolygon( r, z, numRZ );
69 
70  Create( phiStart, phiTotal, rz );
71 
72  // Set original_parameters struct for consistency
73  //
74  //SetOriginalParameters(rz);
75 
76  delete rz;
77 }
78 
79 //
80 // Create
81 //
82 // Generic create routine, called by each constructor after
83 // conversion of arguments
84 //
86  G4double phiTotal,
87  G4ReduciblePolygon *rz )
88 {
89  //
90  // Perform checks of rz values
91  //
92  if (rz->Amin() < 0.0)
93  {
94  std::ostringstream message;
95  message << "Illegal input parameters - " << GetName() << G4endl
96  << " All R values must be >= 0 !";
97  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
98  FatalErrorInArgument, message);
99  }
100 
101  G4double rzArea = rz->Area();
102  if (rzArea < -kCarTolerance)
103  rz->ReverseOrder();
104 
105  else if (rzArea < -kCarTolerance)
106  {
107  std::ostringstream message;
108  message << "Illegal input parameters - " << GetName() << G4endl
109  << " R/Z cross section is zero or near zero: " << rzArea;
110  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
111  FatalErrorInArgument, message);
112  }
113 
115  || (!rz->RemoveRedundantVertices( kCarTolerance )) )
116  {
117  std::ostringstream message;
118  message << "Illegal input parameters - " << GetName() << G4endl
119  << " Too few unique R/Z values !";
120  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
121  FatalErrorInArgument, message);
122  }
123 
124  if (rz->CrossesItself(1/kInfinity))
125  {
126  std::ostringstream message;
127  message << "Illegal input parameters - " << GetName() << G4endl
128  << " R/Z segments cross !";
129  G4Exception("G4GenericPolycone::Create()", "GeomSolids0002",
130  FatalErrorInArgument, message);
131  }
132 
133  numCorner = rz->NumVertices();
134 
135  //
136  // Phi opening? Account for some possible roundoff, and interpret
137  // nonsense value as representing no phi opening
138  //
139  if (phiTotal <= 0 || phiTotal > twopi-1E-10)
140  {
141  phiIsOpen = false;
142  startPhi = 0;
143  endPhi = twopi;
144  }
145  else
146  {
147  phiIsOpen = true;
148 
149  //
150  // Convert phi into our convention
151  //
152  startPhi = phiStart;
153  while( startPhi < 0 ) startPhi += twopi;
154 
155  endPhi = phiStart+phiTotal;
156  while( endPhi < startPhi ) endPhi += twopi;
157  }
158 
159  //
160  // Allocate corner array.
161  //
163 
164  //
165  // Copy corners
166  //
167  G4ReduciblePolygonIterator iterRZ(rz);
168 
169  G4PolyconeSideRZ *next = corners;
170  iterRZ.Begin();
171  do
172  {
173  next->r = iterRZ.GetA();
174  next->z = iterRZ.GetB();
175  } while( ++next, iterRZ.Next() );
176 
177  //
178  // Allocate face pointer array
179  //
181  faces = new G4VCSGface*[numFace];
182 
183  //
184  // Construct conical faces
185  //
186  // But! Don't construct a face if both points are at zero radius!
187  //
188  G4PolyconeSideRZ *corner = corners,
189  *prev = corners + numCorner-1,
190  *nextNext;
191  G4VCSGface **face = faces;
192  do
193  {
194  next = corner+1;
195  if (next >= corners+numCorner) next = corners;
196  nextNext = next+1;
197  if (nextNext >= corners+numCorner) nextNext = corners;
198 
199  if (corner->r < 1/kInfinity && next->r < 1/kInfinity) continue;
200 
201  //
202  // We must decide here if we can dare declare one of our faces
203  // as having a "valid" normal (i.e. allBehind = true). This
204  // is never possible if the face faces "inward" in r.
205  //
206  G4bool allBehind;
207  if (corner->z > next->z)
208  {
209  allBehind = false;
210  }
211  else
212  {
213  //
214  // Otherwise, it is only true if the line passing
215  // through the two points of the segment do not
216  // split the r/z cross section
217  //
218  allBehind = !rz->BisectedBy( corner->r, corner->z,
219  next->r, next->z, kCarTolerance );
220  }
221 
222  *face++ = new G4PolyconeSide( prev, corner, next, nextNext,
223  startPhi, endPhi-startPhi, phiIsOpen, allBehind );
224  } while( prev=corner, corner=next, corner > corners );
225 
226  if (phiIsOpen)
227  {
228  //
229  // Construct phi open edges
230  //
231  *face++ = new G4PolyPhiFace( rz, startPhi, 0, endPhi );
232  *face++ = new G4PolyPhiFace( rz, endPhi, 0, startPhi );
233  }
234 
235  //
236  // We might have dropped a face or two: recalculate numFace
237  //
238  numFace = face-faces;
239 
240  //
241  // Make enclosingCylinder
242  //
244  new G4EnclosingCylinder( rz, phiIsOpen, phiStart, phiTotal );
245 }
246 
247 
248 //
249 // Fake default constructor - sets only member data and allocates memory
250 // for usage restricted to object persistency.
251 //
253  : G4VCSGfaceted(a), startPhi(0.), endPhi(0.), phiIsOpen(false),
254  numCorner(0), corners(0), enclosingCylinder(0)
255 {
256 }
257 
258 
259 //
260 // Destructor
261 //
263 {
264  delete [] corners;
265  delete enclosingCylinder;
266 }
267 
268 
269 //
270 // Copy constructor
271 //
273  : G4VCSGfaceted( source )
274 {
275  CopyStuff( source );
276 }
277 
278 
279 //
280 // Assignment operator
281 //
283 {
284  if (this == &source) return *this;
285 
286  G4VCSGfaceted::operator=( source );
287 
288  delete [] corners;
289  // if (original_parameters) delete original_parameters;
290 
291  delete enclosingCylinder;
292 
293  CopyStuff( source );
294 
295  return *this;
296 }
297 
298 
299 //
300 // CopyStuff
301 //
303 {
304  //
305  // Simple stuff
306  //
307  startPhi = source.startPhi;
308  endPhi = source.endPhi;
309  phiIsOpen = source.phiIsOpen;
310  numCorner = source.numCorner;
311 
312  //
313  // The corner array
314  //
316 
317  G4PolyconeSideRZ *corn = corners,
318  *sourceCorn = source.corners;
319  do
320  {
321  *corn = *sourceCorn;
322  } while( ++sourceCorn, ++corn < corners+numCorner );
323 
324  //
325  // Enclosing cylinder
326  //
328 
330 }
331 
332 
333 //
334 // Reset
335 //
337 {
338 
339  std::ostringstream message;
340  message << "Solid " << GetName() << " built using generic construct."
341  << G4endl << "Not applicable to the generic construct !";
342  G4Exception("G4GenericPolycone::Reset()", "GeomSolids1001",
343  JustWarning, message, "Parameters NOT resetted.");
344  return 1;
345 
346 }
347 
348 
349 //
350 // Inside
351 //
352 // This is an override of G4VCSGfaceted::Inside, created in order
353 // to speed things up by first checking with G4EnclosingCylinder.
354 //
356 {
357  //
358  // Quick test
359  //
360  if (enclosingCylinder->MustBeOutside(p)) return kOutside;
361 
362  //
363  // Long answer
364  //
365  return G4VCSGfaceted::Inside(p);
366 }
367 
368 
369 //
370 // DistanceToIn
371 //
372 // This is an override of G4VCSGfaceted::Inside, created in order
373 // to speed things up by first checking with G4EnclosingCylinder.
374 //
376  const G4ThreeVector &v ) const
377 {
378  //
379  // Quick test
380  //
381  if (enclosingCylinder->ShouldMiss(p,v))
382  return kInfinity;
383 
384  //
385  // Long answer
386  //
387  return G4VCSGfaceted::DistanceToIn( p, v );
388 }
389 
390 
391 //
392 // DistanceToIn
393 //
395 {
396  return G4VCSGfaceted::DistanceToIn(p);
397 }
398 
399 
400 //
401 // ComputeDimensions
402 //
403 /*void G4GenericPolycone::ComputeDimensions( G4VPVParameterisation* p,
404  const G4int n,
405  const G4VPhysicalVolume* pRep )
406 {
407  p->ComputeDimensions(*this,n,pRep);
408 }
409 */
410 //
411 // GetEntityType
412 //
414 {
415  return G4String("G4GenericPolycone");
416 }
417 
418 
419 //
420 // Make a clone of the object
421 //
423 {
424  return new G4GenericPolycone(*this);
425 }
426 
427 //
428 // Stream object contents to an output stream
429 //
430 std::ostream& G4GenericPolycone::StreamInfo( std::ostream& os ) const
431 {
432  G4int oldprc = os.precision(16);
433  os << "-----------------------------------------------------------\n"
434  << " *** Dump for solid - " << GetName() << " ***\n"
435  << " ===================================================\n"
436  << " Solid type: G4GenericPolycone\n"
437  << " Parameters: \n"
438  << " starting phi angle : " << startPhi/degree << " degrees \n"
439  << " ending phi angle : " << endPhi/degree << " degrees \n";
440  G4int i=0;
441 
442  os << " number of RZ points: " << numCorner << "\n"
443  << " RZ values (corners): \n";
444  for (i=0; i<numCorner; i++)
445  {
446  os << " "
447  << corners[i].r << ", " << corners[i].z << "\n";
448  }
449  os << "-----------------------------------------------------------\n";
450  os.precision(oldprc);
451 
452  return os;
453 }
454 
455 
456 
457 //
458 // GetPointOnSurface
459 //
461 {
462  return GetPointOnSurfaceGeneric();
463 
464 }
465 
466 //
467 // CreatePolyhedron
468 //
470 {
471  // The following code prepares for:
472  // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
473  // const double xyz[][3],
474  // const int faces_vec[][4])
475  // Here is an extract from the header file HepPolyhedron.h:
493  const G4int numSide =
494  G4int(G4Polyhedron::GetNumberOfRotationSteps()
495  * (endPhi - startPhi) / twopi) + 1;
496  G4int nNodes;
497  G4int nFaces;
498  typedef G4double double3[3];
499  double3* xyz;
500  typedef G4int int4[4];
501  int4* faces_vec;
502  if (phiIsOpen)
503  {
504  // Triangulate open ends. Simple ear-chopping algorithm...
505  // I'm not sure how robust this algorithm is (J.Allison).
506  //
507  std::vector<G4bool> chopped(numCorner, false);
508  std::vector<G4int*> triQuads;
509  G4int remaining = numCorner;
510  G4int iStarter = 0;
511  while (remaining >= 3)
512  {
513  // Find unchopped corners...
514  //
515  G4int A = -1, B = -1, C = -1;
516  G4int iStepper = iStarter;
517  do
518  {
519  if (A < 0) { A = iStepper; }
520  else if (B < 0) { B = iStepper; }
521  else if (C < 0) { C = iStepper; }
522  do
523  {
524  if (++iStepper >= numCorner) { iStepper = 0; }
525  }
526  while (chopped[iStepper]);
527  }
528  while (C < 0 && iStepper != iStarter);
529 
530  // Check triangle at B is pointing outward (an "ear").
531  // Sign of z cross product determines...
532  //
533  G4double BAr = corners[A].r - corners[B].r;
534  G4double BAz = corners[A].z - corners[B].z;
535  G4double BCr = corners[C].r - corners[B].r;
536  G4double BCz = corners[C].z - corners[B].z;
537  if (BAr * BCz - BAz * BCr < kCarTolerance)
538  {
539  G4int* tq = new G4int[3];
540  tq[0] = A + 1;
541  tq[1] = B + 1;
542  tq[2] = C + 1;
543  triQuads.push_back(tq);
544  chopped[B] = true;
545  --remaining;
546  }
547  else
548  {
549  do
550  {
551  if (++iStarter >= numCorner) { iStarter = 0; }
552  }
553  while (chopped[iStarter]);
554  }
555  }
556  // Transfer to faces...
557  //
558  nNodes = (numSide + 1) * numCorner;
559  nFaces = numSide * numCorner + 2 * triQuads.size();
560  faces_vec = new int4[nFaces];
561  G4int iface = 0;
562  G4int addition = numCorner * numSide;
563  G4int d = numCorner - 1;
564  for (G4int iEnd = 0; iEnd < 2; ++iEnd)
565  {
566  for (size_t i = 0; i < triQuads.size(); ++i)
567  {
568  // Negative for soft/auxiliary/normally invisible edges...
569  //
570  G4int a, b, c;
571  if (iEnd == 0)
572  {
573  a = triQuads[i][0];
574  b = triQuads[i][1];
575  c = triQuads[i][2];
576  }
577  else
578  {
579  a = triQuads[i][0] + addition;
580  b = triQuads[i][2] + addition;
581  c = triQuads[i][1] + addition;
582  }
583  G4int ab = std::abs(b - a);
584  G4int bc = std::abs(c - b);
585  G4int ca = std::abs(a - c);
586  faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a;
587  faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b;
588  faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c;
589  faces_vec[iface][3] = 0;
590  ++iface;
591  }
592  }
593 
594  // Continue with sides...
595 
596  xyz = new double3[nNodes];
597  const G4double dPhi = (endPhi - startPhi) / numSide;
598  G4double phi = startPhi;
599  G4int ixyz = 0;
600  for (G4int iSide = 0; iSide < numSide; ++iSide)
601  {
602  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
603  {
604  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
605  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
606  xyz[ixyz][2] = corners[iCorner].z;
607  if (iSide == 0) // startPhi
608  {
609  if (iCorner < numCorner - 1)
610  {
611  faces_vec[iface][0] = ixyz + 1;
612  faces_vec[iface][1] = -(ixyz + numCorner + 1);
613  faces_vec[iface][2] = ixyz + numCorner + 2;
614  faces_vec[iface][3] = ixyz + 2;
615  }
616  else
617  {
618  faces_vec[iface][0] = ixyz + 1;
619  faces_vec[iface][1] = -(ixyz + numCorner + 1);
620  faces_vec[iface][2] = ixyz + 2;
621  faces_vec[iface][3] = ixyz - numCorner + 2;
622  }
623  }
624  else if (iSide == numSide - 1) // endPhi
625  {
626  if (iCorner < numCorner - 1)
627  {
628  faces_vec[iface][0] = ixyz + 1;
629  faces_vec[iface][1] = ixyz + numCorner + 1;
630  faces_vec[iface][2] = ixyz + numCorner + 2;
631  faces_vec[iface][3] = -(ixyz + 2);
632  }
633  else
634  {
635  faces_vec[iface][0] = ixyz + 1;
636  faces_vec[iface][1] = ixyz + numCorner + 1;
637  faces_vec[iface][2] = ixyz + 2;
638  faces_vec[iface][3] = -(ixyz - numCorner + 2);
639  }
640  }
641  else
642  {
643  if (iCorner < numCorner - 1)
644  {
645  faces_vec[iface][0] = ixyz + 1;
646  faces_vec[iface][1] = -(ixyz + numCorner + 1);
647  faces_vec[iface][2] = ixyz + numCorner + 2;
648  faces_vec[iface][3] = -(ixyz + 2);
649  }
650  else
651  {
652  faces_vec[iface][0] = ixyz + 1;
653  faces_vec[iface][1] = -(ixyz + numCorner + 1);
654  faces_vec[iface][2] = ixyz + 2;
655  faces_vec[iface][3] = -(ixyz - numCorner + 2);
656  }
657  }
658  ++iface;
659  ++ixyz;
660  }
661  phi += dPhi;
662  }
663 
664  // Last corners...
665 
666  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
667  {
668  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
669  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
670  xyz[ixyz][2] = corners[iCorner].z;
671  ++ixyz;
672  }
673  }
674  else // !phiIsOpen - i.e., a complete 360 degrees.
675  {
676  nNodes = numSide * numCorner;
677  nFaces = numSide * numCorner;;
678  xyz = new double3[nNodes];
679  faces_vec = new int4[nFaces];
680  const G4double dPhi = (endPhi - startPhi) / numSide;
681  G4double phi = startPhi;
682  G4int ixyz = 0, iface = 0;
683  for (G4int iSide = 0; iSide < numSide; ++iSide)
684  {
685  for (G4int iCorner = 0; iCorner < numCorner; ++iCorner)
686  {
687  xyz[ixyz][0] = corners[iCorner].r * std::cos(phi);
688  xyz[ixyz][1] = corners[iCorner].r * std::sin(phi);
689  xyz[ixyz][2] = corners[iCorner].z;
690 
691  if (iSide < numSide - 1)
692  {
693  if (iCorner < numCorner - 1)
694  {
695  faces_vec[iface][0] = ixyz + 1;
696  faces_vec[iface][1] = -(ixyz + numCorner + 1);
697  faces_vec[iface][2] = ixyz + numCorner + 2;
698  faces_vec[iface][3] = -(ixyz + 2);
699  }
700  else
701  {
702  faces_vec[iface][0] = ixyz + 1;
703  faces_vec[iface][1] = -(ixyz + numCorner + 1);
704  faces_vec[iface][2] = ixyz + 2;
705  faces_vec[iface][3] = -(ixyz - numCorner + 2);
706  }
707  }
708  else // Last side joins ends...
709  {
710  if (iCorner < numCorner - 1)
711  {
712  faces_vec[iface][0] = ixyz + 1;
713  faces_vec[iface][1] = -(ixyz + numCorner - nFaces + 1);
714  faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
715  faces_vec[iface][3] = -(ixyz + 2);
716  }
717  else
718  {
719  faces_vec[iface][0] = ixyz + 1;
720  faces_vec[iface][1] = -(ixyz - nFaces + numCorner + 1);
721  faces_vec[iface][2] = ixyz - nFaces + 2;
722  faces_vec[iface][3] = -(ixyz - numCorner + 2);
723  }
724  }
725  ++ixyz;
726  ++iface;
727  }
728  phi += dPhi;
729  }
730  }
731  G4Polyhedron* polyhedron = new G4Polyhedron;
732  G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec);
733  delete [] faces_vec;
734  delete [] xyz;
735  if (problem)
736  {
737  std::ostringstream message;
738  message << "Problem creating G4Polyhedron for: " << GetName();
739  G4Exception("G4GenericPolycone::CreatePolyhedron()", "GeomSolids1002",
740  JustWarning, message);
741  delete polyhedron;
742  return 0;
743  }
744  else
745  {
746  return polyhedron;
747  }
748 }
749 
750 #endif
G4String GetName() const
G4bool CrossesItself(G4double tolerance)
G4VCSGfaceted & operator=(const G4VCSGfaceted &source)
G4Polyhedron * fpPolyhedron
static const G4double kInfinity
Definition: geomdefs.hh:42
G4PolyconeSideRZ * corners
CLHEP::Hep3Vector G4ThreeVector
G4double Amin() const
G4double z
Definition: TRTMaterials.hh:39
G4String name
Definition: TRTMaterials.hh:40
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double a
Definition: TRTMaterials.hh:39
G4bool MustBeOutside(const G4ThreeVector &p) const
int G4int
Definition: G4Types.hh:78
G4GenericPolycone(const G4String &name, G4double phiStart, G4double phiTotal, G4int numRZ, const G4double r[], const G4double z[])
G4GeometryType GetEntityType() const
G4Polyhedron * CreatePolyhedron() const
G4bool RemoveDuplicateVertices(G4double tolerance)
G4ThreeVector GetPointOnSurface() const
G4EnclosingCylinder * enclosingCylinder
void Create(G4double phiStart, G4double phiTotal, G4ReduciblePolygon *rz)
G4bool RemoveRedundantVertices(G4double tolerance)
G4ThreeVector GetPointOnSurfaceGeneric() const
G4int NumVertices() const
void CopyStuff(const G4GenericPolycone &source)
bool G4bool
Definition: G4Types.hh:79
std::ostream & StreamInfo(std::ostream &os) const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
static const G4double A[nN]
G4VCSGface ** faces
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4GenericPolycone & operator=(const G4GenericPolycone &source)
static const G4double ab
EInside Inside(const G4ThreeVector &p) const
EInside
Definition: geomdefs.hh:58
static const double degree
Definition: G4SIunits.hh:125
G4bool BisectedBy(G4double a1, G4double b1, G4double a2, G4double b2, G4double tolerance)
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:305
double G4double
Definition: G4Types.hh:76
G4VSolid * Clone() const
virtual EInside Inside(const G4ThreeVector &p) const
G4bool ShouldMiss(const G4ThreeVector &p, const G4ThreeVector &v) const
virtual G4Polyhedron * GetPolyhedron() const