Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // Implementation for G4UTet wrapper class 27 // 28 // 1.11.13 G.Cosmo, CERN 29 // -------------------------------------------------------------------- 30 31 #include "G4Tet.hh" 32 #include "G4UTet.hh" 33 34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G4GEOM_USE_PARTIAL_USOLIDS) ) 35 36 #include "G4AffineTransform.hh" 37 #include "G4VPVParameterisation.hh" 38 #include "G4BoundingEnvelope.hh" 39 40 using namespace CLHEP; 41 42 //////////////////////////////////////////////////////////////////////// 43 // 44 // Constructor - create a tetrahedron 45 // This class is implemented separately from general polyhedra, 46 // because the simplex geometry can be computed very quickly, 47 // which may become important in situations imported from mesh generators, 48 // in which a very large number of G4Tets are created. 49 // A Tet has all of its geometrical information precomputed 50 // 51 G4UTet::G4UTet(const G4String& pName, 52 const G4ThreeVector& anchor, 53 const G4ThreeVector& p1, 54 const G4ThreeVector& p2, 55 const G4ThreeVector& p3, G4bool* degeneracyFlag) 56 : Base_t(pName, U3Vector(anchor.x(),anchor.y(),anchor.z()), 57 U3Vector(p1.x(), p1.y(), p1.z()), 58 U3Vector(p2.x(), p2.y(), p2.z()), 59 U3Vector(p3.x(), p3.y(), p3.z())) 60 { 61 // Check for degeneracy 62 G4bool degenerate = CheckDegeneracy(anchor, p1, p2, p3); 63 if(degeneracyFlag != nullptr) *degeneracyFlag = degenerate; 64 else if (degenerate) 65 { 66 G4Exception("G4UTet::G4UTet()", "GeomSolids0002", FatalException, 67 "Degenerate tetrahedron not allowed."); 68 } 69 70 // Set bounding box 71 for (G4int i = 0; i < 3; ++i) 72 { 73 fBmin[i] = std::min(std::min(std::min(anchor[i], p1[i]), p2[i]), p3[i]); 74 fBmax[i] = std::max(std::max(std::max(anchor[i], p1[i]), p2[i]), p3[i]); 75 } 76 } 77 78 ////////////////////////////////////////////////////////////////////////// 79 // 80 // Fake default constructor - sets only member data and allocates memory 81 // for usage restricted to object persistency. 82 // 83 G4UTet::G4UTet( __void__& a ) 84 : Base_t(a) 85 { 86 } 87 88 ////////////////////////////////////////////////////////////////////////// 89 // 90 // Destructor 91 // 92 G4UTet::~G4UTet() = default; 93 94 /////////////////////////////////////////////////////////////////////////////// 95 // 96 // Copy constructor 97 // 98 G4UTet::G4UTet(const G4UTet& rhs) 99 : Base_t(rhs) 100 { 101 fBmin = rhs.fBmin; 102 fBmax = rhs.fBmax; 103 } 104 105 106 /////////////////////////////////////////////////////////////////////////////// 107 // 108 // Assignment operator 109 // 110 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 case if its minimal 129 // height is less than the degeneracy tolerance 130 // 131 G4bool G4UTet::CheckDegeneracy(const G4ThreeVector& p0, 132 const G4ThreeVector& p1, 133 const G4ThreeVector& p2, 134 const G4ThreeVector& p3) const 135 { 136 G4double hmin = 4. * kCarTolerance; // degeneracy tolerance 137 138 // Calculate volume 139 G4double vol = std::abs((p1 - p0).cross(p2 - p0).dot(p3 - p0)); 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] > ss[k]) k = 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 replication mechanism dimension 159 // computation & modification. 160 // 161 void G4UTet::ComputeDimensions(G4VPVParameterisation*, 162 const G4int, 163 const G4VPhysicalVolume*) 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& anchor, 181 const G4ThreeVector& p1, 182 const G4ThreeVector& p2, 183 const G4ThreeVector& p3, 184 G4bool* degeneracyFlag) 185 { 186 // Check for degeneracy 187 G4bool degenerate = CheckDegeneracy(anchor, p1, p2, p3); 188 if(degeneracyFlag != nullptr) *degeneracyFlag = degenerate; 189 else if (degenerate) 190 { 191 G4Exception("G4UTet::SetVertices()", "GeomSolids0002", FatalException, 192 "Degenerate tetrahedron not allowed."); 193 } 194 195 // Change tetrahedron 196 *this = G4UTet(GetName(), anchor, p1, p2, p3, °enerate); 197 } 198 199 /////////////////////////////////////////////////////////////////////////////// 200 // 201 // Accessors 202 // 203 void G4UTet::GetVertices(G4ThreeVector& anchor, 204 G4ThreeVector& p1, 205 G4ThreeVector& p2, 206 G4ThreeVector& p3) const 207 { 208 std::vector<U3Vector> vec(4); 209 Base_t::GetVertices(vec[0], vec[1], vec[2], vec[3]); 210 anchor = G4ThreeVector(vec[0].x(), vec[0].y(), vec[0].z()); 211 p1 = G4ThreeVector(vec[1].x(), vec[1].y(), vec[1].z()); 212 p2 = G4ThreeVector(vec[2].x(), vec[2].y(), vec[2].z()); 213 p3 = G4ThreeVector(vec[3].x(), vec[3].y(), vec[3].z()); 214 } 215 216 std::vector<G4ThreeVector> G4UTet::GetVertices() const 217 { 218 std::vector<U3Vector> vec(4); 219 Base_t::GetVertices(vec[0], vec[1], vec[2], vec[3]); 220 std::vector<G4ThreeVector> vertices; 221 for (unsigned int i=0; i<4; ++i) 222 { 223 G4ThreeVector v(vec[i].x(), vec[i].y(), vec[i].z()); 224 vertices.push_back(v); 225 } 226 return vertices; 227 } 228 229 //////////////////////////////////////////////////////////////////////// 230 // 231 // Set bounding box 232 // 233 void G4UTet::SetBoundingLimits(const G4ThreeVector& pMin, 234 const G4ThreeVector& pMax) 235 { 236 G4ThreeVector fVertex[4]; 237 GetVertices(fVertex[0], fVertex[1], fVertex[2], fVertex[3]); 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] != 0) 250 { 251 std::ostringstream message; 252 message << "Attempt to set bounding box that does not encapsulate solid: " 253 << GetName() << " !\n" 254 << " Specified bounding box limits:\n" 255 << " pmin: " << pMin << "\n" 256 << " pmax: " << pMax << "\n" 257 << " Tetrahedron vertices:\n" 258 << " anchor " << fVertex[0] << ((iout[0]) != 0 ? " is outside\n" : "\n") 259 << " p1 " << fVertex[1] << ((iout[1]) != 0 ? " is outside\n" : "\n") 260 << " p2 " << fVertex[2] << ((iout[2]) != 0 ? " is outside\n" : "\n") 261 << " p3 " << fVertex[3] << ((iout[3]) != 0 ? " is outside" : ""); 262 G4Exception("G4UTet::SetBoundingLimits()", "GeomSolids0002", 263 FatalException, message); 264 } 265 fBmin = pMin; 266 fBmax = pMax; 267 } 268 269 ////////////////////////////////////////////////////////////////////////// 270 // 271 // Get bounding box 272 273 void G4UTet::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const 274 { 275 pMin = fBmin; 276 pMax = fBmax; 277 } 278 279 ////////////////////////////////////////////////////////////////////////// 280 // 281 // Calculate extent under transform and specified limit 282 283 G4bool 284 G4UTet::CalculateExtent(const EAxis pAxis, 285 const G4VoxelLimits& pVoxelLimit, 286 const G4AffineTransform& pTransform, 287 G4double& pMin, G4double& pMax) const 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 case of complex 3D meshes 297 // 298 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 299 300 #if 0 301 // Precise extent computation (disabled by default for this shape) 302 // 303 G4bool exist; 304 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax)) 305 { 306 return exist = (pMin < pMax) ? true : false; 307 } 308 309 // Set bounding envelope (benv) and calculate extent 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 *> polygons(2); 322 polygons[0] = &anchor; 323 polygons[1] = &base; 324 325 G4BoundingEnvelope benv(bmin,bmax,polygons); 326 return exists = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 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], vec[3]); 338 339 G4double xyz[4][3]; 340 const G4int faces[4][4] = {{1,3,2,0},{1,4,3,0},{1,2,4,0},{2,3,4,0}}; 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