| 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 FastAerosol class
29 // Author: A.Knaian (ara@nklabs.com), N.MacFadden (natemacfadden@gmail.com)
30 // --------------------------------------------------------------------
31
32 #include "FastAerosol.hh"
33
34 #include "G4SystemOfUnits.hh"
35 #include "G4GeometryTolerance.hh"
36
37 #include <numeric> // for summing vectors with accumulate
38
39 // multithreaded safety
40 //#include <atomic>
41 #include "G4AutoLock.hh"
42 namespace
43 {
44 G4Mutex gridMutex = G4MUTEX_INITIALIZER;
45 G4Mutex sphereMutex = G4MUTEX_INITIALIZER;
46 }
47
48 ///////////////////////////////////////////////////////////////////////////////
49 //
50 // Constructor - check inputs
51 // - initialize grid
52 // - cache values
53 //
54 FastAerosol::FastAerosol(const G4String& pName, G4VSolid* pCloud,
55 G4double pR, G4double pMinD, G4double pAvgNumDens,
56 G4double pdR,
57 std::function<G4double (G4ThreeVector)> pDistribution)
58 : fName(pName), fCloud(pCloud), fMinD(pMinD), fDistribution(pDistribution)
59 {
60 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
61
62 // get and set bounding box dimensions
63 G4ThreeVector minBounding, maxBounding;
64 fCloud->BoundingLimits(minBounding, maxBounding);
65 G4ThreeVector halfSizeVec = 0.5*(maxBounding - minBounding);
66
67 G4double pX = halfSizeVec[0];
68 G4double pY = halfSizeVec[1];
69 G4double pZ = halfSizeVec[2];
70
71 if (pX < 2*kCarTolerance ||
72 pY < 2*kCarTolerance ||
73 pZ < 2*kCarTolerance) // limit to thickness of surfaces
74 {
75 std::ostringstream message;
76 message << "Dimensions too small for cloud: "
77 << GetName() << "!" << G4endl
78 << " fDx, fDy, fDz = " << pX << ", " << pY << ", " << pZ;
79 G4Exception("FastAerosol::FastAerosol()", "GeomSolids0002",
80 FatalException, message);
81 }
82 fDx = pX;
83 fDy = pY;
84 fDz = pZ;
85
86
87 // check and set droplet radius
88 if (pR < 0.0)
89 {
90 std::ostringstream message;
91 message << "Invalid radius for cloud: " << GetName() << "!" << G4endl
92 << " Radius must be positive."
93 << " Inputs: pR = " << pR;
94 G4Exception("FastAerosol::FastAerosol()", "GeomSolids0002",
95 FatalErrorInArgument, message);
96 }
97 fR = pR;
98 fR2 = fR*fR;
99
100
101 // check and set droplet radius safety
102 if (pdR < 0.0 || pdR > fR)
103 {
104 std::ostringstream message;
105 message << "Invalid sphericality measure for cloud: " << GetName() << "!" << G4endl
106 << " Radius uncertainty must be between 0 and fR."
107 << " Inputs: pdR = " << pdR
108 << " Inputs: pR = " << pR;
109 G4Exception("FastAerosol::FastAerosol()", "GeomSolids0002",
110 FatalErrorInArgument, message);
111 }
112 fdR = pdR;
113
114
115 // check and set number density
116 if (pAvgNumDens <= 0)
117 {
118 std::ostringstream message;
119 message << "Invalid average number density for cloud: " << GetName() << "!" << G4endl
120 << " pAvgNumDens = " << pAvgNumDens;
121 G4Exception("FastAerosol::FastAerosol()", "GeomSolids0002",
122 FatalException, message);
123 }
124 fAvgNumDens = pAvgNumDens;
125
126
127 // Set collision limit for collsion between equal sized balls with fMinD between them
128 // no droplets will be placed closer than this distance forom each other
129 fCollisionLimit2 = (2*fR + fMinD)*(2*fR + fMinD);
130
131 // set maximum number of droplet that we are allowed to skip before halting with an error
132 fMaxDropCount = (G4int)floor(fAvgNumDens*(fCloud->GetCubicVolume())*(0.01*fMaxDropPercent));
133
134 // initialize grid variables
135 InitializeGrid();
136
137
138 // begin cache of voxelized circles and spheres
139 // see header for more details on these data structures
140 G4AutoLock lockSphere(&sphereMutex); // lock for multithreaded safety. Likely not needed here, but doesn't hurt
141
142 fCircleCollection.push_back({{0, 0}}); // the R=0 circle only has one point {{x1,y1}} = {{0,0}}
143 fSphereCollection.push_back({{{0}}}); // the R=0 sphere only has one point {{{z1}}} = {{{0}}}
144
145 fMaxCircleR = 0;
146 fMaxSphereR = 0;
147
148 MakeSphere(fPreSphereR);
149
150 lockSphere.unlock(); // unlock
151
152 // vector search radius. In terms of voxel width, how far do you search for droplets in vector search
153 // you need to search a larger area if fR is larger than one grid (currently disabled)
154 fVectorSearchRadius = (G4int)ceill(fR/fGridPitch);
155 }
156
157 ///////////////////////////////////////////////////////////////////////////////
158 //
159 // Alternative Constructor
160 //
161 // Same as standard constructor except with a uniform droplet distribution
162 //
163 FastAerosol::FastAerosol(const G4String& pName, G4VSolid* pCloud, G4double pR, G4double pMinD, G4double pNumDens, G4double pdR):
164 FastAerosol(pName, pCloud, pR, pMinD, pNumDens, pdR,
165 [](G4ThreeVector) {return 1.0;})
166 {}
167
168 ///////////////////////////////////////////////////////////////////////////////
169 //
170 // Alternative Constructor
171 //
172 // Same as standard constructor except with a uniform droplet distribution and
173 // assuming no uncertainty in sphericality
174 //
175 FastAerosol::FastAerosol(const G4String& pName, G4VSolid* pCloud, G4double pR, G4double pMinD, G4double pNumDens): FastAerosol(pName, pCloud, pR, pMinD, pNumDens, 0.0,[](G4ThreeVector) {return 1.0;})
176 {}
177
178 ///////////////////////////////////////////////////////////////////////////////
179 //
180 // Initialize grid
181 //
182 // Sets grids to initial values and calculates expected number of droplets
183 // for each voxel
184 //
185 void FastAerosol::InitializeGrid()
186 {
187 // set pitch so, on average, fDropletsPerVoxel droplets are in each voxel
188 fGridPitch = std::pow(fDropletsPerVoxel/fAvgNumDens,1.0/3.0);
189
190 // if a voxel has center farther than this distance from the bulk outside,
191 // we know it is fully contained in the bulk
192 fEdgeDistance = fGridPitch*std::sqrt(3.0)/2.0 + fR;
193
194 // set number of grid cells
195 fNx = (G4int)ceill(2*fDx / fGridPitch);
196 fNy = (G4int)ceill(2*fDy / fGridPitch);
197 fNz = (G4int)ceill(2*fDz / fGridPitch);
198 fNumGridCells = (long) fNx*fNy*fNz;
199 fNxy = fNx*fNy;
200
201 // try... catch... because vectors can only be so long
202 // thus you can't set grid too fine
203 try
204 {
205 std::vector<G4ThreeVector> emptyVoxel{};
206 fGrid.resize(fNumGridCells, emptyVoxel);
207 fGridMean.resize(fNumGridCells, 0);
208 fGridValid = new std::atomic<bool>[fNumGridCells];
209 }
210 catch ( const std::bad_alloc&)
211 {
212 std::ostringstream message;
213 message << "Out of memory! Grid pitch too small for cloud: "
214 << GetName() << "!" << G4endl
215 << " Asked for fNumGridCells = "
216 << fNumGridCells << G4endl
217 << " each with element of size "
218 << sizeof(std::vector<G4ThreeVector>) << G4endl
219 << " each with element of size " << sizeof(bool)
220 << G4endl
221 << " but the max size is " << fGrid.max_size() << "!";
222 G4Exception("FastAerosol::FastAerosol()", "GeomSolids0002",
223 FatalErrorInArgument, message); }
224
225 // initialize grid cache
226 G4ThreeVector voxelCenter;
227 G4bool valid;
228
229 for (int i=0; i<fNumGridCells; i++)
230 {
231 voxelCenter = fGridPitch*(GetIndexCoord(i)+G4ThreeVector(0.5,0.5,0.5)); // center of grid at index i with assuming minimum voxel is at [0,fGridPitch]x[0,fGridPitch]x[0,fGridPitch]
232 voxelCenter -= G4ThreeVector(fDx,fDy,fDz);
233
234 //shift all voxels so that aerosol center is properly at the origin
235 // whether or not the grid is 'finished'
236 // fGridValid[i] is true if we don't plan on populating
237 // more droplets in the voxel
238 // false if otherwise
239 valid = (fCloud->Inside(voxelCenter) != kInside && fCloud->DistanceToIn(voxelCenter) >= fEdgeDistance);
240
241 if (valid) // will not populate the voxel
242 {
243 // have to store 'valid' in this way so that the value is atomic and respects multithreadedness
244 fGridValid[i].store(true, std::memory_order_release);
245 }
246 else // might need to populate the voxel
247 {
248 fGridValid[i].store(false, std::memory_order_release);
249
250 // find the overlap of the voxel with the bulk
251 G4double volScaling=1.0;
252 G4bool edgeVoxel = ( kInside != fCloud->Inside(voxelCenter) || fCloud->DistanceToOut(voxelCenter) < fEdgeDistance );
253 if (edgeVoxel)
254 {
255 volScaling = VoxelOverlap(voxelCenter, 100, 0.0);
256 }
257 // calculates number of droplets based off of voxel center - not ideal
258 fGridMean[i] = std::max(0.0, volScaling*fDistribution(voxelCenter));
259 }
260 }
261
262 // must scale fGridMean[i] so that the desired number densities are actualy achieved
263 // this is because no restrictions are applied to fDistribution
264 G4double tempScaling = (fCloud->GetCubicVolume())*fAvgNumDens/accumulate(fGridMean.begin(), fGridMean.end(), 0.0);
265 for (int i=0; i<fNumGridCells; i++) {fGridMean[i] *= tempScaling;}
266 }
267
268 ///////////////////////////////////////////////////////////////////////////////
269 //
270 // Calculate the overlap of the voxel with the aerosol bulk
271 //
272 // The method is largely copied from G4VSolid::EstimateCubicVolume
273 //
274 G4double FastAerosol::VoxelOverlap(G4ThreeVector voxelCenter, G4int nStat, G4double epsilon)
275 {
276 G4double cenX = voxelCenter.x();
277 G4double cenY = voxelCenter.y();
278 G4double cenZ = voxelCenter.z();
279
280 G4int iInside=0;
281 G4double px,py,pz;
282 G4ThreeVector p;
283 G4bool in;
284
285 if(nStat < 100) nStat = 100;
286 if(epsilon > 0.01) epsilon = 0.01;
287
288 for(G4int i = 0; i < nStat; i++ )
289 {
290 px = cenX+(fGridPitch-epsilon)*(G4UniformRand()-0.5);
291 py = cenY+(fGridPitch-epsilon)*(G4UniformRand()-0.5);
292 pz = cenZ+(fGridPitch-epsilon)*(G4UniformRand()-0.5);
293 p = G4ThreeVector(px,py,pz);
294 in = (fCloud->Inside(p) == kInside) && (fCloud->DistanceToOut(p) >= fR);
295 if(in) iInside++;
296 }
297
298 return (double)iInside/nStat;
299 }
300
301 ///////////////////////////////////////////////////////////////////////////////
302 //
303 // Populate all voxels
304 //
305 // Allows for partially populated clouds.
306 //
307 void FastAerosol::PopulateAllGrids()
308 {
309 unsigned int gi;
310
311 for (int xi=0; xi<fNx; xi++)
312 {
313 for (int yi=0; yi<fNy; yi++)
314 {
315 for (int zi=0; zi<fNz; zi++)
316 {
317 // PopulateGrid takes unsigned arguments
318 // fNx, fNy, and fNz are signed so that some other algebra works
319 // thus we iterate with xi, yi, and zi signed and then cast them
320 PopulateGrid((unsigned)xi,(unsigned)yi,(unsigned)zi,gi);
321 }
322 }
323 }
324 }
325
326 ///////////////////////////////////////////////////////////////////////////////
327 //
328 // Populate a single grid
329 //
330 // If grid is already populated, does nothing.
331 //
332 void FastAerosol::PopulateGrid(unsigned int xi, unsigned int yi, unsigned int zi, unsigned int& gi)
333 {
334 gi = GetGridIndex(xi,yi,zi);
335
336 // Check if this grid needs update
337 bool tmpValid = fGridValid[gi].load(std::memory_order_acquire);
338 // have to do this weirdly because we are using atomic variables so that multithreadedness is safe
339 std::atomic_thread_fence(std::memory_order_acquire);
340
341 if (!tmpValid) // if true then either outside the bulk or in a voxel that is already populated
342 {
343 G4AutoLock lockGrid(&gridMutex);
344 // Check again now that we have the lock - in case grid became valid after first check and before lock
345 tmpValid = fGridValid[gi].load(std::memory_order_acquire);
346
347 if (!tmpValid)
348 {
349 // uniquely set the seed to randomly place droplets
350 // changing global seed gaurantees a totally new batch of seeds
351 fCloudEngine.setSeed((long)gi + (long)fNumGridCells*(long)fSeed);
352
353 // find if the voxel is near the bulk edge. In that case, we need to check whether a placed droplet is inside the bulk
354 // if not an edge voxel, we can place droplets by only considering interference between with other droplets
355 G4ThreeVector voxelCenter = fGridPitch*G4ThreeVector(xi+0.5,yi+0.5,zi+0.5)-G4ThreeVector(fDx,fDy,fDz);
356 G4bool edgeVoxel = ( kInside != fCloud->Inside(voxelCenter) || fCloud->DistanceToOut(voxelCenter) < fEdgeDistance );
357
358 // number of droplets to place
359 unsigned int numDropletsToPlace = CLHEP::RandPoisson::shoot(&fCloudEngine, fGridMean[gi]);
360
361 // actually add the points
362 G4ThreeVector point;
363
364 while (numDropletsToPlace > 0)
365 {
366 // find a new point not overlapping with the others
367 if (FindNewPoint(edgeVoxel, fGridPitch, fGridPitch, fGridPitch, (G4double)xi*fGridPitch-fDx, (G4double)yi*fGridPitch-fDy, (G4double)zi*fGridPitch-fDz, point) )
368 {
369 fGrid[gi].push_back(point);
370 fNumDroplets++;
371 }
372 numDropletsToPlace--;
373 }
374
375 // Memory fence to ensure sequential consistency,
376 // because fGridValid is read without a lock
377 // Voxel data update must be complete before updating fGridValid
378 std::atomic_thread_fence(std::memory_order_release);
379
380 fGridValid[gi].store(true, std::memory_order_release); // set the grid as populated
381 }
382
383 lockGrid.unlock();
384 }
385 }
386
387 ///////////////////////////////////////////////////////////////////////////////
388 //
389 // Find a new droplet position in a box of full (not half) widths dX, dY, and dZ
390 // with minimum corner at (minX, minY, minZ).
391 //
392 // ----------------------------------------------------------------------------
393 //
394 // Ex: FindNewPoint(fGridPitch, fGridPitch, fGridPitch,
395 // (G4double)xi*fGridPitch-fDx, (G4double)yi*fGridPitch-fDy,
396 // (G4double)zi*fGridPitch-fDz);
397 //
398 // 1) CLHEP::RandFlat::shoot(&fCloudEngine)*fGridPitch shoots a point in range
399 // [0, fGridPitch]
400 //
401 // 1.5) Note that the minimum x value of the (xi, yi, zi) grid is
402 // xi*fGridPitch-fDx
403 //
404 // 2) xi*fGridPitch-fDx + (1) adds the minimum x value of the (xi, yi, zi)
405 // voxel
406 //
407 // ----------------------------------------------------------------------------
408 //
409 // Ex: FindNewPoint(2.0*fDx, 2.0*fDy, 2.0*fDz, -fDx, -fDy, -fDz);
410 //
411 // 1) CLHEP::RandFlat::shoot(&fCloudEngine)*2.0*fDx shoots a point in
412 // range [0, 2.0*fDx]
413 //
414 // 2) add -fDx to change range into [-fDx, fDx]
415 //
416 G4bool FastAerosol::FindNewPoint(G4bool edgeVoxel, G4double dX, G4double dY, G4double dZ, G4double minX, G4double minY, G4double minZ, G4ThreeVector &foundPt)
417 {
418 G4int tries = 0;
419 // counter of tries. Give up after fNumNewPointTries (you likely asked for too dense in this case)
420
421 G4double x, y, z;
422
423 G4ThreeVector point;
424 G4bool placedOutside = false;
425 // variable whether or not the point was placed outside. Used for edgeVoxel checking
426
427 // Generate a droplet and check if it overlaps with existing droplets
428 do {
429 tries++;
430 if (tries > fNumNewPointTries)
431 // skip if we tried more than fNumNewPointTries
432 {
433 fNumDropped++;
434 if (fNumDropped < fMaxDropCount)
435 // return error if we skipped more than fMaxDropCount droplets
436 {
437 return false;
438 }
439
440 std::ostringstream message;
441 message << "Threw out too many droplets for cloud: "
442 << GetName() << G4endl
443 << " Tried to place individual droplest "
444 << fNumNewPointTries << " times." << G4endl
445 << " This failed for " << fMaxDropCount
446 << " droplets.";
447 G4Exception("FastAerosol::FindNewPoint()", "GeomSolids0002",
448 FatalErrorInArgument, message);
449 }
450
451 x = minX + CLHEP::RandFlat::shoot(&fCloudEngine)*dX;
452 y = minY + CLHEP::RandFlat::shoot(&fCloudEngine)*dY;
453 z = minZ + CLHEP::RandFlat::shoot(&fCloudEngine)*dZ;
454
455 if (edgeVoxel)
456 {
457 point = G4ThreeVector(x,y,z);
458 placedOutside = (fCloud->Inside(point) != kInside) || (fCloud->DistanceToOut(point) < fR);
459 }
460 }
461 while (CheckCollision(x,y,z) || placedOutside);
462
463 foundPt = G4ThreeVector(x,y,z);
464 return true;
465 }
466
467 ///////////////////////////////////////////////////////////////////////////////
468 //
469 // Get absolutely nearest droplet
470 //
471 // Maximum search radius is maxSearch
472 //
473 // Start in grid of point p, check if any speres are in there
474 // and then go out one more (spherically-shaped) level to the surrounding grids
475 // and so on, until the whole volume has been checked or a droplet has been found
476 // in general as you go out, you always need to check one more level out to make sure that
477 // the one you have (at the closer level) is the actually the nearest one
478 //
479 bool FastAerosol::GetNearestDroplet(const G4ThreeVector &p, G4ThreeVector ¢er, G4double &minRealDistance, G4double maxSearch, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation)
480 {
481 G4double cloudDistance = fCloud->DistanceToIn(p);
482
483 // starting radius/diam of voxel layer
484 G4int searchRad = (G4int)floor(0.5+cloudDistance/fGridPitch); // no reason to search shells totally outside cloud
485
486 // find the voxel containing p
487 int xGrid, yGrid, zGrid;
488 GetGrid(p, xGrid, yGrid, zGrid);
489 // it is OK if the starting point is outside the volume
490
491 // initialize the containers for all candidate closest droplets and their respective distances
492 std::vector<G4ThreeVector> candidates;
493 std::vector<G4double> distances;
494 // initialize distances to indicate that no droplet has yet been found
495 G4double minDistance = kInfinity;
496 G4double maxCheckDistance = kInfinity;
497
498 bool unsafe = true;
499
500 while (unsafe)
501 {
502 // limit maximum search radius
503 if (searchRad*fGridPitch > maxSearch+fGridPitch) { return false; }
504
505 SearchSphere(searchRad, minDistance, candidates, distances, xGrid, yGrid, zGrid, p);
506 maxCheckDistance = minDistance+fdR;
507
508 // theory says that, to safely have searched for centers up to some maxCheckDistance, we must search voxelized spheres at least up to ceil(0.25+maxCheckDistance/fGridPitch)
509 // *** unsure if fully accurate. Calculations for 2D, not 3D... ***
510 unsafe = searchRad < std::ceil(0.25+maxCheckDistance/fGridPitch);
511
512 searchRad++;
513 }
514
515 // delete any collected droplets that are too far out. Check all candidates to see what is closest
516 unsigned int index = 0;
517
518 G4ThreeVector tempCenter;
519 G4double tempDistance;
520
521 minRealDistance = kInfinity;
522
523 while (index < distances.size())
524 {
525 if (distances[index]>maxCheckDistance)
526 {
527 candidates.erase(candidates.begin()+index);
528 distances.erase(distances.begin()+index);
529 }
530 else
531 {
532 tempCenter = candidates[index];
533 tempDistance = droplet->DistanceToIn( rotation(tempCenter).inverse()*(p - tempCenter) );
534
535 if (tempDistance < minRealDistance)
536 {
537 minRealDistance = tempDistance;
538 center = tempCenter;
539 }
540 index++; // move onto next droplet
541 }
542 }
543
544 return true;
545 }
546
547 ///////////////////////////////////////////////////////////////////////////////
548 //
549 // Search sphere of radius R for droplets
550 //
551 // Saves any droplets to "center" if they are closer than minDistance
552 // (updating this minimum found distance at each interval).
553 //
554 // Returns if minDistance<infinity (i.e., if we have yet found a droplet)
555 //
556 void FastAerosol::SearchSphere(G4int searchRad, G4double &minDistance, std::vector<G4ThreeVector> &candidates, std::vector<G4double> &distances, G4int xGrid, G4int yGrid, G4int zGrid, const G4ThreeVector &p)
557 {
558 // search variables
559 G4int xSearch, ySearch, zSearch;
560 // x, y, and z coordinates of the currently searching voxel
561 // in the shell voxel layer variables
562 fSphereType shell; // the shell that we are searching for droplets
563
564 // we pre-calculate spheres up to radius fPreSphereR to speed up calculations
565 // any sphere smaller than that does not need to use locks
566 if (searchRad > fPreSphereR)
567 {
568 // need to consider making a sphere
569 G4AutoLock lockSphere(&sphereMutex);
570 if (searchRad > fMaxSphereR)
571 {
572 // actually need to make a sphere
573 shell = MakeSphere(searchRad);
574 }
575 else
576 {
577 // we have previously made a sphere large enough
578 shell = fSphereCollection[searchRad];
579 }
580 lockSphere.unlock();
581 }
582 else
583 {
584 shell = fSphereCollection[searchRad];
585 }
586
587 // search spherical voxel layer
588 for (int i=0; i<(2*searchRad+1); i++)
589 {
590 xSearch = xGrid+i-searchRad;
591 if (0<=xSearch && xSearch<fNx)
592 {
593 for (int j=0; j<(2*searchRad+1); j++)
594 {
595 ySearch = yGrid+j-searchRad;
596
597 if (0<=ySearch && ySearch<fNy)
598 {
599 for (auto it = shell[i][j].begin(); it != shell[i][j].end(); ++it) {
600 zSearch = *it+zGrid;
601
602 if (0<=zSearch && zSearch<fNz)
603 { GetNearestDropletInsideGrid(minDistance, candidates, distances, (unsigned)xSearch, (unsigned)ySearch, (unsigned)zSearch, p);
604 }
605 }
606 }
607 }
608 }
609 }
610 }
611
612 ///////////////////////////////////////////////////////////////////////////////
613 //
614 // Make voxelized spheres up to radius R for fSphereCollection
615 //
616 // These spheres are used for searching for droplets
617 //
618 // To create them, we first create a helper half-circle "zr" for which each
619 // voxel/point in zr represents a layer of our sphere. The 2nd coordinate of
620 // zr represents the radius of the layer and the 1st coordinate denotes the
621 // layer's displacement from the origin along the z-axis (layers are perp to
622 // z in our convention).
623 //
624 std::vector<std::vector<std::vector<int>>> FastAerosol::MakeSphere(G4int R) {
625 // inductively make smaller spheres since fSphereCollection organizes by index
626 if (fMaxSphereR<(R-1)) { MakeSphere(R-1);}
627
628 std::vector<int> z;
629 // z-coordinates of voxels in (x,y) of sphere
630
631 std::vector<std::vector<int>> y(2*R+1, z);
632 // plane YZ of sphere
633
634 fSphereType x(2*R+1, y); // entire sphere
635
636 x[R][R].push_back(R);// add (0,0,R)
637 x[R][R].push_back(-R);// add (0,0,-R)
638 fCircleType zr = MakeHalfCircle(R);
639 // break sphere into a collection of circles centered at (0,0,z) for different -R<=z<=R
640
641 for (auto it1 = zr.begin(); it1 != zr.end(); ++it1)
642 {
643 std::vector<int> pt1 = *it1;
644 fCircleType circ = MakeCircle(pt1[1]); // make the circle
645
646 for (auto it2 = circ.begin(); it2 != circ.end(); ++it2)
647 {
648 std::vector<int> pt2 = *it2;
649 x[pt2[0]+R][pt2[1]+R].push_back(pt1[0]);
650 }
651 }
652
653 fSphereCollection.push_back(x);
654 fMaxSphereR = R;
655 return x;
656 }
657
658 ///////////////////////////////////////////////////////////////////////////////
659 //
660 // Make voxelized circles up to radius R for fCircleCollection
661 //
662 // These circles are used to create voxelized spheres (these spheres are then
663 // used for searching for droplets). This just ports the calculation to making
664 // a half-circle, adding the points (R,0) and (0,R), and reflecting the half-
665 // circle along the x-axis.
666 //
667 std::vector<std::vector<int>> FastAerosol::MakeCircle(G4int R)
668 {
669 // inductively make smaller circles since fCircleCollection organizes by index
670 if (fMaxCircleR<R)
671 {
672 if (fMaxCircleR<(R-1)) {MakeCircle(R-1);}
673 fCircleType voxels = {{R, 0}, {-R, 0}};
674 // add (R,0) and (-R,0) since half circle excludes them
675 fCircleType halfVoxels = MakeHalfCircle(R);
676
677 // join the voxels, halfVoxels, and -halfVoxels
678 for (auto it = halfVoxels.begin(); it != halfVoxels.end(); ++it)
679 {
680 std::vector<int> pt = *it;
681 voxels.push_back(pt);
682 voxels.push_back({-pt[0], -pt[1]});
683 }
684 fCircleCollection.push_back(voxels);
685 fMaxCircleR++;
686 }
687
688 return fCircleCollection[R];
689 }
690
691 ///////////////////////////////////////////////////////////////////////////////
692 //
693 // Make half circle by a modified midpoint circle algorithm
694 //
695 // Modifies the algorithm so it doesn't have 'holes' when making many circles
696 //
697 // See B. Roget, J. Sitaraman, "Wall distance search algorithim using
698 // voxelized marching spheres," Journal of Computational Physics, Vol. 241,
699 // May 2013, pp. 76-94, https://doi.org/10.1016/j.jcp.2013.01.035
700 //
701 std::vector<std::vector<int>> FastAerosol::MakeHalfCircle(G4int R)
702 {
703 // makes an octant of a voxelized circle in the 1st quadrant starting at y-axis
704 fCircleType voxels = {{0, R}}; // hard code inclusion of (0,R)
705 int x = 1;
706 int y = R;
707 int dx = 4-R;
708 // measure of whether (x+1,y) will be outside the sphere
709 int dxup = 1;
710 // measure of whether (x+1,y+1) will be outside the sphere
711 while (x<y)
712 {
713 // mirror the points to be added to change octant->upper half plane
714 voxels.push_back({ x, y});
715 voxels.push_back({ y, x});
716 voxels.push_back({-x, y});
717 voxels.push_back({-y, x});
718 // if next point outside the sphere, de-increment y
719 if (dx>0)
720 {
721 // if dxup<=0, the layer above just moves to the right
722 // this would create a hole at (x+1,y) unless we add it
723 dxup = dx + 2*y -2*R-1;
724 if ((dxup<=0) && (x+1<y))
725 {
726 voxels.push_back({ 1+x, y});
727 voxels.push_back({ y, 1+x});
728 voxels.push_back({-1-x, y});
729 voxels.push_back({ -y, 1+x});
730 }
731
732 dx += 4-2*y;
733 y--;
734 }
735
736 dx += 1+2*x;
737 x++;
738 }
739
740 // add points on y=+-x
741 if (x==y || (x==y+1 && dxup<=0))
742 {
743 voxels.push_back({x,x});
744 voxels.push_back({-x,x});
745 }
746
747 return voxels;
748 }
749
750 ///////////////////////////////////////////////////////////////////////////////
751 //
752 // Get nearest droplet inside a voxel at (xGrid, yGrid, zGrid)
753 //
754 void FastAerosol::GetNearestDropletInsideGrid(G4double &minDistance, std::vector<G4ThreeVector> &candidates, std::vector<G4double> &distances, unsigned int xGrid, unsigned int yGrid, unsigned int zGrid, const G4ThreeVector &p)
755 {
756 unsigned int gi;
757 PopulateGrid(xGrid, yGrid, zGrid, gi);
758
759 // initialize values
760 G4double foundDistance;
761 // std::vector<G4ThreeVector>::iterator bestPt;
762
763 // find closest droplet
764 for (auto it = fGrid[gi].begin(); it != fGrid[gi].end(); ++it)
765 {
766 foundDistance = std::sqrt((p-*it).mag2());
767
768 if (foundDistance < minDistance+fdR)
769 {
770 minDistance = std::min(minDistance, foundDistance);
771 candidates.push_back(*it);
772 distances.push_back(foundDistance);
773 }
774 }
775 }
776
777 ///////////////////////////////////////////////////////////////////////////////
778 //
779 // Get nearest droplet along vector
780 //
781 // Maximum search radius is fVectorSearchRadius
782 // Maximum distance travelled along vector is maxSearch
783 //
784 // 1) Find the closest voxel along v that is inside the bulk
785 // 2) Search for droplets in a box width (2*fVectorSearchRadius+1) around this voxel
786 // that our ray would intersect
787 // 3) If so, save the droplet that is closest along the ray to p and return true
788 // 4) If not, step 0.99*fGridPitch along v until in a new voxel
789 // 5) If outside bulk, jump until inside bulk and then start again from 2
790 // 6) If inside bulk, search the new voxels in a box centered of width (2*fVectorSearchRadius+1)
791 // centered at our new point
792 // 7) If we find a point, do as in (3). If not, do as in (4)
793 //
794 // This is searching box-shaped regions of voxels, an approach that is only valid
795 // for droplets which fit inside a voxel, which is one reason why the code checks
796 // to make sure that fR < fGridPitch at initialization.
797 //
798 bool FastAerosol::GetNearestDroplet(const G4ThreeVector &p, const G4ThreeVector &normalizedV, G4ThreeVector ¢er, G4double &minDistance, G4double maxSearch, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation)
799 {
800 // get the grid box this initial position is inside
801 int xGrid, yGrid, zGrid;
802 GetGrid(p, xGrid, yGrid, zGrid);
803
804 // initialize some values
805 G4ThreeVector currentP = p;
806 G4double distanceAlongVector = 0.0;
807 G4double deltaDistanceAlongVector;
808
809 int newXGrid = 0; int newYGrid = 0; int newZGrid = 0;
810 int deltaX = 0; int deltaY = 0; int deltaZ = 0;
811 int edgeX = 0 ; int edgeY = 0; int edgeZ = 0;
812 G4bool boxSearch = true;
813
814 G4double distanceToCloud;
815
816 minDistance = kInfinity;
817
818 // actual search
819 while(true)
820 {
821 deltaDistanceAlongVector = 0.0;
822 // Jump gaps
823 distanceToCloud = DistanceToCloud(currentP,normalizedV);
824 if (distanceToCloud == kInfinity)
825 {return(minDistance != kInfinity);}
826 else if (distanceToCloud > fGridPitch)
827 {
828 deltaDistanceAlongVector = distanceToCloud-fGridPitch;
829 distanceAlongVector += deltaDistanceAlongVector;
830 boxSearch = true;// if jumping gap, use box search
831 currentP += deltaDistanceAlongVector*normalizedV;
832 // skip gaps
833 }
834
835 // Search for droplets
836 if (distanceAlongVector > maxSearch)// quit if will jump too far
837 {
838 return(minDistance != kInfinity);
839 }
840 else
841 {
842 if (boxSearch) // we jumped
843 {
844 GetGrid(currentP, xGrid, yGrid, zGrid);
845 GetNearestDropletInsideRegion(minDistance, center, xGrid, yGrid, zGrid,
846 fVectorSearchRadius, fVectorSearchRadius, fVectorSearchRadius, p, normalizedV, droplet, rotation);
847 }
848 else // didn't jump
849 {
850 // Searching endcaps
851 // =================
852 if (deltaX != 0)
853 {
854 GetNearestDropletInsideRegion(minDistance, center,
855 edgeX, yGrid, zGrid, 0, fVectorSearchRadius, fVectorSearchRadius, p, normalizedV,
856 droplet, rotation);}
857
858 if (deltaY != 0)
859 {
860 GetNearestDropletInsideRegion(minDistance, center,
861 xGrid, edgeY, zGrid, fVectorSearchRadius, 0, fVectorSearchRadius, p, normalizedV, droplet, rotation);
862 }
863
864 if (deltaZ != 0)
865 {
866 GetNearestDropletInsideRegion(minDistance, center, xGrid, yGrid, edgeZ, fVectorSearchRadius, fVectorSearchRadius, 0, p, normalizedV, droplet, rotation);
867 }
868
869 // Search bars
870 // ===========
871 // (a bar joins two endcaps)
872 if (deltaX != 0 && deltaY != 0)
873 { GetNearestDropletInsideRegion(minDistance, center, edgeX, edgeY, zGrid,
874 0, 0, fVectorSearchRadius, p, normalizedV,
875 droplet, rotation);
876 }
877
878 if (deltaX != 0 && deltaZ != 0)
879 { GetNearestDropletInsideRegion(minDistance, center, edgeX, yGrid, edgeZ, 0, fVectorSearchRadius, 0, p, normalizedV, droplet, rotation);}
880
881 if (deltaY != 0 && deltaZ != 0)
882 {
883 GetNearestDropletInsideRegion(minDistance, center,
884 xGrid, edgeY, edgeZ,
885 fVectorSearchRadius, 0, 0,
886 p, normalizedV,droplet, rotation);
887 }
888
889 // Search corner
890 // =============
891 if (deltaX != 0 && deltaY != 0 && deltaZ != 0) { GetNearestDropletInsideRegion(minDistance, center, edgeX, edgeY, edgeZ, 0, 0, 0, p, normalizedV, droplet, rotation);}
892
893 // Update position
894 // ================================
895 xGrid = newXGrid; yGrid = newYGrid; zGrid = newZGrid;
896 }
897
898 // Check if we are done
899 // ====================
900 if (distanceAlongVector>minDistance+fdR) {return(true);}
901
902 // walk along the grid
903 // advance by 0.99 fGridPitch (0.99 rather than 1.0 to avoid issues
904 //caused by rounding errors for lines parallel to the grid and based on a grid boundary)
905
906 while (true) {
907 deltaDistanceAlongVector = fGridPitch*0.99;
908 distanceAlongVector += deltaDistanceAlongVector;
909
910 // exit returning false if we have traveled more than MaxVectorFollowingDistance
911 if (distanceAlongVector > maxSearch) { return(false); }
912
913 currentP += deltaDistanceAlongVector*normalizedV;
914 GetGrid(currentP, newXGrid, newYGrid, newZGrid);
915
916 if ((newXGrid != xGrid) || (newYGrid != yGrid) || (newZGrid != zGrid))
917 {
918 deltaX = newXGrid - xGrid;
919 edgeX = xGrid+deltaX*(1+fVectorSearchRadius);
920 deltaY = newYGrid - yGrid;
921 edgeY = yGrid+deltaY*(1+fVectorSearchRadius);
922 deltaZ = newZGrid - zGrid;
923 edgeZ = zGrid+deltaZ*(1+fVectorSearchRadius);
924
925 break;
926 }
927 }
928 boxSearch = false;
929 }
930 }
931
932 return false;
933 }
934
935 ///////////////////////////////////////////////////////////////////////////////
936 //
937 // Get nearest droplet (along a vector normalizedV) inside a voxel centered at
938 // (xGridCenter, yGridCenter, zGridCenter) of width (xWidth, yWidth, zWidth)
939 //
940 // This searches box-shaped regions.
941 //
942 void FastAerosol::GetNearestDropletInsideRegion(G4double &minDistance, G4ThreeVector ¢er, int xGridCenter, int yGridCenter, int zGridCenter, int xWidth, int yWidth, int zWidth, const G4ThreeVector &p, const G4ThreeVector &normalizedV, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation)
943 {
944 for (int xGrid = (xGridCenter-xWidth); xGrid<=(xGridCenter+xWidth); xGrid++)
945 {
946 if (0<=xGrid && xGrid<fNx)
947 {
948 for (int yGrid = (yGridCenter-yWidth); yGrid<=(yGridCenter+yWidth); yGrid++)
949 {
950 if (0<=yGrid && yGrid<fNy)
951 {
952 for (int zGrid = (zGridCenter-zWidth); zGrid<=(zGridCenter+zWidth); zGrid++)
953 {
954 if (0<=zGrid && zGrid<fNz)
955 {
956 GetNearestDropletInsideGrid(minDistance, center, (unsigned)xGrid, (unsigned)yGrid, (unsigned)zGrid, p, normalizedV, droplet, rotation);
957 }
958 }
959 }
960 }
961 }
962 }
963 }
964
965 ///////////////////////////////////////////////////////////////////////////////
966 //
967 // Get nearest droplet inside a voxel at (xGrid, yGrid, zGrid) along a vector
968 //
969 // Return the closest one, as measured along the line
970 //
971 void FastAerosol::GetNearestDropletInsideGrid(G4double &minDistance, G4ThreeVector ¢er, unsigned int xGrid, unsigned int yGrid, unsigned int zGrid, const G4ThreeVector &p, const G4ThreeVector &normalizedV, G4VSolid* droplet, std::function<G4RotationMatrix (G4ThreeVector)> rotation)
972 {
973 unsigned int gi;
974 PopulateGrid(xGrid, yGrid, zGrid, gi);
975
976 G4double foundDistance;
977
978 // find closest droplet
979 for (auto it = fGrid[gi].begin(); it != fGrid[gi].end(); ++it)
980 {
981 // could have the following check to see if the ray pierces the bounding sphere. Currently seems like unnecessary addition
982 /*
983 deltaP = *it-p;
984 foundDistance = normalizedV.dot(deltaP);
985
986 if ((0<=foundDistance) && ((deltaP - normalizedV*foundDistance).mag2() < fR2)) {}
987 */
988
989 G4RotationMatrix irotm = rotation(*it).inverse();
990
991 G4ThreeVector relPos = irotm*(p - *it);
992
993 if (droplet->Inside(relPos) == kInside)
994 {
995 minDistance = 0;
996 center = *it;
997 }
998 else
999 {
1000 foundDistance = droplet->DistanceToIn( relPos , irotm*normalizedV);
1001
1002 if (foundDistance<minDistance)
1003 {
1004 minDistance = foundDistance;
1005 center = *it;
1006 }
1007 }
1008 }
1009 }
1010
1011 ///////////////////////////////////////////////////////////////////////////////
1012 //
1013 // Check if a droplet at the indicated point would collide with an existing droplet
1014 //
1015 // Search neighboring voxels, too. Returns true if there is a collision
1016 //
1017 bool FastAerosol::CheckCollision(G4double x, G4double y, G4double z)
1018 {
1019 G4ThreeVector p(x,y,z);
1020
1021 std::pair<int, int> minMaxXGrid, minMaxYGrid, minMaxZGrid;
1022
1023 int xGrid, yGrid, zGrid;
1024 GetGrid(p, xGrid, yGrid, zGrid);
1025
1026 minMaxXGrid = GetMinMaxSide(xGrid, fNx);
1027 minMaxYGrid = GetMinMaxSide(yGrid, fNy);
1028 minMaxZGrid = GetMinMaxSide(zGrid, fNz);
1029
1030 for (int xi=minMaxXGrid.first; xi <= minMaxXGrid.second; xi++) {
1031 for (int yi = minMaxYGrid.first; yi <= minMaxYGrid.second; yi++) {
1032 for (int zi = minMaxZGrid.first; zi <= minMaxZGrid.second; zi++) {
1033 if (CheckCollisionInsideGrid(x, y, z, (unsigned)xi, (unsigned)yi, (unsigned)zi)) {
1034 fNumCollisions++; // log number of collisions for statistics print
1035 return(true);
1036 }
1037 }
1038 }
1039 }
1040 return(false);
1041 }
1042
1043 ///////////////////////////////////////////////////////////////////////////////
1044 //
1045 // Check if a droplet at the indicated point would collide with an existing
1046 // droplet inside a specific grid
1047 //
1048 // Note that you don't need to lock the mutex since this is only called by code
1049 // that already has the mutex (always called by PopulateGrid).
1050 //
1051 bool FastAerosol::CheckCollisionInsideGrid(G4double x, G4double y, G4double z, unsigned int xi, unsigned int yi, unsigned int zi)
1052 {
1053 std::vector<G4ThreeVector> *thisGrid = &(fGrid[GetGridIndex(xi, yi, zi)]);
1054 unsigned int numel = (unsigned int)thisGrid->size();
1055
1056 for (unsigned int i=0; i < numel; i++)
1057 {
1058 if (CheckCollisionWithDroplet(x, y, z, (*thisGrid)[i]))
1059 return(true);
1060 }
1061
1062 return(false);
1063 }
1064
1065 ///////////////////////////////////////////////////////////////////////////////
1066 //
1067 // Check for collsion with a specific droplet
1068 //
1069 bool FastAerosol::CheckCollisionWithDroplet(G4double x, G4double y, G4double z, G4ThreeVector p )
1070 {
1071 return( std::pow(x-p.x(), 2.0) + std::pow(y-p.y(), 2.0) + std::pow(z-p.z(), 2.0) < fCollisionLimit2 );
1072 }
1073
1074 ///////////////////////////////////////////////////////////////////////////////
1075 //
1076 // Save droplet positions to a file for visualization and analysis
1077 //
1078 void FastAerosol::SaveToFile(const char* filename)
1079 {
1080 G4cout << "Saving droplet positions to " << filename << "..." << G4endl;
1081 std::ofstream file;
1082 file.open(filename);
1083
1084 std::vector<G4ThreeVector> voxel;
1085 G4ThreeVector pt;
1086
1087 for (auto it1 = fGrid.begin(); it1 != fGrid.end(); ++it1)
1088 {
1089 voxel = *it1;
1090
1091 for (auto it2 = voxel.begin(); it2 != voxel.end(); ++it2)
1092 {
1093 pt = *it2;
1094
1095 double x = pt.x();
1096 double y = pt.y();
1097 double z = pt.z();
1098 file << std::setprecision(15) << x/mm << ","
1099 << y/mm << "," << z/mm << "\n";
1100 }
1101 }
1102 file.close();
1103 }
1104