Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4UExtrudedSolid.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 // Implementation of G4UExtrudedSolid wrapper class
 27 //
 28 // 17.11.17 G.Cosmo, CERN
 29 // --------------------------------------------------------------------
 30 
 31 #include "G4ExtrudedSolid.hh"
 32 #include "G4UExtrudedSolid.hh"
 33 
 34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
 35 
 36 #include "G4GeomTools.hh"
 37 #include "G4AffineTransform.hh"
 38 #include "G4BoundingEnvelope.hh"
 39 
 40 ////////////////////////////////////////////////////////////////////////
 41 //
 42 // Constructors
 43 //
 44 G4UExtrudedSolid::G4UExtrudedSolid(const G4String&          name,
 45                                    const std::vector<G4TwoVector>& polygon,
 46                                    const std::vector<ZSection>&    zsections)
 47   : Base_t(name)  // General constructor
 48 {
 49   unsigned int nVertices = polygon.size();
 50   unsigned int nSections = zsections.size();
 51 
 52   auto vertices = new vecgeom::XtruVertex2[nVertices];
 53   auto sections = new vecgeom::XtruSection[nSections];
 54 
 55   for (unsigned int i = 0; i < nVertices; ++i)
 56   {
 57     vertices[i].x = polygon[i].x();
 58     vertices[i].y = polygon[i].y();
 59   }
 60   for (unsigned int i = 0; i < nSections; ++i)
 61   {
 62     sections[i].fOrigin.Set(zsections[i].fOffset.x(),
 63                             zsections[i].fOffset.y(),
 64                             zsections[i].fZ);
 65     sections[i].fScale = zsections[i].fScale;
 66   }
 67   Base_t::Initialize(nVertices, vertices, nSections, sections);
 68   delete[] vertices;
 69   delete[] sections;
 70 }
 71 
 72 
 73 G4UExtrudedSolid::G4UExtrudedSolid(const G4String&          name,
 74                                    const std::vector<G4TwoVector>& polygon,
 75                                    G4double                 halfZ,
 76                                    const G4TwoVector& off1, G4double scale1,
 77                                    const G4TwoVector& off2, G4double scale2)
 78   : Base_t(name)  // Special constructor for 2 sections
 79 {
 80   unsigned int nVertices = polygon.size();
 81   unsigned int nSections = 2;
 82 
 83   auto vertices = new vecgeom::XtruVertex2[nVertices];
 84   auto sections = new vecgeom::XtruSection[nSections];
 85 
 86   for (unsigned int i = 0; i < nVertices; ++i)
 87   {
 88     vertices[i].x = polygon[i].x();
 89     vertices[i].y = polygon[i].y();
 90   }
 91   sections[0].fOrigin.Set(off1.x(), off1.y(), -halfZ);
 92   sections[0].fScale = scale1;
 93   sections[1].fOrigin.Set(off2.x(), off2.y(), halfZ);
 94   sections[1].fScale = scale2;
 95   Base_t::Initialize(nVertices, vertices, nSections, sections);
 96   delete[] vertices;
 97   delete[] sections;
 98 }
 99 
