00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061 #include "G4TriangularFacet.hh"
00062
00063 #include "Randomize.hh"
00064 #include "G4TessellatedGeometryAlgorithms.hh"
00065
00066 using namespace std;
00067
00069
00070
00071
00072
00073
00074
00075 G4TriangularFacet::G4TriangularFacet (const G4ThreeVector &vt0,
00076 const G4ThreeVector &vt1,
00077 const G4ThreeVector &vt2,
00078 G4FacetVertexType vertexType)
00079 : fSqrDist(0.)
00080 {
00081 fVertices = new vector<G4ThreeVector>(3);
00082
00083 SetVertex(0, vt0);
00084 if (vertexType == ABSOLUTE)
00085 {
00086 SetVertex(1, vt1);
00087 SetVertex(2, vt2);
00088 fE1 = vt1 - vt0;
00089 fE2 = vt2 - vt0;
00090 }
00091 else
00092 {
00093 SetVertex(1, vt0 + vt1);
00094 SetVertex(2, vt0 + vt2);
00095 fE1 = vt1;
00096 fE2 = vt2;
00097 }
00098
00099 for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
00100
00101 G4double eMag1 = fE1.mag();
00102 G4double eMag2 = fE2.mag();
00103 G4double eMag3 = (fE2-fE1).mag();
00104
00105 if (eMag1 <= kCarTolerance || eMag2 <= kCarTolerance
00106 || eMag3 <= kCarTolerance)
00107 {
00108 ostringstream message;
00109 message << "Length of sides of facet are too small." << G4endl
00110 << "fVertices[0] = " << GetVertex(0) << G4endl
00111 << "fVertices[1] = " << GetVertex(1) << G4endl
00112 << "fVertices[2] = " << GetVertex(2) << G4endl
00113 << "Side lengths = fVertices[0]->fVertices[1]" << eMag1 << G4endl
00114 << "Side lengths = fVertices[0]->fVertices[2]" << eMag2 << G4endl
00115 << "Side lengths = fVertices[1]->fVertices[2]" << eMag3;
00116 G4Exception("G4TriangularFacet::G4TriangularFacet()",
00117 "GeomSolids1001", JustWarning, message);
00118 fIsDefined = false;
00119 fSurfaceNormal.set(0,0,0);
00120 fA = fB = fC = 0.0;
00121 fDet = 0.0;
00122 fArea = fRadius = 0.0;
00123 }
00124 else
00125 {
00126 fIsDefined = true;
00127 fSurfaceNormal = fE1.cross(fE2).unit();
00128 fA = fE1.mag2();
00129 fB = fE1.dot(fE2);
00130 fC = fE2.mag2();
00131 fDet = fabs(fA*fC - fB*fB);
00132
00133
00134
00135
00136
00137
00138 fArea = 0.5 * (fE1.cross(fE2)).mag();
00139 G4double lambda0 = (fA-fB) * fC / (8.0*fArea*fArea);
00140 G4double lambda1 = (fC-fB) * fA / (8.0*fArea*fArea);
00141 G4ThreeVector p0 = GetVertex(0);
00142 fCircumcentre = p0 + lambda0*fE1 + lambda1*fE2;
00143 G4double radiusSqr = (fCircumcentre-p0).mag2();
00144 fRadius = sqrt(radiusSqr);
00145 }
00146 }
00147
00149
00150 G4TriangularFacet::G4TriangularFacet ()
00151 : fSqrDist(0.)
00152 {
00153 fVertices = new vector<G4ThreeVector>(3);
00154 G4ThreeVector zero(0,0,0);
00155 SetVertex(0, zero);
00156 SetVertex(1, zero);
00157 SetVertex(2, zero);
00158 for (G4int i = 0; i < 3; ++i) fIndices[i] = -1;
00159 fIsDefined = false;
00160 fSurfaceNormal.set(0,0,0);
00161 fA = fB = fC = 0;
00162 fE1 = zero;
00163 fE2 = zero;
00164 fDet = 0.0;
00165 fArea = fRadius = 0.0;
00166 }
00167
00169
00170 G4TriangularFacet::~G4TriangularFacet ()
00171 {
00172 SetVertices(0);
00173 }
00174
00176
00177 void G4TriangularFacet::CopyFrom (const G4TriangularFacet &rhs)
00178 {
00179 char *p = (char *) &rhs;
00180 copy(p, p + sizeof(*this), (char *)this);
00181
00182 if (fIndices[0] < 0 && fVertices)
00183 {
00184 fVertices = new vector<G4ThreeVector>(3);
00185 for (G4int i = 0; i < 3; ++i) (*fVertices)[i] = (*rhs.fVertices)[i];
00186 }
00187 }
00188
00190
00191 G4TriangularFacet::G4TriangularFacet (const G4TriangularFacet &rhs)
00192 : G4VFacet(rhs)
00193 {
00194 CopyFrom(rhs);
00195 }
00196
00198
00199 G4TriangularFacet &
00200 G4TriangularFacet::operator=(const G4TriangularFacet &rhs)
00201 {
00202 SetVertices(0);
00203
00204 if (this != &rhs)
00205 CopyFrom(rhs);
00206
00207 return *this;
00208 }
00209
00211
00212
00213
00214
00215
00216 G4VFacet *G4TriangularFacet::GetClone ()
00217 {
00218 G4TriangularFacet *fc =
00219 new G4TriangularFacet (GetVertex(0), GetVertex(1), GetVertex(2), ABSOLUTE);
00220 return fc;
00221 }
00222
00224
00225
00226
00227
00228
00229
00230 G4TriangularFacet *G4TriangularFacet::GetFlippedFacet ()
00231 {
00232 G4TriangularFacet *flipped =
00233 new G4TriangularFacet (GetVertex(0), GetVertex(1), GetVertex(2), ABSOLUTE);
00234 return flipped;
00235 }
00236
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 G4ThreeVector G4TriangularFacet::Distance (const G4ThreeVector &p)
00252 {
00253 G4ThreeVector D = GetVertex(0) - p;
00254 G4double d = fE1.dot(D);
00255 G4double e = fE2.dot(D);
00256 G4double f = D.mag2();
00257 G4double q = fB*e - fC*d;
00258 G4double t = fB*d - fA*e;
00259 fSqrDist = 0.;
00260
00261 if (q+t <= fDet)
00262 {
00263 if (q < 0.0)
00264 {
00265 if (t < 0.0)
00266 {
00267
00268
00269
00270 if (d < 0.0)
00271 {
00272 t = 0.0;
00273 if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
00274 else {q = -d/fA; fSqrDist = d*q + f;}
00275 }
00276 else
00277 {
00278 q = 0.0;
00279 if (e >= 0.0) {t = 0.0; fSqrDist = f;}
00280 else if (-e >= fC) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
00281 else {t = -e/fC; fSqrDist = e*t + f;}
00282 }
00283 }
00284 else
00285 {
00286
00287
00288
00289 q = 0.0;
00290 if (e >= 0.0) {t = 0.0; fSqrDist = f;}
00291 else if (-e >= fC) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
00292 else {t = -e/fC; fSqrDist = e*t + f;}
00293 }
00294 }
00295 else if (t < 0.0)
00296 {
00297
00298
00299
00300 t = 0.0;
00301 if (d >= 0.0) {q = 0.0; fSqrDist = f;}
00302 else if (-d >= fA) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
00303 else {q = -d/fA; fSqrDist = d*q + f;}
00304 }
00305 else
00306 {
00307
00308
00309
00310 q = q / fDet;
00311 t = t / fDet;
00312 fSqrDist = q*(fA*q + fB*t + 2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
00313 }
00314 }
00315 else
00316 {
00317 if (q < 0.0)
00318 {
00319
00320
00321
00322 G4double tmp0 = fB + d;
00323 G4double tmp1 = fC + e;
00324 if (tmp1 > tmp0)
00325 {
00326 G4double numer = tmp1 - tmp0;
00327 G4double denom = fA - 2.0*fB + fC;
00328 if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
00329 else
00330 {
00331 q = numer/denom;
00332 t = 1.0 - q;
00333 fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
00334 }
00335 }
00336 else
00337 {
00338 q = 0.0;
00339 if (tmp1 <= 0.0) {t = 1.0; fSqrDist = fC + 2.0*e + f;}
00340 else if (e >= 0.0) {t = 0.0; fSqrDist = f;}
00341 else {t = -e/fC; fSqrDist = e*t + f;}
00342 }
00343 }
00344 else if (t < 0.0)
00345 {
00346
00347
00348
00349 G4double tmp0 = fB + e;
00350 G4double tmp1 = fA + d;
00351 if (tmp1 > tmp0)
00352 {
00353 G4double numer = tmp1 - tmp0;
00354 G4double denom = fA - 2.0*fB + fC;
00355 if (numer >= denom) {t = 1.0; q = 0.0; fSqrDist = fC + 2.0*e + f;}
00356 else
00357 {
00358 t = numer/denom;
00359 q = 1.0 - t;
00360 fSqrDist = q*(fA*q + fB*t +2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
00361 }
00362 }
00363 else
00364 {
00365 t = 0.0;
00366 if (tmp1 <= 0.0) {q = 1.0; fSqrDist = fA + 2.0*d + f;}
00367 else if (d >= 0.0) {q = 0.0; fSqrDist = f;}
00368 else {q = -d/fA; fSqrDist = d*q + f;}
00369 }
00370 }
00371 else
00372
00373
00374
00375 {
00376 G4double numer = fC + e - fB - d;
00377 if (numer <= 0.0)
00378 {
00379 q = 0.0;
00380 t = 1.0;
00381 fSqrDist = fC + 2.0*e + f;
00382 }
00383 else
00384 {
00385 G4double denom = fA - 2.0*fB + fC;
00386 if (numer >= denom) {q = 1.0; t = 0.0; fSqrDist = fA + 2.0*d + f;}
00387 else
00388 {
00389 q = numer/denom;
00390 t = 1.0 - q;
00391 fSqrDist = q*(fA*q + fB*t + 2.0*d) + t*(fB*q + fC*t + 2.0*e) + f;
00392 }
00393 }
00394 }
00395 }
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413 if (fSqrDist < 0.0) fSqrDist = 0.;
00414 G4ThreeVector u = D + q*fE1 + t*fE2;
00415 G4double u2 = u.mag2();
00416
00417
00418
00419
00420 if (fSqrDist > u2) fSqrDist = u2;
00421
00422 return u;
00423 }
00424
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435 G4double G4TriangularFacet::Distance (const G4ThreeVector &p,
00436 G4double minDist)
00437 {
00438
00439
00440
00441
00442
00443 G4double dist = kInfinity;
00444 if ((p-fCircumcentre).mag()-fRadius < minDist)
00445 {
00446
00447
00448
00449
00450 dist = Distance(p).mag();
00451 }
00452 return dist;
00453 }
00454
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470 G4double G4TriangularFacet::Distance (const G4ThreeVector &p,
00471 G4double minDist,
00472 const G4bool outgoing)
00473 {
00474
00475
00476
00477
00478
00479 G4double dist = kInfinity;
00480 if ((p-fCircumcentre).mag()-fRadius < minDist)
00481 {
00482
00483
00484
00485
00486 G4ThreeVector v = Distance(p);
00487 G4double dist1 = sqrt(fSqrDist);
00488 G4double dir = v.dot(fSurfaceNormal);
00489 G4bool wrongSide = (dir > 0.0 && !outgoing) || (dir < 0.0 && outgoing);
00490 if (dist1 <= kCarTolerance)
00491 {
00492
00493
00494
00495
00496 if (wrongSide) dist = 0.0;
00497 else dist = dist1;
00498 }
00499 else if (!wrongSide) dist = dist1;
00500 }
00501 return dist;
00502 }
00503
00505
00506
00507
00508
00509
00510
00511 G4double G4TriangularFacet::Extent (const G4ThreeVector axis)
00512 {
00513 G4double ss = GetVertex(0).dot(axis);
00514 G4double sp = GetVertex(1).dot(axis);
00515 if (sp > ss) ss = sp;
00516 sp = GetVertex(2).dot(axis);
00517 if (sp > ss) ss = sp;
00518 return ss;
00519 }
00520
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548 G4bool G4TriangularFacet::Intersect (const G4ThreeVector &p,
00549 const G4ThreeVector &v,
00550 G4bool outgoing,
00551 G4double &distance,
00552 G4double &distFromSurface,
00553 G4ThreeVector &normal)
00554 {
00555
00556
00557
00558
00559
00560 G4double w = v.dot(fSurfaceNormal);
00561 if ((outgoing && w < -dirTolerance) || (!outgoing && w > dirTolerance))
00562 {
00563 distance = kInfinity;
00564 distFromSurface = kInfinity;
00565 normal.set(0,0,0);
00566 return false;
00567 }
00568
00569
00570
00571
00572
00573
00574 const G4ThreeVector &p0 = GetVertex(0);
00575 G4ThreeVector D = p0 - p;
00576 distFromSurface = D.dot(fSurfaceNormal);
00577 G4bool wrongSide = (outgoing && distFromSurface < -0.5*kCarTolerance) ||
00578 (!outgoing && distFromSurface > 0.5*kCarTolerance);
00579 if (wrongSide)
00580 {
00581 distance = kInfinity;
00582 distFromSurface = kInfinity;
00583 normal.set(0,0,0);
00584 return false;
00585 }
00586
00587 wrongSide = (outgoing && distFromSurface < 0.0)
00588 || (!outgoing && distFromSurface > 0.0);
00589 if (wrongSide)
00590 {
00591
00592
00593
00594
00595 G4ThreeVector u = Distance(p);
00596 if (fSqrDist <= kCarTolerance*kCarTolerance)
00597 {
00598
00599
00600
00601
00602
00603 distance = 0.0;
00604 normal = fSurfaceNormal;
00605 return true;
00606 }
00607 else
00608 {
00609
00610
00611
00612
00613
00614 distance = kInfinity;
00615 distFromSurface = kInfinity;
00616 normal.set(0,0,0);
00617 return false;
00618 }
00619 }
00620 if (w < dirTolerance && w > -dirTolerance)
00621 {
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635 G4ThreeVector mu = fE1.unit();
00636 G4ThreeVector nu = fSurfaceNormal.cross(mu);
00637 G4TwoVector pprime(p.dot(mu), p.dot(nu));
00638 G4TwoVector vprime(v.dot(mu), v.dot(nu));
00639 G4TwoVector P0prime(p0.dot(mu), p0.dot(nu));
00640 G4TwoVector E0prime(fE1.mag(), 0.0);
00641 G4TwoVector E1prime(fE2.dot(mu), fE2.dot(nu));
00642 G4TwoVector loc[2];
00643 if (G4TessellatedGeometryAlgorithms::IntersectLineAndTriangle2D(pprime,
00644 vprime, P0prime, E0prime, E1prime, loc))
00645 {
00646
00647
00648
00649
00650
00651 G4double vprimemag = vprime.mag();
00652 G4double s0 = (loc[0] - pprime).mag()/vprimemag;
00653 G4double s1 = (loc[1] - pprime).mag()/vprimemag;
00654 G4double normDist0 = fSurfaceNormal.dot(s0*v) - distFromSurface;
00655 G4double normDist1 = fSurfaceNormal.dot(s1*v) - distFromSurface;
00656
00657 if ((normDist0 < 0.0 && normDist1 < 0.0)
00658 || (normDist0 > 0.0 && normDist1 > 0.0)
00659 || (normDist0 == 0.0 && normDist1 == 0.0) )
00660 {
00661 distance = kInfinity;
00662 distFromSurface = kInfinity;
00663 normal.set(0,0,0);
00664 return false;
00665 }
00666 else
00667 {
00668 G4double dnormDist = normDist1 - normDist0;
00669 if (fabs(dnormDist) < DBL_EPSILON)
00670 {
00671 distance = s0;
00672 normal = fSurfaceNormal;
00673 if (!outgoing) distFromSurface = -distFromSurface;
00674 return true;
00675 }
00676 else
00677 {
00678 distance = s0 - normDist0*(s1-s0)/dnormDist;
00679 normal = fSurfaceNormal;
00680 if (!outgoing) distFromSurface = -distFromSurface;
00681 return true;
00682 }
00683 }
00684 }
00685 else
00686 {
00687 distance = kInfinity;
00688 distFromSurface = kInfinity;
00689 normal.set(0,0,0);
00690 return false;
00691 }
00692 }
00693
00694
00695
00696
00697
00698
00699
00700 distance = distFromSurface / w;
00701 G4ThreeVector pp = p + v*distance;
00702 G4ThreeVector DD = p0 - pp;
00703 G4double d = fE1.dot(DD);
00704 G4double e = fE2.dot(DD);
00705 G4double ss = fB*e - fC*d;
00706 G4double t = fB*d - fA*e;
00707
00708 G4double sTolerance = (fabs(fB)+ fabs(fC) + fabs(d) + fabs(e))*kCarTolerance;
00709 G4double tTolerance = (fabs(fA)+ fabs(fB) + fabs(d) + fabs(e))*kCarTolerance;
00710 G4double detTolerance = (fabs(fA)+ fabs(fC) + 2*fabs(fB) )*kCarTolerance;
00711
00712
00713 if (ss < -sTolerance || t < -tTolerance || ( ss+t - fDet ) > detTolerance)
00714 {
00715
00716
00717
00718 distance = distFromSurface = kInfinity;
00719 normal.set(0,0,0);
00720 return false;
00721 }
00722 else
00723 {
00724
00725
00726
00727 normal = fSurfaceNormal;
00728 if (!outgoing) distFromSurface = -distFromSurface;
00729 return true;
00730 }
00731 }
00732
00734
00735
00736
00737
00738
00739 G4ThreeVector G4TriangularFacet::GetPointOnFace() const
00740 {
00741 G4double alpha = G4RandFlat::shoot(0., 1.);
00742 G4double beta = G4RandFlat::shoot(0., 1.);
00743 G4double lambda1 = alpha*beta;
00744 G4double lambda0 = alpha-lambda1;
00745
00746 return GetVertex(0) + lambda0*fE1 + lambda1*fE2;
00747 }
00748
00750
00751
00752
00753
00754
00755 G4double G4TriangularFacet::GetArea()
00756 {
00757 return fArea;
00758 }
00759
00761
00762 G4GeometryType G4TriangularFacet::GetEntityType () const
00763 {
00764 return "G4TriangularFacet";
00765 }
00766
00768
00769 G4ThreeVector G4TriangularFacet::GetSurfaceNormal () const
00770 {
00771 return fSurfaceNormal;
00772 }
00773
00775
00776 void G4TriangularFacet::SetSurfaceNormal (G4ThreeVector normal)
00777 {
00778 fSurfaceNormal = normal;
00779 }