| 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
27 // --------------------------------------------------------------------
28 // Implementation for FastAerosolSolid class
29 // Author: A.Knaian (ara@nklabs.com), N.MacFadden (natemacfadden@gmail.com)
30 // --------------------------------------------------------------------
31
32 #include "FastAerosolSolid.hh"
33
34 #include "G4SystemOfUnits.hh"
35
36 // calculate extent
37 #include "G4BoundingEnvelope.hh"
38 #include "G4AffineTransform.hh"
39 #include "G4VoxelLimits.hh"
40
41 // visualization
42 #include "G4VGraphicsScene.hh"
43 #include "G4VisExtent.hh"
44
45 // polyhedron
46 #include "G4AutoLock.hh"
47 #include "G4Polyhedron.hh"
48 #include "HepPolyhedronProcessor.h"
49
50 namespace
51 {
52 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
53 }
54
55
56 ///////////////////////////////////////////////////////////////////////////////
57 //
58 // Constructor
59 //
60 FastAerosolSolid::FastAerosolSolid(const G4String& pName, FastAerosol* pCloud, G4VSolid* pDroplet,
61 std::function<G4RotationMatrix (G4ThreeVector)> pRotation)
62 : G4VSolid(pName), fCloud(pCloud), fDroplet(pDroplet), fRotation(pRotation), fRebuildPolyhedron(false), fpPolyhedron(nullptr)
63 {
64 // Get cloud size from fCloud
65 G4ThreeVector cloudPMin, cloudPMax;
66 fCloud->GetBoundingLimits(cloudPMin, cloudPMax);
67
68 fVisDx = cloudPMax.x();
69 fVisDy = cloudPMax.y();
70 fVisDz = cloudPMax.z();
71
72 // Check and set droplet radius
73 G4double pR = fCloud->GetRadius();
74 // would be nice to add a check to make sure pDroplet fits in sphere of radius pR
75 fR = pR;
76
77 fBulk = fCloud->GetBulk();
78
79 farFromCloudDist = fCloud->GetPreSphereR()*fR;
80 }
81
82
83 ///////////////////////////////////////////////////////////////////////////////
84 //
85 // Alternative constructor (constant rotation function)
86 //
87 FastAerosolSolid::FastAerosolSolid(const G4String& pName,
88 FastAerosol* pCloud,
89 G4VSolid* pDroplet):
90 FastAerosolSolid(pName, pCloud, pDroplet,
91 [](G4ThreeVector) {return G4RotationMatrix();})
92 {}
93
94 ///////////////////////////////////////////////////////////////////////////////
95 //
96 // Fake default constructor - sets only member data and allocates memory
97 // for usage restricted to object persistency.
98 //
99 FastAerosolSolid::FastAerosolSolid( __void__& a )
100 : G4VSolid(a), fCloud(nullptr), fDroplet(nullptr),
101 fBulk(nullptr), fR(0.),
102 fVisDx(0.), fVisDy(0.), fVisDz(0.),
103 fCubicVolume(0.), fSurfaceArea(0.),
104 farFromCloudDist(0.),
105 fRotation([](G4ThreeVector) {return G4RotationMatrix();}),
106 fRebuildPolyhedron(false), fpPolyhedron(nullptr)
107 {
108 }
109
110 ///////////////////////////////////////////////////////////////////////////////
111 //
112 // Copy constructor
113 //
114 FastAerosolSolid::FastAerosolSolid(const FastAerosolSolid &rhs)
115 : G4VSolid(rhs), fCloud(rhs.fCloud), fDroplet(rhs.fDroplet),
116 fBulk(rhs.fBulk), fR(rhs.fR),
117 fVisDx(rhs.fVisDx), fVisDy(rhs.fVisDy), fVisDz(rhs.fVisDz),
118 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
119 farFromCloudDist(rhs.farFromCloudDist),
120 fRotation(rhs.fRotation),
121 fRebuildPolyhedron(rhs.fRebuildPolyhedron), fpPolyhedron(rhs.fpPolyhedron)
122 {
123 }
124
125 //////////////////////////////////////////////////////////////////////////
126 //
127 // Assignment operator
128 //
129 FastAerosolSolid &FastAerosolSolid::operator = (const FastAerosolSolid &rhs)
130 {
131 // Check assignment to self
132 //
133 if (this == &rhs)
134 {
135 return *this;
136 }
137
138 // Copy base class data
139 //
140 G4VSolid::operator=(rhs);
141
142 // Copy data
143 //
144 fCloud = rhs.fCloud;
145 fDroplet = rhs.fDroplet;
146 fBulk = rhs.fBulk;
147 fR = rhs.fR;
148 fVisDx = rhs.fVisDx;
149 fVisDy = rhs.fVisDy;
150 fVisDz = rhs.fVisDz;
151 fCubicVolume = rhs.fCubicVolume;
152 fSurfaceArea = rhs.fSurfaceArea;
153 farFromCloudDist = rhs.farFromCloudDist;
154 fRotation = rhs.fRotation;
155 fRebuildPolyhedron = rhs.fRebuildPolyhedron;
156 fpPolyhedron = rhs.fpPolyhedron;
157
158
159 return *this;
160 }
161
162 ///////////////////////////////////////////////////////////////////////////////
163 //
164 // Calculate extent under transform and specified limit
165 //
166 G4bool FastAerosolSolid::CalculateExtent(const EAxis pAxis,
167 const G4VoxelLimits &pVoxelLimit,
168 const G4AffineTransform &pTransform,
169 G4double &pMin, G4double &pMax) const
170 {
171 // Get smallest box to fully contain the cloud of objects, not just the centers
172 //
173 G4ThreeVector bmin, bmax;
174 fCloud->GetBoundingLimits(bmin, bmax);
175
176 // Find extent
177 //
178 G4BoundingEnvelope bbox(bmin, bmax);
179 return bbox.CalculateExtent(pAxis, pVoxelLimit, pTransform, pMin, pMax);
180 }
181
182 ///////////////////////////////////////////////////////////////////////////////
183 //
184 // Return whether point inside/outside/on surface
185 //
186 // This function assumes the cloud has at least 1 droplet
187 //
188 EInside FastAerosolSolid::Inside(const G4ThreeVector &p) const
189 {
190 G4ThreeVector center;
191 G4double closestDistance;
192
193 fCloud->GetNearestDroplet(p, center, closestDistance, fR, fDroplet, fRotation);
194
195 if (closestDistance==0.0)
196 {
197 G4RotationMatrix irotm = fRotation(center).inverse();
198
199 return fDroplet->Inside( irotm*(p - center) );
200 }
201 else
202 {
203 return kOutside;
204 }
205 }
206
207 /////////////////////////////////////////////////////////////////////
208 //
209 // Return unit normal of surface closest to p
210 //
211 // This function assumes the cloud has at least 1 droplet
212 //
213 G4ThreeVector FastAerosolSolid::SurfaceNormal(const G4ThreeVector &p) const
214 {
215 G4ThreeVector center;
216 G4double closestDistance;
217
218 fCloud->GetNearestDroplet(p, center, closestDistance, DBL_MAX, fDroplet, fRotation);
219
220 G4RotationMatrix rotm = fRotation(center);
221
222 return rotm*( fDroplet->SurfaceNormal( rotm.inverse()*(p - center) ) );
223 }
224
225 ///////////////////////////////////////////////////////////////////////////////
226 //
227 // Calculate distance to shape from outside, along normalised vector
228 //
229 // This CANNOT be an underestimate
230 //
231 G4double FastAerosolSolid::DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
232 {
233 G4ThreeVector center;
234 G4double closestDistance;
235
236 if (fCloud->GetNearestDroplet(p, v, center, closestDistance, fStepLim, fDroplet, fRotation)) // if we found a droplet within fStepLim of query
237 {
238 return closestDistance;
239 }
240 else if (fCloud->DistanceToCloud(p,v)<DBL_MAX) // if there is cloud in front of us
241 {
242 return 1.1*fStepLim;
243 }
244 else // flying away from cloud
245 {
246 return kInfinity;
247 }
248 }
249
250 //////////////////////////////////////////////////////////////////////
251 //
252 // Calculate distance (<= actual) to closest surface of shape from outside
253 //
254 // This function assumes the cloud has at least 1 droplet
255 //
256 // This can be an underestimate
257 //
258 G4double FastAerosolSolid::DistanceToIn(const G4ThreeVector &p) const
259 {
260 G4ThreeVector center;
261 G4double closestDistance;
262
263 G4double distanceToCloud = fBulk->DistanceToIn(p);
264
265 if (fBulk->Inside(p)==kOutside && distanceToCloud>=farFromCloudDist)
266 {
267 return distanceToCloud;
268 }
269 else if (fCloud->GetNearestDroplet(p, center, closestDistance, fStepLim, fDroplet, fRotation)) // if we found a droplet within fStepLim of query
270 {
271 return closestDistance;
272 }
273 else
274 {
275 return 1.1*fStepLim;
276 }
277 }
278
279 //////////////////////////////////////////////////////////////////////
280 //
281 // Calculate distance (<= actual) to closest surface of shape from inside
282 //
283 // Despite being a vector distance, we find the absolutely closest
284 // droplet to our point since we assume that p is in a droplet and p
285 // could be past the center
286 //
287 // This CANNOT be an underestimate
288 //
289 G4double FastAerosolSolid::DistanceToOut(const G4ThreeVector &p,
290 const G4ThreeVector &v,
291 const G4bool calcNorm,
292 G4bool *validNorm,
293 G4ThreeVector *n) const
294 {
295 G4ThreeVector center;
296 G4double distanceToIn; // should be 0
297
298 fCloud->GetNearestDroplet(p, center, distanceToIn, fR, fDroplet, fRotation); // if we call this function, must be inside and thus must have a droplet within fR
299
300 G4RotationMatrix rotm = fRotation(center);
301 G4RotationMatrix irotm = rotm.inverse();
302
303 G4ThreeVector relPos = irotm*(p-center);
304
305 if (fDroplet->Inside(relPos) == kOutside) // something went wrong... we should be inside
306 {
307 std::ostringstream message;
308 message << std::setprecision(15) << "The particle at point p = " << p/mm << "mm"
309 << std::setprecision(15) << " called DistanceToOut(p,v) and found the closest droplet to be at center = " << center/mm << "mm"
310 << " but p is outside the droplet!";
311 G4Exception("FastAerosolSolid::DistanceToOut()", "GeomSolids0002",
312 FatalErrorInArgument, message);
313 }
314
315 G4double dist = fDroplet->DistanceToOut(relPos, irotm*v, calcNorm, validNorm, n);
316 *n = rotm*(*n);
317 *validNorm = false; // even if droplet is convex, the aerosol isn't
318
319 return dist;
320 }
321
322 /////////////////////////////////////////////////////////////////////////
323 //
324 // Calculate distance (<=actual) to closest surface of shape from inside
325 //
326 // This can be an underestimate
327 //
328 G4double FastAerosolSolid::DistanceToOut(const G4ThreeVector &p) const
329 {
330 G4ThreeVector center;
331 G4double distanceToIn; // should be 0
332
333 fCloud->GetNearestDroplet(p, center, distanceToIn, fR, fDroplet, fRotation); // if we call this function, must be inside and thus must have a droplet within fR
334
335 G4RotationMatrix irotm = fRotation(center).inverse();
336 G4ThreeVector relPos = irotm*(p-center);
337
338 if (fDroplet->Inside(relPos) == kOutside) // something went wrong... we should be inside
339 {
340 std::ostringstream message;
341 message << "The particle at point p = " << p/mm << "mm"
342 << " called DistanceToOut(p) and found the closest droplet to be at center = " << center/mm << "mm"
343 << " but p is outside the droplet!";
344 G4Exception("FastAerosolSolid::DistanceToOut()", "GeomSolids0002",
345 FatalErrorInArgument, message);
346 }
347
348 return fDroplet->DistanceToOut(relPos);
349 }
350
351 //////////////////////////////////////////////////////////////////////////
352 //
353 // G4EntityType
354 //
355 G4GeometryType FastAerosolSolid::GetEntityType() const
356 {
357 return G4String("FastAerosolSolid");
358 }
359
360 //////////////////////////////////////////////////////////////////////////
361 //
362 // G4EntityType
363 //
364 G4VSolid* FastAerosolSolid::Clone() const
365 {
366 return new FastAerosolSolid(*this);
367 }
368
369 //////////////////////////////////////////////////////////////////////////
370 //
371 // Stream object contents to an output stream
372 //
373 std::ostream &FastAerosolSolid::StreamInfo(std::ostream &os) const
374 {
375 os << "-----------------------------------------------------------\n"
376 << " *** Dump for solid - " << GetName() << " ***\n"
377 << " ===================================================\n"
378 << " Solid type: FastAerosolSolid\n"
379 << " Parameters: \n"
380 << " numDroplets: " << fCloud->GetNumDroplets() << "\n"
381 << " fDroplet type: " << fDroplet->GetName() << "\n"
382 << " fDroplet parameters: \n";
383 fDroplet->StreamInfo(os);
384 os << "-----------------------------------------------------------\n";
385 return os;
386 }
387
388 ////////////////////////////////////////////////////////////////////////////////
389 //
390 // GetPointOnSurface
391 //
392 // Currently hardcoded to look at all droplets, not just the populated ones
393 //
394 G4ThreeVector FastAerosolSolid::GetPointOnSurface() const
395 {
396 G4ThreeVector center;
397 G4double closestDistance;
398
399 G4double fDx = fCloud->GetXHalfLength();
400 G4double fDy = fCloud->GetYHalfLength();
401 G4double fDz = fCloud->GetZHalfLength();
402
403 G4ThreeVector p(2.0*fDx*G4UniformRand(),2.0*fDy*G4UniformRand(),2.0*fDz*G4UniformRand());
404 p -= G4ThreeVector(fDx, fDy, fDz);
405
406 fCloud->GetNearestDroplet(p, center, closestDistance, DBL_MAX, fDroplet, fRotation);
407
408 return(center + fRotation(center)*fDroplet->GetPointOnSurface());
409 }
410
411
412 /////////////////////////////////////////////////////////////////////////////
413 //
414 // Methods for visualisation
415 //
416 void FastAerosolSolid::DescribeYourselfTo (G4VGraphicsScene& scene) const
417 {
418 scene.AddSolid(*this);
419 }
420
421 G4VisExtent FastAerosolSolid::GetExtent() const
422 {
423 return G4VisExtent (-fVisDx, fVisDx, -fVisDy, fVisDy, -fVisDz, fVisDz);
424 }
425
426 G4Polyhedron* FastAerosolSolid::CreatePolyhedron () const
427 {
428 return fBulk->CreatePolyhedron();
429 }
430
431
432 // copied from G4Ellipsoid
433 G4Polyhedron* FastAerosolSolid::GetPolyhedron () const
434 {
435 if (!fpPolyhedron ||
436 fRebuildPolyhedron ||
437 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
438 fpPolyhedron->GetNumberOfRotationSteps())
439 {
440 G4AutoLock l(&polyhedronMutex);
441 delete fpPolyhedron;
442 fpPolyhedron = CreatePolyhedron();
443 fRebuildPolyhedron = false;
444 l.unlock();
445 }
446 return fpPolyhedron;
447 }
448