| Geant4 Cross Reference |
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 G4UTessellatedSolid wrapper class
27 //
28 // 11.01.18 G.Cosmo, CERN
29 // --------------------------------------------------------------------
30
31 #include "G4TessellatedSolid.hh"
32 #include "G4UTessellatedSolid.hh"
33
34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) )
35
36 #include "G4TriangularFacet.hh"
37 #include "G4QuadrangularFacet.hh"
38
39 #include "G4GeomTools.hh"
40 #include "G4AffineTransform.hh"
41 #include "G4BoundingEnvelope.hh"
42
43 ////////////////////////////////////////////////////////////////////////
44 //
45 // Constructors
46 //
47 G4UTessellatedSolid::G4UTessellatedSolid()
48 : Base_t("")
49 {
50 }
51
52 G4UTessellatedSolid::G4UTessellatedSolid(const G4String& name)
53 : Base_t(name)
54 {
55 }
56
57 ////////////////////////////////////////////////////////////////////////
58 //
59 // Fake default constructor - sets only member data and allocates memory
60 // for usage restricted to object persistency.
61 //
62 G4UTessellatedSolid::G4UTessellatedSolid(__void__& a)
63 : Base_t(a)
64 {
65 }
66
67 //////////////////////////////////////////////////////////////////////////
68 //
69 // Destructor
70 //
71 G4UTessellatedSolid::~G4UTessellatedSolid()
72 {
73 std::size_t size = fFacets.size();
74 for (std::size_t i = 0; i < size; ++i) { delete fFacets[i]; }
75 fFacets.clear();
76 }
77
78 //////////////////////////////////////////////////////////////////////////
79 //
80 // Copy constructor
81 //
82 G4UTessellatedSolid::G4UTessellatedSolid(const G4UTessellatedSolid& source)
83 : Base_t(source)
84 {
85 }
86
87 //////////////////////////////////////////////////////////////////////////
88 //
89 // Assignment operator
90 //
91 G4UTessellatedSolid&
92 G4UTessellatedSolid::operator=(const G4UTessellatedSolid& source)
93 {
94 if (this == &source) return *this;
95
96 Base_t::operator=( source );
97
98 return *this;
99 }
100
101 //////////////////////////////////////////////////////////////////////////
102 //
103 // Accessors
104
105 G4bool G4UTessellatedSolid::AddFacet(G4VFacet* aFacet)
106 {
107 // Add a facet to the structure, checking validity.
108 //
109 if (GetSolidClosed())
110 {
111 G4Exception("G4UTessellatedSolid::AddFacet()", "GeomSolids1002",
112 JustWarning, "Attempt to add facets when solid is closed.");
113 return false;
114 }
115 if (!aFacet->IsDefined())
116 {
117 G4Exception("G4UTessellatedSolid::AddFacet()", "GeomSolids1002",
118 JustWarning, "Attempt to add facet not properly defined.");
119 aFacet->StreamInfo(G4cout);
120 return false;
121 }
122 if (aFacet->GetNumberOfVertices() == 3)
123 {
124 auto a3Facet = dynamic_cast<G4TriangularFacet*>(aFacet);
125 return Base_t::AddTriangularFacet(U3Vector(a3Facet->GetVertex(0).x(),
126 a3Facet->GetVertex(0).y(),
127 a3Facet->GetVertex(0).z()),
128 U3Vector(a3Facet->GetVertex(1).x(),
129 a3Facet->GetVertex(1).y(),
130 a3Facet->GetVertex(1).z()),
131 U3Vector(a3Facet->GetVertex(2).x(),
132 a3Facet->GetVertex(2).y(),
133 a3Facet->GetVertex(2).z()),
134 true);
135 }
136 else if (aFacet->GetNumberOfVertices() == 4)
137 {
138 auto a4Facet = dynamic_cast<G4QuadrangularFacet*>(aFacet);
139 return Base_t::AddQuadrilateralFacet(U3Vector(a4Facet->GetVertex(0).x(),
140 a4Facet->GetVertex(0).y(),
141 a4Facet->GetVertex(0).z()),
142 U3Vector(a4Facet->GetVertex(1).x(),
143 a4Facet->GetVertex(1).y(),
144 a4Facet->GetVertex(1).z()),
145 U3Vector(a4Facet->GetVertex(2).x(),
146 a4Facet->GetVertex(2).y(),
147 a4Facet->GetVertex(2).z()),
148 U3Vector(a4Facet->GetVertex(3).x(),
149 a4Facet->GetVertex(3).y(),
150 a4Facet->GetVertex(3).z()),
151 true);
152 }
153 else
154 {
155 G4Exception("G4UTessellatedSolid::AddFacet()", "GeomSolids1002",
156 JustWarning, "Attempt to add facet not properly defined.");
157 aFacet->StreamInfo(G4cout);
158 return false;
159 }
160 }
161
162 G4VFacet* G4UTessellatedSolid::GetFacet(G4int i) const
163 {
164 return fFacets[i];
165 }
166
167 G4int G4UTessellatedSolid::GetNumberOfFacets() const
168 {
169 return GetNFacets();
170 }
171
172 void G4UTessellatedSolid::SetSolidClosed(const G4bool t)
173 {
174 if (t && !Base_t::IsClosed())
175 {
176 Base_t::Close();
177 std::size_t nVertices = fTessellated.fVertices.size();
178 std::size_t nFacets = fTessellated.fFacets.size();
179 for (std::size_t j = 0; j < nVertices; ++j)
180 {
181 U3Vector vt = fTessellated.fVertices[j];
182 fVertexList.emplace_back(vt.x(), vt.y(), vt.z());
183 }
184 for (std::size_t i = 0; i < nFacets; ++i)
185 {
186 vecgeom::TriangleFacet<G4double>* afacet = Base_t::GetFacet(i);
187 std::vector<G4ThreeVector> v;
188 for (const auto & vertex : afacet->fVertices)
189 {
190 v.emplace_back(vertex.x(),
191 vertex.y(),
192 vertex.z());
193 }
194 G4VFacet* facet = new G4TriangularFacet(v[0], v[1], v[2],
195 G4FacetVertexType::ABSOLUTE);
196 facet->SetVertices(&fVertexList);
197 for (G4int k=0; k<3; ++k)
198 {
199 facet->SetVertexIndex(k, afacet->fIndices[k]);
200 }
201 fFacets.push_back(facet);
202 }
203 }
204 }
205
206 G4bool G4UTessellatedSolid::GetSolidClosed() const
207 {
208 return Base_t::IsClosed();
209 }
210
211 void G4UTessellatedSolid::SetMaxVoxels(G4int)
212 {
213 // Not yet implemented !
214 }
215
216 G4double G4UTessellatedSolid::GetMinXExtent() const
217 {
218 U3Vector aMin, aMax;
219 Base_t::Extent(aMin, aMax);
220 return aMin.x();
221 }
222 G4double G4UTessellatedSolid::GetMaxXExtent() const
223 {
224 U3Vector aMin, aMax;
225 Base_t::Extent(aMin, aMax);
226 return aMax.x();
227 }
228 G4double G4UTessellatedSolid::GetMinYExtent() const
229 {
230 U3Vector aMin, aMax;
231 Base_t::Extent(aMin, aMax);
232 return aMin.y();
233 }
234 G4double G4UTessellatedSolid::GetMaxYExtent() const
235 {
236 U3Vector aMin, aMax;
237 Base_t::Extent(aMin, aMax);
238 return aMax.y();
239 }
240 G4double G4UTessellatedSolid::GetMinZExtent() const
241 {
242 U3Vector aMin, aMax;
243 Base_t::Extent(aMin, aMax);
244 return aMin.z();
245 }
246 G4double G4UTessellatedSolid::GetMaxZExtent() const
247 {
248 U3Vector aMin, aMax;
249 Base_t::Extent(aMin, aMax);
250 return aMax.z();
251 }
252
253 G4int G4UTessellatedSolid::AllocatedMemoryWithoutVoxels()
254 {
255 G4int base = sizeof(*this);
256 base += fVertexList.capacity() * sizeof(G4ThreeVector);
257
258 std::size_t limit = fFacets.size();
259 for (std::size_t i = 0; i < limit; ++i)
260 {
261 G4VFacet &facet = *fFacets[i];
262 base += facet.AllocatedMemory();
263 }
264 return base;
265 }
266 G4int G4UTessellatedSolid::AllocatedMemory()
267 {
268 return AllocatedMemoryWithoutVoxels();
269 }
270 void G4UTessellatedSolid::DisplayAllocatedMemory()
271 {
272 G4int without = AllocatedMemoryWithoutVoxels();
273 // G4int with = AllocatedMemory();
274 // G4double ratio = (G4double) with / without;
275 // G4cout << "G4TessellatedSolid - Allocated memory without voxel overhead "
276 // << without << "; with " << with << "; ratio: " << ratio << G4endl;
277 G4cout << "G4TessellatedSolid - Allocated memory without voxel overhead "
278 << without << G4endl;
279 }
280
281
282 ///////////////////////////////////////////////////////////////////////////////
283 //
284 // Get bounding box
285
286 void G4UTessellatedSolid::BoundingLimits(G4ThreeVector& pMin,
287 G4ThreeVector& pMax) const
288 {
289 U3Vector aMin, aMax;
290 Base_t::Extent(aMin, aMax);
291 pMin = G4ThreeVector(aMin.x(), aMin.y(), aMin.z());
292 pMax = G4ThreeVector(aMax.x(), aMax.y(), aMax.z());
293
294 // Check correctness of the bounding box
295 //
296 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
297 {
298 std::ostringstream message;
299 message << "Bad bounding box (min >= max) for solid: "
300 << GetName() << " !"
301 << "\npMin = " << pMin
302 << "\npMax = " << pMax;
303 G4Exception("G4UTessellatedSolid::BoundingLimits()",
304 "GeomMgt0001", JustWarning, message);
305 StreamInfo(G4cout);
306 }
307 }
308
309
310 //////////////////////////////////////////////////////////////////////////////
311 //
312 // Calculate extent under transform and specified limit
313
314 G4bool
315 G4UTessellatedSolid::CalculateExtent(const EAxis pAxis,
316 const G4VoxelLimits& pVoxelLimit,
317 const G4AffineTransform& pTransform,
318 G4double& pMin, G4double& pMax) const
319 {
320 G4ThreeVector bmin, bmax;
321
322 // Check bounding box (bbox)
323 //
324 BoundingLimits(bmin,bmax);
325 G4BoundingEnvelope bbox(bmin,bmax);
326
327 // Use simple bounding-box to help in the case of complex meshes
328 //
329 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
330
331 #if 0
332 // Precise extent computation (disabled by default for this shape)
333 //
334 G4double kCarToleranceHalf = 0.5*kCarTolerance;
335 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
336 {
337 return (pMin < pMax) ? true : false;
338 }
339
340 // The extent is calculated as cumulative extent of the pyramids
341 // formed by facets and the center of the bounding box.
342 //
343 G4double eminlim = pVoxelLimit.GetMinExtent(pAxis);
344 G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis);
345
346 G4ThreeVectorList base;
347 G4ThreeVectorList apex(1);
348 std::vector<const G4ThreeVectorList *> pyramid(2);
349 pyramid[0] = &base;
350 pyramid[1] = &apex;
351 apex[0] = (bmin+bmax)*0.5;
352
353 // main loop along facets
354 pMin = kInfinity;
355 pMax = -kInfinity;
356 for (G4int i=0; i<GetNumberOfFacets(); ++i)
357 {
358 G4VFacet* facet = GetFacet(i);
359 if (std::abs((facet->GetSurfaceNormal()).dot(facet->GetVertex(0)-apex[0]))
360 < kCarToleranceHalf) continue;
361
362 base.resize(3);
363 for (G4int k=0; k<3; ++k) { base[k] = facet->GetVertex(k); }
364 G4double emin,emax;
365 G4BoundingEnvelope benv(pyramid);
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 return (pMin < pMax);
372 #endif
373 }
374
375
376 ///////////////////////////////////////////////////////////////////////////////
377 //
378 // CreatePolyhedron()
379 //
380 G4Polyhedron* G4UTessellatedSolid::CreatePolyhedron () const
381 {
382 auto nVertices = (G4int)fVertexList.size();
383 auto nFacets = (G4int)fFacets.size();
384 auto polyhedron = new G4Polyhedron(nVertices, nFacets);
385 for (auto i = 0; i < nVertices; ++i)
386 {
387 polyhedron->SetVertex(i+1, fVertexList[i]);
388 }
389
390 for (auto i = 0; i < nFacets; ++i)
391 {
392 G4int v[3]; // Only facets with 3 vertices are defined in VecGeom
393 G4VFacet* facet = GetFacet(i);
394 for (auto j = 0; j < 3; ++j) // Retrieve indexing directly from VecGeom
395 {
396 v[j] = facet->GetVertexIndex(j) + 1;
397 }
398 polyhedron->SetFacet(i+1, v[0], v[1], v[2]);
399 }
400 polyhedron->SetReferences();
401
402 return polyhedron;
403 }
404
405 #endif // G4GEOM_USE_USOLIDS
406