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 result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // Implementation of G4UPolyhedra wrapper clas << 27 // 26 // 28 // 31.10.13 G.Cosmo, CERN << 27 // $Id:$ >> 28 // >> 29 // Implementation of G4UPolycone wrapper class 29 // ------------------------------------------- 30 // -------------------------------------------------------------------- 30 31 31 #include "G4Polyhedra.hh" 32 #include "G4Polyhedra.hh" 32 #include "G4UPolyhedra.hh" 33 #include "G4UPolyhedra.hh" 33 << 34 #include "G4Polyhedron.hh" 34 #if ( defined(G4GEOM_USE_USOLIDS) || defined(G << 35 << 36 #include "G4GeomTools.hh" << 37 #include "G4GeometryTolerance.hh" << 38 #include "G4AffineTransform.hh" << 39 #include "G4VPVParameterisation.hh" 35 #include "G4VPVParameterisation.hh" 40 #include "G4BoundingEnvelope.hh" << 41 36 42 using namespace CLHEP; << 37 using CLHEP::twopi; 43 38 44 ////////////////////////////////////////////// 39 //////////////////////////////////////////////////////////////////////// 45 // 40 // 46 // Constructor (GEANT3 style parameters) 41 // Constructor (GEANT3 style parameters) 47 // 42 // 48 // GEANT3 PGON radii are specified in the dist 43 // GEANT3 PGON radii are specified in the distance to the norm of each face. 49 // << 44 // 50 G4UPolyhedra::G4UPolyhedra(const G4String& nam << 45 G4UPolyhedra::G4UPolyhedra(const G4String& name, 51 G4double phiS 46 G4double phiStart, 52 G4double phiT 47 G4double phiTotal, 53 G4int numSide << 48 G4int numSide, 54 G4int numZPla 49 G4int numZPlanes, 55 const G4double zPla 50 const G4double zPlane[], 56 const G4double rInn 51 const G4double rInner[], 57 const G4double rOut 52 const G4double rOuter[] ) 58 : Base_t(name, phiStart, phiTotal, numSide, << 53 : G4USolid(name, new UPolyhedra(name,phiStart, phiTotal, numSide, 59 numZPlanes, zPlane, rInner, rOuter) << 54 numZPlanes, zPlane, rInner, rOuter)) 60 { 55 { 61 fGenericPgon = false; << 62 SetOriginalParameters(); << 63 wrStart = phiStart; << 64 while (wrStart < 0) << 65 { << 66 wrStart += twopi; << 67 } << 68 wrDelta = phiTotal; << 69 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL << 70 { << 71 wrDelta = twopi; << 72 } << 73 wrNumSide = numSide; << 74 G4double convertRad = 1./std::cos(0.5*wrDelt << 75 rzcorners.resize(0); << 76 for (G4int i=0; i<numZPlanes; ++i) << 77 { << 78 G4double z = zPlane[i]; << 79 G4double r = rOuter[i]*convertRad; << 80 rzcorners.emplace_back(r,z); << 81 } << 82 for (G4int i=numZPlanes-1; i>=0; --i) << 83 { << 84 G4double z = zPlane[i]; << 85 G4double r = rInner[i]*convertRad; << 86 rzcorners.emplace_back(r,z); << 87 } << 88 std::vector<G4int> iout; << 89 G4GeomTools::RemoveRedundantVertices(rzcorne << 90 } 56 } 91 57 92 58 93 ////////////////////////////////////////////// 59 //////////////////////////////////////////////////////////////////////// 94 // 60 // 95 // Constructor (generic parameters) 61 // Constructor (generic parameters) 96 // 62 // 97 G4UPolyhedra::G4UPolyhedra(const G4String& nam << 63 G4UPolyhedra::G4UPolyhedra(const G4String& name, 98 G4double phiS 64 G4double phiStart, 99 G4double phiT 65 G4double phiTotal, 100 G4int numS << 66 G4int numSide, 101 G4int numR 67 G4int numRZ, 102 const G4double r[], 68 const G4double r[], 103 const G4double z[] 69 const G4double z[] ) 104 : Base_t(name, phiStart, phiTotal, numSide, << 70 : G4USolid(name, new UPolyhedra(name, phiStart, phiTotal, numSide, 105 { << 71 numRZ, r, z)) 106 fGenericPgon = true; << 72 { 107 SetOriginalParameters(); << 108 wrStart = phiStart; << 109 while (wrStart < 0.) << 110 { << 111 wrStart += twopi; << 112 } << 113 wrDelta = phiTotal; << 114 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL << 115 { << 116 wrDelta = twopi; << 117 } << 118 wrNumSide = numSide; << 119 rzcorners.resize(0); << 120 for (G4int i=0; i<numRZ; ++i) << 121 { << 122 rzcorners.emplace_back(r[i],z[i]); << 123 } << 124 std::vector<G4int> iout; << 125 G4GeomTools::RemoveRedundantVertices(rzcorne << 126 } 73 } 127 74 128 75 129 ////////////////////////////////////////////// 76 //////////////////////////////////////////////////////////////////////// 130 // 77 // 131 // Fake default constructor - sets only member 78 // Fake default constructor - sets only member data and allocates memory 132 // for usage restri 79 // for usage restricted to object persistency. 133 // 80 // 134 G4UPolyhedra::G4UPolyhedra( __void__& a ) 81 G4UPolyhedra::G4UPolyhedra( __void__& a ) 135 : Base_t(a) << 82 : G4USolid(a) 136 { 83 { 137 } 84 } 138 85 139 86 140 ////////////////////////////////////////////// 87 //////////////////////////////////////////////////////////////////////// 141 // 88 // 142 // Destructor 89 // Destructor 143 // 90 // 144 G4UPolyhedra::~G4UPolyhedra() = default; << 91 G4UPolyhedra::~G4UPolyhedra() >> 92 { >> 93 } 145 94 146 95 147 ////////////////////////////////////////////// 96 //////////////////////////////////////////////////////////////////////// 148 // 97 // 149 // Copy constructor 98 // Copy constructor 150 // 99 // 151 G4UPolyhedra::G4UPolyhedra( const G4UPolyhedra << 100 G4UPolyhedra::G4UPolyhedra( const G4UPolyhedra &source ) 152 : Base_t( source ) << 101 : G4USolid( source ) 153 { 102 { 154 fGenericPgon = source.fGenericPgon; << 155 fOriginalParameters = source.fOriginalParame << 156 wrStart = source.wrStart; << 157 wrDelta = source.wrDelta; << 158 wrNumSide = source.wrNumSide; << 159 rzcorners = source.rzcorners; << 160 } 103 } 161 104 162 105 163 ////////////////////////////////////////////// 106 //////////////////////////////////////////////////////////////////////// 164 // 107 // 165 // Assignment operator 108 // Assignment operator 166 // 109 // 167 G4UPolyhedra& G4UPolyhedra::operator=( const G << 110 G4UPolyhedra& G4UPolyhedra::operator=( const G4UPolyhedra &source ) 168 { 111 { 169 if (this == &source) return *this; 112 if (this == &source) return *this; 170 113 171 Base_t::operator=( source ); << 114 G4USolid::operator=( source ); 172 fGenericPgon = source.fGenericPgon; << 115 173 fOriginalParameters = source.fOriginalParame << 174 wrStart = source.wrStart; << 175 wrDelta = source.wrDelta; << 176 wrNumSide = source.wrNumSide; << 177 rzcorners = source.rzcorners; << 178 << 179 return *this; 116 return *this; 180 } 117 } 181 118 182 119 183 ////////////////////////////////////////////// 120 //////////////////////////////////////////////////////////////////////// 184 // 121 // 185 // Accessors & modifiers << 186 // << 187 G4int G4UPolyhedra::GetNumSide() const << 188 { << 189 return wrNumSide; << 190 } << 191 G4double G4UPolyhedra::GetStartPhi() const << 192 { << 193 return wrStart; << 194 } << 195 G4double G4UPolyhedra::GetEndPhi() const << 196 { << 197 return (wrStart + wrDelta); << 198 } << 199 G4double G4UPolyhedra::GetSinStartPhi() const << 200 { << 201 G4double phi = GetStartPhi(); << 202 return std::sin(phi); << 203 } << 204 G4double G4UPolyhedra::GetCosStartPhi() const << 205 { << 206 G4double phi = GetStartPhi(); << 207 return std::cos(phi); << 208 } << 209 G4double G4UPolyhedra::GetSinEndPhi() const << 210 { << 211 G4double phi = GetEndPhi(); << 212 return std::sin(phi); << 213 } << 214 G4double G4UPolyhedra::GetCosEndPhi() const << 215 { << 216 G4double phi = GetEndPhi(); << 217 return std::cos(phi); << 218 } << 219 G4bool G4UPolyhedra::IsOpen() const << 220 { << 221 return (wrDelta < twopi); << 222 } << 223 G4bool G4UPolyhedra::IsGeneric() const << 224 { << 225 return fGenericPgon; << 226 } << 227 G4int G4UPolyhedra::GetNumRZCorner() const << 228 { << 229 return rzcorners.size(); << 230 } << 231 G4PolyhedraSideRZ G4UPolyhedra::GetCorner(G4in << 232 { << 233 G4TwoVector rz = rzcorners.at(index); << 234 G4PolyhedraSideRZ psiderz = { rz.x(), rz.y() << 235 << 236 return psiderz; << 237 } << 238 G4PolyhedraHistorical* G4UPolyhedra::GetOrigin << 239 { << 240 return new G4PolyhedraHistorical(fOriginalPa << 241 } << 242 void G4UPolyhedra::SetOriginalParameters() << 243 { << 244 G4double startPhi = GetPhiStart(); << 245 G4double deltaPhi = GetPhiDelta(); << 246 G4int numPlanes = GetZSegmentCount() + 1; << 247 G4int numSides = GetSideCount(); << 248 << 249 fOriginalParameters.Start_angle = startPhi << 250 fOriginalParameters.Opening_angle = deltaPhi << 251 fOriginalParameters.Num_z_planes = numPlane << 252 fOriginalParameters.numSide = numSides << 253 << 254 delete [] fOriginalParameters.Z_values; << 255 delete [] fOriginalParameters.Rmin; << 256 delete [] fOriginalParameters.Rmax; << 257 fOriginalParameters.Z_values = new G4double[ << 258 fOriginalParameters.Rmin = new G4double[numP << 259 fOriginalParameters.Rmax = new G4double[numP << 260 << 261 G4double convertRad = fGenericPgon << 262 ? 1.0 : std::cos(0.5*del << 263 for (G4int i=0; i<numPlanes; ++i) << 264 { << 265 fOriginalParameters.Z_values[i] = GetZPlan << 266 fOriginalParameters.Rmax[i] = GetRMax( << 267 fOriginalParameters.Rmin[i] = GetRMin( << 268 } << 269 } << 270 void G4UPolyhedra::SetOriginalParameters(G4Pol << 271 { << 272 fOriginalParameters = *pars; << 273 fRebuildPolyhedron = true; << 274 Reset(); << 275 } << 276 << 277 G4bool G4UPolyhedra::Reset() << 278 { << 279 if (fGenericPgon) << 280 { << 281 std::ostringstream message; << 282 message << "Solid " << GetName() << " buil << 283 << G4endl << "Not applicable to th << 284 G4Exception("G4UPolyhedra::Reset()", "Geom << 285 JustWarning, message, "Paramet << 286 return true; // error code set << 287 } << 288 << 289 // << 290 // Rebuild polyhedra based on original param << 291 // << 292 wrStart = fOriginalParameters.Start_angle; << 293 while (wrStart < 0.) << 294 { << 295 wrStart += twopi; << 296 } << 297 wrDelta = fOriginalParameters.Opening_angle; << 298 if (wrDelta <= 0. || wrDelta >= twopi*(1-DBL << 299 { << 300 wrDelta = twopi; << 301 } << 302 wrNumSide = fOriginalParameters.numSide; << 303 rzcorners.resize(0); << 304 for (G4int i=0; i<fOriginalParameters.Num_z_ << 305 { << 306 G4double z = fOriginalParameters.Z_values[ << 307 G4double r = fOriginalParameters.Rmax[i]; << 308 rzcorners.emplace_back(r,z); << 309 } << 310 for (G4int i=fOriginalParameters.Num_z_plane << 311 { << 312 G4double z = fOriginalParameters.Z_values[ << 313 G4double r = fOriginalParameters.Rmin[i]; << 314 rzcorners.emplace_back(r,z); << 315 } << 316 std::vector<G4int> iout; << 317 G4GeomTools::RemoveRedundantVertices(rzcorne << 318 << 319 return false; // error code unset << 320 } << 321 << 322 << 323 ////////////////////////////////////////////// << 324 // << 325 // Dispatch to parameterisation for replicatio 122 // Dispatch to parameterisation for replication mechanism dimension 326 // computation & modification. 123 // computation & modification. 327 // 124 // 328 void G4UPolyhedra::ComputeDimensions(G4VPVPara 125 void G4UPolyhedra::ComputeDimensions(G4VPVParameterisation* p, 329 const G4i 126 const G4int n, 330 const G4V 127 const G4VPhysicalVolume* pRep) 331 { 128 { 332 p->ComputeDimensions(*(G4Polyhedra*)this,n,p 129 p->ComputeDimensions(*(G4Polyhedra*)this,n,pRep); 333 } 130 } 334 131 335 132 336 ////////////////////////////////////////////// << 133 //////////////////////////////////////////////////////////////////////// 337 // 134 // 338 // Make a clone of the object << 135 // CreatePolyhedron 339 << 340 G4VSolid* G4UPolyhedra::Clone() const << 341 { << 342 return new G4UPolyhedra(*this); << 343 } << 344 << 345 << 346 ////////////////////////////////////////////// << 347 // 136 // 348 // Get bounding box << 137 G4Polyhedron* G4UPolyhedra::CreatePolyhedron() const 349 << 350 void G4UPolyhedra::BoundingLimits(G4ThreeVecto << 351 G4ThreeVecto << 352 { 138 { 353 static G4bool checkBBox = true; << 139 if (!IsGeneric()) 354 static G4bool checkPhi = true; << 355 << 356 G4double rmin = kInfinity, rmax = -kInfinity << 357 G4double zmin = kInfinity, zmax = -kInfinity << 358 for (G4int i=0; i<GetNumRZCorner(); ++i) << 359 { << 360 G4PolyhedraSideRZ corner = GetCorner(i); << 361 if (corner.r < rmin) rmin = corner.r; << 362 if (corner.r > rmax) rmax = corner.r; << 363 if (corner.z < zmin) zmin = corner.z; << 364 if (corner.z > zmax) zmax = corner.z; << 365 } << 366 << 367 G4double sphi = GetStartPhi(); << 368 G4double ephi = GetEndPhi(); << 369 G4double dphi = IsOpen() ? ephi-sphi : tw << 370 G4int ksteps = GetNumSide(); << 371 G4double astep = dphi/ksteps; << 372 G4double sinStep = std::sin(astep); << 373 G4double cosStep = std::cos(astep); << 374 << 375 G4double sinCur = GetSinStartPhi(); << 376 G4double cosCur = GetCosStartPhi(); << 377 if (!IsOpen()) rmin = 0.; << 378 G4double xmin = rmin*cosCur, xmax = xmin; << 379 G4double ymin = rmin*sinCur, ymax = ymin; << 380 for (G4int k=0; k<ksteps+1; ++k) << 381 { 140 { 382 G4double x = rmax*cosCur; << 141 G4PolyhedraHistorical* original_parameters = GetOriginalParameters(); 383 if (x < xmin) xmin = x; << 142 G4PolyhedronPgon* 384 if (x > xmax) xmax = x; << 143 polyhedron = new G4PolyhedronPgon( GetOriginalParameters()->Start_angle, 385 G4double y = rmax*sinCur; << 144 GetOriginalParameters()->Opening_angle, 386 if (y < ymin) ymin = y; << 145 GetOriginalParameters()->numSide, 387 if (y > ymax) ymax = y; << 146 GetOriginalParameters()->Num_z_planes, 388 if (rmin > 0.) << 147 GetOriginalParameters()->Z_values, >> 148 GetOriginalParameters()->Rmin, >> 149 GetOriginalParameters()->Rmax); >> 150 delete original_parameters; // delete local copy >> 151 return polyhedron; >> 152 } >> 153 else >> 154 { >> 155 // The following code prepares for: >> 156 // HepPolyhedron::createPolyhedron(int Nnodes, int Nfaces, >> 157 // const double xyz[][3], >> 158 // const int faces_vec[][4]) >> 159 // Here is an extract from the header file HepPolyhedron.h: >> 160 /** >> 161 * Creates user defined polyhedron. >> 162 * This function allows to the user to define arbitrary polyhedron. >> 163 * The faces of the polyhedron should be either triangles or planar >> 164 * quadrilateral. Nodes of a face are defined by indexes pointing to >> 165 * the elements in the xyz array. Numeration of the elements in the >> 166 * array starts from 1 (like in fortran). The indexes can be positive >> 167 * or negative. Negative sign means that the corresponding edge is >> 168 * invisible. The normal of the face should be directed to exterior >> 169 * of the polyhedron. >> 170 * >> 171 * @param Nnodes number of nodes >> 172 * @param Nfaces number of faces >> 173 * @param xyz nodes >> 174 * @param faces_vec faces (quadrilaterals or triangles) >> 175 * @return status of the operation - is non-zero in case of problem >> 176 */ >> 177 G4int nNodes; >> 178 G4int nFaces; >> 179 typedef G4double double3[3]; >> 180 double3* xyz; >> 181 typedef G4int int4[4]; >> 182 int4* faces_vec; >> 183 if (IsOpen()) 389 { 184 { 390 G4double xx = rmin*cosCur; << 185 // Triangulate open ends. Simple ear-chopping algorithm... 391 if (xx < xmin) xmin = xx; << 186 // I'm not sure how robust this algorithm is (J.Allison). 392 if (xx > xmax) xmax = xx; << 187 // 393 G4double yy = rmin*sinCur; << 188 std::vector<G4bool> chopped(GetNumRZCorner(), false); 394 if (yy < ymin) ymin = yy; << 189 std::vector<G4int*> triQuads; 395 if (yy > ymax) ymax = yy; << 190 G4int remaining = GetNumRZCorner(); >> 191 G4int iStarter = 0; >> 192 while (remaining >= 3) >> 193 { >> 194 // Find unchopped corners... >> 195 // >> 196 G4int A = -1, B = -1, C = -1; >> 197 G4int iStepper = iStarter; >> 198 do >> 199 { >> 200 if (A < 0) { A = iStepper; } >> 201 else if (B < 0) { B = iStepper; } >> 202 else if (C < 0) { C = iStepper; } >> 203 do >> 204 { >> 205 if (++iStepper >= GetNumRZCorner()) iStepper = 0; >> 206 } >> 207 while (chopped[iStepper]); >> 208 } >> 209 while (C < 0 && iStepper != iStarter); >> 210 >> 211 // Check triangle at B is pointing outward (an "ear"). >> 212 // Sign of z cross product determines... >> 213 >> 214 G4double BAr = GetCorner(A).r - GetCorner(B).r; >> 215 G4double BAz = GetCorner(A).z - GetCorner(B).z; >> 216 G4double BCr = GetCorner(C).r - GetCorner(B).r; >> 217 G4double BCz = GetCorner(C).z - GetCorner(B).z; >> 218 if (BAr * BCz - BAz * BCr < kCarTolerance) >> 219 { >> 220 G4int* tq = new G4int[3]; >> 221 tq[0] = A + 1; >> 222 tq[1] = B + 1; >> 223 tq[2] = C + 1; >> 224 triQuads.push_back(tq); >> 225 chopped[B] = true; >> 226 --remaining; >> 227 } >> 228 else >> 229 { >> 230 do >> 231 { >> 232 if (++iStarter >= GetNumRZCorner()) { iStarter = 0; } >> 233 } >> 234 while (chopped[iStarter]); >> 235 } >> 236 } >> 237 >> 238 // Transfer to faces... >> 239 G4int numSide=GetNumSide(); >> 240 nNodes = (numSide + 1) * GetNumRZCorner(); >> 241 nFaces = numSide * GetNumRZCorner() + 2 * triQuads.size(); >> 242 faces_vec = new int4[nFaces]; >> 243 G4int iface = 0; >> 244 G4int addition = GetNumRZCorner() * numSide; >> 245 G4int d = GetNumRZCorner() - 1; >> 246 for (G4int iEnd = 0; iEnd < 2; ++iEnd) >> 247 { >> 248 for (size_t i = 0; i < triQuads.size(); ++i) >> 249 { >> 250 // Negative for soft/auxiliary/normally invisible edges... >> 251 // >> 252 G4int a, b, c; >> 253 if (iEnd == 0) >> 254 { >> 255 a = triQuads[i][0]; >> 256 b = triQuads[i][1]; >> 257 c = triQuads[i][2]; >> 258 } >> 259 else >> 260 { >> 261 a = triQuads[i][0] + addition; >> 262 b = triQuads[i][2] + addition; >> 263 c = triQuads[i][1] + addition; >> 264 } >> 265 G4int ab = std::abs(b - a); >> 266 G4int bc = std::abs(c - b); >> 267 G4int ca = std::abs(a - c); >> 268 faces_vec[iface][0] = (ab == 1 || ab == d)? a: -a; >> 269 faces_vec[iface][1] = (bc == 1 || bc == d)? b: -b; >> 270 faces_vec[iface][2] = (ca == 1 || ca == d)? c: -c; >> 271 faces_vec[iface][3] = 0; >> 272 ++iface; >> 273 } >> 274 } >> 275 >> 276 // Continue with sides... >> 277 >> 278 xyz = new double3[nNodes]; >> 279 const G4double dPhi = (GetEndPhi() - GetStartPhi()) / numSide; >> 280 G4double phi = GetStartPhi(); >> 281 G4int ixyz = 0; >> 282 for (G4int iSide = 0; iSide < numSide; ++iSide) >> 283 { >> 284 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner) >> 285 { >> 286 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi); >> 287 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi); >> 288 xyz[ixyz][2] = GetCorner(iCorner).z; >> 289 if (iCorner < GetNumRZCorner() - 1) >> 290 { >> 291 faces_vec[iface][0] = ixyz + 1; >> 292 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1; >> 293 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2; >> 294 faces_vec[iface][3] = ixyz + 2; >> 295 } >> 296 else >> 297 { >> 298 faces_vec[iface][0] = ixyz + 1; >> 299 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1; >> 300 faces_vec[iface][2] = ixyz + 2; >> 301 faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2; >> 302 } >> 303 ++iface; >> 304 ++ixyz; >> 305 } >> 306 phi += dPhi; >> 307 } >> 308 >> 309 // Last GetCorner... >> 310 >> 311 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner) >> 312 { >> 313 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi); >> 314 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi); >> 315 xyz[ixyz][2] = GetCorner(iCorner).z; >> 316 ++ixyz; >> 317 } 396 } 318 } 397 G4double sinTmp = sinCur; << 319 else // !phiIsOpen - i.e., a complete 360 degrees. 398 sinCur = sinCur*cosStep + cosCur*sinStep; << 399 cosCur = cosCur*cosStep - sinTmp*sinStep; << 400 } << 401 pMin.set(xmin,ymin,zmin); << 402 pMax.set(xmax,ymax,zmax); << 403 << 404 // Check correctness of the bounding box << 405 // << 406 if (pMin.x() >= pMax.x() || pMin.y() >= pMax << 407 { << 408 std::ostringstream message; << 409 message << "Bad bounding box (min >= max) << 410 << GetName() << " !" << 411 << "\npMin = " << pMin << 412 << "\npMax = " << pMax; << 413 G4Exception("G4UPolyhedra::BoundingLimits( << 414 JustWarning, message); << 415 StreamInfo(G4cout); << 416 } << 417 << 418 // Check consistency of bounding boxes << 419 // << 420 if (checkBBox) << 421 { << 422 U3Vector vmin, vmax; << 423 Extent(vmin,vmax); << 424 if (std::abs(pMin.x()-vmin.x()) > kCarTole << 425 std::abs(pMin.y()-vmin.y()) > kCarTole << 426 std::abs(pMin.z()-vmin.z()) > kCarTole << 427 std::abs(pMax.x()-vmax.x()) > kCarTole << 428 std::abs(pMax.y()-vmax.y()) > kCarTole << 429 std::abs(pMax.z()-vmax.z()) > kCarTole << 430 { 320 { 431 std::ostringstream message; << 321 nNodes = GetNumSide() * GetNumRZCorner(); 432 message << "Inconsistency in bounding bo << 322 nFaces = GetNumSide() * GetNumRZCorner();; 433 << GetName() << " !" << 323 xyz = new double3[nNodes]; 434 << "\nBBox min: wrapper = " << p << 324 faces_vec = new int4[nFaces]; 435 << "\nBBox max: wrapper = " << p << 325 // const G4double dPhi = (endPhi - startPhi) / numSide; 436 G4Exception("G4UPolyhedra::BoundingLimit << 326 const G4double dPhi = twopi / GetNumSide(); // !phiIsOpen endPhi-startPhi = 360 degrees. 437 JustWarning, message); << 327 G4double phi = GetStartPhi(); 438 checkBBox = false; << 328 G4int ixyz = 0, iface = 0; >> 329 for (G4int iSide = 0; iSide < GetNumSide(); ++iSide) >> 330 { >> 331 for (G4int iCorner = 0; iCorner < GetNumRZCorner(); ++iCorner) >> 332 { >> 333 xyz[ixyz][0] = GetCorner(iCorner).r * std::cos(phi); >> 334 xyz[ixyz][1] = GetCorner(iCorner).r * std::sin(phi); >> 335 xyz[ixyz][2] = GetCorner(iCorner).z; >> 336 if (iSide < GetNumSide() - 1) >> 337 { >> 338 if (iCorner < GetNumRZCorner() - 1) >> 339 { >> 340 faces_vec[iface][0] = ixyz + 1; >> 341 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1; >> 342 faces_vec[iface][2] = ixyz + GetNumRZCorner() + 2; >> 343 faces_vec[iface][3] = ixyz + 2; >> 344 } >> 345 else >> 346 { >> 347 faces_vec[iface][0] = ixyz + 1; >> 348 faces_vec[iface][1] = ixyz + GetNumRZCorner() + 1; >> 349 faces_vec[iface][2] = ixyz + 2; >> 350 faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2; >> 351 } >> 352 } >> 353 else // Last side joins ends... >> 354 { >> 355 if (iCorner < GetNumRZCorner() - 1) >> 356 { >> 357 faces_vec[iface][0] = ixyz + 1; >> 358 faces_vec[iface][1] = ixyz + GetNumRZCorner() - nFaces + 1; >> 359 faces_vec[iface][2] = ixyz + GetNumRZCorner() - nFaces + 2; >> 360 faces_vec[iface][3] = ixyz + 2; >> 361 } >> 362 else >> 363 { >> 364 faces_vec[iface][0] = ixyz + 1; >> 365 faces_vec[iface][1] = ixyz - nFaces + GetNumRZCorner() + 1; >> 366 faces_vec[iface][2] = ixyz - nFaces + 2; >> 367 faces_vec[iface][3] = ixyz - GetNumRZCorner() + 2; >> 368 } >> 369 } >> 370 ++ixyz; >> 371 ++iface; >> 372 } >> 373 phi += dPhi; >> 374 } 439 } 375 } 440 } << 376 G4Polyhedron* polyhedron = new G4Polyhedron; 441 << 377 G4int problem = polyhedron->createPolyhedron(nNodes, nFaces, xyz, faces_vec); 442 // Check consistency of angles << 378 delete [] faces_vec; 443 // << 379 delete [] xyz; 444 if (checkPhi) << 380 if (problem) 445 { << 446 if (GetStartPhi() != GetPhiStart() || << 447 GetEndPhi() != GetPhiEnd() || << 448 GetNumSide() != GetSideCount() || << 449 IsOpen() != (Base_t::GetPhiDelta( << 450 { 381 { 451 std::ostringstream message; 382 std::ostringstream message; 452 message << "Inconsistency in Phi angles << 383 message << "Problem creating G4Polyhedron for: " << GetName(); 453 << GetName() << " !" << 384 G4Exception("G4Polyhedra::CreatePolyhedron()", "GeomSolids1002", 454 << "\nPhi start : wrapper = " < << 455 << " solid = " << GetPhiStar << 456 << "\nPhi end : wrapper = " < << 457 << " solid = " << GetPhiEnd( << 458 << "\nPhi # sides: wrapper = " < << 459 << " solid = " << GetSideCou << 460 << "\nPhi is open: wrapper = " < << 461 << " solid = " << 462 << ((Base_t::GetPhiDelta() < two << 463 G4Exception("G4UPolyhedra::BoundingLimit << 464 JustWarning, message); 385 JustWarning, message); 465 checkPhi = false; << 386 delete polyhedron; >> 387 return 0; 466 } 388 } 467 } << 389 else 468 } << 469 << 470 ////////////////////////////////////////////// << 471 // << 472 // Calculate extent under transform and specif << 473 << 474 G4bool << 475 G4UPolyhedra::CalculateExtent(const EAxis pAxi << 476 const G4VoxelLim << 477 const G4AffineTr << 478 G4double& << 479 { << 480 G4ThreeVector bmin, bmax; << 481 G4bool exist; << 482 << 483 // Check bounding box (bbox) << 484 // << 485 BoundingLimits(bmin,bmax); << 486 G4BoundingEnvelope bbox(bmin,bmax); << 487 #ifdef G4BBOX_EXTENT << 488 if (true) return bbox.CalculateExtent(pAxis, << 489 #endif << 490 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVox << 491 { << 492 return exist = pMin < pMax; << 493 } << 494 << 495 // To find the extent, RZ contour of the pol << 496 // in triangles. The extent is calculated as << 497 // all sub-polycones formed by rotation of t << 498 // << 499 G4TwoVectorList contourRZ; << 500 G4TwoVectorList triangles; << 501 std::vector<G4int> iout; << 502 G4double eminlim = pVoxelLimit.GetMinExtent( << 503 G4double emaxlim = pVoxelLimit.GetMaxExtent( << 504 << 505 // get RZ contour, ensure anticlockwise orde << 506 for (G4int i=0; i<GetNumRZCorner(); ++i) << 507 { << 508 G4PolyhedraSideRZ corner = GetCorner(i); << 509 contourRZ.emplace_back(corner.r,corner.z); << 510 } << 511 G4GeomTools::RemoveRedundantVertices(contour << 512 G4double area = G4GeomTools::PolygonArea(con << 513 if (area < 0.) std::reverse(contourRZ.begin( << 514 << 515 // triangulate RZ countour << 516 if (!G4GeomTools::TriangulatePolygon(contour << 517 { << 518 std::ostringstream message; << 519 message << "Triangulation of RZ contour ha << 520 << GetName() << " !" << 521 << "\nExtent has been calculated u << 522 G4Exception("G4UPolyhedra::CalculateExtent << 523 "GeomMgt1002",JustWarning,mess << 524 return bbox.CalculateExtent(pAxis,pVoxelLi << 525 } << 526 << 527 // set trigonometric values << 528 G4double sphi = GetStartPhi(); << 529 G4double ephi = GetEndPhi(); << 530 G4double dphi = IsOpen() ? ephi-sphi : t << 531 G4int ksteps = GetNumSide(); << 532 G4double astep = dphi/ksteps; << 533 G4double sinStep = std::sin(astep); << 534 G4double cosStep = std::cos(astep); << 535 G4double sinStart = GetSinStartPhi(); << 536 G4double cosStart = GetCosStartPhi(); << 537 << 538 // allocate vector lists << 539 std::vector<const G4ThreeVectorList *> polyg << 540 polygons.resize(ksteps+1); << 541 for (G4int k=0; k<ksteps+1; ++k) << 542 { << 543 polygons[k] = new G4ThreeVectorList(3); << 544 } << 545 << 546 // main loop along triangles << 547 pMin = kInfinity; << 548 pMax = -kInfinity; << 549 G4int ntria = triangles.size()/3; << 550 for (G4int i=0; i<ntria; ++i) << 551 { << 552 G4double sinCur = sinStart; << 553 G4double cosCur = cosStart; << 554 G4int i3 = i*3; << 555 for (G4int k=0; k<ksteps+1; ++k) // rotate << 556 { 390 { 557 auto ptr = const_cast<G4ThreeVectorList* << 391 return polyhedron; 558 auto iter = ptr->begin(); << 559 iter->set(triangles[i3+0].x()*cosCur, << 560 triangles[i3+0].x()*sinCur, << 561 triangles[i3+0].y()); << 562 iter++; << 563 iter->set(triangles[i3+1].x()*cosCur, << 564 triangles[i3+1].x()*sinCur, << 565 triangles[i3+1].y()); << 566 iter++; << 567 iter->set(triangles[i3+2].x()*cosCur, << 568 triangles[i3+2].x()*sinCur, << 569 triangles[i3+2].y()); << 570 << 571 G4double sinTmp = sinCur; << 572 sinCur = sinCur*cosStep + cosCur*sinStep << 573 cosCur = cosCur*cosStep - sinTmp*sinStep << 574 } 392 } 575 << 576 // set sub-envelope and adjust extent << 577 G4double emin,emax; << 578 G4BoundingEnvelope benv(polygons); << 579 if (!benv.CalculateExtent(pAxis,pVoxelLimi << 580 if (emin < pMin) pMin = emin; << 581 if (emax > pMax) pMax = emax; << 582 if (eminlim > pMin && emaxlim < pMax) brea << 583 } 393 } 584 // free memory << 585 for (G4int k=0; k<ksteps+1; ++k) { delete po << 586 return (pMin < pMax); << 587 } << 588 << 589 << 590 ////////////////////////////////////////////// << 591 // << 592 // CreatePolyhedron << 593 // << 594 G4Polyhedron* G4UPolyhedra::CreatePolyhedron() << 595 { << 596 return new G4PolyhedronPgon(wrStart, wrDelta << 597 } 394 } 598 << 599 #endif // G4GEOM_USE_USOLIDS << 600 395