100 ////////////////////////////////////////////////////////////////////////
101 //
102 // Fake default constructor - sets only member data and allocates memory
103 //                            for usage restricted to object persistency.
104 //
105 G4UExtrudedSolid::G4UExtrudedSolid(__void__& a)
106   : Base_t(a)
107 {
108 }
109 
110 
111 //////////////////////////////////////////////////////////////////////////
112 //
113 // Destructor
114 //
115 G4UExtrudedSolid::~G4UExtrudedSolid() = default;
116 
117 
118 //////////////////////////////////////////////////////////////////////////
119 //
120 // Copy constructor
121 //
122 G4UExtrudedSolid::G4UExtrudedSolid(const G4UExtrudedSolid &source)
123   : Base_t(source)
124 {
125 }
126 
127 
128 //////////////////////////////////////////////////////////////////////////
129 //
130 // Assignment operator
131 //
132 G4UExtrudedSolid&
133 G4UExtrudedSolid::operator=(const G4UExtrudedSolid &source)
134 {
135   if (this == &source) return *this;
136   
137   Base_t::operator=( source );
138   
139   return *this;
140 }
141 
142 
143 //////////////////////////////////////////////////////////////////////////
144 //
145 // Accessors
146 
147 G4int G4UExtrudedSolid::GetNofVertices() const
148 {
149   return Base_t::GetNVertices();
150 }
151 
152 G4TwoVector G4UExtrudedSolid::GetVertex(G4int i) const
153 {
154   G4double xx, yy;
155   Base_t::GetVertex(i, xx, yy);
156   return { xx, yy };
157 }
158 
159 std::vector<G4TwoVector> G4UExtrudedSolid::GetPolygon() const
160 {
161   std::vector<G4TwoVector> pol;
162   for (unsigned int i = 0; i < Base_t::GetNVertices(); ++i)
163   {
164     pol.push_back(GetVertex(i));
165   }
166   return pol;
167 }
168 
169 G4int G4UExtrudedSolid::GetNofZSections() const
170 {
171   return Base_t::GetNSections();
172 }
173 
174 G4UExtrudedSolid::ZSection G4UExtrudedSolid::GetZSection(G4int i) const
175 {
176   vecgeom::XtruSection sect = Base_t::GetSection(i);
177   return { sect.fOrigin[2],
178            G4TwoVector(sect.fOrigin[0], sect.fOrigin[1]),
179            sect.fScale };
180 }
181 
182 std::vector<G4UExtrudedSolid::ZSection> G4UExtrudedSolid::GetZSections() const
183 {
184   std::vector<G4UExtrudedSolid::ZSection> sections;
185   for (unsigned int i = 0; i < Base_t::GetNSections(); ++i)
186   {
187     vecgeom::XtruSection sect = Base_t::GetSection(i);
188     sections.emplace_back(sect.fOrigin[2],
189                           G4TwoVector(sect.fOrigin[0], sect.fOrigin[1]),
190                           sect.fScale);
191   }
192   return sections;
193 }
194 
195 
196 ///////////////////////////////////////////////////////////////////////////////
197 //
198 // Get bounding box
199 
200 void G4UExtrudedSolid::BoundingLimits(G4ThreeVector& pMin,
201                                       G4ThreeVector& pMax) const
202 {
203   static G4bool checkBBox = true;
204 
205   G4double xmin0 = kInfinity, xmax0 = -kInfinity;
206   G4double ymin0 = kInfinity, ymax0 = -kInfinity;
207 
208   for (G4int i=0; i<GetNofVertices(); ++i)
209   {
210     G4TwoVector vertex = GetVertex(i);
211     G4double x = vertex.x();
212     if (x < xmin0) xmin0 = x;
213     if (x > xmax0) xmax0 = x;
214     G4double y = vertex.y();
215     if (y < ymin0) ymin0 = y;
216     if (y > ymax0) ymax0 = y;
217   }
218 
219   G4double xmin = kInfinity, xmax = -kInfinity;
220   G4double ymin = kInfinity, ymax = -kInfinity;
221 
222   G4int nsect = GetNofZSections();
223   for (G4int i=0; i<nsect; ++i)
224   {
225     ZSection zsect = GetZSection(i);
226     G4double dx    = zsect.fOffset.x();
227     G4double dy    = zsect.fOffset.y();
228     G4double scale = zsect.fScale;
229     xmin = std::min(xmin,xmin0*scale+dx);
230     xmax = std::max(xmax,xmax0*scale+dx);
231     ymin = std::min(ymin,ymin0*scale+dy);
232     ymax = std::max(ymax,ymax0*scale+dy);
233   }
234 
235   G4double zmin = GetZSection(0).fZ;
236   G4double zmax = GetZSection(nsect-1).fZ;
237 
238   pMin.set(xmin,ymin,zmin);
239   pMax.set(xmax,ymax,zmax);
240 
241   // Check correctness of the bounding box
242   //
243   if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
244   {
245     std::ostringstream message;
246     message << "Bad bounding box (min >= max) for solid: "
247             << GetName() << " !"
248             << "\npMin = " << pMin
249             << "\npMax = " << pMax;
250     G4Exception("G4UExtrudedSolid::BoundingLimits()", "GeomMgt0001",
251                 JustWarning, message);
252     StreamInfo(G4cout);
253   }
254 
255   // Check consistency of bounding boxes
256   //
257   if (checkBBox)
258   {
259     U3Vector vmin, vmax;
260     Base_t::Extent(vmin,vmax);
261     if (std::abs(pMin.x()-vmin.x()) > kCarTolerance ||
262         std::abs(pMin.y()-vmin.y()) > kCarTolerance ||
263         std::abs(pMin.z()-vmin.z()) > kCarTolerance ||
264         std::abs(pMax.x()-vmax.x()) > kCarTolerance ||
265         std::abs(pMax.y()-vmax.y()) > kCarTolerance ||
266         std::abs(pMax.z()-vmax.z()) > kCarTolerance)
267     {
268       std::ostringstream message;
269       message << "Inconsistency in bounding boxes for solid: "
270               << GetName() << " !"
271               << "\nBBox min: wrapper = " << pMin << " solid = " << vmin
272               << "\nBBox max: wrapper = " << pMax << " solid = " << vmax;
273       G4Exception("G4UExtrudedSolid::BoundingLimits()", "GeomMgt0001",
274                   JustWarning, message);
275       checkBBox = false;
276     }
277   }
278 }
279 
280 
281 //////////////////////////////////////////////////////////////////////////////
282 //
283 // Calculate extent under transform and specified limit
284 
285 G4bool
286 G4UExtrudedSolid::CalculateExtent(const EAxis pAxis,
287                                   const G4VoxelLimits& pVoxelLimit,
288                                   const G4AffineTransform& pTransform,
289                                         G4double& pMin, G4double& pMax) const
290 {
291   G4ThreeVector bmin, bmax;
292   G4bool exist;
293 
294   // Check bounding box (bbox)
295   //
296   BoundingLimits(bmin,bmax);
297   G4BoundingEnvelope bbox(bmin,bmax);
298 #ifdef G4BBOX_EXTENT
299   return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
300 #endif
301   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
302   {
303     return exist = pMin < pMax;
304   }
305 
306   // To find the extent, the base polygon is subdivided in triangles.
307   // The extent is calculated as cumulative extent of the parts
308   // formed by extrusion of the triangles
309   //
310   G4TwoVectorList basePolygon = GetPolygon();
311   G4TwoVectorList triangles;
312   G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
313   G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
314 
315   // triangulate the base polygon
316   if (!G4GeomTools::TriangulatePolygon(basePolygon,triangles))
317   {
318     std::ostringstream message;
319     message << "Triangulation of the base polygon has failed for solid: "
320             << GetName() << " !"
321             << "\nExtent has been calculated using boundary box";
322     G4Exception("G4UExtrudedSolid::CalculateExtent()",
323                 "GeomMgt1002",JustWarning,message);
324     return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
325   }
326 
327   // allocate vector lists
328   G4int nsect = GetNofZSections();
329   std::vector<const G4ThreeVectorList *> polygons;
330   polygons.resize(nsect);
331   for (G4int k=0; k<nsect; ++k) { polygons[k] = new G4ThreeVectorList(3); }
332 
333   // main loop along triangles
334   pMin =  kInfinity;
335   pMax = -kInfinity;
336   G4int ntria = triangles.size()/3;
337   for (G4int i=0; i<ntria; ++i)
338   {
339     G4int i3 = i*3;
340     for (G4int k=0; k<nsect; ++k) // extrude triangle
341     {
342       ZSection zsect = GetZSection(k);
343       G4double z     = zsect.fZ;
344       G4double dx    = zsect.fOffset.x();
345       G4double dy    = zsect.fOffset.y();
346       G4double scale = zsect.fScale;
347 
348       auto ptr = const_cast<G4ThreeVectorList*>(polygons[k]);
349       auto iter = ptr->begin();
350       G4double x0 = triangles[i3+0].x()*scale+dx;
351       G4double y0 = triangles[i3+0].y()*scale+dy;
352       iter->set(x0,y0,z);
353       iter++;
354       G4double x1 = triangles[i3+1].x()*scale+dx;
355       G4double y1 = triangles[i3+1].y()*scale+dy;
356       iter->set(x1,y1,z);
357       iter++;
358       G4double x2 = triangles[i3+2].x()*scale+dx;
359       G4double y2 = triangles[i3+2].y()*scale+dy;
360       iter->set(x2,y2,z);
361     }
362 
363     // set sub-envelope and adjust extent
364     G4double emin,emax;
365     G4BoundingEnvelope benv(polygons);
366     if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue;
367     if (emin < pMin) pMin = emin;
368     if (emax > pMax) pMax = emax;
369     if (eminlim > pMin && emaxlim < pMax) break; // max possible extent
370   }
371   // free memory
372   for (G4int k=0; k<nsect; ++k) { delete polygons[k]; polygons[k]=nullptr;}
373   return (pMin < pMax);
374 }
375 
376 
377 ///////////////////////////////////////////////////////////////////////////////
378 //
379 // CreatePolyhedron()
380 //
381 G4Polyhedron* G4UExtrudedSolid::CreatePolyhedron () const
382 {
383   unsigned int nFacets = Base_t::GetStruct().fTslHelper.fFacets.size();
384   unsigned int nVertices = Base_t::GetStruct().fTslHelper.fVertices.size();
385 
386   auto polyhedron = new G4Polyhedron(nVertices, nFacets);
387 
388   // Copy vertices
389   for (unsigned int i = 0; i < nVertices; ++i)
390   {
391     U3Vector v = Base_t::GetStruct().fTslHelper.fVertices[i];
392     polyhedron->SetVertex(i+1, G4ThreeVector(v.x(), v.y(), v.z()));
393   }
394 
395   // Copy facets
396   for (unsigned int i = 0; i < nFacets; ++i)
397   {
398     // Facets are only triangular in VecGeom
399     G4int i1 = Base_t::GetStruct().fTslHelper.fFacets[i]->fIndices[0] + 1;
400     G4int i2 = Base_t::GetStruct().fTslHelper.fFacets[i]->fIndices[1] + 1;
401     G4int i3 = Base_t::GetStruct().fTslHelper.fFacets[i]->fIndices[2] + 1;
402     polyhedron->SetFacet(i+1, i1, i2, i3);
403   }
404   polyhedron->SetReferences();
405 
406   return polyhedron;
407 }
408 
409 #endif  // G4GEOM_USE_USOLIDS
410