fastAerosol/src/FastAerosolSolid.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026
0027 // --------------------------------------------------------------------
0028 // Implementation for FastAerosolSolid class
0029 // Author: A.Knaian (ara@nklabs.com), N.MacFadden (natemacfadden@gmail.com)
0030 // --------------------------------------------------------------------
0031
0032 #include "FastAerosolSolid.hh"
0033
0034 #include "G4SystemOfUnits.hh"
0035
0036 // calculate extent
0037 #include "G4BoundingEnvelope.hh"
0038 #include "G4AffineTransform.hh"
0039 #include "G4VoxelLimits.hh"
0040
0041 // visualization
0042 #include "G4VGraphicsScene.hh"
0043 #include "G4VisExtent.hh"
0044
0045 // polyhedron
0046 #include "G4AutoLock.hh"
0047 #include "G4Polyhedron.hh"
0048 #include "HepPolyhedronProcessor.h"
0049
0050 namespace
0051 {
0052  G4Mutex polyhedronMutex = G4MUTEX_INITIALIZER;
0053 }
0054
0055
0056 ///////////////////////////////////////////////////////////////////////////////
0057 //
0058 // Constructor
0059 //
0060 FastAerosolSolid::FastAerosolSolid(const G4String& pName,                                FastAerosol* pCloud,                                        G4VSolid* pDroplet,
0061                                          std::function<G4RotationMatrix (G4ThreeVector)> pRotation)
0062     : G4VSolid(pName), fCloud(pCloud), fDroplet(pDroplet), fRotation(pRotation), fRebuildPolyhedron(false), fpPolyhedron(nullptr)
0063 {
0064     // Get cloud size from fCloud
0065     G4ThreeVector cloudPMin, cloudPMax;
0066     fCloud->GetBoundingLimits(cloudPMin, cloudPMax);
0067
0068     fVisDx = cloudPMax.x();
0069     fVisDy = cloudPMax.y();
0070     fVisDz = cloudPMax.z();
0071
0072     // Check and set droplet radius
0073     G4double pR = fCloud->GetRadius();
0074     // would be nice to add a check to make sure pDroplet fits in sphere of radius pR
0075     fR = pR;
0076
0077     fBulk = fCloud->GetBulk();
0078
0079     farFromCloudDist = fCloud->GetPreSphereR()*fR;
0080 }
0081
0082
0083 ///////////////////////////////////////////////////////////////////////////////
0084 //
0085 // Alternative constructor (constant rotation function)
0086 //
0087 FastAerosolSolid::FastAerosolSolid(const G4String& pName,
0088                                          FastAerosol* pCloud,
0089                                          G4VSolid* pDroplet):
0090     FastAerosolSolid(pName, pCloud, pDroplet,
0091                      [](G4ThreeVector) {return G4RotationMatrix();})
0092 {}
0093
0094 ///////////////////////////////////////////////////////////////////////////////
0095 //
0096 // Fake default constructor - sets only member data and allocates memory
0097 //                            for usage restricted to object persistency.
0098 //
0099 FastAerosolSolid::FastAerosolSolid( __void__& a )
0100     : G4VSolid(a), fCloud(nullptr), fDroplet(nullptr),
0101     fBulk(nullptr), fR(0.),
0102     fVisDx(0.), fVisDy(0.), fVisDz(0.),
0103     fCubicVolume(0.), fSurfaceArea(0.),
0104     farFromCloudDist(0.),
0105     fRotation([](G4ThreeVector) {return G4RotationMatrix();}),
0106     fRebuildPolyhedron(false), fpPolyhedron(nullptr)
0107 {
0108 }
0109
0110 ///////////////////////////////////////////////////////////////////////////////
0111 //
0112 // Copy constructor
0113 //
0114 FastAerosolSolid::FastAerosolSolid(const FastAerosolSolid &rhs)
0115     : G4VSolid(rhs), fCloud(rhs.fCloud), fDroplet(rhs.fDroplet),
0116     fBulk(rhs.fBulk), fR(rhs.fR),
0117     fVisDx(rhs.fVisDx), fVisDy(rhs.fVisDy), fVisDz(rhs.fVisDz),
0118     fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
0119     farFromCloudDist(rhs.farFromCloudDist),
0120     fRotation(rhs.fRotation),
0121     fRebuildPolyhedron(rhs.fRebuildPolyhedron), fpPolyhedron(rhs.fpPolyhedron)
0122 {
0123 }
0124
0125 //////////////////////////////////////////////////////////////////////////
0126 //
0127 // Assignment operator
0128 //
0129 FastAerosolSolid &FastAerosolSolid::operator = (const FastAerosolSolid &rhs)
0130 {
0131     // Check assignment to self
0132     //
0133     if (this == &rhs)
0134     {
0135         return *this;
0136     }
0137
0138     // Copy base class data
0139     //
0140     G4VSolid::operator=(rhs);
0141
0142     // Copy data
0143     //
0144     fCloud = rhs.fCloud;
0145     fDroplet = rhs.fDroplet;
0146     fBulk = rhs.fBulk;
0147     fR = rhs.fR;
0148     fVisDx = rhs.fVisDx;
0149     fVisDy = rhs.fVisDy;
0150     fVisDz = rhs.fVisDz;
0151     fCubicVolume = rhs.fCubicVolume;
0152     fSurfaceArea = rhs.fSurfaceArea;
0153     farFromCloudDist = rhs.farFromCloudDist;
0154     fRotation = rhs.fRotation;
0155     fRebuildPolyhedron = rhs.fRebuildPolyhedron;
0156     fpPolyhedron = rhs.fpPolyhedron;
0157
0158
0159     return *this;
0160 }
0161
0162 ///////////////////////////////////////////////////////////////////////////////
0163 //
0164 // Calculate extent under transform and specified limit
0165 //
0166 G4bool FastAerosolSolid::CalculateExtent(const EAxis pAxis,
0167                                             const G4VoxelLimits &pVoxelLimit,
0168                                             const G4AffineTransform &pTransform,
0169                                             G4double &pMin, G4double &pMax) const
0170 {
0171     // Get smallest box to fully contain the cloud of objects, not just the centers
0172     //
0173     G4ThreeVector bmin, bmax;
0174     fCloud->GetBoundingLimits(bmin, bmax);
0175
0176     // Find extent
0177     //
0178     G4BoundingEnvelope bbox(bmin, bmax);
0179     return bbox.CalculateExtent(pAxis, pVoxelLimit, pTransform, pMin, pMax);
0180 }
0181
0182 ///////////////////////////////////////////////////////////////////////////////
0183 //
0184 // Return whether point inside/outside/on surface
0185 //
0186 // This function assumes the cloud has at least 1 droplet
0187 //
0188 EInside FastAerosolSolid::Inside(const G4ThreeVector &p) const
0189 {
0190     G4ThreeVector center;
0191     G4double closestDistance;
0192
0193     fCloud->GetNearestDroplet(p, center, closestDistance, fR, fDroplet, fRotation);
0194
0195     if (closestDistance==0.0)
0196     {
0197         G4RotationMatrix irotm = fRotation(center).inverse();
0198
0199         return fDroplet->Inside( irotm*(p - center) );
0200     }
0201     else
0202     {
0203         return kOutside;
0204     }
0205 }
0206
0207 /////////////////////////////////////////////////////////////////////
0208 //
0209 // Return unit normal of surface closest to p
0210 //
0211 // This function assumes the cloud has at least 1 droplet
0212 //
0213 G4ThreeVector FastAerosolSolid::SurfaceNormal(const G4ThreeVector &p) const
0214 {
0215     G4ThreeVector center;
0216     G4double closestDistance;
0217
0218     fCloud->GetNearestDroplet(p, center, closestDistance, DBL_MAX, fDroplet, fRotation);
0219
0220     G4RotationMatrix rotm = fRotation(center);
0221
0222     return rotm*( fDroplet->SurfaceNormal( rotm.inverse()*(p - center) ) );
0223 }
0224
0225 ///////////////////////////////////////////////////////////////////////////////
0226 //
0227 // Calculate distance to shape from outside, along normalised vector
0228 //
0229 // This CANNOT be an underestimate
0230 //
0231 G4double FastAerosolSolid::DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
0232 {
0233     G4ThreeVector center;
0234     G4double closestDistance;
0235
0236     if (fCloud->GetNearestDroplet(p, v, center, closestDistance, fStepLim, fDroplet, fRotation))    // if we found a droplet within fStepLim of query
0237     {
0238         return closestDistance;
0239     }
0240     else if (fCloud->DistanceToCloud(p,v)<DBL_MAX)          // if there is cloud in front of us
0241     {
0242         return 1.1*fStepLim;
0243     }
0244     else                                                    // flying away from cloud
0245     {
0246         return kInfinity;
0247     }
0248 }
0249
0250 //////////////////////////////////////////////////////////////////////
0251 //
0252 // Calculate distance (<= actual) to closest surface of shape from outside
0253 //
0254 // This function assumes the cloud has at least 1 droplet
0255 //
0256 // This can be an underestimate
0257 //
0258 G4double FastAerosolSolid::DistanceToIn(const G4ThreeVector &p) const
0259 {
0260     G4ThreeVector center;
0261     G4double closestDistance;
0262
0263     G4double distanceToCloud = fBulk->DistanceToIn(p);
0264
0265     if (fBulk->Inside(p)==kOutside && distanceToCloud>=farFromCloudDist)
0266     {
0267         return distanceToCloud;
0268     }
0269     else if (fCloud->GetNearestDroplet(p, center, closestDistance, fStepLim, fDroplet, fRotation))      // if we found a droplet within fStepLim of query
0270     {
0271         return closestDistance;
0272     }
0273     else
0274     {
0275         return 1.1*fStepLim;
0276     }
0277 }
0278
0279 //////////////////////////////////////////////////////////////////////
0280 //
0281 // Calculate distance (<= actual) to closest surface of shape from inside
0282 //
0283 // Despite being a vector distance, we find the absolutely closest
0284 // droplet to our point since we assume that p is in a droplet and p
0285 // could be past the center
0286 //
0287 // This CANNOT be an underestimate
0288 //
0289 G4double FastAerosolSolid::DistanceToOut(const G4ThreeVector &p,
0290                                          const G4ThreeVector &v,
0291                                          const G4bool calcNorm,
0292                                          G4bool *validNorm,
0293                                          G4ThreeVector *n) const
0294 {
0295     G4ThreeVector center;
0296     G4double distanceToIn; // should be 0
0297
0298     fCloud->GetNearestDroplet(p, center, distanceToIn, fR, fDroplet, fRotation);    // if we call this function, must be inside and thus must have a droplet within fR
0299
0300     G4RotationMatrix rotm = fRotation(center);
0301     G4RotationMatrix irotm = rotm.inverse();
0302
0303     G4ThreeVector relPos = irotm*(p-center);
0304
0305     if (fDroplet->Inside(relPos) == kOutside) // something went wrong... we should be inside
0306     {
0307         std::ostringstream message;
0308         message << std::setprecision(15) << "The particle at point p = " << p/mm << "mm"
0309                 << std::setprecision(15) << " called DistanceToOut(p,v) and found the closest droplet to be at center = " << center/mm << "mm"
0310                 << " but p is outside the droplet!";
0311         G4Exception("FastAerosolSolid::DistanceToOut()", "GeomSolids0002",
0312                     FatalErrorInArgument, message);
0313     }
0314
0315     G4double dist = fDroplet->DistanceToOut(relPos, irotm*v, calcNorm, validNorm, n);
0316     *n = rotm*(*n);
0317     *validNorm = false; // even if droplet is convex, the aerosol isn't
0318
0319     return dist;
0320 }
0321
0322 /////////////////////////////////////////////////////////////////////////
0323 //
0324 // Calculate distance (<=actual) to closest surface of shape from inside
0325 //
0326 // This can be an underestimate
0327 //
0328 G4double FastAerosolSolid::DistanceToOut(const G4ThreeVector &p) const
0329 {
0330     G4ThreeVector center;
0331     G4double distanceToIn; // should be 0
0332
0333     fCloud->GetNearestDroplet(p, center, distanceToIn, fR, fDroplet, fRotation);    // if we call this function, must be inside and thus must have a droplet within fR
0334
0335     G4RotationMatrix irotm = fRotation(center).inverse();
0336     G4ThreeVector relPos = irotm*(p-center);
0337
0338     if (fDroplet->Inside(relPos) == kOutside) // something went wrong... we should be inside
0339     {
0340         std::ostringstream message;
0341         message << "The particle at point p = " << p/mm << "mm"
0342                 << " called DistanceToOut(p) and found the closest droplet to be at center = " << center/mm << "mm"
0343                 << " but p is outside the droplet!";
0344         G4Exception("FastAerosolSolid::DistanceToOut()", "GeomSolids0002",
0345                     FatalErrorInArgument, message);
0346     }
0347
0348     return fDroplet->DistanceToOut(relPos);
0349 }
0350
0351 //////////////////////////////////////////////////////////////////////////
0352 //
0353 // G4EntityType
0354 //
0355 G4GeometryType FastAerosolSolid::GetEntityType() const
0356 {
0357     return G4String("FastAerosolSolid");
0358 }
0359
0360 //////////////////////////////////////////////////////////////////////////
0361 //
0362 // G4EntityType
0363 //
0364 G4VSolid* FastAerosolSolid::Clone() const
0365 {
0366     return new FastAerosolSolid(*this);
0367 }
0368
0369 //////////////////////////////////////////////////////////////////////////
0370 //
0371 // Stream object contents to an output stream
0372 //
0373 std::ostream &FastAerosolSolid::StreamInfo(std::ostream &os) const
0374 {
0375     os << "-----------------------------------------------------------\n"
0376        << "    *** Dump for solid - " << GetName() << " ***\n"
0377        << "    ===================================================\n"
0378        << " Solid type: FastAerosolSolid\n"
0379        << " Parameters: \n"
0380        << "    numDroplets: " << fCloud->GetNumDroplets() << "\n"
0381        << "    fDroplet type: " << fDroplet->GetName() << "\n"
0382        << "    fDroplet parameters: \n";
0383        fDroplet->StreamInfo(os);
0384       os  << "-----------------------------------------------------------\n";
0385     return os;
0386 }
0387
0388 ////////////////////////////////////////////////////////////////////////////////
0389 //
0390 // GetPointOnSurface
0391 //
0392 // Currently hardcoded to look at all droplets, not just the populated ones
0393 //
0394 G4ThreeVector FastAerosolSolid::GetPointOnSurface() const
0395 {
0396     G4ThreeVector center;
0397     G4double closestDistance;
0398
0399     G4double fDx = fCloud->GetXHalfLength();
0400     G4double fDy = fCloud->GetYHalfLength();
0401     G4double fDz = fCloud->GetZHalfLength();
0402
0403     G4ThreeVector p(2.0*fDx*G4UniformRand(),2.0*fDy*G4UniformRand(),2.0*fDz*G4UniformRand());
0404     p -= G4ThreeVector(fDx, fDy, fDz);
0405
0406     fCloud->GetNearestDroplet(p, center, closestDistance, DBL_MAX, fDroplet, fRotation);
0407
0408     return(center + fRotation(center)*fDroplet->GetPointOnSurface());
0409 }
0410
0411
0412 /////////////////////////////////////////////////////////////////////////////
0413 //
0414 // Methods for visualisation
0415 //
0416 void FastAerosolSolid::DescribeYourselfTo (G4VGraphicsScene& scene) const
0417 {
0418     scene.AddSolid(*this);
0419 }
0420
0421 G4VisExtent FastAerosolSolid::GetExtent() const
0422 {
0423     return G4VisExtent (-fVisDx, fVisDx, -fVisDy, fVisDy, -fVisDz, fVisDz);
0424 }
0425
0426 G4Polyhedron* FastAerosolSolid::CreatePolyhedron () const
0427 {
0428     return fBulk->CreatePolyhedron();
0429 }
0430
0431
0432 // copied from G4Ellipsoid
0433 G4Polyhedron* FastAerosolSolid::GetPolyhedron () const
0434 {
0435     if (!fpPolyhedron ||
0436         fRebuildPolyhedron ||
0437         fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
0438         fpPolyhedron->GetNumberOfRotationSteps())
0439     {
0440         G4AutoLock l(&polyhedronMutex);
0441         delete fpPolyhedron;
0442         fpPolyhedron = CreatePolyhedron();
0443         fRebuildPolyhedron = false;
0444         l.unlock();
0445     }
0446     return fpPolyhedron;
0447 }