[831] | 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 intellectual property of the * |
---|
| 19 | // * Vanderbilt University Free Electron Laser Center * |
---|
| 20 | // * Vanderbilt University, Nashville, TN, USA * |
---|
| 21 | // * Development supported by: * |
---|
| 22 | // * United States MFEL program under grant FA9550-04-1-0045 * |
---|
| 23 | // * and NASA under contract number NNG04CT05P * |
---|
| 24 | // * Written by Marcus H. Mendenhall and Robert A. Weller. * |
---|
| 25 | // * * |
---|
| 26 | // * Contributed to the Geant4 Core, January, 2005. * |
---|
| 27 | // * * |
---|
| 28 | // ******************************************************************** |
---|
| 29 | // |
---|
| 30 | // $Id: G4Tet.cc,v 1.11 2006/11/13 08:58:03 gcosmo Exp $ |
---|
[1058] | 31 | // GEANT4 tag $Name: geant4-09-02-ref-02 $ |
---|
[831] | 32 | // |
---|
| 33 | // class G4Tet |
---|
| 34 | // |
---|
| 35 | // Implementation for G4Tet class |
---|
| 36 | // |
---|
| 37 | // History: |
---|
| 38 | // |
---|
| 39 | // 20040903 - Marcus Mendenhall, created G4Tet |
---|
| 40 | // 20041101 - Marcus Mendenhall, optimized constant dot products with |
---|
| 41 | // fCdotNijk values |
---|
| 42 | // 20041101 - MHM removed tracking error by clipping DistanceToOut to 0 |
---|
| 43 | // for surface cases |
---|
| 44 | // 20041101 - MHM many speed optimizations in if statements |
---|
| 45 | // 20041101 - MHM changed vdotn comparisons to 1e-12 instead of 0.0 to |
---|
| 46 | // avoid nearly-parallel problems |
---|
| 47 | // 20041102 - MHM Added extra distance into solid to DistanceToIn(p,v) |
---|
| 48 | // hit testing |
---|
| 49 | // 20041102 - MHM added ability to check for degeneracy without throwing |
---|
| 50 | // G4Exception |
---|
| 51 | // 20041103 - MHM removed many unused variables from class |
---|
| 52 | // 20040803 - Dionysios Anninos, added GetPointOnSurface() method |
---|
| 53 | // 20061112 - MHM added code for G4VSolid GetSurfaceArea() |
---|
| 54 | // |
---|
| 55 | // -------------------------------------------------------------------- |
---|
| 56 | |
---|
| 57 | #include "G4Tet.hh" |
---|
| 58 | |
---|
| 59 | const char G4Tet::CVSVers[]="$Id: G4Tet.cc,v 1.11 2006/11/13 08:58:03 gcosmo Exp $"; |
---|
| 60 | |
---|
| 61 | #include "G4VoxelLimits.hh" |
---|
| 62 | #include "G4AffineTransform.hh" |
---|
| 63 | |
---|
| 64 | #include "G4VPVParameterisation.hh" |
---|
| 65 | |
---|
| 66 | #include "Randomize.hh" |
---|
| 67 | |
---|
| 68 | #include "G4VGraphicsScene.hh" |
---|
| 69 | #include "G4Polyhedron.hh" |
---|
| 70 | #include "G4NURBS.hh" |
---|
| 71 | #include "G4NURBSbox.hh" |
---|
| 72 | #include "G4VisExtent.hh" |
---|
| 73 | |
---|
| 74 | #include "G4ThreeVector.hh" |
---|
| 75 | |
---|
| 76 | #include <cmath> |
---|
| 77 | |
---|
| 78 | using namespace CLHEP; |
---|
| 79 | |
---|
| 80 | //////////////////////////////////////////////////////////////////////// |
---|
| 81 | // |
---|
| 82 | // Constructor - create a tetrahedron |
---|
| 83 | // This class is implemented separately from general polyhedra, |
---|
| 84 | // because the simplex geometry can be computed very quickly, |
---|
| 85 | // which may become important in situations imported from mesh generators, |
---|
| 86 | // in which a very large number of G4Tets are created. |
---|
| 87 | // A Tet has all of its geometrical information precomputed |
---|
| 88 | |
---|
| 89 | G4Tet::G4Tet(const G4String& pName, |
---|
| 90 | G4ThreeVector anchor, |
---|
| 91 | G4ThreeVector p2, |
---|
| 92 | G4ThreeVector p3, |
---|
| 93 | G4ThreeVector p4, G4bool *degeneracyFlag) |
---|
| 94 | : G4VSolid(pName), fpPolyhedron(0), warningFlag(0) |
---|
| 95 | { |
---|
| 96 | // fV<x><y> is vector from vertex <y> to vertex <x> |
---|
| 97 | // |
---|
| 98 | G4ThreeVector fV21=p2-anchor; |
---|
| 99 | G4ThreeVector fV31=p3-anchor; |
---|
| 100 | G4ThreeVector fV41=p4-anchor; |
---|
| 101 | |
---|
| 102 | // make sure this is a correctly oriented set of points for the tetrahedron |
---|
| 103 | // |
---|
| 104 | G4double signed_vol=fV21.cross(fV31).dot(fV41); |
---|
| 105 | if(signed_vol<0.0) |
---|
| 106 | { |
---|
| 107 | G4ThreeVector temp(p4); |
---|
| 108 | p4=p3; |
---|
| 109 | p3=temp; |
---|
| 110 | temp=fV41; |
---|
| 111 | fV41=fV31; |
---|
| 112 | fV31=temp; |
---|
| 113 | } |
---|
| 114 | fCubicVolume = std::abs(signed_vol) / 6.; |
---|
| 115 | |
---|
| 116 | G4ThreeVector fV24=p2-p4; |
---|
| 117 | G4ThreeVector fV43=p4-p3; |
---|
| 118 | G4ThreeVector fV32=p3-p2; |
---|
| 119 | |
---|
| 120 | fXMin=std::min(std::min(std::min(anchor.x(), p2.x()),p3.x()),p4.x()); |
---|
| 121 | fXMax=std::max(std::max(std::max(anchor.x(), p2.x()),p3.x()),p4.x()); |
---|
| 122 | fYMin=std::min(std::min(std::min(anchor.y(), p2.y()),p3.y()),p4.y()); |
---|
| 123 | fYMax=std::max(std::max(std::max(anchor.y(), p2.y()),p3.y()),p4.y()); |
---|
| 124 | fZMin=std::min(std::min(std::min(anchor.z(), p2.z()),p3.z()),p4.z()); |
---|
| 125 | fZMax=std::max(std::max(std::max(anchor.z(), p2.z()),p3.z()),p4.z()); |
---|
| 126 | |
---|
| 127 | fDx=(fXMax-fXMin)*0.5; fDy=(fYMax-fYMin)*0.5; fDz=(fZMax-fZMin)*0.5; |
---|
| 128 | |
---|
| 129 | fMiddle=G4ThreeVector(fXMax+fXMin, fYMax+fYMin, fZMax+fZMin)*0.5; |
---|
| 130 | fMaxSize=std::max(std::max(std::max((anchor-fMiddle).mag(), |
---|
| 131 | (p2-fMiddle).mag()), |
---|
| 132 | (p3-fMiddle).mag()), |
---|
| 133 | (p4-fMiddle).mag()); |
---|
| 134 | |
---|
| 135 | G4bool degenerate=std::abs(signed_vol) < 1e-9*fMaxSize*fMaxSize*fMaxSize; |
---|
| 136 | |
---|
| 137 | if(degeneracyFlag) *degeneracyFlag=degenerate; |
---|
| 138 | else if (degenerate) |
---|
| 139 | { |
---|
| 140 | G4Exception("G4Tet::G4Tet()", "InvalidSetup", FatalException, |
---|
| 141 | "Degenerate tetrahedron not allowed."); |
---|
| 142 | } |
---|
| 143 | |
---|
| 144 | fTol=1e-9*(std::abs(fXMin)+std::abs(fXMax)+std::abs(fYMin) |
---|
| 145 | +std::abs(fYMax)+std::abs(fZMin)+std::abs(fZMax)); |
---|
| 146 | //fTol=kCarTolerance; |
---|
| 147 | |
---|
| 148 | fAnchor=anchor; |
---|
| 149 | fP2=p2; |
---|
| 150 | fP3=p3; |
---|
| 151 | fP4=p4; |
---|
| 152 | |
---|
| 153 | G4ThreeVector fCenter123=(anchor+p2+p3)*(1.0/3.0); // face center |
---|
| 154 | G4ThreeVector fCenter134=(anchor+p4+p3)*(1.0/3.0); |
---|
| 155 | G4ThreeVector fCenter142=(anchor+p4+p2)*(1.0/3.0); |
---|
| 156 | G4ThreeVector fCenter234=(p2+p3+p4)*(1.0/3.0); |
---|
| 157 | |
---|
| 158 | // compute area of each triangular face by cross product |
---|
| 159 | // and sum for total surface area |
---|
| 160 | |
---|
| 161 | G4ThreeVector normal123=fV31.cross(fV21); |
---|
| 162 | G4ThreeVector normal134=fV41.cross(fV31); |
---|
| 163 | G4ThreeVector normal142=fV21.cross(fV41); |
---|
| 164 | G4ThreeVector normal234=fV32.cross(fV43); |
---|
| 165 | |
---|
| 166 | fSurfaceArea=( |
---|
| 167 | normal123.mag()+ |
---|
| 168 | normal134.mag()+ |
---|
| 169 | normal142.mag()+ |
---|
| 170 | normal234.mag() |
---|
| 171 | )/2.0; |
---|
| 172 | |
---|
| 173 | fNormal123=normal123.unit(); |
---|
| 174 | fNormal134=normal134.unit(); |
---|
| 175 | fNormal142=normal142.unit(); |
---|
| 176 | fNormal234=normal234.unit(); |
---|
| 177 | |
---|
| 178 | fCdotN123=fCenter123.dot(fNormal123); |
---|
| 179 | fCdotN134=fCenter134.dot(fNormal134); |
---|
| 180 | fCdotN142=fCenter142.dot(fNormal142); |
---|
| 181 | fCdotN234=fCenter234.dot(fNormal234); |
---|
| 182 | } |
---|
| 183 | |
---|
| 184 | ////////////////////////////////////////////////////////////////////////// |
---|
| 185 | // |
---|
| 186 | // Fake default constructor - sets only member data and allocates memory |
---|
| 187 | // for usage restricted to object persistency. |
---|
| 188 | // |
---|
| 189 | G4Tet::G4Tet( __void__& a ) |
---|
| 190 | : G4VSolid(a), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0), |
---|
| 191 | fAnchor(0,0,0), fP2(0,0,0), fP3(0,0,0), fP4(0,0,0), fMiddle(0,0,0), |
---|
| 192 | fNormal123(0,0,0), fNormal142(0,0,0), fNormal134(0,0,0), |
---|
| 193 | fNormal234(0,0,0), warningFlag(0), |
---|
| 194 | fCdotN123(0.), fCdotN142(0.), fCdotN134(0.), fCdotN234(0.), |
---|
| 195 | fXMin(0.), fXMax(0.), fYMin(0.), fYMax(0.), fZMin(0.), fZMax(0.), |
---|
| 196 | fDx(0.), fDy(0.), fDz(0.), fTol(0.), fMaxSize(0.) |
---|
| 197 | { |
---|
| 198 | } |
---|
| 199 | |
---|
| 200 | ////////////////////////////////////////////////////////////////////////// |
---|
| 201 | // |
---|
| 202 | // Destructor |
---|
| 203 | |
---|
| 204 | G4Tet::~G4Tet() |
---|
| 205 | { |
---|
| 206 | delete fpPolyhedron; |
---|
| 207 | } |
---|
| 208 | |
---|
| 209 | ////////////////////////////////////////////////////////////////////////// |
---|
| 210 | // |
---|
| 211 | // CheckDegeneracy |
---|
| 212 | |
---|
| 213 | G4bool G4Tet::CheckDegeneracy( G4ThreeVector anchor, |
---|
| 214 | G4ThreeVector p2, |
---|
| 215 | G4ThreeVector p3, |
---|
| 216 | G4ThreeVector p4 ) |
---|
| 217 | { |
---|
| 218 | G4bool result; |
---|
| 219 | G4Tet *object=new G4Tet("temp",anchor,p2,p3,p4,&result); |
---|
| 220 | delete object; |
---|
| 221 | return result; |
---|
| 222 | } |
---|
| 223 | |
---|
| 224 | ////////////////////////////////////////////////////////////////////////// |
---|
| 225 | // |
---|
| 226 | // Dispatch to parameterisation for replication mechanism dimension |
---|
| 227 | // computation & modification. |
---|
| 228 | |
---|
| 229 | void G4Tet::ComputeDimensions(G4VPVParameterisation* , |
---|
| 230 | const G4int , |
---|
| 231 | const G4VPhysicalVolume* ) |
---|
| 232 | { |
---|
| 233 | } |
---|
| 234 | |
---|
| 235 | ////////////////////////////////////////////////////////////////////////// |
---|
| 236 | // |
---|
| 237 | // Calculate extent under transform and specified limit |
---|
| 238 | |
---|
| 239 | G4bool G4Tet::CalculateExtent(const EAxis pAxis, |
---|
| 240 | const G4VoxelLimits& pVoxelLimit, |
---|
| 241 | const G4AffineTransform& pTransform, |
---|
| 242 | G4double& pMin, G4double& pMax) const |
---|
| 243 | { |
---|
| 244 | G4double xMin,xMax; |
---|
| 245 | G4double yMin,yMax; |
---|
| 246 | G4double zMin,zMax; |
---|
| 247 | |
---|
| 248 | if (pTransform.IsRotated()) |
---|
| 249 | { |
---|
| 250 | G4ThreeVector pp0=pTransform.TransformPoint(fAnchor); |
---|
| 251 | G4ThreeVector pp1=pTransform.TransformPoint(fP2); |
---|
| 252 | G4ThreeVector pp2=pTransform.TransformPoint(fP3); |
---|
| 253 | G4ThreeVector pp3=pTransform.TransformPoint(fP4); |
---|
| 254 | |
---|
| 255 | xMin = std::min(std::min(std::min(pp0.x(), pp1.x()),pp2.x()),pp3.x()); |
---|
| 256 | xMax = std::max(std::max(std::max(pp0.x(), pp1.x()),pp2.x()),pp3.x()); |
---|
| 257 | yMin = std::min(std::min(std::min(pp0.y(), pp1.y()),pp2.y()),pp3.y()); |
---|
| 258 | yMax = std::max(std::max(std::max(pp0.y(), pp1.y()),pp2.y()),pp3.y()); |
---|
| 259 | zMin = std::min(std::min(std::min(pp0.z(), pp1.z()),pp2.z()),pp3.z()); |
---|
| 260 | zMax = std::max(std::max(std::max(pp0.z(), pp1.z()),pp2.z()),pp3.z()); |
---|
| 261 | |
---|
| 262 | } |
---|
| 263 | else |
---|
| 264 | { |
---|
| 265 | G4double xoffset = pTransform.NetTranslation().x() ; |
---|
| 266 | xMin = xoffset + fXMin; |
---|
| 267 | xMax = xoffset + fXMax; |
---|
| 268 | G4double yoffset = pTransform.NetTranslation().y() ; |
---|
| 269 | yMin = yoffset + fYMin; |
---|
| 270 | yMax = yoffset + fYMax; |
---|
| 271 | G4double zoffset = pTransform.NetTranslation().z() ; |
---|
| 272 | zMin = zoffset + fZMin; |
---|
| 273 | zMax = zoffset + fZMax; |
---|
| 274 | } |
---|
| 275 | |
---|
| 276 | if (pVoxelLimit.IsXLimited()) |
---|
| 277 | { |
---|
| 278 | if ( (xMin > pVoxelLimit.GetMaxXExtent()+fTol) || |
---|
| 279 | (xMax < pVoxelLimit.GetMinXExtent()-fTol) ) { return false; } |
---|
| 280 | else |
---|
| 281 | { |
---|
| 282 | xMin = std::max(xMin, pVoxelLimit.GetMinXExtent()); |
---|
| 283 | xMax = std::min(xMax, pVoxelLimit.GetMaxXExtent()); |
---|
| 284 | } |
---|
| 285 | } |
---|
| 286 | |
---|
| 287 | if (pVoxelLimit.IsYLimited()) |
---|
| 288 | { |
---|
| 289 | if ( (yMin > pVoxelLimit.GetMaxYExtent()+fTol) || |
---|
| 290 | (yMax < pVoxelLimit.GetMinYExtent()-fTol) ) { return false; } |
---|
| 291 | else |
---|
| 292 | { |
---|
| 293 | yMin = std::max(yMin, pVoxelLimit.GetMinYExtent()); |
---|
| 294 | yMax = std::min(yMax, pVoxelLimit.GetMaxYExtent()); |
---|
| 295 | } |
---|
| 296 | } |
---|
| 297 | |
---|
| 298 | if (pVoxelLimit.IsZLimited()) |
---|
| 299 | { |
---|
| 300 | if ( (zMin > pVoxelLimit.GetMaxZExtent()+fTol) || |
---|
| 301 | (zMax < pVoxelLimit.GetMinZExtent()-fTol) ) { return false; } |
---|
| 302 | else |
---|
| 303 | { |
---|
| 304 | zMin = std::max(zMin, pVoxelLimit.GetMinZExtent()); |
---|
| 305 | zMax = std::min(zMax, pVoxelLimit.GetMaxZExtent()); |
---|
| 306 | } |
---|
| 307 | } |
---|
| 308 | |
---|
| 309 | switch (pAxis) |
---|
| 310 | { |
---|
| 311 | case kXAxis: |
---|
| 312 | pMin=xMin; |
---|
| 313 | pMax=xMax; |
---|
| 314 | break; |
---|
| 315 | case kYAxis: |
---|
| 316 | pMin=yMin; |
---|
| 317 | pMax=yMax; |
---|
| 318 | break; |
---|
| 319 | case kZAxis: |
---|
| 320 | pMin=zMin; |
---|
| 321 | pMax=zMax; |
---|
| 322 | break; |
---|
| 323 | default: |
---|
| 324 | break; |
---|
| 325 | } |
---|
| 326 | |
---|
| 327 | return true; |
---|
| 328 | } |
---|
| 329 | |
---|
| 330 | ///////////////////////////////////////////////////////////////////////// |
---|
| 331 | // |
---|
| 332 | // Return whether point inside/outside/on surface, using tolerance |
---|
| 333 | |
---|
| 334 | EInside G4Tet::Inside(const G4ThreeVector& p) const |
---|
| 335 | { |
---|
| 336 | G4double r123, r134, r142, r234; |
---|
| 337 | |
---|
| 338 | // this is written to allow if-statement truncation so the outside test |
---|
| 339 | // (where most of the world is) can fail very quickly and efficiently |
---|
| 340 | |
---|
| 341 | if ( (r123=p.dot(fNormal123)-fCdotN123) > fTol || |
---|
| 342 | (r134=p.dot(fNormal134)-fCdotN134) > fTol || |
---|
| 343 | (r142=p.dot(fNormal142)-fCdotN142) > fTol || |
---|
| 344 | (r234=p.dot(fNormal234)-fCdotN234) > fTol ) |
---|
| 345 | { |
---|
| 346 | return kOutside; // at least one is out! |
---|
| 347 | } |
---|
| 348 | else if( (r123 < -fTol)&&(r134 < -fTol)&&(r142 < -fTol)&&(r234 < -fTol) ) |
---|
| 349 | { |
---|
| 350 | return kInside; // all are definitively inside |
---|
| 351 | } |
---|
| 352 | else |
---|
| 353 | { |
---|
| 354 | return kSurface; // too close to tell |
---|
| 355 | } |
---|
| 356 | } |
---|
| 357 | |
---|
| 358 | /////////////////////////////////////////////////////////////////////// |
---|
| 359 | // |
---|
| 360 | // Calculate side nearest to p, and return normal |
---|
| 361 | // If two sides are equidistant, normal of first side (x/y/z) |
---|
| 362 | // encountered returned. |
---|
| 363 | // This assumes that we are looking from the inside! |
---|
| 364 | |
---|
| 365 | G4ThreeVector G4Tet::SurfaceNormal( const G4ThreeVector& p) const |
---|
| 366 | { |
---|
| 367 | G4double r123=std::abs(p.dot(fNormal123)-fCdotN123); |
---|
| 368 | G4double r134=std::abs(p.dot(fNormal134)-fCdotN134); |
---|
| 369 | G4double r142=std::abs(p.dot(fNormal142)-fCdotN142); |
---|
| 370 | G4double r234=std::abs(p.dot(fNormal234)-fCdotN234); |
---|
| 371 | |
---|
| 372 | if( (r123<=r134) && (r123<=r142) && (r123<=r234) ) { return fNormal123; } |
---|
| 373 | else if ( (r134<=r142) && (r134<=r234) ) { return fNormal134; } |
---|
| 374 | else if (r142 <= r234) { return fNormal142; } |
---|
| 375 | return fNormal234; |
---|
| 376 | } |
---|
| 377 | |
---|
| 378 | /////////////////////////////////////////////////////////////////////////// |
---|
| 379 | // |
---|
| 380 | // Calculate distance to box from an outside point |
---|
| 381 | // - return kInfinity if no intersection. |
---|
| 382 | // All this is very unrolled, for speed. |
---|
| 383 | |
---|
| 384 | G4double G4Tet::DistanceToIn(const G4ThreeVector& p, |
---|
| 385 | const G4ThreeVector& v) const |
---|
| 386 | { |
---|
| 387 | G4ThreeVector vu(v.unit()), hp; |
---|
| 388 | G4double vdotn, t, tmin=kInfinity; |
---|
| 389 | |
---|
| 390 | G4double extraDistance=10.0*fTol; // a little ways into the solid |
---|
| 391 | |
---|
| 392 | vdotn=-vu.dot(fNormal123); |
---|
| 393 | if(vdotn > 1e-12) |
---|
| 394 | { // this is a candidate face, since it is pointing at us |
---|
| 395 | t=(p.dot(fNormal123)-fCdotN123)/vdotn; // # distance to intersection |
---|
| 396 | if( (t>=-fTol) && (t<tmin) ) |
---|
| 397 | { // if not true, we're going away from this face or it's not close |
---|
| 398 | hp=p+vu*(t+extraDistance); // a little beyond point of intersection |
---|
| 399 | if ( ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) && |
---|
| 400 | ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) && |
---|
| 401 | ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) ) |
---|
| 402 | { |
---|
| 403 | tmin=t; |
---|
| 404 | } |
---|
| 405 | } |
---|
| 406 | } |
---|
| 407 | |
---|
| 408 | vdotn=-vu.dot(fNormal134); |
---|
| 409 | if(vdotn > 1e-12) |
---|
| 410 | { // # this is a candidate face, since it is pointing at us |
---|
| 411 | t=(p.dot(fNormal134)-fCdotN134)/vdotn; // # distance to intersection |
---|
| 412 | if( (t>=-fTol) && (t<tmin) ) |
---|
| 413 | { // if not true, we're going away from this face |
---|
| 414 | hp=p+vu*(t+extraDistance); // a little beyond point of intersection |
---|
| 415 | if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) && |
---|
| 416 | ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) && |
---|
| 417 | ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) ) |
---|
| 418 | { |
---|
| 419 | tmin=t; |
---|
| 420 | } |
---|
| 421 | } |
---|
| 422 | } |
---|
| 423 | |
---|
| 424 | vdotn=-vu.dot(fNormal142); |
---|
| 425 | if(vdotn > 1e-12) |
---|
| 426 | { // # this is a candidate face, since it is pointing at us |
---|
| 427 | t=(p.dot(fNormal142)-fCdotN142)/vdotn; // # distance to intersection |
---|
| 428 | if( (t>=-fTol) && (t<tmin) ) |
---|
| 429 | { // if not true, we're going away from this face |
---|
| 430 | hp=p+vu*(t+extraDistance); // a little beyond point of intersection |
---|
| 431 | if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) && |
---|
| 432 | ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) && |
---|
| 433 | ( hp.dot(fNormal234)-fCdotN234 < 0.0 ) ) |
---|
| 434 | { |
---|
| 435 | tmin=t; |
---|
| 436 | } |
---|
| 437 | } |
---|
| 438 | } |
---|
| 439 | |
---|
| 440 | vdotn=-vu.dot(fNormal234); |
---|
| 441 | if(vdotn > 1e-12) |
---|
| 442 | { // # this is a candidate face, since it is pointing at us |
---|
| 443 | t=(p.dot(fNormal234)-fCdotN234)/vdotn; // # distance to intersection |
---|
| 444 | if( (t>=-fTol) && (t<tmin) ) |
---|
| 445 | { // if not true, we're going away from this face |
---|
| 446 | hp=p+vu*(t+extraDistance); // a little beyond point of intersection |
---|
| 447 | if ( ( hp.dot(fNormal123)-fCdotN123 < 0.0 ) && |
---|
| 448 | ( hp.dot(fNormal134)-fCdotN134 < 0.0 ) && |
---|
| 449 | ( hp.dot(fNormal142)-fCdotN142 < 0.0 ) ) |
---|
| 450 | { |
---|
| 451 | tmin=t; |
---|
| 452 | } |
---|
| 453 | } |
---|
| 454 | } |
---|
| 455 | |
---|
| 456 | return std::max(0.0,tmin); |
---|
| 457 | } |
---|
| 458 | |
---|
| 459 | ////////////////////////////////////////////////////////////////////////// |
---|
| 460 | // |
---|
| 461 | // Approximate distance to tet. |
---|
| 462 | // returns distance to sphere centered on bounding box |
---|
| 463 | // - If inside return 0 |
---|
| 464 | |
---|
| 465 | G4double G4Tet::DistanceToIn(const G4ThreeVector& p) const |
---|
| 466 | { |
---|
| 467 | G4double dd=(p-fMiddle).mag() - fMaxSize - fTol; |
---|
| 468 | return std::max(0.0, dd); |
---|
| 469 | } |
---|
| 470 | |
---|
| 471 | ///////////////////////////////////////////////////////////////////////// |
---|
| 472 | // |
---|
| 473 | // Calcluate distance to surface of box from inside |
---|
| 474 | // by calculating distances to box's x/y/z planes. |
---|
| 475 | // Smallest distance is exact distance to exiting. |
---|
| 476 | |
---|
| 477 | G4double G4Tet::DistanceToOut( const G4ThreeVector& p,const G4ThreeVector& v, |
---|
| 478 | const G4bool calcNorm, |
---|
| 479 | G4bool *validNorm, G4ThreeVector *n) const |
---|
| 480 | { |
---|
| 481 | G4ThreeVector vu(v.unit()); |
---|
| 482 | G4double t1=kInfinity,t2=kInfinity,t3=kInfinity,t4=kInfinity, vdotn, tt; |
---|
| 483 | |
---|
| 484 | vdotn=vu.dot(fNormal123); |
---|
| 485 | if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate |
---|
| 486 | { |
---|
| 487 | t1=(fCdotN123-p.dot(fNormal123))/vdotn; // # distance to intersection |
---|
| 488 | } |
---|
| 489 | |
---|
| 490 | vdotn=vu.dot(fNormal134); |
---|
| 491 | if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate |
---|
| 492 | { |
---|
| 493 | t2=(fCdotN134-p.dot(fNormal134))/vdotn; // # distance to intersection |
---|
| 494 | } |
---|
| 495 | |
---|
| 496 | vdotn=vu.dot(fNormal142); |
---|
| 497 | if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate |
---|
| 498 | { |
---|
| 499 | t3=(fCdotN142-p.dot(fNormal142))/vdotn; // # distance to intersection |
---|
| 500 | } |
---|
| 501 | |
---|
| 502 | vdotn=vu.dot(fNormal234); |
---|
| 503 | if(vdotn > 1e-12) // #we're heading towards this face, so it is a candidate |
---|
| 504 | { |
---|
| 505 | t4=(fCdotN234-p.dot(fNormal234))/vdotn; // # distance to intersection |
---|
| 506 | } |
---|
| 507 | |
---|
| 508 | tt=std::min(std::min(std::min(t1,t2),t3),t4); |
---|
| 509 | |
---|
| 510 | if (warningFlag && (tt == kInfinity || tt < -fTol)) |
---|
| 511 | { |
---|
| 512 | DumpInfo(); |
---|
| 513 | G4cout << "p = " << p / mm << "mm" << G4endl; |
---|
| 514 | G4cout << "v = " << v << G4endl; |
---|
| 515 | G4cout << "t1, t2, t3, t4 (mm) " |
---|
| 516 | << t1/mm << ", " << t2/mm << ", " << t3/mm << ", " << t4/mm |
---|
| 517 | << G4endl << G4endl; |
---|
| 518 | G4Exception("G4Tet::DistanceToOut(p,v,...)", "Notification", JustWarning, |
---|
| 519 | "No good intersection found or already outside!?" ); |
---|
| 520 | if(validNorm) |
---|
| 521 | { |
---|
| 522 | *validNorm=false; // flag normal as meaningless |
---|
| 523 | } |
---|
| 524 | } |
---|
| 525 | else if(calcNorm && n) |
---|
| 526 | { |
---|
| 527 | static G4ThreeVector normal; |
---|
| 528 | if(tt==t1) { normal=fNormal123; } |
---|
| 529 | else if (tt==t2) { normal=fNormal134; } |
---|
| 530 | else if (tt==t3) { normal=fNormal142; } |
---|
| 531 | else if (tt==t4) { normal=fNormal234; } |
---|
| 532 | n=&normal; |
---|
| 533 | if(validNorm) { *validNorm=true; } |
---|
| 534 | } |
---|
| 535 | |
---|
| 536 | return std::max(tt,0.0); // avoid tt<0.0 by a tiny bit |
---|
| 537 | // if we are right on a face |
---|
| 538 | } |
---|
| 539 | |
---|
| 540 | //////////////////////////////////////////////////////////////////////////// |
---|
| 541 | // |
---|
| 542 | // Calculate exact shortest distance to any boundary from inside |
---|
| 543 | // - If outside return 0 |
---|
| 544 | |
---|
| 545 | G4double G4Tet::DistanceToOut(const G4ThreeVector& p) const |
---|
| 546 | { |
---|
| 547 | G4double t1,t2,t3,t4; |
---|
| 548 | t1=fCdotN123-p.dot(fNormal123); // distance to plane, positive if inside |
---|
| 549 | t2=fCdotN134-p.dot(fNormal134); // distance to plane |
---|
| 550 | t3=fCdotN142-p.dot(fNormal142); // distance to plane |
---|
| 551 | t4=fCdotN234-p.dot(fNormal234); // distance to plane |
---|
| 552 | |
---|
| 553 | // if any one of these is negative, we are outside, |
---|
| 554 | // so return zero in that case |
---|
| 555 | |
---|
| 556 | G4double tmin=std::min(std::min(std::min(t1,t2),t3),t4); |
---|
| 557 | return (tmin < fTol)? 0:tmin; |
---|
| 558 | } |
---|
| 559 | |
---|
| 560 | //////////////////////////////////////////////////////////////////////// |
---|
| 561 | // |
---|
| 562 | // Create a List containing the transformed vertices |
---|
| 563 | // Note: Caller has deletion responsibility |
---|
| 564 | |
---|
| 565 | G4ThreeVectorList* |
---|
| 566 | G4Tet::CreateRotatedVertices(const G4AffineTransform& pTransform) const |
---|
| 567 | { |
---|
| 568 | G4ThreeVectorList* vertices = new G4ThreeVectorList(); |
---|
| 569 | vertices->reserve(4); |
---|
| 570 | |
---|
| 571 | if (vertices) |
---|
| 572 | { |
---|
| 573 | G4ThreeVector vertex0(fAnchor); |
---|
| 574 | G4ThreeVector vertex1(fP2); |
---|
| 575 | G4ThreeVector vertex2(fP3); |
---|
| 576 | G4ThreeVector vertex3(fP4); |
---|
| 577 | |
---|
| 578 | vertices->push_back(pTransform.TransformPoint(vertex0)); |
---|
| 579 | vertices->push_back(pTransform.TransformPoint(vertex1)); |
---|
| 580 | vertices->push_back(pTransform.TransformPoint(vertex2)); |
---|
| 581 | vertices->push_back(pTransform.TransformPoint(vertex3)); |
---|
| 582 | } |
---|
| 583 | else |
---|
| 584 | { |
---|
| 585 | DumpInfo(); |
---|
| 586 | G4Exception("G4Tet::CreateRotatedVertices()", |
---|
| 587 | "FatalError", FatalException, |
---|
| 588 | "Error in allocation of vertices. Out of memory !"); |
---|
| 589 | } |
---|
| 590 | return vertices; |
---|
| 591 | } |
---|
| 592 | |
---|
| 593 | ////////////////////////////////////////////////////////////////////////// |
---|
| 594 | // |
---|
| 595 | // GetEntityType |
---|
| 596 | |
---|
| 597 | G4GeometryType G4Tet::GetEntityType() const |
---|
| 598 | { |
---|
| 599 | return G4String("G4Tet"); |
---|
| 600 | } |
---|
| 601 | |
---|
| 602 | ////////////////////////////////////////////////////////////////////////// |
---|
| 603 | // |
---|
| 604 | // Stream object contents to an output stream |
---|
| 605 | |
---|
| 606 | std::ostream& G4Tet::StreamInfo(std::ostream& os) const |
---|
| 607 | { |
---|
| 608 | os << "-----------------------------------------------------------\n" |
---|
| 609 | << " *** Dump for solid - " << GetName() << " ***\n" |
---|
| 610 | << " ===================================================\n" |
---|
| 611 | << " Solid type: G4Tet\n" |
---|
| 612 | << " Parameters: \n" |
---|
| 613 | << " anchor: " << fAnchor/mm << " mm \n" |
---|
| 614 | << " p2: " << fP2/mm << " mm \n" |
---|
| 615 | << " p3: " << fP3/mm << " mm \n" |
---|
| 616 | << " p4: " << fP4/mm << " mm \n" |
---|
| 617 | << " normal123: " << fNormal123 << " \n" |
---|
| 618 | << " normal134: " << fNormal134 << " \n" |
---|
| 619 | << " normal142: " << fNormal142 << " \n" |
---|
| 620 | << " normal234: " << fNormal234 << " \n" |
---|
| 621 | << "-----------------------------------------------------------\n"; |
---|
| 622 | |
---|
| 623 | return os; |
---|
| 624 | } |
---|
| 625 | |
---|
| 626 | |
---|
| 627 | //////////////////////////////////////////////////////////////////////// |
---|
| 628 | // |
---|
| 629 | // GetPointOnFace |
---|
| 630 | // |
---|
| 631 | // Auxiliary method for get point on surface |
---|
| 632 | |
---|
| 633 | G4ThreeVector G4Tet::GetPointOnFace(G4ThreeVector p1, G4ThreeVector p2, |
---|
| 634 | G4ThreeVector p3, G4double& area) const |
---|
| 635 | { |
---|
| 636 | G4double lambda1,lambda2; |
---|
| 637 | G4ThreeVector v, w; |
---|
| 638 | |
---|
| 639 | v = p3 - p1; |
---|
| 640 | w = p1 - p2; |
---|
| 641 | |
---|
| 642 | lambda1 = RandFlat::shoot(0.,1.); |
---|
| 643 | lambda2 = RandFlat::shoot(0.,lambda1); |
---|
| 644 | |
---|
| 645 | area = 0.5*(v.cross(w)).mag(); |
---|
| 646 | |
---|
| 647 | return (p2 + lambda1*w + lambda2*v); |
---|
| 648 | } |
---|
| 649 | |
---|
| 650 | //////////////////////////////////////////////////////////////////////////// |
---|
| 651 | // |
---|
| 652 | // GetPointOnSurface |
---|
| 653 | |
---|
| 654 | G4ThreeVector G4Tet::GetPointOnSurface() const |
---|
| 655 | { |
---|
| 656 | G4double chose,aOne,aTwo,aThree,aFour; |
---|
| 657 | G4ThreeVector p1, p2, p3, p4; |
---|
| 658 | |
---|
| 659 | p1 = GetPointOnFace(fAnchor,fP2,fP3,aOne); |
---|
| 660 | p2 = GetPointOnFace(fAnchor,fP4,fP3,aTwo); |
---|
| 661 | p3 = GetPointOnFace(fAnchor,fP4,fP2,aThree); |
---|
| 662 | p4 = GetPointOnFace(fP4,fP3,fP2,aFour); |
---|
| 663 | |
---|
| 664 | chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour); |
---|
| 665 | if( (chose>=0.) && (chose <aOne) ) {return p1;} |
---|
| 666 | else if( (chose>=aOne) && (chose < aOne+aTwo) ) {return p2;} |
---|
| 667 | else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) ) {return p3;} |
---|
| 668 | return p4; |
---|
| 669 | } |
---|
| 670 | |
---|
| 671 | //////////////////////////////////////////////////////////////////////// |
---|
| 672 | // |
---|
| 673 | // GetVertices |
---|
| 674 | |
---|
| 675 | std::vector<G4ThreeVector> G4Tet::GetVertices() const |
---|
| 676 | { |
---|
| 677 | std::vector<G4ThreeVector> vertices(4); |
---|
| 678 | vertices[0] = fAnchor; |
---|
| 679 | vertices[1] = fP2; |
---|
| 680 | vertices[2] = fP3; |
---|
| 681 | vertices[3] = fP4; |
---|
| 682 | |
---|
| 683 | return vertices; |
---|
| 684 | } |
---|
| 685 | |
---|
| 686 | //////////////////////////////////////////////////////////////////////// |
---|
| 687 | // |
---|
| 688 | // GetCubicVolume |
---|
| 689 | |
---|
| 690 | G4double G4Tet::GetCubicVolume() |
---|
| 691 | { |
---|
| 692 | return fCubicVolume; |
---|
| 693 | } |
---|
| 694 | |
---|
| 695 | //////////////////////////////////////////////////////////////////////// |
---|
| 696 | // |
---|
| 697 | // GetSurfaceArea |
---|
| 698 | |
---|
| 699 | G4double G4Tet::GetSurfaceArea() |
---|
| 700 | { |
---|
| 701 | return fSurfaceArea; |
---|
| 702 | } |
---|
| 703 | |
---|
| 704 | // Methods for visualisation |
---|
| 705 | |
---|
| 706 | //////////////////////////////////////////////////////////////////////// |
---|
| 707 | // |
---|
| 708 | // DescribeYourselfTo |
---|
| 709 | |
---|
| 710 | void G4Tet::DescribeYourselfTo (G4VGraphicsScene& scene) const |
---|
| 711 | { |
---|
| 712 | scene.AddSolid (*this); |
---|
| 713 | } |
---|
| 714 | |
---|
| 715 | //////////////////////////////////////////////////////////////////////// |
---|
| 716 | // |
---|
| 717 | // GetExtent |
---|
| 718 | |
---|
| 719 | G4VisExtent G4Tet::GetExtent() const |
---|
| 720 | { |
---|
| 721 | return G4VisExtent (fXMin, fXMax, fYMin, fYMax, fZMin, fZMax); |
---|
| 722 | } |
---|
| 723 | |
---|
| 724 | //////////////////////////////////////////////////////////////////////// |
---|
| 725 | // |
---|
| 726 | // CreatePolyhedron |
---|
| 727 | |
---|
| 728 | G4Polyhedron* G4Tet::CreatePolyhedron () const |
---|
| 729 | { |
---|
| 730 | G4Polyhedron *ph=new G4Polyhedron; |
---|
| 731 | G4double xyz[4][3]; |
---|
| 732 | static G4int faces[4][4]={{1,3,2,0},{1,4,3,0},{1,2,4,0},{2,3,4,0}}; |
---|
| 733 | xyz[0][0]=fAnchor.x(); xyz[0][1]=fAnchor.y(); xyz[0][2]=fAnchor.z(); |
---|
| 734 | xyz[1][0]=fP2.x(); xyz[1][1]=fP2.y(); xyz[1][2]=fP2.z(); |
---|
| 735 | xyz[2][0]=fP3.x(); xyz[2][1]=fP3.y(); xyz[2][2]=fP3.z(); |
---|
| 736 | xyz[3][0]=fP4.x(); xyz[3][1]=fP4.y(); xyz[3][2]=fP4.z(); |
---|
| 737 | |
---|
| 738 | ph->createPolyhedron(4,4,xyz,faces); |
---|
| 739 | |
---|
| 740 | return ph; |
---|
| 741 | } |
---|
| 742 | |
---|
| 743 | //////////////////////////////////////////////////////////////////////// |
---|
| 744 | // |
---|
| 745 | // CreateNURBS |
---|
| 746 | |
---|
| 747 | G4NURBS* G4Tet::CreateNURBS () const |
---|
| 748 | { |
---|
| 749 | return new G4NURBSbox (fDx, fDy, fDz); |
---|
| 750 | } |
---|
| 751 | |
---|
| 752 | //////////////////////////////////////////////////////////////////////// |
---|
| 753 | // |
---|
| 754 | // GetPolyhedron |
---|
| 755 | |
---|
| 756 | G4Polyhedron* G4Tet::GetPolyhedron () const |
---|
| 757 | { |
---|
| 758 | if (!fpPolyhedron || |
---|
| 759 | fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != |
---|
| 760 | fpPolyhedron->GetNumberOfRotationSteps()) |
---|
| 761 | { |
---|
| 762 | delete fpPolyhedron; |
---|
| 763 | fpPolyhedron = CreatePolyhedron(); |
---|
| 764 | } |
---|
| 765 | return fpPolyhedron; |
---|
| 766 | } |
---|