Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result << 18 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 19 // * Vanderbilt University Free Electron Laser Center * 20 // * By using, copying, modifying or distri << 20 // * Vanderbilt University, Nashville, TN, USA * 21 // * any work based on the software) you ag << 21 // * Development supported by: * 22 // * use in resulting scientific publicati << 22 // * United States MFEL program under grant FA9550-04-1-0045 * 23 // * acceptance of all terms of the Geant4 Sof << 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 // * * 24 // ******************************************* 28 // ******************************************************************** 25 // 29 // 26 // Implementation for G4UTet wrapper class << 30 // $Id:$ 27 // 31 // 28 // 1.11.13 G.Cosmo, CERN << 32 // >> 33 // Implementation for G4UTet wrapper class 29 // ------------------------------------------- 34 // -------------------------------------------------------------------- 30 35 31 #include "G4Tet.hh" 36 #include "G4Tet.hh" 32 #include "G4UTet.hh" 37 #include "G4UTet.hh" 33 38 34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G << 35 << 36 #include "G4AffineTransform.hh" << 37 #include "G4VPVParameterisation.hh" << 38 #include "G4BoundingEnvelope.hh" << 39 << 40 using namespace CLHEP; << 41 << 42 ////////////////////////////////////////////// 39 //////////////////////////////////////////////////////////////////////// 43 // 40 // 44 // Constructor - create a tetrahedron 41 // Constructor - create a tetrahedron 45 // This class is implemented separately from g 42 // This class is implemented separately from general polyhedra, 46 // because the simplex geometry can be compute 43 // because the simplex geometry can be computed very quickly, 47 // which may become important in situations im 44 // which may become important in situations imported from mesh generators, 48 // in which a very large number of G4Tets are 45 // in which a very large number of G4Tets are created. 49 // A Tet has all of its geometrical informatio 46 // A Tet has all of its geometrical information precomputed 50 // 47 // 51 G4UTet::G4UTet(const G4String& pName, 48 G4UTet::G4UTet(const G4String& pName, 52 const G4ThreeVector& anchor, << 49 G4ThreeVector anchor, 53 const G4ThreeVector& p1, << 50 G4ThreeVector p2, 54 const G4ThreeVector& p2, << 51 G4ThreeVector p3, 55 const G4ThreeVector& p3, G4bool << 52 G4ThreeVector p4, G4bool* degeneracyFlag) 56 : Base_t(pName, U3Vector(anchor.x(),anchor.y << 53 : G4USolid(pName, new UTet(pName, 57 U3Vector(p1.x(), p1.y(), p1. << 54 UVector3(anchor.x(),anchor.y(),anchor.z()), 58 U3Vector(p2.x(), p2.y(), p2. << 55 UVector3(p2.x(), p2.y(), p2.z()), 59 U3Vector(p3.x(), p3.y(), p3. << 56 UVector3(p3.x(), p3.y(), p3.z()), >> 57 UVector3(p4.x(), p4.y(), p4.z()), >> 58 degeneracyFlag)) 60 { 59 { 61 // Check for degeneracy << 62 G4bool degenerate = CheckDegeneracy(anchor, << 63 if(degeneracyFlag != nullptr) *degeneracyFla << 64 else if (degenerate) << 65 { << 66 G4Exception("G4UTet::G4UTet()", "GeomSolid << 67 "Degenerate tetrahedron not al << 68 } << 69 << 70 // Set bounding box << 71 for (G4int i = 0; i < 3; ++i) << 72 { << 73 fBmin[i] = std::min(std::min(std::min(anch << 74 fBmax[i] = std::max(std::max(std::max(anch << 75 } << 76 } 60 } 77 61 78 ////////////////////////////////////////////// 62 ////////////////////////////////////////////////////////////////////////// 79 // 63 // 80 // Fake default constructor - sets only member 64 // Fake default constructor - sets only member data and allocates memory 81 // for usage restri 65 // for usage restricted to object persistency. 82 // 66 // 83 G4UTet::G4UTet( __void__& a ) 67 G4UTet::G4UTet( __void__& a ) 84 : Base_t(a) << 68 : G4USolid(a) 85 { 69 { 86 } 70 } 87 71 88 ////////////////////////////////////////////// 72 ////////////////////////////////////////////////////////////////////////// 89 // 73 // 90 // Destructor 74 // Destructor 91 // 75 // 92 G4UTet::~G4UTet() = default; << 76 G4UTet::~G4UTet() >> 77 { >> 78 } 93 79 94 ////////////////////////////////////////////// 80 /////////////////////////////////////////////////////////////////////////////// 95 // 81 // 96 // Copy constructor 82 // Copy constructor 97 // 83 // 98 G4UTet::G4UTet(const G4UTet& rhs) 84 G4UTet::G4UTet(const G4UTet& rhs) 99 : Base_t(rhs) << 85 : G4USolid(rhs) 100 { 86 { 101 fBmin = rhs.fBmin; << 102 fBmax = rhs.fBmax; << 103 } 87 } 104 88 105 89 106 ////////////////////////////////////////////// 90 /////////////////////////////////////////////////////////////////////////////// 107 // 91 // 108 // Assignment operator 92 // Assignment operator 109 // 93 // 110 G4UTet& G4UTet::operator = (const G4UTet& rhs) << 94 G4UTet& G4UTet::operator = (const G4UTet& rhs) 111 { << 112 // Check assignment to self << 113 if (this == &rhs) { return *this; } << 114 << 115 // Copy base class data << 116 Base_t::operator=(rhs); << 117 << 118 // Copy bounding box << 119 fBmin = rhs.fBmin; << 120 fBmax = rhs.fBmax; << 121 << 122 return *this; << 123 } << 124 << 125 ////////////////////////////////////////////// << 126 // << 127 // Return true if tetrahedron is degenerate << 128 // Tetrahedron is concidered as degenerate in << 129 // height is less than the degeneracy toleranc << 130 // << 131 G4bool G4UTet::CheckDegeneracy(const G4ThreeVe << 132 const G4ThreeVe << 133 const G4ThreeVe << 134 const G4ThreeVe << 135 { << 136 G4double hmin = 4. * kCarTolerance; // degen << 137 << 138 // Calculate volume << 139 G4double vol = std::abs((p1 - p0).cross(p2 - << 140 << 141 // Calculate face areas squared << 142 G4double ss[4]; << 143 ss[0] = ((p1 - p0).cross(p2 - p0)).mag2(); << 144 ss[1] = ((p2 - p0).cross(p3 - p0)).mag2(); << 145 ss[2] = ((p3 - p0).cross(p1 - p0)).mag2(); << 146 ss[3] = ((p2 - p1).cross(p3 - p1)).mag2(); << 147 << 148 // Find face with max area << 149 G4int k = 0; << 150 for (G4int i = 1; i < 4; ++i) { if (ss[i] > << 151 << 152 // Check: vol^2 / s^2 <= hmin^2 << 153 return (vol*vol <= ss[k]*hmin*hmin); << 154 } << 155 << 156 ////////////////////////////////////////////// << 157 // << 158 // Dispatch to parameterisation for replicatio << 159 // computation & modification. << 160 // << 161 void G4UTet::ComputeDimensions(G4VPVParameteri << 162 const G4int, << 163 const G4VPhysic << 164 { << 165 } << 166 << 167 ////////////////////////////////////////////// << 168 // << 169 // Make a clone of the object << 170 // << 171 G4VSolid* G4UTet::Clone() const << 172 { << 173 return new G4UTet(*this); << 174 } << 175 << 176 ////////////////////////////////////////////// << 177 // << 178 // Modifier << 179 // << 180 void G4UTet::SetVertices(const G4ThreeVector& << 181 const G4ThreeVector& << 182 const G4ThreeVector& << 183 const G4ThreeVector& << 184 G4bool* degeneracyFla << 185 { << 186 // Check for degeneracy << 187 G4bool degenerate = CheckDegeneracy(anchor, << 188 if(degeneracyFlag != nullptr) *degeneracyFla << 189 else if (degenerate) << 190 { << 191 G4Exception("G4UTet::SetVertices()", "Geom << 192 "Degenerate tetrahedron not al << 193 } << 194 << 195 // Change tetrahedron << 196 *this = G4UTet(GetName(), anchor, p1, p2, p3 << 197 } << 198 << 199 ////////////////////////////////////////////// << 200 // << 201 // Accessors << 202 // << 203 void G4UTet::GetVertices(G4ThreeVector& anchor << 204 G4ThreeVector& p1, << 205 G4ThreeVector& p2, << 206 G4ThreeVector& p3) co << 207 { 95 { 208 std::vector<U3Vector> vec(4); << 96 // Check assignment to self 209 Base_t::GetVertices(vec[0], vec[1], vec[2], << 97 // 210 anchor = G4ThreeVector(vec[0].x(), vec[0].y( << 98 if (this == &rhs) { return *this; } 211 p1 = G4ThreeVector(vec[1].x(), vec[1].y(), v << 212 p2 = G4ThreeVector(vec[2].x(), vec[2].y(), v << 213 p3 = G4ThreeVector(vec[3].x(), vec[3].y(), v << 214 } << 215 99 216 std::vector<G4ThreeVector> G4UTet::GetVertices << 100 // Copy base class data 217 { << 101 // 218 std::vector<U3Vector> vec(4); << 102 G4USolid::operator=(rhs); 219 Base_t::GetVertices(vec[0], vec[1], vec[2], << 220 std::vector<G4ThreeVector> vertices; << 221 for (unsigned int i=0; i<4; ++i) << 222 { << 223 G4ThreeVector v(vec[i].x(), vec[i].y(), ve << 224 vertices.push_back(v); << 225 } << 226 return vertices; << 227 } << 228 103 229 ////////////////////////////////////////////// << 104 return *this; 230 // << 231 // Set bounding box << 232 // << 233 void G4UTet::SetBoundingLimits(const G4ThreeVe << 234 const G4ThreeVe << 235 { << 236 G4ThreeVector fVertex[4]; << 237 GetVertices(fVertex[0], fVertex[1], fVertex[ << 238 << 239 G4int iout[4] = { 0, 0, 0, 0 }; << 240 for (G4int i = 0; i < 4; ++i) << 241 { << 242 iout[i] = (G4int)(fVertex[i].x() < pMin.x( << 243 fVertex[i].y() < pMin.y( << 244 fVertex[i].z() < pMin.z( << 245 fVertex[i].x() > pMax.x( << 246 fVertex[i].y() > pMax.y( << 247 fVertex[i].z() > pMax.z( << 248 } << 249 if (iout[0] + iout[1] + iout[2] + iout[3] != << 250 { << 251 std::ostringstream message; << 252 message << "Attempt to set bounding box th << 253 << GetName() << " !\n" << 254 << " Specified bounding box limit << 255 << " pmin: " << pMin << "\n" << 256 << " pmax: " << pMax << "\n" << 257 << " Tetrahedron vertices:\n" << 258 << " anchor " << fVertex[0] << << 259 << " p1 " << fVertex[1] << << 260 << " p2 " << fVertex[2] << << 261 << " p3 " << fVertex[3] << << 262 G4Exception("G4UTet::SetBoundingLimits()", << 263 FatalException, message); << 264 } << 265 fBmin = pMin; << 266 fBmax = pMax; << 267 } << 268 << 269 ////////////////////////////////////////////// << 270 // << 271 // Get bounding box << 272 << 273 void G4UTet::BoundingLimits(G4ThreeVector& pMi << 274 { << 275 pMin = fBmin; << 276 pMax = fBmax; << 277 } 105 } 278 << 279 ////////////////////////////////////////////// << 280 // << 281 // Calculate extent under transform and specif << 282 << 283 G4bool << 284 G4UTet::CalculateExtent(const EAxis pAxis, << 285 const G4VoxelLimits& p << 286 const G4AffineTransfor << 287 G4double& pMin, << 288 { << 289 G4ThreeVector bmin, bmax; << 290 << 291 // Check bounding box (bbox) << 292 // << 293 BoundingLimits(bmin,bmax); << 294 G4BoundingEnvelope bbox(bmin,bmax); << 295 << 296 // Use simple bounding-box to help in the ca << 297 // << 298 return bbox.CalculateExtent(pAxis,pVoxelLimi << 299 << 300 #if 0 << 301 // Precise extent computation (disabled by d << 302 // << 303 G4bool exist; << 304 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 305 { << 306 return exist = (pMin < pMax) ? true : fals << 307 } << 308 << 309 // Set bounding envelope (benv) and calculat << 310 // << 311 std::vector<G4ThreeVector> vec = GetVertices << 312 << 313 G4ThreeVectorList anchor(1); << 314 anchor[0] = vec[0]; << 315 << 316 G4ThreeVectorList base(3); << 317 base[0] = vec[1]; << 318 base[1] = vec[2]; << 319 base[2] = vec[3]; << 320 << 321 std::vector<const G4ThreeVectorList *> polyg << 322 polygons[0] = &anchor; << 323 polygons[1] = &base; << 324 << 325 G4BoundingEnvelope benv(bmin,bmax,polygons); << 326 return exists = benv.CalculateExtent(pAxis,p << 327 #endif << 328 } << 329 << 330 ////////////////////////////////////////////// << 331 // << 332 // CreatePolyhedron << 333 // << 334 G4Polyhedron* G4UTet::CreatePolyhedron() const << 335 { << 336 std::vector<U3Vector> vec(4); << 337 Base_t::GetVertices(vec[0], vec[1], vec[2], << 338 << 339 G4double xyz[4][3]; << 340 const G4int faces[4][4] = {{1,3,2,0},{1,4,3, << 341 for (unsigned int i=0; i<4; ++i) << 342 { << 343 xyz[i][0] = vec[i].x(); << 344 xyz[i][1] = vec[i].y(); << 345 xyz[i][2] = vec[i].z(); << 346 } << 347 << 348 auto ph = new G4Polyhedron; << 349 ph->createPolyhedron(4,4,xyz,faces); << 350 return ph; << 351 } << 352 << 353 #endif // G4GEOM_USE_USOLIDS << 354 106