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 // G4ExtrudedSolid implementation 27 // 28 // Author: Ivana Hrivnacova, IPN Orsay 29 // 30 // CHANGE HISTORY 31 // -------------- 32 // 33 // 31.10.2017 E.Tcherniaev: added implementation for a non-convex 34 // right prism 35 // 08.09.2017 E.Tcherniaev: added implementation for a convex 36 // right prism 37 // 21.10.2016 E.Tcherniaev: reimplemented CalculateExtent(), 38 // used G4GeomTools::PolygonArea() to calculate area, 39 // replaced IsConvex() with G4GeomTools::IsConvex() 40 // 02.03.2016 E.Tcherniaev: added CheckPolygon() to remove 41 // collinear and coincident points from polygon 42 // -------------------------------------------------------------------- 43 44 #include "G4ExtrudedSolid.hh" 45 46 #if !defined(G4GEOM_USE_UEXTRUDEDSOLID) 47 48 #include <set> 49 #include <algorithm> 50 #include <cmath> 51 #include <iomanip> 52 53 #include "G4GeomTools.hh" 54 #include "G4VoxelLimits.hh" 55 #include "G4AffineTransform.hh" 56 #include "G4BoundingEnvelope.hh" 57 58 #include "G4GeometryTolerance.hh" 59 #include "G4PhysicalConstants.hh" 60 #include "G4SystemOfUnits.hh" 61 #include "G4TriangularFacet.hh" 62 #include "G4QuadrangularFacet.hh" 63 64 //_____________________________________________________________________________ 65 66 G4ExtrudedSolid::G4ExtrudedSolid( const G4String& pName, 67 const std::vector<G4TwoVector>& polygon, 68 const std::vector<ZSection>& zsections) 69 : G4TessellatedSolid(pName), 70 fNv(polygon.size()), 71 fNz(zsections.size()), 72 fIsConvex(false), 73 fGeometryType("G4ExtrudedSolid"), 74 fSolidType(0) 75 { 76 // General constructor 77 78 // First check input parameters 79 80 if (fNv < 3) 81 { 82 std::ostringstream message; 83 message << "Number of vertices in polygon < 3 - " << pName; 84 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002", 85 FatalErrorInArgument, message); 86 } 87 88 if (fNz < 2) 89 { 90 std::ostringstream message; 91 message << "Number of z-sides < 2 - " << pName; 92 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002", 93 FatalErrorInArgument, message); 94 } 95 96 for ( std::size_t i=0; i<fNz-1; ++i ) 97 { 98 if ( zsections[i].fZ > zsections[i+1].fZ ) 99 { 100 std::ostringstream message; 101 message << "Z-sections have to be ordered by z value (z0 < z1 < z2...) - " 102 << pName; 103 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002", 104 FatalErrorInArgument, message); 105 } 106 if ( std::fabs( zsections[i+1].fZ - zsections[i].fZ ) < kCarToleranceHalf ) 107 { 108 std::ostringstream message; 109 message << "Z-sections with the same z position are not supported - " 110 << pName; 111 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0001", 112 FatalException, message); 113 } 114 } 115 116 // Copy polygon 117 // 118 fPolygon = polygon; 119 120 // Remove collinear and coincident vertices, if any 121 // 122 std::vector<G4int> removedVertices; 123 G4GeomTools::RemoveRedundantVertices(fPolygon,removedVertices, 124 2*kCarTolerance); 125 if (!removedVertices.empty()) 126 { 127 std::size_t nremoved = removedVertices.size(); 128 std::ostringstream message; 129 message << "The following "<< nremoved 130 << " vertices have been removed from polygon in " << pName 131 << "\nas collinear or coincident with other vertices: " 132 << removedVertices[0]; 133 for (std::size_t i=1; i<nremoved; ++i) message << ", " << removedVertices[i]; 134 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001", 135 JustWarning, message); 136 } 137 138 fNv = fPolygon.size(); 139 if (fNv < 3) 140 { 141 std::ostringstream message; 142 message << "Number of vertices in polygon after removal < 3 - " << pName; 143 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002", 144 FatalErrorInArgument, message); 145 } 146 147 // Check if polygon vertices are defined clockwise 148 // (the area is positive if polygon vertices are defined anti-clockwise) 149 // 150 if (G4GeomTools::PolygonArea(fPolygon) > 0.) 151 { 152 // Polygon vertices are defined anti-clockwise, we revert them 153 // G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001", 154 // JustWarning, 155 // "Polygon vertices defined anti-clockwise, reverting polygon"); 156 std::reverse(fPolygon.begin(),fPolygon.end()); 157 } 158 159 // Copy z-sections 160 // 161 fZSections = zsections; 162 163 G4bool result = MakeFacets(); 164 if (!result) 165 { 166 std::ostringstream message; 167 message << "Making facets failed - " << pName; 168 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0003", 169 FatalException, message); 170 } 171 fIsConvex = G4GeomTools::IsConvex(fPolygon); 172 173 ComputeProjectionParameters(); 174 175 // Check if the solid is a right prism, if so then set lateral planes 176 // 177 if ((fNz == 2) 178 && (fZSections[0].fScale == 1) && (fZSections[1].fScale == 1) 179 && (fZSections[0].fOffset == G4TwoVector(0,0)) 180 && (fZSections[1].fOffset == G4TwoVector(0,0))) 181 { 182 fSolidType = (fIsConvex) ? 1 : 2; // 1 - convex, 2 - non-convex right prism 183 ComputeLateralPlanes(); 184 } 185 } 186 187 //_____________________________________________________________________________ 188 189 G4ExtrudedSolid::G4ExtrudedSolid( const G4String& pName, 190 const std::vector<G4TwoVector>& polygon, 191 G4double dz, 192 const G4TwoVector& off1, G4double scale1, 193 const G4TwoVector& off2, G4double scale2 ) 194 : G4TessellatedSolid(pName), 195 fNv(polygon.size()), 196 fNz(2), 197 fGeometryType("G4ExtrudedSolid") 198 { 199 // Special constructor for solid with 2 z-sections 200 201 // First check input parameters 202 // 203 if (fNv < 3) 204 { 205 std::ostringstream message; 206 message << "Number of vertices in polygon < 3 - " << pName; 207 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002", 208 FatalErrorInArgument, message); 209 } 210 211 // Copy polygon 212 // 213 fPolygon = polygon; 214 215 // Remove collinear and coincident vertices, if any 216 // 217 std::vector<G4int> removedVertices; 218 G4GeomTools::RemoveRedundantVertices(fPolygon,removedVertices, 219 2*kCarTolerance); 220 if (!removedVertices.empty()) 221 { 222 std::size_t nremoved = removedVertices.size(); 223 std::ostringstream message; 224 message << "The following "<< nremoved 225 << " vertices have been removed from polygon in " << pName 226 << "\nas collinear or coincident with other vertices: " 227 << removedVertices[0]; 228 for (std::size_t i=1; i<nremoved; ++i) message << ", " << removedVertices[i]; 229 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001", 230 JustWarning, message); 231 } 232 233 fNv = fPolygon.size(); 234 if (fNv < 3) 235 { 236 std::ostringstream message; 237 message << "Number of vertices in polygon after removal < 3 - " << pName; 238 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0002", 239 FatalErrorInArgument, message); 240 } 241 242 // Check if polygon vertices are defined clockwise 243 // (the area is positive if polygon vertices are defined anti-clockwise) 244 // 245 if (G4GeomTools::PolygonArea(fPolygon) > 0.) 246 { 247 // Polygon vertices are defined anti-clockwise, we revert them 248 // G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids1001", 249 // JustWarning, 250 // "Polygon vertices defined anti-clockwise, reverting polygon"); 251 std::reverse(fPolygon.begin(),fPolygon.end()); 252 } 253 254 // Copy z-sections 255 // 256 fZSections.emplace_back(-dz, off1, scale1); 257 fZSections.emplace_back( dz, off2, scale2); 258 259 G4bool result = MakeFacets(); 260 if (!result) 261 { 262 std::ostringstream message; 263 message << "Making facets failed - " << pName; 264 G4Exception("G4ExtrudedSolid::G4ExtrudedSolid()", "GeomSolids0003", 265 FatalException, message); 266 } 267 fIsConvex = G4GeomTools::IsConvex(fPolygon); 268 269 ComputeProjectionParameters(); 270 271 // Check if the solid is a right prism, if so then set lateral planes 272 // 273 if ((scale1 == 1) && (scale2 == 1) 274 && (off1 == G4TwoVector(0,0)) && (off2 == G4TwoVector(0,0))) 275 { 276 fSolidType = (fIsConvex) ? 1 : 2; // 1 - convex, 2 - non-convex right prism 277 ComputeLateralPlanes(); 278 } 279 } 280 281 //_____________________________________________________________________________ 282 283 G4ExtrudedSolid::G4ExtrudedSolid( __void__& a ) 284 : G4TessellatedSolid(a), fGeometryType("G4ExtrudedSolid") 285 { 286 // Fake default constructor - sets only member data and allocates memory 287 // for usage restricted to object persistency. 288 } 289 290 //_____________________________________________________________________________ 291 292 G4ExtrudedSolid::G4ExtrudedSolid(const G4ExtrudedSolid&) = default; 293 294 //_____________________________________________________________________________ 295 296 G4ExtrudedSolid& G4ExtrudedSolid::operator = (const G4ExtrudedSolid& rhs) 297 { 298 // Check assignment to self 299 // 300 if (this == &rhs) { return *this; } 301 302 // Copy base class data 303 // 304 G4TessellatedSolid::operator=(rhs); 305 306 // Copy data 307 // 308 fNv = rhs.fNv; fNz = rhs.fNz; 309 fPolygon = rhs.fPolygon; fZSections = rhs.fZSections; 310 fTriangles = rhs.fTriangles; fIsConvex = rhs.fIsConvex; 311 fGeometryType = rhs.fGeometryType; 312 fSolidType = rhs.fSolidType; fPlanes = rhs.fPlanes; 313 fLines = rhs.fLines; fLengths = rhs.fLengths; 314 fKScales = rhs.fKScales; fScale0s = rhs.fScale0s; 315 fKOffsets = rhs.fKOffsets; fOffset0s = rhs.fOffset0s; 316 317 return *this; 318 } 319 320 //_____________________________________________________________________________ 321 322 G4ExtrudedSolid::~G4ExtrudedSolid() 323 { 324 // Destructor 325 } 326 327 //_____________________________________________________________________________ 328 329 void G4ExtrudedSolid::ComputeProjectionParameters() 330 { 331 // Compute parameters for point projections p(z) 332 // to the polygon scale & offset: 333 // scale(z) = k*z + scale0 334 // offset(z) = l*z + offset0 335 // p(z) = scale(z)*p0 + offset(z) 336 // p0 = (p(z) - offset(z))/scale(z); 337 // 338 339 for (std::size_t iz=0; iz<fNz-1; ++iz) 340 { 341 G4double z1 = fZSections[iz].fZ; 342 G4double z2 = fZSections[iz+1].fZ; 343 G4double scale1 = fZSections[iz].fScale; 344 G4double scale2 = fZSections[iz+1].fScale; 345 G4TwoVector off1 = fZSections[iz].fOffset; 346 G4TwoVector off2 = fZSections[iz+1].fOffset; 347 348 G4double kscale = (scale2 - scale1)/(z2 - z1); 349 G4double scale0 = scale2 - kscale*(z2 - z1)/2.0; 350 G4TwoVector koff = (off2 - off1)/(z2 - z1); 351 G4TwoVector off0 = off2 - koff*(z2 - z1)/2.0; 352 353 fKScales.push_back(kscale); 354 fScale0s.push_back(scale0); 355 fKOffsets.push_back(koff); 356 fOffset0s.push_back(off0); 357 } 358 } 359 360 //_____________________________________________________________________________ 361 362 void G4ExtrudedSolid::ComputeLateralPlanes() 363 { 364 // Compute lateral planes: a*x + b*y + c*z + d = 0 365 // 366 std::size_t Nv = fPolygon.size(); 367 fPlanes.resize(Nv); 368 for (std::size_t i=0, k=Nv-1; i<Nv; k=i++) 369 { 370 G4TwoVector norm = (fPolygon[i] - fPolygon[k]).unit(); 371 fPlanes[i].a = -norm.y(); 372 fPlanes[i].b = norm.x(); 373 fPlanes[i].c = 0; 374 fPlanes[i].d = norm.y()*fPolygon[i].x() - norm.x()*fPolygon[i].y(); 375 } 376 377 // Compute edge equations: x = k*y + m 378 // and edge lengths 379 // 380 fLines.resize(Nv); 381 fLengths.resize(Nv); 382 for (std::size_t i=0, k=Nv-1; i<Nv; k=i++) 383 { 384 if (fPolygon[k].y() == fPolygon[i].y()) 385 { 386 fLines[i].k = 0; 387 fLines[i].m = fPolygon[i].x(); 388 } 389 else 390 { 391 G4double ctg = (fPolygon[k].x()-fPolygon[i].x())/(fPolygon[k].y()-fPolygon[i].y()); 392 fLines[i].k = ctg; 393 fLines[i].m = fPolygon[i].x() - ctg*fPolygon[i].y(); 394 } 395 fLengths[i] = (fPolygon[i] - fPolygon[k]).mag(); 396 } 397 } 398 399 //_____________________________________________________________________________ 400 401 G4ThreeVector G4ExtrudedSolid::GetVertex(G4int iz, G4int ind) const 402 { 403 // Shift and scale vertices 404 405 return { fPolygon[ind].x() * fZSections[iz].fScale 406 + fZSections[iz].fOffset.x(), 407 fPolygon[ind].y() * fZSections[iz].fScale 408 + fZSections[iz].fOffset.y(), 409 fZSections[iz].fZ }; 410 } 411 412 //_____________________________________________________________________________ 413 414 G4TwoVector G4ExtrudedSolid::ProjectPoint(const G4ThreeVector& point) const 415 { 416 // Project point in the polygon scale 417 // scale(z) = k*z + scale0 418 // offset(z) = l*z + offset0 419 // p(z) = scale(z)*p0 + offset(z) 420 // p0 = (p(z) - offset(z))/scale(z); 421 422 // Select projection (z-segment of the solid) according to p.z() 423 // 424 std::size_t iz = 0; 425 while ( point.z() > fZSections[iz+1].fZ && iz < fNz-2 ) { ++iz; } 426 // Loop checking, 13.08.2015, G.Cosmo 427 428 G4double z0 = ( fZSections[iz+1].fZ + fZSections[iz].fZ )/2.0; 429 G4TwoVector p2(point.x(), point.y()); 430 G4double pscale = fKScales[iz]*(point.z()-z0) + fScale0s[iz]; 431 G4TwoVector poffset = fKOffsets[iz]*(point.z()-z0) + fOffset0s[iz]; 432 433 // G4cout << point << " projected to " 434 // << iz << "-th z-segment polygon as " 435 // << (p2 - poffset)/pscale << G4endl; 436 437 // pscale is always >0 as it is an interpolation between two 438 // positive scale values 439 // 440 return (p2 - poffset)/pscale; 441 } 442 443 //_____________________________________________________________________________ 444 445 G4bool G4ExtrudedSolid::IsSameLine(const G4TwoVector& p, 446 const G4TwoVector& l1, 447 const G4TwoVector& l2) const 448 { 449 // Return true if p is on the line through l1, l2 450 451 if ( l1.x() == l2.x() ) 452 { 453 return std::fabs(p.x() - l1.x()) < kCarToleranceHalf; 454 } 455 G4double slope= ((l2.y() - l1.y())/(l2.x() - l1.x())); 456 G4double predy= l1.y() + slope *(p.x() - l1.x()); 457 G4double dy= p.y() - predy; 458 459 // Calculate perpendicular distance 460 // 461 // G4double perpD= std::fabs(dy) / std::sqrt( 1 + slope * slope ); 462 // G4bool simpleComp= (perpD<kCarToleranceHalf); 463 464 // Check perpendicular distance vs tolerance 'directly' 465 // 466 G4bool squareComp = (dy*dy < (1+slope*slope) 467 * kCarToleranceHalf * kCarToleranceHalf); 468 469 // return simpleComp; 470 return squareComp; 471 } 472 473 //_____________________________________________________________________________ 474 475 G4bool G4ExtrudedSolid::IsSameLineSegment(const G4TwoVector& p, 476 const G4TwoVector& l1, 477 const G4TwoVector& l2) const 478 { 479 // Return true if p is on the line through l1, l2 and lies between 480 // l1 and l2 481 482 if ( p.x() < std::min(l1.x(), l2.x()) - kCarToleranceHalf || 483 p.x() > std::max(l1.x(), l2.x()) + kCarToleranceHalf || 484 p.y() < std::min(l1.y(), l2.y()) - kCarToleranceHalf || 485 p.y() > std::max(l1.y(), l2.y()) + kCarToleranceHalf ) 486 { 487 return false; 488 } 489 490 return IsSameLine(p, l1, l2); 491 } 492 493 //_____________________________________________________________________________ 494 495 G4bool G4ExtrudedSolid::IsSameSide(const G4TwoVector& p1, 496 const G4TwoVector& p2, 497 const G4TwoVector& l1, 498 const G4TwoVector& l2) const 499 { 500 // Return true if p1 and p2 are on the same side of the line through l1, l2 501 502 return ( (p1.x() - l1.x()) * (l2.y() - l1.y()) 503 - (l2.x() - l1.x()) * (p1.y() - l1.y()) ) 504 * ( (p2.x() - l1.x()) * (l2.y() - l1.y()) 505 - (l2.x() - l1.x()) * (p2.y() - l1.y()) ) > 0; 506 } 507 508 //_____________________________________________________________________________ 509 510 G4bool G4ExtrudedSolid::IsPointInside(const G4TwoVector& a, 511 const G4TwoVector& b, 512 const G4TwoVector& c, 513 const G4TwoVector& p) const 514 { 515 // Return true if p is inside of triangle abc or on its edges, 516 // else returns false 517 518 // Check extent first 519 // 520 if ( ( p.x() < a.x() && p.x() < b.x() && p.x() < c.x() ) || 521 ( p.x() > a.x() && p.x() > b.x() && p.x() > c.x() ) || 522 ( p.y() < a.y() && p.y() < b.y() && p.y() < c.y() ) || 523 ( p.y() > a.y() && p.y() > b.y() && p.y() > c.y() ) ) return false; 524 525 G4bool inside 526 = IsSameSide(p, a, b, c) 527 && IsSameSide(p, b, a, c) 528 && IsSameSide(p, c, a, b); 529 530 G4bool onEdge 531 = IsSameLineSegment(p, a, b) 532 || IsSameLineSegment(p, b, c) 533 || IsSameLineSegment(p, c, a); 534 535 return inside || onEdge; 536 } 537 538 //_____________________________________________________________________________ 539 540 G4double 541 G4ExtrudedSolid::GetAngle(const G4TwoVector& po, 542 const G4TwoVector& pa, 543 const G4TwoVector& pb) const 544 { 545 // Return the angle of the vertex in po 546 547 G4TwoVector t1 = pa - po; 548 G4TwoVector t2 = pb - po; 549 550 G4double result = (std::atan2(t1.y(), t1.x()) - std::atan2(t2.y(), t2.x())); 551 552 if ( result < 0 ) result += 2*pi; 553 554 return result; 555 } 556 557 //_____________________________________________________________________________ 558 559 G4VFacet* 560 G4ExtrudedSolid::MakeDownFacet(G4int ind1, G4int ind2, G4int ind3) const 561 { 562 // Create a triangular facet from the polygon points given by indices 563 // forming the down side ( the normal goes in -z) 564 565 std::vector<G4ThreeVector> vertices; 566 vertices.push_back(GetVertex(0, ind1)); 567 vertices.push_back(GetVertex(0, ind2)); 568 vertices.push_back(GetVertex(0, ind3)); 569 570 // first vertex most left 571 // 572 G4ThreeVector cross 573 = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]); 574 575 if ( cross.z() > 0.0 ) 576 { 577 // vertices ordered clock wise has to be reordered 578 579 // G4cout << "G4ExtrudedSolid::MakeDownFacet: reordering vertices " 580 // << ind1 << ", " << ind2 << ", " << ind3 << G4endl; 581 582 G4ThreeVector tmp = vertices[1]; 583 vertices[1] = vertices[2]; 584 vertices[2] = tmp; 585 } 586 587 return new G4TriangularFacet(vertices[0], vertices[1], 588 vertices[2], ABSOLUTE); 589 } 590 591 //_____________________________________________________________________________ 592 593 G4VFacet* 594 G4ExtrudedSolid::MakeUpFacet(G4int ind1, G4int ind2, G4int ind3) const 595 { 596 // Creates a triangular facet from the polygon points given by indices 597 // forming the upper side ( z>0 ) 598 599 std::vector<G4ThreeVector> vertices; 600 vertices.push_back(GetVertex((G4int)fNz-1, ind1)); 601 vertices.push_back(GetVertex((G4int)fNz-1, ind2)); 602 vertices.push_back(GetVertex((G4int)fNz-1, ind3)); 603 604 // first vertex most left 605 // 606 G4ThreeVector cross 607 = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]); 608 609 if ( cross.z() < 0.0 ) 610 { 611 // vertices ordered clock wise has to be reordered 612 613 // G4cout << "G4ExtrudedSolid::MakeUpFacet: reordering vertices " 614 // << ind1 << ", " << ind2 << ", " << ind3 << G4endl; 615 616 G4ThreeVector tmp = vertices[1]; 617 vertices[1] = vertices[2]; 618 vertices[2] = tmp; 619 } 620 621 return new G4TriangularFacet(vertices[0], vertices[1], 622 vertices[2], ABSOLUTE); 623 } 624 625 //_____________________________________________________________________________ 626 627 G4bool G4ExtrudedSolid::AddGeneralPolygonFacets() 628 { 629 // Decompose polygonal sides in triangular facets 630 631 typedef std::pair < G4TwoVector, G4int > Vertex; 632 633 static const G4double kAngTolerance = 634 G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 635 636 // Fill one more vector 637 // 638 std::vector< Vertex > verticesToBeDone; 639 for ( G4int i=0; i<(G4int)fNv; ++i ) 640 { 641 verticesToBeDone.emplace_back(fPolygon[i], i); 642 } 643 std::vector< Vertex > ears; 644 645 auto c1 = verticesToBeDone.begin(); 646 auto c2 = c1+1; 647 auto c3 = c1+2; 648 while ( verticesToBeDone.size()>2 ) // Loop checking, 13.08.2015, G.Cosmo 649 { 650 651 // G4cout << "Looking at triangle : " 652 // << c1->second << " " << c2->second 653 // << " " << c3->second << G4endl; 654 //G4cout << "Looking at triangle : " 655 // << c1->first << " " << c2->first 656 // << " " << c3->first << G4endl; 657 658 // skip concave vertices 659 // 660 G4double angle = GetAngle(c2->first, c3->first, c1->first); 661 662 //G4cout << "angle " << angle << G4endl; 663 664 std::size_t counter = 0; 665 while ( angle >= (pi-kAngTolerance) ) // Loop checking, 13.08.2015, G.Cosmo 666 { 667 // G4cout << "Skipping concave vertex " << c2->second << G4endl; 668 669 // try next three consecutive vertices 670 // 671 c1 = c2; 672 c2 = c3; 673 ++c3; 674 if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); } 675 676 //G4cout << "Looking at triangle : " 677 // << c1->first << " " << c2->first 678 // << " " << c3->first << G4endl; 679 680 angle = GetAngle(c2->first, c3->first, c1->first); 681 //G4cout << "angle " << angle << G4endl; 682 683 ++counter; 684 685 if ( counter > fNv ) 686 { 687 G4Exception("G4ExtrudedSolid::AddGeneralPolygonFacets", 688 "GeomSolids0003", FatalException, 689 "Triangularisation has failed."); 690 break; 691 } 692 } 693 694 G4bool good = true; 695 for ( auto it=verticesToBeDone.cbegin(); it!=verticesToBeDone.cend(); ++it ) 696 { 697 // skip vertices of tested triangle 698 // 699 if ( it == c1 || it == c2 || it == c3 ) { continue; } 700 701 if ( IsPointInside(c1->first, c2->first, c3->first, it->first) ) 702 { 703 // G4cout << "Point " << it->second << " is inside" << G4endl; 704 good = false; 705 706 // try next three consecutive vertices 707 // 708 c1 = c2; 709 c2 = c3; 710 ++c3; 711 if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); } 712 break; 713 } 714 // else 715 // { G4cout << "Point " << it->second << " is outside" << G4endl; } 716 } 717 if ( good ) 718 { 719 // all points are outside triangle, we can make a facet 720 721 // G4cout << "Found triangle : " 722 // << c1->second << " " << c2->second 723 // << " " << c3->second << G4endl; 724 725 G4bool result; 726 result = AddFacet( MakeDownFacet(c1->second, c2->second, c3->second) ); 727 if ( ! result ) { return false; } 728 729 result = AddFacet( MakeUpFacet(c1->second, c2->second, c3->second) ); 730 if ( ! result ) { return false; } 731 732 std::vector<G4int> triangle(3); 733 triangle[0] = c1->second; 734 triangle[1] = c2->second; 735 triangle[2] = c3->second; 736 fTriangles.push_back(std::move(triangle)); 737 738 // remove the ear point from verticesToBeDone 739 // 740 verticesToBeDone.erase(c2); 741 c1 = verticesToBeDone.begin(); 742 c2 = c1+1; 743 c3 = c1+2; 744 } 745 } 746 return true; 747 } 748 749 //_____________________________________________________________________________ 750 751 G4bool G4ExtrudedSolid::MakeFacets() 752 { 753 // Define facets 754 755 G4bool good; 756 757 // Decomposition of polygonal sides in the facets 758 // 759 if ( fNv == 3 ) 760 { 761 good = AddFacet( new G4TriangularFacet( GetVertex(0, 0), GetVertex(0, 1), 762 GetVertex(0, 2), ABSOLUTE) ); 763 if ( ! good ) { return false; } 764 765 good = AddFacet( new G4TriangularFacet( GetVertex((G4int)fNz-1, 2), 766 GetVertex((G4int)fNz-1, 1), 767 GetVertex((G4int)fNz-1, 0), 768 ABSOLUTE) ); 769 if ( ! good ) { return false; } 770 771 std::vector<G4int> triangle(3); 772 triangle[0] = 0; 773 triangle[1] = 1; 774 triangle[2] = 2; 775 fTriangles.push_back(std::move(triangle)); 776 } 777 778 else if ( fNv == 4 ) 779 { 780 good = AddFacet( new G4QuadrangularFacet( GetVertex(0, 0),GetVertex(0, 1), 781 GetVertex(0, 2),GetVertex(0, 3), 782 ABSOLUTE) ); 783 if ( ! good ) { return false; } 784 785 good = AddFacet( new G4QuadrangularFacet( GetVertex((G4int)fNz-1, 3), 786 GetVertex((G4int)fNz-1, 2), 787 GetVertex((G4int)fNz-1, 1), 788 GetVertex((G4int)fNz-1, 0), 789 ABSOLUTE) ); 790 if ( ! good ) { return false; } 791 792 std::vector<G4int> triangle1(3); 793 triangle1[0] = 0; 794 triangle1[1] = 1; 795 triangle1[2] = 2; 796 fTriangles.push_back(std::move(triangle1)); 797 798 std::vector<G4int> triangle2(3); 799 triangle2[0] = 0; 800 triangle2[1] = 2; 801 triangle2[2] = 3; 802 fTriangles.push_back(std::move(triangle2)); 803 } 804 else 805 { 806 good = AddGeneralPolygonFacets(); 807 if ( ! good ) { return false; } 808 } 809 810 // The quadrangular sides 811 // 812 for ( G4int iz = 0; iz < (G4int)fNz-1; ++iz ) 813 { 814 for ( G4int i = 0; i < (G4int)fNv; ++i ) 815 { 816 G4int j = (i+1) % fNv; 817 good = AddFacet( new G4QuadrangularFacet 818 ( GetVertex(iz, j), GetVertex(iz, i), 819 GetVertex(iz+1, i), GetVertex(iz+1, j), ABSOLUTE) ); 820 if ( ! good ) { return false; } 821 } 822 } 823 824 SetSolidClosed(true); 825 826 return good; 827 } 828 829 //_____________________________________________________________________________ 830 831 G4GeometryType G4ExtrudedSolid::GetEntityType () const 832 { 833 // Return entity type 834 835 return fGeometryType; 836 } 837 838 //_____________________________________________________________________________ 839 840 G4bool G4ExtrudedSolid::IsFaceted () const 841 { 842 return true; 843 } 844 845 //_____________________________________________________________________________ 846 847 G4VSolid* G4ExtrudedSolid::Clone() const 848 { 849 return new G4ExtrudedSolid(*this); 850 } 851 852 //_____________________________________________________________________________ 853 854 EInside G4ExtrudedSolid::Inside(const G4ThreeVector &p) const 855 { 856 switch (fSolidType) 857 { 858 case 1: // convex right prism 859 { 860 G4double dist = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ); 861 if (dist > kCarToleranceHalf) { return kOutside; } 862 863 std::size_t np = fPlanes.size(); 864 for (std::size_t i=0; i<np; ++i) 865 { 866 G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d; 867 if (dd > dist) { dist = dd; } 868 } 869 if (dist > kCarToleranceHalf) { return kOutside; } 870 return (dist > -kCarToleranceHalf) ? kSurface : kInside; 871 } 872 case 2: // non-convex right prism 873 { 874 G4double distz = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ); 875 if (distz > kCarToleranceHalf) { return kOutside; } 876 877 G4bool in = PointInPolygon(p); 878 if (distz > -kCarToleranceHalf && in) { return kSurface; } 879 880 G4double dd = DistanceToPolygonSqr(p) - kCarToleranceHalf*kCarToleranceHalf; 881 if (in) 882 { 883 return (dd >= 0) ? kInside : kSurface; 884 } 885 else 886 { 887 return (dd > 0) ? kOutside : kSurface; 888 } 889 } 890 } 891 892 // Override the base class function as it fails in case of concave polygon. 893 // Project the point in the original polygon scale and check if it is inside 894 // for each triangle. 895 896 // Check first if outside extent 897 // 898 if ( p.x() < GetMinXExtent() - kCarToleranceHalf || 899 p.x() > GetMaxXExtent() + kCarToleranceHalf || 900 p.y() < GetMinYExtent() - kCarToleranceHalf || 901 p.y() > GetMaxYExtent() + kCarToleranceHalf || 902 p.z() < GetMinZExtent() - kCarToleranceHalf || 903 p.z() > GetMaxZExtent() + kCarToleranceHalf ) 904 { 905 // G4cout << "G4ExtrudedSolid::Outside extent: " << p << G4endl; 906 return kOutside; 907 } 908 909 // Project point p(z) to the polygon scale p0 910 // 911 G4TwoVector pscaled = ProjectPoint(p); 912 913 // Check if on surface of polygon 914 // 915 for ( G4int i=0; i<(G4int)fNv; ++i ) 916 { 917 G4int j = (i+1) % fNv; 918 if ( IsSameLineSegment(pscaled, fPolygon[i], fPolygon[j]) ) 919 { 920 // G4cout << "G4ExtrudedSolid::Inside return Surface (on polygon) " 921 // << G4endl; 922 923 return kSurface; 924 } 925 } 926 927 // Now check if inside triangles 928 // 929 auto it = fTriangles.cbegin(); 930 G4bool inside = false; 931 do // Loop checking, 13.08.2015, G.Cosmo 932 { 933 if ( IsPointInside(fPolygon[(*it)[0]], fPolygon[(*it)[1]], 934 fPolygon[(*it)[2]], pscaled) ) { inside = true; } 935 ++it; 936 } while ( (!inside) && (it != fTriangles.cend()) ); 937 938 if ( inside ) 939 { 940 // Check if on surface of z sides 941 // 942 if ( std::fabs( p.z() - fZSections[0].fZ ) < kCarToleranceHalf || 943 std::fabs( p.z() - fZSections[fNz-1].fZ ) < kCarToleranceHalf ) 944 { 945 // G4cout << "G4ExtrudedSolid::Inside return Surface (on z side)" 946 // << G4endl; 947 948 return kSurface; 949 } 950 951 // G4cout << "G4ExtrudedSolid::Inside return Inside" << G4endl; 952 953 return kInside; 954 } 955 956 // G4cout << "G4ExtrudedSolid::Inside return Outside " << G4endl; 957 958 return kOutside; 959 } 960 961 //_____________________________________________________________________________ 962 963 G4ThreeVector G4ExtrudedSolid::SurfaceNormal(const G4ThreeVector& p) const 964 { 965 G4int nsurf = 0; 966 G4double nx = 0., ny = 0., nz = 0.; 967 switch (fSolidType) 968 { 969 case 1: // convex right prism 970 { 971 if (std::abs(p.z() - fZSections[0].fZ) <= kCarToleranceHalf) 972 { 973 nz = -1; ++nsurf; 974 } 975 if (std::abs(p.z() - fZSections[1].fZ) <= kCarToleranceHalf) 976 { 977 nz = 1; ++nsurf; 978 } 979 for (std::size_t i=0; i<fNv; ++i) 980 { 981 G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d; 982 if (std::abs(dd) > kCarToleranceHalf) continue; 983 nx += fPlanes[i].a; 984 ny += fPlanes[i].b; 985 ++nsurf; 986 } 987 break; 988 } 989 case 2: // non-convex right prism 990 { 991 if (std::abs(p.z() - fZSections[0].fZ) <= kCarToleranceHalf) 992 { 993 nz = -1; ++nsurf; 994 } 995 if (std::abs(p.z() - fZSections[1].fZ) <= kCarToleranceHalf) 996 { 997 nz = 1; ++nsurf; 998 } 999 1000 G4double sqrCarToleranceHalf = kCarToleranceHalf*kCarToleranceHalf; 1001 for (std::size_t i=0, k=fNv-1; i<fNv; k=i++) 1002 { 1003 G4double ix = p.x() - fPolygon[i].x(); 1004 G4double iy = p.y() - fPolygon[i].y(); 1005 G4double u = fPlanes[i].a*iy - fPlanes[i].b*ix; 1006 if (u < 0) 1007 { 1008 if (ix*ix + iy*iy > sqrCarToleranceHalf) continue; 1009 } 1010 else if (u > fLengths[i]) 1011 { 1012 G4double kx = p.x() - fPolygon[k].x(); 1013 G4double ky = p.y() - fPolygon[k].y(); 1014 if (kx*kx + ky*ky > sqrCarToleranceHalf) continue; 1015 } 1016 else 1017 { 1018 G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d; 1019 if (dd*dd > sqrCarToleranceHalf) continue; 1020 } 1021 nx += fPlanes[i].a; 1022 ny += fPlanes[i].b; 1023 ++nsurf; 1024 } 1025 break; 1026 } 1027 default: 1028 { 1029 return G4TessellatedSolid::SurfaceNormal(p); 1030 } 1031 } 1032 1033 // Return normal (right prism) 1034 // 1035 if (nsurf == 1) 1036 { 1037 return { nx,ny,nz }; 1038 } 1039 else if (nsurf != 0) // edge or corner 1040 { 1041 return G4ThreeVector(nx,ny,nz).unit(); 1042 } 1043 else 1044 { 1045 // Point is not on the surface, compute approximate normal 1046 // 1047 #ifdef G4CSGDEBUG 1048 std::ostringstream message; 1049 G4long oldprc = message.precision(16); 1050 message << "Point p is not on surface (!?) of solid: " 1051 << GetName() << G4endl; 1052 message << "Position:\n"; 1053 message << " p.x() = " << p.x()/mm << " mm\n"; 1054 message << " p.y() = " << p.y()/mm << " mm\n"; 1055 message << " p.z() = " << p.z()/mm << " mm"; 1056 G4cout.precision(oldprc) ; 1057 G4Exception("G4TesselatedSolid::SurfaceNormal(p)", "GeomSolids1002", 1058 JustWarning, message ); 1059 DumpInfo(); 1060 #endif 1061 return ApproxSurfaceNormal(p); 1062 } 1063 } 1064 1065 //_____________________________________________________________________________ 1066 1067 G4ThreeVector G4ExtrudedSolid::ApproxSurfaceNormal(const G4ThreeVector& p) const 1068 { 1069 // This method is valid only for right prisms and 1070 // normally should not be called 1071 1072 if (fSolidType == 1 || fSolidType == 2) 1073 { 1074 // Find distances to z-planes 1075 // 1076 G4double dz0 = fZSections[0].fZ - p.z(); 1077 G4double dz1 = p.z() - fZSections[1].fZ; 1078 G4double ddz0 = dz0*dz0; 1079 G4double ddz1 = dz1*dz1; 1080 1081 // Find nearest lateral side and distance to it 1082 // 1083 std::size_t iside = 0; 1084 G4double dd = DBL_MAX; 1085 for (std::size_t i=0, k=fNv-1; i<fNv; k=i++) 1086 { 1087 G4double ix = p.x() - fPolygon[i].x(); 1088 G4double iy = p.y() - fPolygon[i].y(); 1089 G4double u = fPlanes[i].a*iy - fPlanes[i].b*ix; 1090 if (u < 0) 1091 { 1092 G4double tmp = ix*ix + iy*iy; 1093 if (tmp < dd) { dd = tmp; iside = i; } 1094 } 1095 else if (u > fLengths[i]) 1096 { 1097 G4double kx = p.x() - fPolygon[k].x(); 1098 G4double ky = p.y() - fPolygon[k].y(); 1099 G4double tmp = kx*kx + ky*ky; 1100 if (tmp < dd) { dd = tmp; iside = i; } 1101 } 1102 else 1103 { 1104 G4double tmp = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d; 1105 tmp *= tmp; 1106 if (tmp < dd) { dd = tmp; iside = i; } 1107 } 1108 } 1109 1110 // Find region 1111 // 1112 // 3 | 1 | 3 1113 // ----+-------+---- 1114 // 2 | 0 | 2 1115 // ----+-------+---- 1116 // 3 | 1 | 3 1117 // 1118 G4int iregion = 0; 1119 if (std::max(dz0,dz1) > 0) iregion = 1; 1120 1121 G4bool in = PointInPolygon(p); 1122 if (!in) iregion += 2; 1123 1124 // Return normal 1125 // 1126 switch (iregion) 1127 { 1128 case 0: 1129 { 1130 if (ddz0 <= ddz1 && ddz0 <= dd) return {0, 0,-1}; 1131 if (ddz1 <= ddz0 && ddz1 <= dd) return {0, 0, 1}; 1132 return { fPlanes[iside].a, fPlanes[iside].b, 0 }; 1133 } 1134 case 1: 1135 { 1136 return { 0, 0, (G4double)((dz0 > dz1) ? -1 : 1) }; 1137 } 1138 case 2: 1139 { 1140 return { fPlanes[iside].a, fPlanes[iside].b, 0 }; 1141 } 1142 case 3: 1143 { 1144 G4double dzmax = std::max(dz0,dz1); 1145 if (dzmax*dzmax > dd) return { 0, 0, (G4double)((dz0 > dz1) ? -1 : 1) }; 1146 return { fPlanes[iside].a, fPlanes[iside].b, 0 }; 1147 } 1148 } 1149 } 1150 return {0,0,0}; 1151 } 1152 1153 //_____________________________________________________________________________ 1154 1155 G4double G4ExtrudedSolid::DistanceToIn(const G4ThreeVector& p, 1156 const G4ThreeVector& v) const 1157 { 1158 G4double z0 = fZSections[0].fZ; 1159 G4double z1 = fZSections[fNz-1].fZ; 1160 if ((p.z() <= z0 + kCarToleranceHalf) && v.z() <= 0) return kInfinity; 1161 if ((p.z() >= z1 - kCarToleranceHalf) && v.z() >= 0) return kInfinity; 1162 1163 switch (fSolidType) 1164 { 1165 case 1: // convex right prism 1166 { 1167 // Intersection with Z planes 1168 // 1169 G4double dz = (z1 - z0)*0.5; 1170 G4double pz = p.z() - dz - z0; 1171 1172 G4double invz = (v.z() == 0) ? DBL_MAX : -1./v.z(); 1173 G4double ddz = (invz < 0) ? dz : -dz; 1174 G4double tzmin = (pz + ddz)*invz; 1175 G4double tzmax = (pz - ddz)*invz; 1176 1177 // Intersection with lateral planes 1178 // 1179 std::size_t np = fPlanes.size(); 1180 G4double txmin = tzmin, txmax = tzmax; 1181 for (std::size_t i=0; i<np; ++i) 1182 { 1183 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y(); 1184 G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].d; 1185 if (dist >= -kCarToleranceHalf) 1186 { 1187 if (cosa >= 0) { return kInfinity; } 1188 G4double tmp = -dist/cosa; 1189 if (txmin < tmp) { txmin = tmp; } 1190 } 1191 else if (cosa > 0) 1192 { 1193 G4double tmp = -dist/cosa; 1194 if (txmax > tmp) { txmax = tmp; } 1195 } 1196 } 1197 1198 // Find distance 1199 // 1200 G4double tmin = txmin, tmax = txmax; 1201 if (tmax <= tmin + kCarToleranceHalf) // touch or no hit 1202 { 1203 return kInfinity; 1204 } 1205 return (tmin < kCarToleranceHalf) ? 0. : tmin; 1206 } 1207 case 2: // non-convex right prism 1208 { 1209 } 1210 } 1211 return G4TessellatedSolid::DistanceToIn(p,v); 1212 } 1213 1214 //_____________________________________________________________________________ 1215 1216 G4double G4ExtrudedSolid::DistanceToIn (const G4ThreeVector& p) const 1217 { 1218 switch (fSolidType) 1219 { 1220 case 1: // convex right prism 1221 { 1222 G4double dist = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ); 1223 std::size_t np = fPlanes.size(); 1224 for (std::size_t i=0; i<np; ++i) 1225 { 1226 G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d; 1227 if (dd > dist) dist = dd; 1228 } 1229 return (dist > 0) ? dist : 0.; 1230 } 1231 case 2: // non-convex right prism 1232 { 1233 G4bool in = PointInPolygon(p); 1234 if (in) 1235 { 1236 G4double distz= std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ); 1237 return (distz > 0) ? distz : 0; 1238 } 1239 else 1240 { 1241 G4double distz= std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ); 1242 G4double dd = DistanceToPolygonSqr(p); 1243 if (distz > 0) dd += distz*distz; 1244 return std::sqrt(dd); 1245 } 1246 } 1247 } 1248 1249 // General case: use tessellated solid 1250 return G4TessellatedSolid::DistanceToIn(p); 1251 } 1252 1253 //_____________________________________________________________________________ 1254 1255 G4double G4ExtrudedSolid::DistanceToOut (const G4ThreeVector &p, 1256 const G4ThreeVector &v, 1257 const G4bool calcNorm, 1258 G4bool* validNorm, 1259 G4ThreeVector* n) const 1260 { 1261 G4bool getnorm = calcNorm; 1262 if (getnorm) *validNorm = true; 1263 1264 G4double z0 = fZSections[0].fZ; 1265 G4double z1 = fZSections[fNz-1].fZ; 1266 if ((p.z() <= z0 + kCarToleranceHalf) && v.z() < 0) 1267 { 1268 if (getnorm) n->set(0,0,-1); 1269 return 0; 1270 } 1271 if ((p.z() >= z1 - kCarToleranceHalf) && v.z() > 0) 1272 { 1273 if (getnorm) n->set(0,0,1); 1274 return 0; 1275 } 1276 1277 switch (fSolidType) 1278 { 1279 case 1: // convex right prism 1280 { 1281 // Intersection with Z planes 1282 // 1283 G4double dz = (z1 - z0)*0.5; 1284 G4double pz = p.z() - 0.5 * (z0 + z1); 1285 1286 G4double vz = v.z(); 1287 G4double tmax = (vz == 0) ? DBL_MAX : (std::copysign(dz,vz) - pz)/vz; 1288 G4int iside = (vz < 0) ? -4 : -2; // little trick: (-4+3)=-1, (-2+3)=+1 1289 1290 // Intersection with lateral planes 1291 // 1292 std::size_t np = fPlanes.size(); 1293 for (std::size_t i=0; i<np; ++i) 1294 { 1295 G4double cosa = fPlanes[i].a*v.x()+fPlanes[i].b*v.y(); 1296 if (cosa > 0) 1297 { 1298 G4double dist = fPlanes[i].a*p.x()+fPlanes[i].b*p.y()+fPlanes[i].d; 1299 if (dist >= -kCarToleranceHalf) 1300 { 1301 if (getnorm) n->set(fPlanes[i].a, fPlanes[i].b, fPlanes[i].c); 1302 return 0; 1303 } 1304 G4double tmp = -dist/cosa; 1305 if (tmax > tmp) { tmax = tmp; iside = (G4int)i; } 1306 } 1307 } 1308 1309 // Set normal, if required, and return distance 1310 // 1311 if (getnorm) 1312 { 1313 if (iside < 0) 1314 { n->set(0, 0, iside + 3); } // (-4+3)=-1, (-2+3)=+1 1315 else 1316 { n->set(fPlanes[iside].a, fPlanes[iside].b, fPlanes[iside].c); } 1317 } 1318 return tmax; 1319 } 1320 case 2: // non-convex right prism 1321 { 1322 } 1323 } 1324 1325 // Override the base class function to redefine validNorm 1326 // (the solid can be concave) 1327 1328 G4double distOut = 1329 G4TessellatedSolid::DistanceToOut(p, v, calcNorm, validNorm, n); 1330 if (validNorm != nullptr) { *validNorm = fIsConvex; } 1331 1332 return distOut; 1333 } 1334 1335 //_____________________________________________________________________________ 1336 1337 G4double G4ExtrudedSolid::DistanceToOut(const G4ThreeVector& p) const 1338 { 1339 switch (fSolidType) 1340 { 1341 case 1: // convex right prism 1342 { 1343 G4double dist = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ); 1344 std::size_t np = fPlanes.size(); 1345 for (std::size_t i=0; i<np; ++i) 1346 { 1347 G4double dd = fPlanes[i].a*p.x() + fPlanes[i].b*p.y() + fPlanes[i].d; 1348 if (dd > dist) dist = dd; 1349 } 1350 return (dist < 0) ? -dist : 0.; 1351 } 1352 case 2: // non-convex right prism 1353 { 1354 G4double distz = std::max(fZSections[0].fZ-p.z(),p.z()-fZSections[1].fZ); 1355 G4bool in = PointInPolygon(p); 1356 if (distz >= 0 || (!in)) return 0; // point is outside 1357 return std::min(-distz,std::sqrt(DistanceToPolygonSqr(p))); 1358 } 1359 } 1360 1361 // General case: use tessellated solid 1362 return G4TessellatedSolid::DistanceToOut(p); 1363 } 1364 1365 //_____________________________________________________________________________ 1366 // Get bounding box 1367 1368 void G4ExtrudedSolid::BoundingLimits(G4ThreeVector& pMin, 1369 G4ThreeVector& pMax) const 1370 { 1371 G4double xmin0 = kInfinity, xmax0 = -kInfinity; 1372 G4double ymin0 = kInfinity, ymax0 = -kInfinity; 1373 1374 for (G4int i=0; i<GetNofVertices(); ++i) 1375 { 1376 G4double x = fPolygon[i].x(); 1377 if (x < xmin0) xmin0 = x; 1378 if (x > xmax0) xmax0 = x; 1379 G4double y = fPolygon[i].y(); 1380 if (y < ymin0) ymin0 = y; 1381 if (y > ymax0) ymax0 = y; 1382 } 1383 1384 G4double xmin = kInfinity, xmax = -kInfinity; 1385 G4double ymin = kInfinity, ymax = -kInfinity; 1386 1387 G4int nsect = GetNofZSections(); 1388 for (G4int i=0; i<nsect; ++i) 1389 { 1390 ZSection zsect = GetZSection(i); 1391 G4double dx = zsect.fOffset.x(); 1392 G4double dy = zsect.fOffset.y(); 1393 G4double scale = zsect.fScale; 1394 xmin = std::min(xmin,xmin0*scale+dx); 1395 xmax = std::max(xmax,xmax0*scale+dx); 1396 ymin = std::min(ymin,ymin0*scale+dy); 1397 ymax = std::max(ymax,ymax0*scale+dy); 1398 } 1399 1400 G4double zmin = GetZSection(0).fZ; 1401 G4double zmax = GetZSection(nsect-1).fZ; 1402 1403 pMin.set(xmin,ymin,zmin); 1404 pMax.set(xmax,ymax,zmax); 1405 1406 // Check correctness of the bounding box 1407 // 1408 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z()) 1409 { 1410 std::ostringstream message; 1411 message << "Bad bounding box (min >= max) for solid: " 1412 << GetName() << " !" 1413 << "\npMin = " << pMin 1414 << "\npMax = " << pMax; 1415 G4Exception("G4ExtrudedSolid::BoundingLimits()", 1416 "GeomMgt0001", JustWarning, message); 1417 DumpInfo(); 1418 } 1419 } 1420 1421 //_____________________________________________________________________________ 1422 // Calculate extent under transform and specified limit 1423 1424 G4bool 1425 G4ExtrudedSolid::CalculateExtent(const EAxis pAxis, 1426 const G4VoxelLimits& pVoxelLimit, 1427 const G4AffineTransform& pTransform, 1428 G4double& pMin, G4double& pMax) const 1429 { 1430 G4ThreeVector bmin, bmax; 1431 G4bool exist; 1432 1433 // Check bounding box (bbox) 1434 // 1435 BoundingLimits(bmin,bmax); 1436 G4BoundingEnvelope bbox(bmin,bmax); 1437 #ifdef G4BBOX_EXTENT 1438 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 1439 #endif 1440 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax)) 1441 { 1442 return exist = pMin < pMax; 1443 } 1444 1445 // To find the extent, the base polygon is subdivided in triangles. 1446 // The extent is calculated as cumulative extent of the parts 1447 // formed by extrusion of the triangles 1448 // 1449 G4TwoVectorList triangles; 1450 G4double eminlim = pVoxelLimit.GetMinExtent(pAxis); 1451 G4double emaxlim = pVoxelLimit.GetMaxExtent(pAxis); 1452 1453 // triangulate the base polygon 1454 if (!G4GeomTools::TriangulatePolygon(fPolygon,triangles)) 1455 { 1456 std::ostringstream message; 1457 message << "Triangulation of the base polygon has failed for solid: " 1458 << GetName() << " !" 1459 << "\nExtent has been calculated using boundary box"; 1460 G4Exception("G4ExtrudedSolid::CalculateExtent()", 1461 "GeomMgt1002",JustWarning,message); 1462 return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax); 1463 } 1464 1465 // allocate vector lists 1466 G4int nsect = GetNofZSections(); 1467 std::vector<const G4ThreeVectorList *> polygons; 1468 polygons.resize(nsect); 1469 for (G4int k=0; k<nsect; ++k) { polygons[k] = new G4ThreeVectorList(3); } 1470 1471 // main loop along triangles 1472 pMin = kInfinity; 1473 pMax = -kInfinity; 1474 G4int ntria = (G4int)triangles.size()/3; 1475 for (G4int i=0; i<ntria; ++i) 1476 { 1477 G4int i3 = i*3; 1478 for (G4int k=0; k<nsect; ++k) // extrude triangle 1479 { 1480 ZSection zsect = GetZSection(k); 1481 G4double z = zsect.fZ; 1482 G4double dx = zsect.fOffset.x(); 1483 G4double dy = zsect.fOffset.y(); 1484 G4double scale = zsect.fScale; 1485 1486 auto ptr = const_cast<G4ThreeVectorList*>(polygons[k]); 1487 auto iter = ptr->begin(); 1488 G4double x0 = triangles[i3+0].x()*scale+dx; 1489 G4double y0 = triangles[i3+0].y()*scale+dy; 1490 iter->set(x0,y0,z); 1491 iter++; 1492 G4double x1 = triangles[i3+1].x()*scale+dx; 1493 G4double y1 = triangles[i3+1].y()*scale+dy; 1494 iter->set(x1,y1,z); 1495 iter++; 1496 G4double x2 = triangles[i3+2].x()*scale+dx; 1497 G4double y2 = triangles[i3+2].y()*scale+dy; 1498 iter->set(x2,y2,z); 1499 } 1500 1501 // set sub-envelope and adjust extent 1502 G4double emin,emax; 1503 G4BoundingEnvelope benv(polygons); 1504 if (!benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax)) continue; 1505 if (emin < pMin) pMin = emin; 1506 if (emax > pMax) pMax = emax; 1507 if (eminlim > pMin && emaxlim < pMax) break; // max possible extent 1508 } 1509 // free memory 1510 for (G4int k=0; k<nsect; ++k) { delete polygons[k]; polygons[k]=nullptr;} 1511 return (pMin < pMax); 1512 } 1513 1514 //_____________________________________________________________________________ 1515 1516 std::ostream& G4ExtrudedSolid::StreamInfo(std::ostream &os) const 1517 { 1518 G4long oldprc = os.precision(16); 1519 os << "-----------------------------------------------------------\n" 1520 << " *** Dump for solid - " << GetName() << " ***\n" 1521 << " ===================================================\n" 1522 << " Solid geometry type: " << fGeometryType << G4endl; 1523 1524 if ( fIsConvex) 1525 { os << " Convex polygon; list of vertices:" << G4endl; } 1526 else 1527 { os << " Concave polygon; list of vertices:" << G4endl; } 1528 1529 for ( std::size_t i=0; i<fNv; ++i ) 1530 { 1531 os << std::setw(5) << "#" << i 1532 << " vx = " << fPolygon[i].x()/mm << " mm" 1533 << " vy = " << fPolygon[i].y()/mm << " mm" << G4endl; 1534 } 1535 1536 os << " Sections:" << G4endl; 1537 for ( std::size_t iz=0; iz<fNz; ++iz ) 1538 { 1539 os << " z = " << fZSections[iz].fZ/mm << " mm " 1540 << " x0= " << fZSections[iz].fOffset.x()/mm << " mm " 1541 << " y0= " << fZSections[iz].fOffset.y()/mm << " mm " 1542 << " scale= " << fZSections[iz].fScale << G4endl; 1543 } 1544 1545 /* 1546 // Triangles (for debugging) 1547 os << G4endl; 1548 os << " Triangles:" << G4endl; 1549 os << " Triangle # vertex1 vertex2 vertex3" << G4endl; 1550 1551 G4int counter = 0; 1552 std::vector< std::vector<G4int> >::const_iterator it; 1553 for ( it = fTriangles.begin(); it != fTriangles.end(); it++ ) { 1554 std::vector<G4int> triangle = *it; 1555 os << std::setw(10) << counter++ 1556 << std::setw(10) << triangle[0] << std::setw(10) << triangle[1] 1557 << std::setw(10) << triangle[2] 1558 << G4endl; 1559 } 1560 */ 1561 os.precision(oldprc); 1562 1563 return os; 1564 } 1565 1566 #endif 1567