Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/management/src/G4VSolid.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // G4VSolid implementation for solid base class
 27 //
 28 // 10.10.18 E.Tcherniaev, more robust EstimateSurfaceArea() based on distance
 29 // 30.06.95 P.Kent, Created.
 30 // --------------------------------------------------------------------
 31 
 32 #include "G4VSolid.hh"
 33 #include "G4SolidStore.hh"
 34 #include "globals.hh"
 35 #include "G4QuickRand.hh"
 36 #include "G4GeometryTolerance.hh"
 37 
 38 #include "G4VoxelLimits.hh"
 39 #include "G4AffineTransform.hh"
 40 #include "G4VisExtent.hh"
 41 
 42 //////////////////////////////////////////////////////////////////////////
 43 //
 44 // Streaming operator dumping solid contents
 45 
 46 std::ostream& operator<< ( std::ostream& os, const G4VSolid& e )
 47 {
 48     return e.StreamInfo(os);
 49 }
 50 
 51 //////////////////////////////////////////////////////////////////////////
 52 //
 53 // Constructor
 54 //  - Copies name
 55 //  - Add ourselves to solid Store
 56 
 57 G4VSolid::G4VSolid(const G4String& name)
 58   : fshapeName(name)
 59 {
 60     kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
 61 
 62     // Register to store
 63     //
 64     G4SolidStore::GetInstance()->Register(this);
 65 }
 66 
 67 //////////////////////////////////////////////////////////////////////////
 68 //
 69 // Copy constructor
 70 //
 71 
 72 G4VSolid::G4VSolid(const G4VSolid& rhs)
 73   : kCarTolerance(rhs.kCarTolerance), fshapeName(rhs.fshapeName)
 74 {
 75     // Register to store
 76     //
 77     G4SolidStore::GetInstance()->Register(this);
 78 }
 79 
 80 //////////////////////////////////////////////////////////////////////////
 81 //
 82 // Fake default constructor - sets only member data and allocates memory
 83 //                            for usage restricted to object persistency.
 84 //
 85 G4VSolid::G4VSolid( __void__& )
 86   : fshapeName("")
 87 {
 88     // Register to store
 89     //
 90     G4SolidStore::GetInstance()->Register(this);
 91 }
 92 
 93 //////////////////////////////////////////////////////////////////////////
 94 //
 95 // Destructor (virtual)
 96 // - Remove ourselves from solid Store
 97 
 98 G4VSolid::~G4VSolid()
 99 {
100     G4SolidStore::GetInstance()->DeRegister(this);
101 }
102 
103 //////////////////////////////////////////////////////////////////////////
104 //
105 // Assignment operator
106 
107 G4VSolid& G4VSolid::operator = (const G4VSolid& rhs)
108 {
109    // Check assignment to self
110    //
111    if (this == &rhs)  { return *this; }
112 
113    // Copy data
114    //
115    kCarTolerance = rhs.kCarTolerance;
116    fshapeName = rhs.fshapeName;
117 
118    return *this;
119 }
120 
121 
122 
123 //////////////////////////////////////////////////////////////////////////
124 //
125 // Set solid name and notify store of the change
126 
127 void G4VSolid::SetName(const G4String& name)
128 {
129   fshapeName = name;
130   G4SolidStore::GetInstance()->SetMapValid(false);
131 }
132 
133 //////////////////////////////////////////////////////////////////////////
134 //
135 // Throw exception if ComputeDimensions called for illegal derived class
136 
137 void G4VSolid::ComputeDimensions(G4VPVParameterisation*,
138                                  const G4int,
139                                  const G4VPhysicalVolume*)
140 {
141     std::ostringstream message;
142     message << "Illegal call to G4VSolid::ComputeDimensions()" << G4endl
143             << "Method not overloaded by derived class !";
144     G4Exception("G4VSolid::ComputeDimensions()", "GeomMgt0003",
145                 FatalException, message);
146 }
147 
148 //////////////////////////////////////////////////////////////////////////
149 //
150 // Throw exception (warning) for solids not implementing the method
151 
152 G4ThreeVector G4VSolid::GetPointOnSurface() const
153 {
154     std::ostringstream message;
155     message << "Not implemented for solid: "
156             << GetEntityType() << " !" << G4endl
157             << "Returning origin.";
158     G4Exception("G4VSolid::GetPointOnSurface()", "GeomMgt1001",
159                 JustWarning, message);
160     return {0,0,0};
161 }
162 
163 //////////////////////////////////////////////////////////////////////////
164 //
165 // Returns total number of constituents that was used for construction
166 // of the solid. For non-Boolean solids the return value is one.
167 
168 G4int G4VSolid::GetNumOfConstituents() const
169 { return 1; }
170 
171 //////////////////////////////////////////////////////////////////////////
172 //
173 // Returns true if the solid has only planar faces, false otherwise.
174 
175 G4bool G4VSolid::IsFaceted() const
176 { return false; }
177 
178 //////////////////////////////////////////////////////////////////////////
179 //
180 // Dummy implementations ...
181 
182 const G4VSolid* G4VSolid::GetConstituentSolid(G4int) const
183 { return nullptr; }
184 
185 G4VSolid* G4VSolid::GetConstituentSolid(G4int)
186 { return nullptr; }
187 
188 const G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() const
189 { return nullptr; }
190 
191 G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr()
192 { return nullptr; }
193 
194 ////////////////////////////////////////////////////////////////
195 //
196 // Returns an estimation of the solid volume in internal units.
197 // The number of statistics and error accuracy is fixed.
198 // This method may be overloaded by derived classes to compute the
199 // exact geometrical quantity for solids where this is possible.
200 // or anyway to cache the computed value.
201 // This implementation does NOT cache the computed value.
202 
203 G4double G4VSolid::GetCubicVolume()
204 {
205   G4int cubVolStatistics = 1000000;
206   G4double cubVolEpsilon = 0.001;
207   return EstimateCubicVolume(cubVolStatistics, cubVolEpsilon);
208 }
209 
210 ////////////////////////////////////////////////////////////////
211 //
212 // Calculate cubic volume based on Inside() method.
213 // Accuracy is limited by the second argument or the statistics
214 // expressed by the first argument.
215 // Implementation is courtesy of Vasiliki Despoina Mitsou,
216 // University of Athens.
217 
218 G4double G4VSolid::EstimateCubicVolume(G4int nStat, G4double epsilon) const
219 {
220   G4int iInside=0;
221   G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,volume,halfepsilon;
222   G4ThreeVector p;
223   EInside in;
224 
225   // values needed for CalculateExtent signature
226 
227   G4VoxelLimits limit;                // Unlimited
228   G4AffineTransform origin;
229 
230   // min max extents of pSolid along X,Y,Z
231 
232   CalculateExtent(kXAxis,limit,origin,minX,maxX);
233   CalculateExtent(kYAxis,limit,origin,minY,maxY);
234   CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
235 
236   // limits
237 
238   if(nStat < 100)    nStat   = 100;
239   if(epsilon > 0.01) epsilon = 0.01;
240   halfepsilon = 0.5*epsilon;
241 
242   for(auto i = 0; i < nStat; ++i )
243   {
244     px = minX-halfepsilon+(maxX-minX+epsilon)*G4QuickRand();
245     py = minY-halfepsilon+(maxY-minY+epsilon)*G4QuickRand();
246     pz = minZ-halfepsilon+(maxZ-minZ+epsilon)*G4QuickRand();
247     p  = G4ThreeVector(px,py,pz);
248     in = Inside(p);
249     if(in != kOutside) ++iInside;
250   }
251   volume = (maxX-minX+epsilon)*(maxY-minY+epsilon)
252          * (maxZ-minZ+epsilon)*iInside/nStat;
253   return volume;
254 }
255 
256 ////////////////////////////////////////////////////////////////
257 //
258 // Returns an estimation of the solid surface area in internal units.
259 // The number of statistics and error accuracy is fixed.
260 // This method may be overloaded by derived classes to compute the
261 // exact geometrical quantity for solids where this is possible.
262 // or anyway to cache the computed value.
263 // This implementation does NOT cache the computed value.
264 
265 G4double G4VSolid::GetSurfaceArea()
266 {
267   G4int stat = 1000000;
268   G4double ell = -1.;
269   return EstimateSurfaceArea(stat,ell);
270 }
271 
272 //////////////////////////////////////////////////////////////////////////
273 //
274 // Calculate surface area by estimating volume of a thin shell
275 // surrounding the surface using Monte-Carlo method.
276 // Input parameters:
277 //    nstat - statistics (number of random points)
278 //    eps   - shell thinkness
279 
280 G4double G4VSolid::EstimateSurfaceArea(G4int nstat, G4double ell) const
281 {
282   static const G4double s2 = 1./std::sqrt(2.);
283   static const G4double s3 = 1./std::sqrt(3.);
284   static const G4ThreeVector directions[64] =
285   {
286     G4ThreeVector(  0,  0,  0), G4ThreeVector( -1,  0,  0), // (  ,  ,  ) ( -,  ,  )
287     G4ThreeVector(  1,  0,  0), G4ThreeVector( -1,  0,  0), // ( +,  ,  ) (-+,  ,  )
288     G4ThreeVector(  0, -1,  0), G4ThreeVector(-s2,-s2,  0), // (  , -,  ) ( -, -,  )
289     G4ThreeVector( s2, -s2, 0), G4ThreeVector(  0, -1,  0), // ( +, -,  ) (-+, -,  )
290 
291     G4ThreeVector(  0,  1,  0), G4ThreeVector( -s2, s2, 0), // (  , +,  ) ( -, +,  )
292     G4ThreeVector( s2, s2,  0), G4ThreeVector(  0,  1,  0), // ( +, +,  ) (-+, +,  )
293     G4ThreeVector(  0, -1,  0), G4ThreeVector( -1,  0,  0), // (  ,-+,  ) ( -,-+,  )
294     G4ThreeVector(  1,  0,  0), G4ThreeVector( -1,  0,  0), // ( +,-+,  ) (-+,-+,  )
295 
296     G4ThreeVector(  0,  0, -1), G4ThreeVector(-s2,  0,-s2), // (  ,  , -) ( -,  , -)
297     G4ThreeVector( s2,  0,-s2), G4ThreeVector(  0,  0, -1), // ( +,  , -) (-+,  , -)
298     G4ThreeVector(  0,-s2,-s2), G4ThreeVector(-s3,-s3,-s3), // (  , -, -) ( -, -, -)
299     G4ThreeVector( s3,-s3,-s3), G4ThreeVector(  0,-s2,-s2), // ( +, -, -) (-+, -, -)
300 
301     G4ThreeVector(  0, s2,-s2), G4ThreeVector(-s3, s3,-s3), // (  , +, -) ( -, +, -)
302     G4ThreeVector( s3, s3,-s3), G4ThreeVector(  0, s2,-s2), // ( +, +, -) (-+, +, -)
303     G4ThreeVector(  0,  0, -1), G4ThreeVector(-s2,  0,-s2), // (  ,-+, -) ( -,-+, -)
304     G4ThreeVector( s2,  0,-s2), G4ThreeVector(  0,  0, -1), // ( +,-+, -) (-+,-+, -)
305 
306     G4ThreeVector(  0,  0,  1), G4ThreeVector(-s2,  0, s2), // (  ,  , +) ( -,  , +)
307     G4ThreeVector( s2,  0, s2), G4ThreeVector(  0,  0,  1), // ( +,  , +) (-+,  , +)
308     G4ThreeVector(  0,-s2, s2), G4ThreeVector(-s3,-s3, s3), // (  , -, +) ( -, -, +)
309     G4ThreeVector( s3,-s3, s3), G4ThreeVector(  0,-s2, s2), // ( +, -, +) (-+, -, +)
310 
311     G4ThreeVector(  0, s2, s2), G4ThreeVector(-s3, s3, s3), // (  , +, +) ( -, +, +)
312     G4ThreeVector( s3, s3, s3), G4ThreeVector(  0, s2, s2), // ( +, +, +) (-+, +, +)
313     G4ThreeVector(  0,  0,  1), G4ThreeVector(-s2,  0, s2), // (  ,-+, +) ( -,-+, +)
314     G4ThreeVector( s2,  0, s2), G4ThreeVector(  0,  0,  1), // ( +,-+, +) (-+,-+, +)
315 
316     G4ThreeVector(  0,  0, -1), G4ThreeVector( -1,  0,  0), // (  ,  ,-+) ( -,  ,-+)
317     G4ThreeVector(  1,  0,  0), G4ThreeVector( -1,  0,  0), // ( +,  ,-+) (-+,  ,-+)
318     G4ThreeVector(  0, -1,  0), G4ThreeVector(-s2,-s2,  0), // (  , -,-+) ( -, -,-+)
319     G4ThreeVector( s2, -s2, 0), G4ThreeVector(  0, -1,  0), // ( +, -,-+) (-+, -,-+)
320 
321     G4ThreeVector(  0,  1,  0), G4ThreeVector( -s2, s2, 0), // (  , +,-+) ( -, +,-+)
322     G4ThreeVector( s2, s2,  0), G4ThreeVector(  0,  1,  0), // ( +, +,-+) (-+, +,-+)
323     G4ThreeVector(  0, -1,  0), G4ThreeVector( -1,  0,  0), // (  ,-+,-+) ( -,-+,-+)
324     G4ThreeVector(  1,  0,  0), G4ThreeVector( -1,  0,  0), // ( +,-+,-+) (-+,-+,-+)
325   };
326 
327   G4ThreeVector bmin, bmax;
328   BoundingLimits(bmin, bmax);
329 
330   G4double dX = bmax.x() - bmin.x();
331   G4double dY = bmax.y() - bmin.y();
332   G4double dZ = bmax.z() - bmin.z();
333 
334   // Define statistics and shell thickness
335   //
336   G4int npoints = (nstat < 1000) ? 1000 : nstat;
337   G4double coeff = 0.5 / std::cbrt(G4double(npoints));
338   G4double eps = (ell > 0) ? ell : coeff * std::min(std::min(dX, dY), dZ);
339   G4double del = 1.8 * eps; // shold be more than sqrt(3.)
340 
341   G4double minX = bmin.x() - eps;
342   G4double minY = bmin.y() - eps;
343   G4double minZ = bmin.z() - eps;
344 
345   G4double dd = 2. * eps;
346   dX += dd;
347   dY += dd;
348   dZ += dd;
349 
350   // Calculate surface area
351   //
352   G4int icount = 0;
353   for(auto i = 0; i < npoints; ++i)
354   {
355     G4double px = minX + dX*G4QuickRand();
356     G4double py = minY + dY*G4QuickRand();
357     G4double pz = minZ + dZ*G4QuickRand();
358     G4ThreeVector p  = G4ThreeVector(px, py, pz);
359     EInside in = Inside(p);
360     G4double dist = 0;
361     if (in == kInside)
362     {
363       if (DistanceToOut(p) >= eps) continue;
364       G4int icase = 0;
365       if (Inside(G4ThreeVector(px-del, py, pz)) != kInside) icase += 1;
366       if (Inside(G4ThreeVector(px+del, py, pz)) != kInside) icase += 2;
367       if (Inside(G4ThreeVector(px, py-del, pz)) != kInside) icase += 4;
368       if (Inside(G4ThreeVector(px, py+del, pz)) != kInside) icase += 8;
369       if (Inside(G4ThreeVector(px, py, pz-del)) != kInside) icase += 16;
370       if (Inside(G4ThreeVector(px, py, pz+del)) != kInside) icase += 32;
371       if (icase == 0) continue;
372       G4ThreeVector v = directions[icase];
373       dist = DistanceToOut(p, v);
374       G4ThreeVector n = SurfaceNormal(p + v*dist);
375       dist *= v.dot(n);
376     }
377     else if (in == kOutside)
378     {
379       if (DistanceToIn(p) >= eps) continue;
380       G4int icase = 0;
381       if (Inside(G4ThreeVector(px-del, py, pz)) != kOutside) icase += 1;
382       if (Inside(G4ThreeVector(px+del, py, pz)) != kOutside) icase += 2;
383       if (Inside(G4ThreeVector(px, py-del, pz)) != kOutside) icase += 4;
384       if (Inside(G4ThreeVector(px, py+del, pz)) != kOutside) icase += 8;
385       if (Inside(G4ThreeVector(px, py, pz-del)) != kOutside) icase += 16;
386       if (Inside(G4ThreeVector(px, py, pz+del)) != kOutside) icase += 32;
387       if (icase == 0) continue;
388       G4ThreeVector v = directions[icase];
389       dist = DistanceToIn(p, v);
390       if (dist == kInfinity) continue;
391       G4ThreeVector n = SurfaceNormal(p + v*dist);
392       dist *= -(v.dot(n));
393     }
394     if (dist < eps) ++icount;
395   }
396   return dX*dY*dZ*icount/npoints/dd;
397 }
398 
399 ///////////////////////////////////////////////////////////////////////////
400 //
401 // Returns a pointer of a dynamically allocated copy of the solid.
402 // Returns NULL pointer with warning in case the concrete solid does not
403 // implement this method. The caller has responsibility for ownership.
404 //
405 
406 G4VSolid* G4VSolid::Clone() const
407 {
408   std::ostringstream message;
409   message << "Clone() method not implemented for type: "
410           << GetEntityType() << "!" << G4endl
411           << "Returning NULL pointer!";
412   G4Exception("G4VSolid::Clone()", "GeomMgt1001", JustWarning, message);
413   return nullptr;
414 }
415 
416 ///////////////////////////////////////////////////////////////////////////
417 //
418 // Calculate the maximum and minimum extents of the polygon described
419 // by the vertices: pSectionIndex->pSectionIndex+1->
420 //                   pSectionIndex+2->pSectionIndex+3->pSectionIndex
421 // in the List pVertices
422 //
423 // If the minimum is <pMin pMin is set to the new minimum
424 // If the maximum is >pMax pMax is set to the new maximum
425 //
426 // No modifications are made to pVertices
427 //
428 
429 void G4VSolid::ClipCrossSection(       G4ThreeVectorList* pVertices,
430                                  const G4int pSectionIndex,
431                                  const G4VoxelLimits& pVoxelLimit,
432                                  const EAxis pAxis,
433                                        G4double& pMin, G4double& pMax) const
434 {
435 
436   G4ThreeVectorList polygon;
437   polygon.reserve(4);
438   polygon.push_back((*pVertices)[pSectionIndex]);
439   polygon.push_back((*pVertices)[pSectionIndex+1]);
440   polygon.push_back((*pVertices)[pSectionIndex+2]);
441   polygon.push_back((*pVertices)[pSectionIndex+3]);
442   CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
443   return;
444 }
445 
446 //////////////////////////////////////////////////////////////////////////////////
447 //
448 // Calculate the maximum and minimum extents of the polygons
449 // joining the CrossSections at pSectionIndex->pSectionIndex+3 and
450 //                              pSectionIndex+4->pSectionIndex7
451 //
452 // in the List pVertices, within the boundaries of the voxel limits pVoxelLimit
453 //
454 // If the minimum is <pMin pMin is set to the new minimum
455 // If the maximum is >pMax pMax is set to the new maximum
456 //
457 // No modifications are made to pVertices
458 
459 void G4VSolid::ClipBetweenSections(      G4ThreeVectorList* pVertices,
460                                    const G4int pSectionIndex,
461                                    const G4VoxelLimits& pVoxelLimit,
462                                    const EAxis pAxis,
463                                          G4double& pMin, G4double& pMax) const
464 {
465   G4ThreeVectorList polygon;
466   polygon.reserve(4);
467   polygon.push_back((*pVertices)[pSectionIndex]);
468   polygon.push_back((*pVertices)[pSectionIndex+4]);
469   polygon.push_back((*pVertices)[pSectionIndex+5]);
470   polygon.push_back((*pVertices)[pSectionIndex+1]);
471   CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
472   polygon.clear();
473 
474   polygon.push_back((*pVertices)[pSectionIndex+1]);
475   polygon.push_back((*pVertices)[pSectionIndex+5]);
476   polygon.push_back((*pVertices)[pSectionIndex+6]);
477   polygon.push_back((*pVertices)[pSectionIndex+2]);
478   CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
479   polygon.clear();
480 
481   polygon.push_back((*pVertices)[pSectionIndex+2]);
482   polygon.push_back((*pVertices)[pSectionIndex+6]);
483   polygon.push_back((*pVertices)[pSectionIndex+7]);
484   polygon.push_back((*pVertices)[pSectionIndex+3]);
485   CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
486   polygon.clear();
487 
488   polygon.push_back((*pVertices)[pSectionIndex+3]);
489   polygon.push_back((*pVertices)[pSectionIndex+7]);
490   polygon.push_back((*pVertices)[pSectionIndex+4]);
491   polygon.push_back((*pVertices)[pSectionIndex]);
492   CalculateClippedPolygonExtent(polygon,pVoxelLimit,pAxis,pMin,pMax);
493   return;
494 }
495 
496 
497 ///////////////////////////////////////////////////////////////////////////////
498 //
499 // Calculate the maximum and minimum extents of the convex polygon pPolygon
500 // along the axis pAxis, within the limits pVoxelLimit
501 //
502 
503 void
504 G4VSolid::CalculateClippedPolygonExtent(G4ThreeVectorList& pPolygon,
505                                   const G4VoxelLimits& pVoxelLimit,
506                                   const EAxis pAxis,
507                                         G4double& pMin,
508                                         G4double& pMax) const
509 {
510   G4int noLeft,i;
511   G4double component;
512 
513   ClipPolygon(pPolygon,pVoxelLimit,pAxis);
514   noLeft = (G4int)pPolygon.size();
515 
516   if ( noLeft != 0 )
517   {
518     for (i=0; i<noLeft; ++i)
519     {
520       component = pPolygon[i].operator()(pAxis);
521 
522       if (component < pMin)
523       {
524         pMin = component;
525       }
526       if (component > pMax)
527       {
528         pMax = component;
529       }
530     }
531   }
532 }
533 
534 /////////////////////////////////////////////////////////////////////////////
535 //
536 // Clip the convex polygon described by the vertices at
537 // pSectionIndex ->pSectionIndex+3 within pVertices to the limits pVoxelLimit
538 //
539 // Set pMin to the smallest
540 //
541 // Calculate the extent of the polygon along pAxis, when clipped to the
542 // limits pVoxelLimit. If the polygon exists after clippin, set pMin to
543 // the polygon's minimum extent along the axis if <pMin, and set pMax to
544 // the polygon's maximum extent along the axis if >pMax.
545 //
546 // The polygon is described by a set of vectors, where each vector represents
547 // a vertex, so that the polygon is described by the vertex sequence:
548 //   0th->1st 1st->2nd 2nd->... nth->0th
549 //
550 // Modifications to the polygon are made
551 //
552 // NOTE: Execessive copying during clipping
553 
554 void G4VSolid::ClipPolygon(      G4ThreeVectorList& pPolygon,
555                            const G4VoxelLimits& pVoxelLimit,
556                            const EAxis                        ) const
557 {
558   G4ThreeVectorList outputPolygon;
559 
560   if ( pVoxelLimit.IsLimited() )
561   {
562     if (pVoxelLimit.IsXLimited() ) // && pAxis != kXAxis)
563     {
564       G4VoxelLimits simpleLimit1;
565       simpleLimit1.AddLimit(kXAxis,pVoxelLimit.GetMinXExtent(),kInfinity);
566       ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
567 
568       pPolygon.clear();
569 
570       if ( outputPolygon.empty() )  return;
571 
572       G4VoxelLimits simpleLimit2;
573       simpleLimit2.AddLimit(kXAxis,-kInfinity,pVoxelLimit.GetMaxXExtent());
574       ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
575 
576       if ( pPolygon.empty() )       return;
577       else                          outputPolygon.clear();
578     }
579     if ( pVoxelLimit.IsYLimited() ) // && pAxis != kYAxis)
580     {
581       G4VoxelLimits simpleLimit1;
582       simpleLimit1.AddLimit(kYAxis,pVoxelLimit.GetMinYExtent(),kInfinity);
583       ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
584 
585       // Must always clear pPolygon - for clip to simpleLimit2 and in case of
586       // early exit
587 
588       pPolygon.clear();
589 
590       if ( outputPolygon.empty() )  return;
591 
592       G4VoxelLimits simpleLimit2;
593       simpleLimit2.AddLimit(kYAxis,-kInfinity,pVoxelLimit.GetMaxYExtent());
594       ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
595 
596       if ( pPolygon.empty() )       return;
597       else                          outputPolygon.clear();
598     }
599     if ( pVoxelLimit.IsZLimited() ) // && pAxis != kZAxis)
600     {
601       G4VoxelLimits simpleLimit1;
602       simpleLimit1.AddLimit(kZAxis,pVoxelLimit.GetMinZExtent(),kInfinity);
603       ClipPolygonToSimpleLimits(pPolygon,outputPolygon,simpleLimit1);
604 
605       // Must always clear pPolygon - for clip to simpleLimit2 and in case of
606       // early exit
607 
608       pPolygon.clear();
609 
610       if ( outputPolygon.empty() )  return;
611 
612       G4VoxelLimits simpleLimit2;
613       simpleLimit2.AddLimit(kZAxis,-kInfinity,pVoxelLimit.GetMaxZExtent());
614       ClipPolygonToSimpleLimits(outputPolygon,pPolygon,simpleLimit2);
615 
616       // Return after final clip - no cleanup
617     }
618   }
619 }
620 
621 ////////////////////////////////////////////////////////////////////////////
622 //
623 // pVoxelLimits must be only limited along one axis, and either the maximum
624 // along the axis must be +kInfinity, or the minimum -kInfinity
625 
626 void
627 G4VSolid::ClipPolygonToSimpleLimits( G4ThreeVectorList& pPolygon,
628                                      G4ThreeVectorList& outputPolygon,
629                                const G4VoxelLimits& pVoxelLimit       ) const
630 {
631   G4int i;
632   auto  noVertices = (G4int)pPolygon.size();
633   G4ThreeVector vEnd,vStart;
634 
635   for (i = 0 ; i < noVertices ; ++i )
636   {
637     vStart = pPolygon[i];
638     if ( i == noVertices-1 )    vEnd = pPolygon[0];
639     else                        vEnd = pPolygon[i+1];
640 
641     if ( pVoxelLimit.Inside(vStart) )
642     {
643       if (pVoxelLimit.Inside(vEnd))
644       {
645         // vStart and vEnd inside -> output end point
646         //
647         outputPolygon.push_back(vEnd);
648       }
649       else
650       {
651         // vStart inside, vEnd outside -> output crossing point
652         //
653         pVoxelLimit.ClipToLimits(vStart,vEnd);
654         outputPolygon.push_back(vEnd);
655       }
656     }
657     else
658     {
659       if (pVoxelLimit.Inside(vEnd))
660       {
661         // vStart outside, vEnd inside -> output inside section
662         //
663         pVoxelLimit.ClipToLimits(vStart,vEnd);
664         outputPolygon.push_back(vStart);
665         outputPolygon.push_back(vEnd);
666       }
667       else  // Both point outside -> no output
668       {
669         // outputPolygon.push_back(vStart);
670         // outputPolygon.push_back(vEnd);
671       }
672     }
673   }
674 }
675 
676 //////////////////////////////////////////////////////////////////////////
677 //
678 // Throw exception (warning) for solids not implementing the method
679 
680 void G4VSolid::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
681 {
682   std::ostringstream message;
683   message << "Not implemented for solid: "
684           << GetEntityType() << " !"
685           << "\nReturning infinite boundinx box.";
686   G4Exception("G4VSolid::BoundingLimits()", "GeomMgt1001",
687               JustWarning, message);
688 
689   pMin.set(-kInfinity,-kInfinity,-kInfinity);
690   pMax.set( kInfinity, kInfinity, kInfinity);
691 }
692 
693 //////////////////////////////////////////////////////////////////////////
694 //
695 // Get G4VisExtent - bounding box for graphics
696 
697 G4VisExtent G4VSolid::GetExtent () const
698 {
699   G4VisExtent extent;
700   G4VoxelLimits voxelLimits;  // Defaults to "infinite" limits.
701   G4AffineTransform affineTransform;
702   G4double vmin, vmax;
703   CalculateExtent(kXAxis,voxelLimits,affineTransform,vmin,vmax);
704   extent.SetXmin (vmin);
705   extent.SetXmax (vmax);
706   CalculateExtent(kYAxis,voxelLimits,affineTransform,vmin,vmax);
707   extent.SetYmin (vmin);
708   extent.SetYmax (vmax);
709   CalculateExtent(kZAxis,voxelLimits,affineTransform,vmin,vmax);
710   extent.SetZmin (vmin);
711   extent.SetZmax (vmax);
712   return extent;
713 }
714 
715 G4Polyhedron* G4VSolid::CreatePolyhedron () const
716 {
717   return nullptr;
718 }
719 
720 G4Polyhedron* G4VSolid::GetPolyhedron () const
721 {
722   return nullptr;
723 }
724