[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 result of the scientific and * |
---|
| 19 | // * technical work of the GEANT4 collaboration and of QinetiQ Ltd, * |
---|
| 20 | // * subject to DEFCON 705 IPR conditions. * |
---|
| 21 | // * By using, copying, modifying or distributing the software (or * |
---|
| 22 | // * any work based on the software) you agree to acknowledge its * |
---|
| 23 | // * use in resulting scientific publications, and indicate your * |
---|
| 24 | // * acceptance of all terms of the Geant4 Software license. * |
---|
| 25 | // ******************************************************************** |
---|
| 26 | // |
---|
[1228] | 27 | // $Id: G4TessellatedSolid.cc,v 1.19 2009/04/27 08:06:27 gcosmo Exp $ |
---|
| 28 | // GEANT4 tag $Name: geant4-09-03 $ |
---|
[831] | 29 | // |
---|
| 30 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 31 | // |
---|
| 32 | // MODULE: G4TessellatedSolid.cc |
---|
| 33 | // |
---|
| 34 | // Date: 15/06/2005 |
---|
| 35 | // Author: P R Truscott |
---|
| 36 | // Organisation: QinetiQ Ltd, UK |
---|
| 37 | // Customer: UK Ministry of Defence : RAO CRP TD Electronic Systems |
---|
| 38 | // Contract: C/MAT/N03517 |
---|
| 39 | // |
---|
| 40 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 41 | // |
---|
| 42 | // CHANGE HISTORY |
---|
| 43 | // -------------- |
---|
| 44 | // |
---|
| 45 | // 14 November 2007 P R Truscott, QinetiQ & Stan Seibert, U Texas |
---|
| 46 | // Bug fixes to CalculateExtent |
---|
| 47 | // |
---|
| 48 | // 17 September 2007, P R Truscott, QinetiQ Ltd & Richard Holmberg |
---|
| 49 | // Updated extensively prior to this date to deal with |
---|
| 50 | // concaved tessellated surfaces, based on the algorithm |
---|
| 51 | // of Richard Holmberg. This had been slightly modified |
---|
| 52 | // to determine with inside the geometry by projecting |
---|
| 53 | // random rays from the point provided. Now random rays |
---|
| 54 | // are predefined rather than making use of random |
---|
| 55 | // number generator at run-time. |
---|
| 56 | // |
---|
| 57 | // 22 November 2005, F Lei |
---|
| 58 | // - Changed ::DescribeYourselfTo(), line 464 |
---|
| 59 | // - added GetPolyHedron() |
---|
| 60 | // |
---|
| 61 | // 31 October 2004, P R Truscott, QinetiQ Ltd, UK |
---|
| 62 | // - Created. |
---|
| 63 | // |
---|
| 64 | // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
| 65 | |
---|
| 66 | #include "G4TessellatedSolid.hh" |
---|
| 67 | #include "G4PolyhedronArbitrary.hh" |
---|
| 68 | #include "globals.hh" |
---|
| 69 | #include "Randomize.hh" |
---|
| 70 | |
---|
| 71 | #include <iostream> |
---|
| 72 | |
---|
| 73 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 74 | // |
---|
| 75 | // Standard contructor has blank name and defines no facets. |
---|
| 76 | // |
---|
| 77 | G4TessellatedSolid::G4TessellatedSolid () |
---|
| 78 | : G4VSolid("dummy"), fpPolyhedron(0), cubicVolume(0.), surfaceArea(0.) |
---|
| 79 | { |
---|
| 80 | dirTolerance = 1.0E-14; |
---|
| 81 | |
---|
| 82 | geometryType = "G4TessellatedSolid"; |
---|
| 83 | facets.clear(); |
---|
| 84 | solidClosed = false; |
---|
| 85 | |
---|
| 86 | xMinExtent = kInfinity; |
---|
| 87 | xMaxExtent = -kInfinity; |
---|
| 88 | yMinExtent = kInfinity; |
---|
| 89 | yMaxExtent = -kInfinity; |
---|
| 90 | zMinExtent = kInfinity; |
---|
| 91 | zMaxExtent = -kInfinity; |
---|
| 92 | |
---|
| 93 | SetRandomVectorSet(); |
---|
| 94 | } |
---|
| 95 | |
---|
| 96 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 97 | // |
---|
| 98 | // Alternative constructor. Simple define name and geometry type - no facets |
---|
| 99 | // to detine. |
---|
| 100 | // |
---|
| 101 | G4TessellatedSolid::G4TessellatedSolid (const G4String &name) |
---|
| 102 | : G4VSolid(name), fpPolyhedron(0), cubicVolume(0.), surfaceArea(0.) |
---|
| 103 | { |
---|
| 104 | dirTolerance = 1.0E-14; |
---|
| 105 | |
---|
| 106 | geometryType = "G4TessellatedSolid"; |
---|
| 107 | facets.clear(); |
---|
| 108 | solidClosed = false; |
---|
| 109 | |
---|
| 110 | xMinExtent = kInfinity; |
---|
| 111 | xMaxExtent = -kInfinity; |
---|
| 112 | yMinExtent = kInfinity; |
---|
| 113 | yMaxExtent = -kInfinity; |
---|
| 114 | zMinExtent = kInfinity; |
---|
| 115 | zMaxExtent = -kInfinity; |
---|
| 116 | |
---|
| 117 | SetRandomVectorSet(); |
---|
| 118 | } |
---|
| 119 | |
---|
| 120 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 121 | // |
---|
| 122 | // Fake default constructor - sets only member data and allocates memory |
---|
| 123 | // for usage restricted to object persistency. |
---|
| 124 | // |
---|
| 125 | G4TessellatedSolid::G4TessellatedSolid( __void__& a ) |
---|
| 126 | : G4VSolid(a), fpPolyhedron(0), facets(0), |
---|
| 127 | geometryType("G4TessellatedSolid"), cubicVolume(0.), surfaceArea(0.), |
---|
| 128 | vertexList(), xMinExtent(0.), xMaxExtent(0.), |
---|
| 129 | yMinExtent(0.), yMaxExtent(0.), zMinExtent(0.), zMaxExtent(0.), |
---|
| 130 | solidClosed(false) |
---|
| 131 | { |
---|
| 132 | SetRandomVectorSet(); |
---|
| 133 | } |
---|
| 134 | |
---|
| 135 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 136 | // |
---|
| 137 | // Destructor. |
---|
| 138 | // |
---|
| 139 | G4TessellatedSolid::~G4TessellatedSolid () |
---|
| 140 | { |
---|
| 141 | DeleteObjects (); |
---|
| 142 | } |
---|
| 143 | |
---|
| 144 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 145 | // |
---|
| 146 | // Define copy constructor. |
---|
| 147 | // |
---|
| 148 | G4TessellatedSolid::G4TessellatedSolid (const G4TessellatedSolid &s) |
---|
| 149 | : G4VSolid(s) |
---|
| 150 | { |
---|
| 151 | if (&s == this) { return; } |
---|
| 152 | |
---|
| 153 | dirTolerance = 1.0E-14; |
---|
| 154 | |
---|
| 155 | geometryType = "G4TessellatedSolid"; |
---|
| 156 | facets.clear(); |
---|
| 157 | solidClosed = false; |
---|
| 158 | |
---|
| 159 | xMinExtent = kInfinity; |
---|
| 160 | xMaxExtent = -kInfinity; |
---|
| 161 | yMinExtent = kInfinity; |
---|
| 162 | yMaxExtent = -kInfinity; |
---|
| 163 | zMinExtent = kInfinity; |
---|
| 164 | zMaxExtent = -kInfinity; |
---|
| 165 | |
---|
| 166 | SetRandomVectorSet(); |
---|
| 167 | |
---|
| 168 | CopyObjects (s); |
---|
| 169 | } |
---|
| 170 | |
---|
| 171 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 172 | // |
---|
| 173 | // Define assignment operator. |
---|
| 174 | // |
---|
| 175 | const G4TessellatedSolid & |
---|
| 176 | G4TessellatedSolid::operator= (const G4TessellatedSolid &s) |
---|
| 177 | { |
---|
| 178 | if (&s == this) { return *this; } |
---|
| 179 | |
---|
| 180 | DeleteObjects (); |
---|
| 181 | CopyObjects (s); |
---|
| 182 | |
---|
| 183 | return *this; |
---|
| 184 | } |
---|
| 185 | |
---|
| 186 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 187 | // |
---|
| 188 | void G4TessellatedSolid::DeleteObjects () |
---|
| 189 | { |
---|
| 190 | for (std::vector<G4VFacet *>::iterator f=facets.begin(); f!=facets.end(); f++) |
---|
| 191 | { |
---|
| 192 | delete *f; |
---|
| 193 | } |
---|
| 194 | facets.clear(); |
---|
| 195 | } |
---|
| 196 | |
---|
| 197 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 198 | // |
---|
| 199 | void G4TessellatedSolid::CopyObjects (const G4TessellatedSolid &s) |
---|
| 200 | { |
---|
| 201 | size_t n = s.GetNumberOfFacets(); |
---|
| 202 | for (size_t i=0; i<n; i++) |
---|
| 203 | { |
---|
| 204 | G4VFacet *facetClone = (s.GetFacet(i))->GetClone(); |
---|
| 205 | AddFacet(facetClone); |
---|
| 206 | } |
---|
| 207 | |
---|
| 208 | if ( s.GetSolidClosed() ) { SetSolidClosed(true); } |
---|
| 209 | |
---|
| 210 | // cubicVolume = s.GetCubicVolume(); |
---|
| 211 | } |
---|
| 212 | |
---|
| 213 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 214 | // |
---|
| 215 | // Add a facet to the facet list. Note that you can add, but you cannot |
---|
| 216 | // delete. |
---|
| 217 | // |
---|
| 218 | G4bool G4TessellatedSolid::AddFacet (G4VFacet *aFacet) |
---|
| 219 | { |
---|
| 220 | // Add the facet to the vector. |
---|
| 221 | |
---|
| 222 | if (solidClosed) |
---|
| 223 | { |
---|
| 224 | G4Exception("G4TessellatedSolid::AddFacet()", "InvalidSetup", |
---|
| 225 | JustWarning, "Attempt to add facets when solid is closed."); |
---|
| 226 | return false; |
---|
| 227 | } |
---|
| 228 | else if (aFacet->IsDefined()) |
---|
| 229 | { |
---|
| 230 | if (facets.size() == 0) |
---|
| 231 | { |
---|
| 232 | facets.push_back(aFacet); |
---|
| 233 | } |
---|
| 234 | else |
---|
| 235 | { |
---|
| 236 | G4bool found = false; |
---|
| 237 | FacetI it = facets.begin(); |
---|
| 238 | do |
---|
| 239 | { |
---|
| 240 | found = (**it == *aFacet); |
---|
| 241 | } while (!found && ++it!=facets.end()); |
---|
| 242 | |
---|
| 243 | if (found) |
---|
| 244 | { |
---|
| 245 | delete *it; |
---|
| 246 | facets.erase(it); |
---|
| 247 | } |
---|
| 248 | else |
---|
| 249 | { |
---|
| 250 | facets.push_back(aFacet); |
---|
| 251 | } |
---|
| 252 | } |
---|
| 253 | |
---|
| 254 | return true; |
---|
| 255 | } |
---|
| 256 | else |
---|
| 257 | { |
---|
| 258 | G4Exception("G4TessellatedSolid::AddFacet()", "InvalidSetup", |
---|
| 259 | JustWarning, "Attempt to add facet not properly defined."); |
---|
| 260 | G4cerr << "Facet attributes:" << G4endl; |
---|
| 261 | aFacet->StreamInfo(G4cerr); |
---|
| 262 | G4cerr << G4endl; |
---|
| 263 | |
---|
| 264 | return false; |
---|
| 265 | } |
---|
| 266 | } |
---|
| 267 | |
---|
| 268 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 269 | // |
---|
| 270 | void G4TessellatedSolid::SetSolidClosed (const G4bool t) |
---|
| 271 | { |
---|
| 272 | if (t) |
---|
| 273 | { |
---|
| 274 | vertexList.clear(); |
---|
| 275 | for (FacetCI it=facets.begin(); it!=facets.end(); it++) |
---|
| 276 | { |
---|
| 277 | size_t m = vertexList.size(); |
---|
| 278 | G4ThreeVector p(0.0,0.0,0.0); |
---|
| 279 | for (size_t i=0; i<(*it)->GetNumberOfVertices(); i++) |
---|
| 280 | { |
---|
| 281 | p = (*it)->GetVertex(i); |
---|
| 282 | G4bool found = false; |
---|
| 283 | size_t j = 0; |
---|
| 284 | while (j < m && !found) |
---|
| 285 | { |
---|
| 286 | G4ThreeVector q = vertexList[j]; |
---|
| 287 | found = (q-p).mag() < 0.5*kCarTolerance; |
---|
| 288 | if (!found) j++; |
---|
| 289 | } |
---|
| 290 | |
---|
| 291 | if (!found) |
---|
| 292 | { |
---|
| 293 | vertexList.push_back(p); |
---|
| 294 | (*it)->SetVertexIndex(i,vertexList.size()-1); |
---|
| 295 | } |
---|
| 296 | else |
---|
| 297 | { |
---|
| 298 | (*it)->SetVertexIndex(i,j); |
---|
| 299 | } |
---|
| 300 | } |
---|
| 301 | } |
---|
| 302 | // |
---|
| 303 | // Now update the maximum x, y and z limits of the volume. |
---|
| 304 | // |
---|
| 305 | for (size_t i=0; i<vertexList.size(); i++) |
---|
| 306 | { |
---|
| 307 | G4ThreeVector p = vertexList[i]; |
---|
| 308 | G4double x = p.x(); |
---|
| 309 | G4double y = p.y(); |
---|
| 310 | G4double z = p.z(); |
---|
| 311 | |
---|
| 312 | if (i > 0) |
---|
| 313 | { |
---|
| 314 | if (x > xMaxExtent) xMaxExtent = x; |
---|
| 315 | if (x < xMinExtent) xMinExtent = x; |
---|
| 316 | if (y > yMaxExtent) yMaxExtent = y; |
---|
| 317 | if (y < yMinExtent) yMinExtent = y; |
---|
| 318 | if (z > zMaxExtent) zMaxExtent = z; |
---|
| 319 | if (z < zMinExtent) zMinExtent = z; |
---|
| 320 | } |
---|
| 321 | else |
---|
| 322 | { |
---|
| 323 | xMaxExtent = x; |
---|
| 324 | xMinExtent = x; |
---|
| 325 | yMaxExtent = y; |
---|
| 326 | yMinExtent = y; |
---|
| 327 | zMaxExtent = z; |
---|
| 328 | zMinExtent = z; |
---|
| 329 | } |
---|
| 330 | } |
---|
| 331 | // |
---|
| 332 | // |
---|
| 333 | // Compute extremeFacets, i.e. find those facets that have surface |
---|
| 334 | // planes that bound the volume. |
---|
| 335 | // Note that this is going to reject concaved surfaces as being extreme. Also |
---|
| 336 | // note that if the vertex is on the facet, displacement is zero, so IsInside |
---|
| 337 | // returns true. So will this work?? Need non-equality |
---|
| 338 | // "G4bool inside = displacement < 0.0;" |
---|
| 339 | // or |
---|
| 340 | // "G4bool inside = displacement <= -0.5*kCarTolerance" |
---|
| 341 | // (Notes from PT 13/08/2007). |
---|
| 342 | // |
---|
| 343 | for (FacetCI it=facets.begin(); it!=facets.end(); it++) |
---|
| 344 | { |
---|
| 345 | G4bool isExtreme = true; |
---|
| 346 | for (size_t i=0; i<vertexList.size(); i++) |
---|
| 347 | { |
---|
| 348 | if (!(*it)->IsInside(vertexList[i])) |
---|
| 349 | { |
---|
| 350 | isExtreme = false; |
---|
| 351 | break; |
---|
| 352 | } |
---|
| 353 | } |
---|
| 354 | if (isExtreme) |
---|
| 355 | extremeFacets.insert(*it); |
---|
| 356 | } |
---|
| 357 | solidClosed = true; |
---|
| 358 | } |
---|
| 359 | else |
---|
| 360 | { |
---|
| 361 | solidClosed = false; |
---|
| 362 | } |
---|
| 363 | } |
---|
| 364 | |
---|
| 365 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 366 | // |
---|
| 367 | // GetSolidClosed |
---|
| 368 | // |
---|
| 369 | // Used to determine whether the solid is closed to adding further facets. |
---|
| 370 | // |
---|
| 371 | G4bool G4TessellatedSolid::GetSolidClosed () const |
---|
| 372 | {return solidClosed;} |
---|
| 373 | |
---|
| 374 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 375 | // |
---|
| 376 | // operator+= |
---|
| 377 | // |
---|
| 378 | // This operator allows the user to add two tessellated solids together, so |
---|
| 379 | // that the solid on the left then includes all of the facets in the solid |
---|
| 380 | // on the right. Note that copies of the facets are generated, rather than |
---|
| 381 | // using the original facet set of the solid on the right. |
---|
| 382 | // |
---|
| 383 | const G4TessellatedSolid &G4TessellatedSolid::operator+= |
---|
| 384 | (const G4TessellatedSolid &right) |
---|
| 385 | { |
---|
| 386 | for (size_t i=0; i<right.GetNumberOfFacets(); i++) |
---|
| 387 | AddFacet(right.GetFacet(i)->GetClone()); |
---|
| 388 | return *this; |
---|
| 389 | } |
---|
| 390 | |
---|
| 391 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 392 | // |
---|
| 393 | // GetFacet |
---|
| 394 | // |
---|
| 395 | // Access pointer to facet in solid, indexed by integer i. |
---|
| 396 | // |
---|
| 397 | G4VFacet *G4TessellatedSolid::GetFacet (size_t i) const |
---|
| 398 | { |
---|
| 399 | return facets[i]; |
---|
| 400 | } |
---|
| 401 | |
---|
| 402 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 403 | // |
---|
| 404 | // GetNumberOfFacets |
---|
| 405 | // |
---|
| 406 | size_t G4TessellatedSolid::GetNumberOfFacets () const |
---|
| 407 | { |
---|
| 408 | return facets.size(); |
---|
| 409 | } |
---|
| 410 | |
---|
| 411 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 412 | // |
---|
| 413 | // EInside G4TessellatedSolid::Inside (const G4ThreeVector &p) const |
---|
| 414 | // |
---|
| 415 | // This method must return: |
---|
| 416 | // * kOutside if the point at offset p is outside the shape |
---|
| 417 | // boundaries plus kCarTolerance/2, |
---|
| 418 | // * kSurface if the point is <= kCarTolerance/2 from a surface, or |
---|
| 419 | // * kInside otherwise. |
---|
| 420 | // |
---|
| 421 | EInside G4TessellatedSolid::Inside (const G4ThreeVector &p) const |
---|
| 422 | { |
---|
| 423 | // |
---|
| 424 | // First the simple test - check if we're outside of the X-Y-Z extremes |
---|
| 425 | // of the tessellated solid. |
---|
| 426 | // |
---|
| 427 | if ( p.x() < xMinExtent - kCarTolerance || |
---|
| 428 | p.x() > xMaxExtent + kCarTolerance || |
---|
| 429 | p.y() < yMinExtent - kCarTolerance || |
---|
| 430 | p.y() > yMaxExtent + kCarTolerance || |
---|
| 431 | p.z() < zMinExtent - kCarTolerance || |
---|
| 432 | p.z() > zMaxExtent + kCarTolerance ) |
---|
| 433 | { |
---|
| 434 | return kOutside; |
---|
| 435 | } |
---|
| 436 | |
---|
| 437 | G4double minDist = kInfinity; |
---|
| 438 | // |
---|
| 439 | // |
---|
| 440 | // Check if we are close to a surface |
---|
| 441 | // |
---|
| 442 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 443 | { |
---|
| 444 | G4double dist = (*f)->Distance(p,minDist); |
---|
| 445 | if (dist < minDist) minDist = dist; |
---|
| 446 | if (dist <= 0.5*kCarTolerance) |
---|
| 447 | { |
---|
| 448 | return kSurface; |
---|
| 449 | } |
---|
| 450 | } |
---|
| 451 | // |
---|
| 452 | // |
---|
| 453 | // The following is something of an adaptation of the method implemented by |
---|
| 454 | // Rickard Holmberg augmented with information from Schneider & Eberly, |
---|
| 455 | // "Geometric Tools for Computer Graphics," pp700-701, 2003. In essence, we're |
---|
| 456 | // trying to determine whether we're inside the volume by projecting a few rays |
---|
| 457 | // and determining if the first surface crossed is has a normal vector between |
---|
| 458 | // 0 to pi/2 (out-going) or pi/2 to pi (in-going). We should also avoid rays |
---|
| 459 | // which are nearly within the plane of the tessellated surface, and therefore |
---|
| 460 | // produce rays randomly. For the moment, this is a bit over-engineered |
---|
| 461 | // (belt-braces-and-ducttape). |
---|
| 462 | // |
---|
| 463 | #if G4SPECSDEBUG |
---|
| 464 | G4int nTry = 7; |
---|
| 465 | #else |
---|
| 466 | G4int nTry = 3; |
---|
| 467 | #endif |
---|
| 468 | G4double distOut = kInfinity; |
---|
| 469 | G4double distIn = kInfinity; |
---|
| 470 | G4double distO = 0.0; |
---|
| 471 | G4double distI = 0.0; |
---|
| 472 | G4double distFromSurfaceO = 0.0; |
---|
| 473 | G4double distFromSurfaceI = 0.0; |
---|
| 474 | G4ThreeVector normalO(0.0,0.0,0.0); |
---|
| 475 | G4ThreeVector normalI(0.0,0.0,0.0); |
---|
| 476 | G4bool crossingO = false; |
---|
| 477 | G4bool crossingI = false; |
---|
| 478 | EInside location = kOutside; |
---|
| 479 | EInside locationprime = kOutside; |
---|
| 480 | G4int m = 0; |
---|
| 481 | |
---|
| 482 | for (G4int i=0; i<nTry; i++) |
---|
| 483 | { |
---|
| 484 | G4bool nearParallel = false; |
---|
| 485 | do |
---|
| 486 | { |
---|
| 487 | // |
---|
| 488 | // |
---|
| 489 | // We loop until we find direction where the vector is not nearly parallel |
---|
| 490 | // to the surface of any facet since this causes ambiguities. The usual |
---|
| 491 | // case is that the angles should be sufficiently different, but there are 20 |
---|
| 492 | // random directions to select from - hopefully sufficient. |
---|
| 493 | // |
---|
| 494 | distOut = kInfinity; |
---|
| 495 | distIn = kInfinity; |
---|
| 496 | G4ThreeVector v = randir[m]; |
---|
| 497 | m++; |
---|
| 498 | FacetCI f = facets.begin(); |
---|
| 499 | do |
---|
| 500 | { |
---|
| 501 | // |
---|
| 502 | // |
---|
| 503 | // Here we loop through the facets to find out if there is an intersection |
---|
| 504 | // between the ray and that facet. The test if performed separately whether |
---|
| 505 | // the ray is entering the facet or exiting. |
---|
| 506 | // |
---|
| 507 | crossingO = ((*f)->Intersect(p,v,true,distO,distFromSurfaceO,normalO)); |
---|
| 508 | crossingI = ((*f)->Intersect(p,v,false,distI,distFromSurfaceI,normalI)); |
---|
| 509 | if (crossingO || crossingI) |
---|
| 510 | { |
---|
| 511 | nearParallel = (crossingO && std::abs(normalO.dot(v))<dirTolerance) || |
---|
| 512 | (crossingI && std::abs(normalI.dot(v))<dirTolerance); |
---|
| 513 | if (!nearParallel) |
---|
| 514 | { |
---|
| 515 | if (crossingO && distO > 0.0 && distO < distOut) distOut = distO; |
---|
| 516 | if (crossingI && distI > 0.0 && distI < distIn) distIn = distI; |
---|
| 517 | } |
---|
| 518 | } |
---|
| 519 | } while (!nearParallel && ++f!=facets.end()); |
---|
| 520 | } while (nearParallel && m!=maxTries); |
---|
| 521 | |
---|
[1228] | 522 | #ifdef G4VERBOSE |
---|
[831] | 523 | if (m == maxTries) |
---|
| 524 | { |
---|
| 525 | // |
---|
| 526 | // |
---|
| 527 | // We've run out of random vector directions. If nTries is set sufficiently |
---|
| 528 | // low (nTries <= 0.5*maxTries) then this would indicate that there is |
---|
| 529 | // something wrong with geometry. |
---|
| 530 | // |
---|
[1228] | 531 | G4cout.precision(16) ; |
---|
| 532 | G4cout << G4endl ; |
---|
| 533 | G4cout << "Solid name = " << GetName() << G4endl; |
---|
| 534 | G4cout << "Geometry Type = " << geometryType << G4endl; |
---|
| 535 | G4cout << "Number of facets = " << facets.size() << G4endl; |
---|
| 536 | G4cout << "Position:" << G4endl << G4endl ; |
---|
| 537 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
| 538 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
| 539 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
[831] | 540 | G4Exception("G4TessellatedSolid::Inside()", |
---|
[1228] | 541 | "UnknownInsideOutside-MaxTries", JustWarning, |
---|
[831] | 542 | "Cannot determine whether point is inside or outside volume!"); |
---|
| 543 | } |
---|
[1228] | 544 | #endif |
---|
[831] | 545 | // |
---|
| 546 | // |
---|
| 547 | // In the next if-then-elseif string the logic is as follows: |
---|
| 548 | // (1) You don't hit anything so cannot be inside volume, provided volume |
---|
| 549 | // constructed correctly! |
---|
| 550 | // (2) Distance to inside (ie. nearest facet such that you enter facet) is |
---|
| 551 | // shorter than distance to outside (nearest facet such that you exit |
---|
| 552 | // facet) - on condition of safety distance - therefore we're outside. |
---|
| 553 | // (3) Distance to outside is shorter than distance to inside therefore we're |
---|
| 554 | // inside. |
---|
| 555 | // |
---|
| 556 | if (distIn == kInfinity && distOut == kInfinity) |
---|
| 557 | locationprime = kOutside; |
---|
| 558 | else if (distIn <= distOut - kCarTolerance*0.5) |
---|
| 559 | locationprime = kOutside; |
---|
| 560 | else if (distOut <= distIn - kCarTolerance*0.5) |
---|
| 561 | locationprime = kInside; |
---|
| 562 | |
---|
[1228] | 563 | if (i == 0) { location = locationprime; } |
---|
| 564 | #ifdef G4VERBOSE |
---|
[831] | 565 | else if (locationprime != location) |
---|
| 566 | { |
---|
| 567 | // |
---|
| 568 | // |
---|
| 569 | // Different ray directions result in different answer. Seems like the |
---|
| 570 | // geometry is not constructed correctly. |
---|
| 571 | // |
---|
[1228] | 572 | G4cout.precision(16) ; |
---|
| 573 | G4cout << G4endl ; |
---|
| 574 | G4cout << "Solid name = " << GetName() << G4endl; |
---|
| 575 | G4cout << "Geometry Type = " << geometryType << G4endl; |
---|
| 576 | G4cout << "Number of facets = " << facets.size() << G4endl; |
---|
| 577 | G4cout << "Position:" << G4endl << G4endl ; |
---|
| 578 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
| 579 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
| 580 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
[831] | 581 | G4Exception("G4TessellatedSolid::Inside()", |
---|
[1228] | 582 | "UnknownInsideOutside", JustWarning, |
---|
[831] | 583 | "Cannot determine whether point is inside or outside volume!"); |
---|
| 584 | } |
---|
[1228] | 585 | #endif |
---|
[831] | 586 | } |
---|
| 587 | |
---|
| 588 | return location; |
---|
| 589 | } |
---|
| 590 | |
---|
| 591 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 592 | // |
---|
| 593 | // G4ThreeVector G4TessellatedSolid::SurfaceNormal (const G4ThreeVector &p) const |
---|
| 594 | // |
---|
| 595 | // Return the outwards pointing unit normal of the shape for the |
---|
| 596 | // surface closest to the point at offset p. |
---|
| 597 | |
---|
| 598 | G4ThreeVector G4TessellatedSolid::SurfaceNormal (const G4ThreeVector &p) const |
---|
| 599 | { |
---|
| 600 | FacetCI minFacet; |
---|
| 601 | G4double minDist = kInfinity; |
---|
| 602 | G4double dist = 0.0; |
---|
| 603 | G4ThreeVector normal; |
---|
| 604 | |
---|
| 605 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 606 | { |
---|
| 607 | dist = (*f)->Distance(p,minDist); |
---|
| 608 | if (dist < minDist) |
---|
| 609 | { |
---|
| 610 | minDist = dist; |
---|
| 611 | minFacet = f; |
---|
| 612 | } |
---|
| 613 | } |
---|
| 614 | |
---|
| 615 | if (minDist != kInfinity) |
---|
| 616 | { |
---|
| 617 | normal = (*minFacet)->GetSurfaceNormal(); |
---|
| 618 | } |
---|
| 619 | else |
---|
| 620 | { |
---|
| 621 | #ifdef G4VERBOSE |
---|
| 622 | G4cout << "WARNING - G4TessellatedSolid::SurfaceNormal(p)" << G4endl |
---|
| 623 | << " No facets found for point: " << p << " !" << G4endl |
---|
| 624 | << " Returning approximated value for normal." << G4endl; |
---|
| 625 | G4Exception("G4TessellatedSolid::SurfaceNormal(p)", "Notification", |
---|
| 626 | JustWarning, "Point p is not on surface !?" ); |
---|
| 627 | #endif |
---|
| 628 | normal = (p.z()>0 ? G4ThreeVector(0,0,1) : G4ThreeVector(0,0,-1)); |
---|
| 629 | } |
---|
| 630 | |
---|
| 631 | return normal; |
---|
| 632 | } |
---|
| 633 | |
---|
| 634 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 635 | // |
---|
| 636 | // G4double DistanceToIn(const G4ThreeVector& p, const G4ThreeVector& v) |
---|
| 637 | // |
---|
| 638 | // Return the distance along the normalised vector v to the shape, |
---|
| 639 | // from the point at offset p. If there is no intersection, return |
---|
| 640 | // kInfinity. The first intersection resulting from âleavingâ a |
---|
| 641 | // surface/volume is discarded. Hence, this is tolerant of points on |
---|
| 642 | // surface of shape. |
---|
| 643 | |
---|
| 644 | G4double G4TessellatedSolid::DistanceToIn (const G4ThreeVector &p, |
---|
| 645 | const G4ThreeVector &v) const |
---|
| 646 | { |
---|
| 647 | G4double minDist = kInfinity; |
---|
| 648 | G4double dist = 0.0; |
---|
| 649 | G4double distFromSurface = 0.0; |
---|
| 650 | G4ThreeVector normal(0.0,0.0,0.0); |
---|
| 651 | |
---|
| 652 | #if G4SPECSDEBUG |
---|
| 653 | if ( Inside(p) == kInside ) |
---|
| 654 | { |
---|
| 655 | G4cout.precision(16) ; |
---|
| 656 | G4cout << G4endl ; |
---|
| 657 | // DumpInfo(); |
---|
| 658 | G4cout << "Position:" << G4endl << G4endl ; |
---|
| 659 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
| 660 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
| 661 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
| 662 | G4cout << "DistanceToOut(p) == " << DistanceToOut(p) << G4endl; |
---|
| 663 | G4Exception("G4TriangularFacet::DistanceToIn(p,v)", "Notification", JustWarning, |
---|
| 664 | "Point p is already inside!?" ); |
---|
| 665 | } |
---|
| 666 | #endif |
---|
| 667 | |
---|
| 668 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 669 | { |
---|
| 670 | if ((*f)->Intersect(p,v,false,dist,distFromSurface,normal)) |
---|
| 671 | { |
---|
| 672 | if (distFromSurface > 0.5*kCarTolerance && dist >= 0.0 && dist < minDist) |
---|
| 673 | { |
---|
| 674 | minDist = dist; |
---|
| 675 | } |
---|
| 676 | } |
---|
| 677 | } |
---|
| 678 | |
---|
| 679 | return minDist; |
---|
| 680 | } |
---|
| 681 | |
---|
| 682 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 683 | // |
---|
| 684 | // G4double DistanceToIn(const G4ThreeVector& p) |
---|
| 685 | // |
---|
| 686 | // Calculate distance to nearest surface of shape from an outside point p. The |
---|
| 687 | // distance can be an underestimate. |
---|
| 688 | |
---|
| 689 | G4double G4TessellatedSolid::DistanceToIn (const G4ThreeVector &p) const |
---|
| 690 | { |
---|
| 691 | G4double minDist = kInfinity; |
---|
| 692 | G4double dist = 0.0; |
---|
| 693 | |
---|
| 694 | #if G4SPECSDEBUG |
---|
| 695 | if ( Inside(p) == kInside ) |
---|
| 696 | { |
---|
| 697 | G4cout.precision(16) ; |
---|
| 698 | G4cout << G4endl ; |
---|
| 699 | // DumpInfo(); |
---|
| 700 | G4cout << "Position:" << G4endl << G4endl ; |
---|
| 701 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
| 702 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
| 703 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
| 704 | G4cout << "DistanceToOut(p) == " << DistanceToOut(p) << G4endl; |
---|
| 705 | G4Exception("G4TriangularFacet::DistanceToIn(p)", "Notification", JustWarning, |
---|
| 706 | "Point p is already inside!?" ); |
---|
| 707 | } |
---|
| 708 | #endif |
---|
| 709 | |
---|
| 710 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 711 | { |
---|
| 712 | dist = (*f)->Distance(p,minDist,false); |
---|
| 713 | if (dist < minDist) { minDist = dist; } |
---|
| 714 | } |
---|
| 715 | |
---|
| 716 | return minDist; |
---|
| 717 | } |
---|
| 718 | |
---|
| 719 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 720 | // |
---|
| 721 | // G4double DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v, |
---|
| 722 | // const G4bool calcNorm=false, |
---|
| 723 | // G4bool *validNorm=0, G4ThreeVector *n=0); |
---|
| 724 | // |
---|
| 725 | // Return distance along the normalised vector v to the shape, from a |
---|
| 726 | // point at an offset p inside or on the surface of the |
---|
| 727 | // shape. Intersections with surfaces, when the point is not greater |
---|
| 728 | // than kCarTolerance/2 from a surface, must be ignored. |
---|
| 729 | // If calcNorm is true, then it must also set validNorm to either |
---|
| 730 | // * true, if the solid lies entirely behind or on the exiting |
---|
| 731 | // surface. Then it must set n to the outwards normal vector |
---|
| 732 | // (the Magnitude of the vector is not defined). |
---|
| 733 | // * false, if the solid does not lie entirely behind or on the |
---|
| 734 | // exiting surface. |
---|
| 735 | // If calcNorm is false, then validNorm and n are unused. |
---|
| 736 | |
---|
| 737 | G4double G4TessellatedSolid::DistanceToOut (const G4ThreeVector &p, |
---|
| 738 | const G4ThreeVector &v, const G4bool calcNorm, |
---|
| 739 | G4bool *validNorm, G4ThreeVector *n) const |
---|
| 740 | { |
---|
| 741 | G4double minDist = kInfinity; |
---|
| 742 | G4double dist = 0.0; |
---|
| 743 | G4double distFromSurface = 0.0; |
---|
| 744 | G4ThreeVector normal(0.0,0.0,0.0); |
---|
| 745 | G4ThreeVector minNormal(0.0,0.0,0.0); |
---|
| 746 | |
---|
| 747 | #if G4SPECSDEBUG |
---|
| 748 | if ( Inside(p) == kOutside ) |
---|
| 749 | { |
---|
| 750 | G4cout.precision(16) ; |
---|
| 751 | G4cout << G4endl ; |
---|
| 752 | // DumpInfo(); |
---|
| 753 | G4cout << "Position:" << G4endl << G4endl ; |
---|
| 754 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
| 755 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
| 756 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
| 757 | G4cout << "DistanceToIn(p) == " << DistanceToIn(p) << G4endl; |
---|
| 758 | G4Exception("G4TriangularFacet::DistanceToOut(p)", "Notification", JustWarning, |
---|
| 759 | "Point p is already outside !?" ); |
---|
| 760 | } |
---|
| 761 | #endif |
---|
| 762 | |
---|
| 763 | G4bool isExtreme = false; |
---|
| 764 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 765 | { |
---|
| 766 | if ((*f)->Intersect(p,v,true,dist,distFromSurface,normal)) |
---|
| 767 | { |
---|
| 768 | if (distFromSurface > 0.0 && distFromSurface <= 0.5*kCarTolerance && |
---|
| 769 | (*f)->Distance(p,kCarTolerance) <= 0.5*kCarTolerance) |
---|
| 770 | { |
---|
| 771 | // We are on a surface. Return zero. |
---|
| 772 | if (calcNorm) { |
---|
| 773 | *validNorm = extremeFacets.count(*f); |
---|
| 774 | *n = SurfaceNormal(p); |
---|
| 775 | } |
---|
| 776 | return 0.0; |
---|
| 777 | } |
---|
| 778 | if (dist >= 0.0 && dist < minDist) |
---|
| 779 | { |
---|
| 780 | minDist = dist; |
---|
| 781 | minNormal = normal; |
---|
| 782 | isExtreme = extremeFacets.count(*f); |
---|
| 783 | } |
---|
| 784 | } |
---|
| 785 | } |
---|
| 786 | |
---|
| 787 | if (minDist < kInfinity) |
---|
| 788 | { |
---|
| 789 | if (calcNorm) |
---|
| 790 | { |
---|
| 791 | *validNorm = isExtreme; |
---|
| 792 | *n = minNormal; |
---|
| 793 | } |
---|
| 794 | return minDist; |
---|
| 795 | } |
---|
| 796 | else |
---|
| 797 | { |
---|
| 798 | // No intersection found |
---|
| 799 | if (calcNorm) |
---|
| 800 | { |
---|
| 801 | *validNorm = false; |
---|
| 802 | *n = SurfaceNormal(p); |
---|
| 803 | } |
---|
| 804 | return 0.0; |
---|
| 805 | } |
---|
| 806 | } |
---|
| 807 | |
---|
| 808 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 809 | // |
---|
| 810 | // G4double DistanceToOut(const G4ThreeVector& p) |
---|
| 811 | // |
---|
| 812 | // Calculate distance to nearest surface of shape from an inside |
---|
| 813 | // point. The distance can be an underestimate. |
---|
| 814 | |
---|
| 815 | G4double G4TessellatedSolid::DistanceToOut (const G4ThreeVector &p) const |
---|
| 816 | { |
---|
| 817 | G4double minDist = kInfinity; |
---|
| 818 | G4double dist = 0.0; |
---|
| 819 | |
---|
| 820 | #if G4SPECSDEBUG |
---|
| 821 | if ( Inside(p) == kOutside ) |
---|
| 822 | { |
---|
| 823 | G4cout.precision(16) ; |
---|
| 824 | G4cout << G4endl ; |
---|
| 825 | // DumpInfo(); |
---|
| 826 | G4cout << "Position:" << G4endl << G4endl ; |
---|
| 827 | G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; |
---|
| 828 | G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; |
---|
| 829 | G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; |
---|
| 830 | G4cout << "DistanceToIn(p) == " << DistanceToIn(p) << G4endl; |
---|
| 831 | G4Exception("G4TriangularFacet::DistanceToOut(p)", "Notification", JustWarning, |
---|
| 832 | "Point p is already outside !?" ); |
---|
| 833 | } |
---|
| 834 | #endif |
---|
| 835 | |
---|
| 836 | for (FacetCI f=facets.begin(); f!=facets.end(); f++) |
---|
| 837 | { |
---|
| 838 | dist = (*f)->Distance(p,minDist,true); |
---|
| 839 | if (dist < minDist) minDist = dist; |
---|
| 840 | } |
---|
| 841 | |
---|
| 842 | return minDist; |
---|
| 843 | } |
---|
| 844 | |
---|
| 845 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 846 | // |
---|
| 847 | // G4GeometryType GetEntityType() const; |
---|
| 848 | // |
---|
| 849 | // Provide identification of the class of an object (required for |
---|
| 850 | // persistency and STEP interface). |
---|
| 851 | // |
---|
| 852 | G4GeometryType G4TessellatedSolid::GetEntityType () const |
---|
| 853 | { |
---|
| 854 | return geometryType; |
---|
| 855 | } |
---|
| 856 | |
---|
| 857 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 858 | // |
---|
| 859 | void G4TessellatedSolid::DescribeYourselfTo (G4VGraphicsScene& scene) const |
---|
| 860 | { |
---|
| 861 | scene.AddSolid (*this); |
---|
| 862 | } |
---|
| 863 | |
---|
| 864 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 865 | // |
---|
| 866 | // Dispatch to parameterisation for replication mechanism dimension |
---|
| 867 | // computation & modification. |
---|
| 868 | // |
---|
| 869 | //void G4TessellatedSolid::ComputeDimensions (G4VPVParameterisation* p, |
---|
| 870 | // const G4int n, const G4VPhysicalVolume* pRep) const |
---|
| 871 | //{ |
---|
| 872 | // G4VSolid *ptr = 0; |
---|
| 873 | // ptr = *this; |
---|
| 874 | // p->ComputeDimensions(ptr,n,pRep); |
---|
| 875 | //} |
---|
| 876 | |
---|
| 877 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 878 | // |
---|
| 879 | std::ostream &G4TessellatedSolid::StreamInfo(std::ostream &os) const |
---|
| 880 | { |
---|
| 881 | os << G4endl; |
---|
| 882 | os << "Geometry Type = " << geometryType << G4endl; |
---|
| 883 | os << "Number of facets = " << facets.size() << G4endl; |
---|
| 884 | |
---|
| 885 | for (FacetCI f = facets.begin(); f != facets.end(); f++) |
---|
| 886 | { |
---|
| 887 | os << "FACET # = " << f-facets.begin()+1 << G4endl; |
---|
| 888 | (*f)->StreamInfo(os); |
---|
| 889 | } |
---|
| 890 | os <<G4endl; |
---|
| 891 | |
---|
| 892 | return os; |
---|
| 893 | } |
---|
| 894 | |
---|
| 895 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 896 | // |
---|
| 897 | G4Polyhedron *G4TessellatedSolid::CreatePolyhedron () const |
---|
| 898 | { |
---|
| 899 | size_t nVertices = vertexList.size(); |
---|
| 900 | size_t nFacets = facets.size(); |
---|
| 901 | G4PolyhedronArbitrary *polyhedron = |
---|
| 902 | new G4PolyhedronArbitrary (nVertices, nFacets); |
---|
| 903 | for (G4ThreeVectorList::const_iterator v = vertexList.begin(); |
---|
| 904 | v!=vertexList.end(); v++) polyhedron->AddVertex(*v); |
---|
| 905 | |
---|
| 906 | for (FacetCI f=facets.begin(); f != facets.end(); f++) |
---|
| 907 | { |
---|
| 908 | size_t v[4]; |
---|
| 909 | for (size_t j=0; j<4; j++) |
---|
| 910 | { |
---|
| 911 | size_t i = (*f)->GetVertexIndex(j); |
---|
| 912 | if (i == 999999999) v[j] = 0; |
---|
| 913 | else v[j] = i+1; |
---|
| 914 | } |
---|
| 915 | if ((*f)->GetEntityType() == "G4RectangularFacet") |
---|
| 916 | { |
---|
| 917 | size_t i = v[3]; |
---|
| 918 | v[3] = v[2]; |
---|
| 919 | v[2] = i; |
---|
| 920 | } |
---|
| 921 | polyhedron->AddFacet(v[0],v[1],v[2],v[3]); |
---|
| 922 | } |
---|
| 923 | |
---|
| 924 | return (G4Polyhedron*) polyhedron; |
---|
| 925 | } |
---|
| 926 | |
---|
| 927 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 928 | // |
---|
| 929 | G4NURBS *G4TessellatedSolid::CreateNURBS () const |
---|
| 930 | { |
---|
| 931 | return 0; |
---|
| 932 | } |
---|
| 933 | |
---|
| 934 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 935 | // |
---|
| 936 | // GetPolyhedron |
---|
| 937 | // |
---|
| 938 | G4Polyhedron* G4TessellatedSolid::GetPolyhedron () const |
---|
| 939 | { |
---|
| 940 | if (!fpPolyhedron || |
---|
| 941 | fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != |
---|
| 942 | fpPolyhedron->GetNumberOfRotationSteps()) |
---|
| 943 | { |
---|
| 944 | delete fpPolyhedron; |
---|
| 945 | fpPolyhedron = CreatePolyhedron(); |
---|
| 946 | } |
---|
| 947 | return fpPolyhedron; |
---|
| 948 | } |
---|
| 949 | |
---|
| 950 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 951 | // |
---|
| 952 | // CalculateExtent |
---|
| 953 | // |
---|
| 954 | // Based on correction provided by Stan Seibert, University of Texas. |
---|
| 955 | // |
---|
| 956 | G4bool |
---|
| 957 | G4TessellatedSolid::CalculateExtent(const EAxis pAxis, |
---|
| 958 | const G4VoxelLimits& pVoxelLimit, |
---|
| 959 | const G4AffineTransform& pTransform, |
---|
| 960 | G4double& pMin, G4double& pMax) const |
---|
| 961 | { |
---|
| 962 | G4ThreeVectorList transVertexList(vertexList); |
---|
| 963 | |
---|
| 964 | // Put solid into transformed frame |
---|
| 965 | for (size_t i=0; i<vertexList.size(); i++) |
---|
| 966 | { pTransform.ApplyPointTransform(transVertexList[i]); } |
---|
| 967 | |
---|
| 968 | // Find min and max extent in each dimension |
---|
| 969 | G4ThreeVector minExtent(kInfinity, kInfinity, kInfinity); |
---|
| 970 | G4ThreeVector maxExtent(-kInfinity, -kInfinity, -kInfinity); |
---|
| 971 | for (size_t i=0; i<transVertexList.size(); i++) |
---|
| 972 | { |
---|
| 973 | for (G4int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; axis++) |
---|
| 974 | { |
---|
| 975 | G4double coordinate = transVertexList[i][axis]; |
---|
| 976 | if (coordinate < minExtent[axis]) |
---|
| 977 | { minExtent[axis] = coordinate; } |
---|
| 978 | if (coordinate > maxExtent[axis]) |
---|
| 979 | { maxExtent[axis] = coordinate; } |
---|
| 980 | } |
---|
| 981 | } |
---|
| 982 | |
---|
| 983 | // Check for containment and clamp to voxel boundaries |
---|
| 984 | for (G4int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; axis++) |
---|
| 985 | { |
---|
| 986 | EAxis geomAxis = kXAxis; // G4 geom classes use different index type |
---|
| 987 | switch(axis) |
---|
| 988 | { |
---|
| 989 | case G4ThreeVector::X: geomAxis = kXAxis; break; |
---|
| 990 | case G4ThreeVector::Y: geomAxis = kYAxis; break; |
---|
| 991 | case G4ThreeVector::Z: geomAxis = kZAxis; break; |
---|
| 992 | } |
---|
| 993 | G4bool isLimited = pVoxelLimit.IsLimited(geomAxis); |
---|
| 994 | G4double voxelMinExtent = pVoxelLimit.GetMinExtent(geomAxis); |
---|
| 995 | G4double voxelMaxExtent = pVoxelLimit.GetMaxExtent(geomAxis); |
---|
| 996 | |
---|
| 997 | if (isLimited) |
---|
| 998 | { |
---|
| 999 | if ( minExtent[axis] > voxelMaxExtent+kCarTolerance || |
---|
| 1000 | maxExtent[axis] < voxelMinExtent-kCarTolerance ) |
---|
| 1001 | { |
---|
| 1002 | return false ; |
---|
| 1003 | } |
---|
| 1004 | else |
---|
| 1005 | { |
---|
| 1006 | if (minExtent[axis] < voxelMinExtent) |
---|
| 1007 | { |
---|
| 1008 | minExtent[axis] = voxelMinExtent ; |
---|
| 1009 | } |
---|
| 1010 | if (maxExtent[axis] > voxelMaxExtent) |
---|
| 1011 | { |
---|
| 1012 | maxExtent[axis] = voxelMaxExtent; |
---|
| 1013 | } |
---|
| 1014 | } |
---|
| 1015 | } |
---|
| 1016 | } |
---|
| 1017 | |
---|
| 1018 | // Convert pAxis into G4ThreeVector index |
---|
| 1019 | G4int vecAxis=0; |
---|
| 1020 | switch(pAxis) |
---|
| 1021 | { |
---|
| 1022 | case kXAxis: vecAxis = G4ThreeVector::X; break; |
---|
| 1023 | case kYAxis: vecAxis = G4ThreeVector::Y; break; |
---|
| 1024 | case kZAxis: vecAxis = G4ThreeVector::Z; break; |
---|
| 1025 | default: break; |
---|
| 1026 | } |
---|
| 1027 | |
---|
| 1028 | pMin = minExtent[vecAxis] - kCarTolerance; |
---|
| 1029 | pMax = maxExtent[vecAxis] + kCarTolerance; |
---|
| 1030 | |
---|
| 1031 | return true; |
---|
| 1032 | } |
---|
| 1033 | |
---|
| 1034 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1035 | // |
---|
| 1036 | G4double G4TessellatedSolid::GetMinXExtent () const |
---|
| 1037 | {return xMinExtent;} |
---|
| 1038 | |
---|
| 1039 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1040 | // |
---|
| 1041 | G4double G4TessellatedSolid::GetMaxXExtent () const |
---|
| 1042 | {return xMaxExtent;} |
---|
| 1043 | |
---|
| 1044 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1045 | // |
---|
| 1046 | G4double G4TessellatedSolid::GetMinYExtent () const |
---|
| 1047 | {return yMinExtent;} |
---|
| 1048 | |
---|
| 1049 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1050 | // |
---|
| 1051 | G4double G4TessellatedSolid::GetMaxYExtent () const |
---|
| 1052 | {return yMaxExtent;} |
---|
| 1053 | |
---|
| 1054 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1055 | // |
---|
| 1056 | G4double G4TessellatedSolid::GetMinZExtent () const |
---|
| 1057 | {return zMinExtent;} |
---|
| 1058 | |
---|
| 1059 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1060 | // |
---|
| 1061 | G4double G4TessellatedSolid::GetMaxZExtent () const |
---|
| 1062 | {return zMaxExtent;} |
---|
| 1063 | |
---|
| 1064 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1065 | // |
---|
| 1066 | G4VisExtent G4TessellatedSolid::GetExtent () const |
---|
| 1067 | { |
---|
| 1068 | return G4VisExtent (xMinExtent, xMaxExtent, yMinExtent, yMaxExtent, |
---|
| 1069 | zMinExtent, zMaxExtent); |
---|
| 1070 | } |
---|
| 1071 | |
---|
| 1072 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1073 | // |
---|
| 1074 | G4double G4TessellatedSolid::GetCubicVolume () |
---|
| 1075 | { |
---|
| 1076 | if(cubicVolume != 0.) {;} |
---|
| 1077 | else { cubicVolume = G4VSolid::GetCubicVolume(); } |
---|
| 1078 | return cubicVolume; |
---|
| 1079 | } |
---|
| 1080 | |
---|
| 1081 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1082 | // |
---|
| 1083 | G4double G4TessellatedSolid::GetSurfaceArea () |
---|
| 1084 | { |
---|
| 1085 | if(surfaceArea != 0.) { return surfaceArea; } |
---|
| 1086 | |
---|
| 1087 | for (FacetI f=facets.begin(); f!=facets.end(); f++) |
---|
| 1088 | { |
---|
| 1089 | surfaceArea += (*f)->GetArea(); |
---|
| 1090 | } |
---|
| 1091 | return surfaceArea; |
---|
| 1092 | } |
---|
| 1093 | |
---|
| 1094 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1095 | // |
---|
| 1096 | G4ThreeVector G4TessellatedSolid::GetPointOnSurface() const |
---|
| 1097 | { |
---|
| 1098 | // Select randomly a facet and return a random point on it |
---|
| 1099 | |
---|
| 1100 | G4int i = CLHEP::RandFlat::shootInt(facets.size()); |
---|
| 1101 | return facets[i]->GetPointOnFace(); |
---|
| 1102 | } |
---|
| 1103 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1104 | // |
---|
| 1105 | // SetRandomVectorSet |
---|
| 1106 | // |
---|
| 1107 | // This is a set of predefined random vectors (if that isn't a contradition |
---|
| 1108 | // in terms!) used to generate rays from a user-defined point. The member |
---|
| 1109 | // function Inside uses these to determine whether the point is inside or |
---|
| 1110 | // outside of the tessellated solid. All vectors should be unit vectors. |
---|
| 1111 | // |
---|
| 1112 | void G4TessellatedSolid::SetRandomVectorSet() |
---|
| 1113 | { |
---|
| 1114 | randir[0] = G4ThreeVector(-0.9577428892113370, 0.2732676269591740, 0.0897405271949221); |
---|
| 1115 | randir[1] = G4ThreeVector(-0.8331264504940770,-0.5162067214954600,-0.1985722492445700); |
---|
| 1116 | randir[2] = G4ThreeVector(-0.1516671651108820, 0.9666292616127460, 0.2064580868390110); |
---|
| 1117 | randir[3] = G4ThreeVector( 0.6570250350323190,-0.6944539025883300, 0.2933460081893360); |
---|
| 1118 | randir[4] = G4ThreeVector(-0.4820456281280320,-0.6331060000098690,-0.6056474264406270); |
---|
| 1119 | randir[5] = G4ThreeVector( 0.7629032554236800, 0.1016854697539910,-0.6384658864065180); |
---|
| 1120 | randir[6] = G4ThreeVector( 0.7689540409061150, 0.5034929891988220, 0.3939600142169160); |
---|
| 1121 | randir[7] = G4ThreeVector( 0.5765188359255740, 0.5997271636278330,-0.5549354566343150); |
---|
| 1122 | randir[8] = G4ThreeVector( 0.6660632777862070,-0.6362809868288380, 0.3892379937580790); |
---|
| 1123 | randir[9] = G4ThreeVector( 0.3824415020414780, 0.6541792713761380,-0.6525243125110690); |
---|
| 1124 | randir[10] = G4ThreeVector(-0.5107726564526760, 0.6020905056811610, 0.6136760679616570); |
---|
| 1125 | randir[11] = G4ThreeVector( 0.7459135439578050, 0.6618796061649330, 0.0743530220183488); |
---|
| 1126 | randir[12] = G4ThreeVector( 0.1536405855311580, 0.8117477913978260,-0.5634359711967240); |
---|
| 1127 | randir[13] = G4ThreeVector( 0.0744395301705579,-0.8707110101772920,-0.4861286795736560); |
---|
| 1128 | randir[14] = G4ThreeVector(-0.1665874645185400, 0.6018553940549240,-0.7810369397872780); |
---|
| 1129 | randir[15] = G4ThreeVector( 0.7766902003633100, 0.6014617505959970,-0.1870724331097450); |
---|
| 1130 | randir[16] = G4ThreeVector(-0.8710128685847430,-0.1434320216603030,-0.4698551243971010); |
---|
| 1131 | randir[17] = G4ThreeVector( 0.8901082092766820,-0.4388411398893870, 0.1229871120030100); |
---|
| 1132 | randir[18] = G4ThreeVector(-0.6430417431544370,-0.3295938228697690, 0.6912779675984150); |
---|
| 1133 | randir[19] = G4ThreeVector( 0.6331124368380410, 0.6306211461665000, 0.4488714875425340); |
---|
| 1134 | |
---|
| 1135 | maxTries = 20; |
---|
| 1136 | } |
---|