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 #include "G4PolyPhiFace.hh"
00042 #include "G4ClippablePolygon.hh"
00043 #include "G4ReduciblePolygon.hh"
00044 #include "G4AffineTransform.hh"
00045 #include "G4SolidExtentList.hh"
00046 #include "G4GeometryTolerance.hh"
00047
00048 #include "Randomize.hh"
00049 #include "G4TwoVector.hh"
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061 G4PolyPhiFace::G4PolyPhiFace( const G4ReduciblePolygon *rz,
00062 G4double phi,
00063 G4double deltaPhi,
00064 G4double phiOther )
00065 : fSurfaceArea(0.), triangles(0)
00066 {
00067 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00068
00069 numEdges = rz->NumVertices();
00070
00071 rMin = rz->Amin();
00072 rMax = rz->Amax();
00073 zMin = rz->Bmin();
00074 zMax = rz->Bmax();
00075
00076
00077
00078
00079 G4bool start = (phiOther > phi);
00080
00081
00082
00083
00084 radial = G4ThreeVector( std::cos(phi), std::sin(phi), 0.0 );
00085
00086
00087
00088
00089 G4double zSign = start ? 1 : -1;
00090 normal = G4ThreeVector( zSign*radial.y(), -zSign*radial.x(), 0 );
00091
00092
00093
00094
00095 allBehind = (zSign*(std::cos(phiOther)*radial.y() - std::sin(phiOther)*radial.x()) < 0);
00096
00097
00098
00099
00100 G4double midPhi = phi + (start ? +0.5 : -0.5)*deltaPhi;
00101 G4double cosMid = std::cos(midPhi),
00102 sinMid = std::sin(midPhi);
00103
00104
00105
00106
00107 corners = new G4PolyPhiFaceVertex[numEdges];
00108
00109
00110
00111 G4ReduciblePolygonIterator iterRZ(rz);
00112
00113 G4PolyPhiFaceVertex *corn = corners;
00114 G4PolyPhiFaceVertex *helper=corners;
00115
00116 iterRZ.Begin();
00117 do
00118 {
00119 corn->r = iterRZ.GetA();
00120 corn->z = iterRZ.GetB();
00121 corn->x = corn->r*radial.x();
00122 corn->y = corn->r*radial.y();
00123
00124
00125
00126 if( corn == corners )
00127 { corn->prev = corners+numEdges-1;}
00128 else
00129 { corn->prev = helper; }
00130
00131
00132
00133 if( corn < corners+numEdges-1 )
00134 { corn->next = corn+1;}
00135 else
00136 { corn->next = corners; }
00137
00138 helper = corn;
00139 } while( ++corn, iterRZ.Next() );
00140
00141
00142
00143
00144 edges = new G4PolyPhiFaceEdge[numEdges];
00145
00146
00147
00148
00149 G4double rFact = std::cos(0.5*deltaPhi);
00150 G4double rFactNormalize = 1.0/std::sqrt(1.0+rFact*rFact);
00151
00152 G4PolyPhiFaceVertex *prev = corners+numEdges-1,
00153 *here = corners;
00154 G4PolyPhiFaceEdge *edge = edges;
00155 do
00156 {
00157 G4ThreeVector sideNorm;
00158
00159 edge->v0 = prev;
00160 edge->v1 = here;
00161
00162 G4double dr = here->r - prev->r,
00163 dz = here->z - prev->z;
00164
00165 edge->length = std::sqrt( dr*dr + dz*dz );
00166
00167 edge->tr = dr/edge->length;
00168 edge->tz = dz/edge->length;
00169
00170 if ((here->r < DBL_MIN) && (prev->r < DBL_MIN))
00171 {
00172
00173
00174
00175
00176
00177 G4double zSignOther = start ? -1 : 1;
00178 sideNorm = G4ThreeVector( zSignOther*std::sin(phiOther),
00179 -zSignOther*std::cos(phiOther), 0 );
00180 }
00181 else
00182 {
00183 sideNorm = G4ThreeVector( edge->tz*cosMid,
00184 edge->tz*sinMid,
00185 -edge->tr*rFact );
00186 sideNorm *= rFactNormalize;
00187 }
00188 sideNorm += normal;
00189
00190 edge->norm3D = sideNorm.unit();
00191 } while( edge++, prev=here, ++here < corners+numEdges );
00192
00193
00194
00195
00196 G4PolyPhiFaceEdge *prevEdge = edges+numEdges-1;
00197 edge = edges;
00198 do
00199 {
00200
00201
00202
00203 G4double rPart = prevEdge->tr + edge->tr;
00204 G4double zPart = prevEdge->tz + edge->tz;
00205 G4double norm = std::sqrt( rPart*rPart + zPart*zPart );
00206 G4double rNorm = +zPart/norm;
00207 G4double zNorm = -rPart/norm;
00208
00209 edge->v0->rNorm = rNorm;
00210 edge->v0->zNorm = zNorm;
00211
00212
00213
00214
00215
00216
00217
00218 G4ThreeVector xyVector;
00219
00220 if (edge->v0->r < DBL_MIN)
00221 {
00222
00223
00224
00225
00226
00227
00228
00229 G4double zSignOther = start ? -1 : 1;
00230 G4ThreeVector normalOther( zSignOther*std::sin(phiOther),
00231 -zSignOther*std::cos(phiOther), 0 );
00232
00233 xyVector = - normal - normalOther;
00234 }
00235 else
00236 {
00237
00238
00239
00240
00241
00242 xyVector = G4ThreeVector( cosMid, sinMid, 0 );
00243 if (rNorm < 0)
00244 xyVector -= normal;
00245 else
00246 xyVector += normal;
00247 }
00248
00249
00250
00251
00252 edge->v0->norm3D = rNorm*xyVector.unit() + G4ThreeVector( 0, 0, zNorm );
00253 } while( prevEdge=edge, ++edge < edges+numEdges );
00254
00255
00256
00257
00258 G4double rAve = 0.5*(rMax-rMin),
00259 zAve = 0.5*(zMax-zMin);
00260 surface = G4ThreeVector( rAve*radial.x(), rAve*radial.y(), zAve );
00261 }
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272 void G4PolyPhiFace::Diagnose( G4VSolid *owner )
00273 {
00274 G4PolyPhiFaceVertex *corner = corners;
00275 do
00276 {
00277 G4ThreeVector test(corner->x, corner->y, corner->z);
00278 test -= 1E-6*corner->norm3D;
00279
00280 if (owner->Inside(test) != kInside)
00281 G4Exception( "G4PolyPhiFace::Diagnose()", "GeomSolids0002",
00282 FatalException, "Bad vertex normal found." );
00283 } while( ++corner < corners+numEdges );
00284 }
00285
00286
00287
00288
00289
00290
00291 G4PolyPhiFace::G4PolyPhiFace( __void__&)
00292 : numEdges(0), edges(0), corners(0), rMin(0.), rMax(0.), zMin(0.), zMax(0.),
00293 allBehind(false), kCarTolerance(0.), fSurfaceArea(0.), triangles(0)
00294 {
00295 }
00296
00297
00298
00299
00300
00301 G4PolyPhiFace::~G4PolyPhiFace()
00302 {
00303 delete [] edges;
00304 delete [] corners;
00305 }
00306
00307
00308
00309
00310
00311 G4PolyPhiFace::G4PolyPhiFace( const G4PolyPhiFace &source )
00312 : G4VCSGface()
00313 {
00314 CopyStuff( source );
00315 }
00316
00317
00318
00319
00320
00321 G4PolyPhiFace& G4PolyPhiFace::operator=( const G4PolyPhiFace &source )
00322 {
00323 if (this == &source) { return *this; }
00324
00325 delete [] edges;
00326 delete [] corners;
00327
00328 CopyStuff( source );
00329
00330 return *this;
00331 }
00332
00333
00334
00335
00336
00337 void G4PolyPhiFace::CopyStuff( const G4PolyPhiFace &source )
00338 {
00339
00340
00341
00342 numEdges = source.numEdges;
00343 normal = source.normal;
00344 radial = source.radial;
00345 surface = source.surface;
00346 rMin = source.rMin;
00347 rMax = source.rMax;
00348 zMin = source.zMin;
00349 zMax = source.zMax;
00350 allBehind = source.allBehind;
00351 triangles = 0;
00352
00353 kCarTolerance = source.kCarTolerance;
00354 fSurfaceArea = source.fSurfaceArea;
00355
00356
00357
00358
00359 corners = new G4PolyPhiFaceVertex[numEdges];
00360 G4PolyPhiFaceVertex *corn = corners,
00361 *sourceCorn = source.corners;
00362 do
00363 {
00364 *corn = *sourceCorn;
00365 } while( ++sourceCorn, ++corn < corners+numEdges );
00366
00367
00368
00369
00370 edges = new G4PolyPhiFaceEdge[numEdges];
00371
00372 G4PolyPhiFaceVertex *prev = corners+numEdges-1,
00373 *here = corners;
00374 G4PolyPhiFaceEdge *edge = edges,
00375 *sourceEdge = source.edges;
00376 do
00377 {
00378 *edge = *sourceEdge;
00379 edge->v0 = prev;
00380 edge->v1 = here;
00381 } while( ++sourceEdge, ++edge, prev=here, ++here < corners+numEdges );
00382 }
00383
00384
00385
00386
00387
00388 G4bool G4PolyPhiFace::Intersect( const G4ThreeVector &p,
00389 const G4ThreeVector &v,
00390 G4bool outgoing,
00391 G4double surfTolerance,
00392 G4double &distance,
00393 G4double &distFromSurface,
00394 G4ThreeVector &aNormal,
00395 G4bool &isAllBehind )
00396 {
00397 G4double normSign = outgoing ? +1 : -1;
00398
00399
00400
00401
00402 isAllBehind = allBehind;
00403 aNormal = normal;
00404
00405
00406
00407
00408
00409 G4double dotProd = normSign*normal.dot(v);
00410
00411 if (dotProd <= 0) return false;
00412
00413
00414
00415
00416
00417 G4ThreeVector ps = p - surface;
00418 distFromSurface = -normSign*ps.dot(normal);
00419
00420 if (distFromSurface < -surfTolerance) return false;
00421
00422
00423
00424
00425
00426 distance = distFromSurface/dotProd;
00427
00428
00429
00430
00431 G4ThreeVector ip = p + distance*v;
00432
00433 G4double r = radial.dot(ip);
00434
00435
00436
00437
00438 return InsideEdgesExact( r, ip.z(), normSign, p, v );
00439 }
00440
00441
00442
00443
00444
00445 G4double G4PolyPhiFace::Distance( const G4ThreeVector &p, G4bool outgoing )
00446 {
00447 G4double normSign = outgoing ? +1 : -1;
00448
00449
00450
00451 G4ThreeVector ps = p - surface;
00452 G4double distPhi = -normSign*normal.dot(ps);
00453
00454 if (distPhi < -0.5*kCarTolerance)
00455 return kInfinity;
00456 else if (distPhi < 0)
00457 distPhi = 0.0;
00458
00459
00460
00461
00462 G4double r = radial.dot(p);
00463
00464
00465
00466
00467 G4double distRZ2;
00468
00469 if (InsideEdges( r, p.z(), &distRZ2, 0 ))
00470 {
00471
00472
00473
00474 return distPhi;
00475 }
00476 else
00477 {
00478
00479
00480
00481 return std::sqrt( distPhi*distPhi + distRZ2 );
00482 }
00483 }
00484
00485
00486
00487
00488
00489 EInside G4PolyPhiFace::Inside( const G4ThreeVector &p,
00490 G4double tolerance,
00491 G4double *bestDistance )
00492 {
00493
00494
00495
00496
00497 G4ThreeVector ps = p - surface;
00498 G4double distPhi = normal.dot(ps);
00499
00500
00501
00502
00503 G4double r = radial.dot(p);
00504
00505
00506
00507
00508 G4double distRZ2;
00509 G4PolyPhiFaceVertex *base3Dnorm;
00510 G4ThreeVector *head3Dnorm;
00511
00512 if (InsideEdges( r, p.z(), &distRZ2, &base3Dnorm, &head3Dnorm ))
00513 {
00514
00515
00516
00517 *bestDistance = std::fabs(distPhi);
00518
00519
00520
00521
00522 if (distPhi < -tolerance) return kInside;
00523 if (distPhi < tolerance) return kSurface;
00524 return kOutside;
00525 }
00526 else
00527 {
00528
00529
00530
00531
00532 *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
00533
00534
00535
00536
00537 G4ThreeVector cc( base3Dnorm->r*radial.x(),
00538 base3Dnorm->r*radial.y(),
00539 base3Dnorm->z );
00540 cc = p - cc;
00541 G4double normDist = head3Dnorm->dot(cc);
00542 if ( distRZ2 > tolerance*tolerance )
00543 {
00544
00545
00546
00547 return normDist < 0 ? kInside : kOutside;
00548 }
00549
00550 if (normDist < -tolerance) return kInside;
00551 if (normDist < tolerance) return kSurface;
00552 return kOutside;
00553 }
00554 }
00555
00556
00557
00558
00559
00560
00561
00562
00563 G4ThreeVector G4PolyPhiFace::Normal( const G4ThreeVector &p,
00564 G4double *bestDistance )
00565 {
00566
00567
00568
00569
00570 G4double distPhi = normal.dot(p);
00571
00572
00573
00574
00575 G4double r = radial.dot(p);
00576
00577
00578
00579
00580 G4double distRZ2;
00581
00582 if (InsideEdges( r, p.z(), &distRZ2, 0 ))
00583 {
00584
00585
00586
00587 *bestDistance = std::fabs(distPhi);
00588 }
00589 else
00590 {
00591
00592
00593
00594 *bestDistance = std::sqrt( distPhi*distPhi + distRZ2 );
00595 }
00596
00597 return normal;
00598 }
00599
00600
00601
00602
00603
00604
00605
00606 G4double G4PolyPhiFace::Extent( const G4ThreeVector axis )
00607 {
00608 G4double max = -kInfinity;
00609
00610 G4PolyPhiFaceVertex *corner = corners;
00611 do
00612 {
00613 G4double here = axis.x()*corner->r*radial.x()
00614 + axis.y()*corner->r*radial.y()
00615 + axis.z()*corner->z;
00616 if (here > max) max = here;
00617 } while( ++corner < corners + numEdges );
00618
00619 return max;
00620 }
00621
00622
00623
00624
00625
00626
00627
00628 void G4PolyPhiFace::CalculateExtent( const EAxis axis,
00629 const G4VoxelLimits &voxelLimit,
00630 const G4AffineTransform &transform,
00631 G4SolidExtentList &extentList )
00632 {
00633
00634
00635
00636
00637
00638 G4ClippablePolygon polygon;
00639
00640 G4PolyPhiFaceVertex *corner = corners;
00641 do
00642 {
00643 G4ThreeVector point( 0, 0, corner->z );
00644 point += radial*corner->r;
00645
00646 polygon.AddVertexInOrder( transform.TransformPoint( point ) );
00647 } while( ++corner < corners + numEdges );
00648
00649
00650
00651
00652 if (polygon.PartialClip( voxelLimit, axis ))
00653 {
00654
00655
00656
00657 polygon.SetNormal( transform.TransformAxis(normal) );
00658 extentList.AddSurface( polygon );
00659 }
00660 }
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683 G4bool G4PolyPhiFace::InsideEdgesExact( G4double r, G4double z,
00684 G4double normSign,
00685 const G4ThreeVector &p,
00686 const G4ThreeVector &v )
00687 {
00688
00689
00690
00691 if ( (r < rMin-kCarTolerance)
00692 || (r > rMax+kCarTolerance) ) return false;
00693
00694 if ( (z < zMin-kCarTolerance)
00695 || (z > zMax+kCarTolerance) ) return false;
00696
00697
00698
00699
00700 G4double qx = p.x() + v.x(),
00701 qy = p.y() + v.y(),
00702 qz = p.z() + v.z();
00703
00704 G4int answer = 0;
00705 G4PolyPhiFaceVertex *corn = corners,
00706 *prev = corners+numEdges-1;
00707
00708 G4double cornZ, prevZ;
00709
00710 prevZ = ExactZOrder( z, qx, qy, qz, v, normSign, prev );
00711 do
00712 {
00713
00714
00715
00716 cornZ = ExactZOrder( z, qx, qy, qz, v, normSign, corn );
00717
00718 if (cornZ < 0)
00719 {
00720 if (prevZ < 0) continue;
00721 }
00722 else if (cornZ > 0)
00723 {
00724 if (prevZ > 0) continue;
00725 }
00726 else
00727 {
00728
00729
00730
00731
00732
00733 if (prevZ == 0) continue;
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743 G4PolyPhiFaceVertex *next = corn;
00744 G4double nextZ;
00745 do
00746 {
00747 next++;
00748 if (next == corners+numEdges) next = corners;
00749
00750 nextZ = ExactZOrder( z, qx, qy, qz, v, normSign, next );
00751 } while( nextZ == 0 );
00752
00753
00754
00755
00756 if (nextZ*prevZ < 0) continue;
00757 }
00758
00759
00760
00761
00762
00763
00764
00765 G4ThreeVector qa( qx - prev->x, qy - prev->y, qz - prev->z ),
00766 qb( qx - corn->x, qy - corn->y, qz - corn->z );
00767
00768 G4double aboveOrBelow = normSign*qa.cross(qb).dot(v);
00769
00770 if (aboveOrBelow > 0)
00771 answer++;
00772 else if (aboveOrBelow < 0)
00773 answer--;
00774 else
00775 {
00776
00777
00778
00779
00780
00781 return true;
00782 }
00783 } while( prevZ = cornZ, prev=corn, ++corn < corners+numEdges );
00784
00785
00786
00787
00788
00789
00790
00791 return answer!=0;
00792 }
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822 G4bool G4PolyPhiFace::InsideEdges( G4double r, G4double z )
00823 {
00824
00825
00826
00827 if ( r < rMin || r > rMax ) return false;
00828 if ( z < zMin || z > zMax ) return false;
00829
00830
00831
00832
00833 G4double notUsed;
00834
00835 return InsideEdges( r, z, ¬Used, 0 );
00836 }
00837
00838
00839
00840
00841
00842
00843
00844 G4bool G4PolyPhiFace::InsideEdges( G4double r, G4double z,
00845 G4double *bestDist2,
00846 G4PolyPhiFaceVertex **base3Dnorm,
00847 G4ThreeVector **head3Dnorm )
00848 {
00849 G4double bestDistance2 = kInfinity;
00850 G4bool answer = 0;
00851
00852 G4PolyPhiFaceEdge *edge = edges;
00853 do
00854 {
00855 G4PolyPhiFaceVertex *testMe;
00856
00857
00858
00859 G4double dr = (r-edge->v0->r), dz = (z-edge->v0->z);
00860
00861 G4double distOut = dr*edge->tz - dz*edge->tr;
00862 G4double distance2 = distOut*distOut;
00863 if (distance2 > bestDistance2) continue;
00864
00865
00866
00867
00868 G4double q = dr*edge->tr + dz*edge->tz;
00869
00870
00871
00872
00873 if (q < 0)
00874 {
00875 distance2 += q*q;
00876 testMe = edge->v0;
00877 }
00878 else if (q > edge->length)
00879 {
00880 G4double s2 = q-edge->length;
00881 distance2 += s2*s2;
00882 testMe = edge->v1;
00883 }
00884 else
00885 {
00886 testMe = 0;
00887 }
00888
00889
00890
00891
00892 if (distance2 < bestDistance2)
00893 {
00894 bestDistance2 = distance2;
00895 if (testMe)
00896 {
00897 G4double distNorm = dr*testMe->rNorm + dz*testMe->zNorm;
00898 answer = (distNorm <= 0);
00899 if (base3Dnorm)
00900 {
00901 *base3Dnorm = testMe;
00902 *head3Dnorm = &testMe->norm3D;
00903 }
00904 }
00905 else
00906 {
00907 answer = (distOut <= 0);
00908 if (base3Dnorm)
00909 {
00910 *base3Dnorm = edge->v0;
00911 *head3Dnorm = &edge->norm3D;
00912 }
00913 }
00914 }
00915 } while( ++edge < edges + numEdges );
00916
00917 *bestDist2 = bestDistance2;
00918 return answer;
00919 }
00920
00921
00922
00923
00924
00925 G4double G4PolyPhiFace::SurfaceTriangle( G4ThreeVector p1,
00926 G4ThreeVector p2,
00927 G4ThreeVector p3,
00928 G4ThreeVector *p4 )
00929 {
00930 G4ThreeVector v, w;
00931
00932 v = p3 - p1;
00933 w = p1 - p2;
00934 G4double lambda1 = G4UniformRand();
00935 G4double lambda2 = lambda1*G4UniformRand();
00936
00937 *p4=p2 + lambda1*w + lambda2*v;
00938 return 0.5*(v.cross(w)).mag();
00939 }
00940
00941
00942
00943
00944 G4double G4PolyPhiFace::SurfaceArea()
00945 {
00946 if ( fSurfaceArea==0. ) { Triangulate(); }
00947 return fSurfaceArea;
00948 }
00949
00950
00951
00952
00953 G4ThreeVector G4PolyPhiFace::GetPointOnFace()
00954 {
00955 Triangulate();
00956 return surface_point;
00957 }
00958
00959
00960
00961
00962
00963
00964
00965
00966 G4double G4PolyPhiFace::Area2( G4TwoVector a,
00967 G4TwoVector b,
00968 G4TwoVector c )
00969 {
00970 return ((b.x()-a.x())*(c.y()-a.y())-
00971 (c.x()-a.x())*(b.y()-a.y()));
00972 }
00973
00974
00975
00976
00977 G4bool G4PolyPhiFace::Left( G4TwoVector a,
00978 G4TwoVector b,
00979 G4TwoVector c )
00980 {
00981 return Area2(a,b,c)>0;
00982 }
00983
00984
00985
00986
00987 G4bool G4PolyPhiFace::LeftOn( G4TwoVector a,
00988 G4TwoVector b,
00989 G4TwoVector c )
00990 {
00991 return Area2(a,b,c)>=0;
00992 }
00993
00994
00995
00996
00997 G4bool G4PolyPhiFace::Collinear( G4TwoVector a,
00998 G4TwoVector b,
00999 G4TwoVector c )
01000 {
01001 return Area2(a,b,c)==0;
01002 }
01003
01004
01005
01006
01007
01008 G4bool G4PolyPhiFace::IntersectProp( G4TwoVector a,
01009 G4TwoVector b,
01010 G4TwoVector c, G4TwoVector d )
01011 {
01012 if( Collinear(a,b,c) || Collinear(a,b,d)||
01013 Collinear(c,d,a) || Collinear(c,d,b) ) { return false; }
01014
01015 G4bool Positive;
01016 Positive = !(Left(a,b,c))^!(Left(a,b,d));
01017 return Positive && (!Left(c,d,a)^!Left(c,d,b));
01018 }
01019
01020
01021
01022
01023
01024 G4bool G4PolyPhiFace::Between( G4TwoVector a, G4TwoVector b, G4TwoVector c )
01025 {
01026 if( !Collinear(a,b,c) ) { return false; }
01027
01028 if(a.x()!=b.x())
01029 {
01030 return ((a.x()<=c.x())&&(c.x()<=b.x()))||
01031 ((a.x()>=c.x())&&(c.x()>=b.x()));
01032 }
01033 else
01034 {
01035 return ((a.y()<=c.y())&&(c.y()<=b.y()))||
01036 ((a.y()>=c.y())&&(c.y()>=b.y()));
01037 }
01038 }
01039
01040
01041
01042
01043
01044 G4bool G4PolyPhiFace::Intersect( G4TwoVector a,
01045 G4TwoVector b,
01046 G4TwoVector c, G4TwoVector d )
01047 {
01048 if( IntersectProp(a,b,c,d) )
01049 { return true; }
01050 else if( Between(a,b,c)||
01051 Between(a,b,d)||
01052 Between(c,d,a)||
01053 Between(c,d,b) )
01054 { return true; }
01055 else
01056 { return false; }
01057 }
01058
01059
01060
01061
01062
01063 G4bool G4PolyPhiFace::Diagonalie( G4PolyPhiFaceVertex *a,
01064 G4PolyPhiFaceVertex *b )
01065 {
01066 G4PolyPhiFaceVertex *corner = triangles;
01067 G4PolyPhiFaceVertex *corner_next=triangles;
01068
01069
01070 do
01071 {
01072 corner_next=corner->next;
01073
01074
01075
01076 if( (corner!=a)&&(corner_next!=a)
01077 &&(corner!=b)&&(corner_next!=b) )
01078 {
01079 G4TwoVector rz1,rz2,rz3,rz4;
01080 rz1 = G4TwoVector(a->r,a->z);
01081 rz2 = G4TwoVector(b->r,b->z);
01082 rz3 = G4TwoVector(corner->r,corner->z);
01083 rz4 = G4TwoVector(corner_next->r,corner_next->z);
01084 if( Intersect(rz1,rz2,rz3,rz4) ) { return false; }
01085 }
01086 corner=corner->next;
01087
01088 } while( corner != triangles );
01089
01090 return true;
01091 }
01092
01093
01094
01095
01096
01097 G4bool G4PolyPhiFace::InCone( G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b )
01098 {
01099
01100
01101 G4PolyPhiFaceVertex *a0,*a1;
01102 a1=a->next;
01103 a0=a->prev;
01104
01105 G4TwoVector arz,arz0,arz1,brz;
01106 arz=G4TwoVector(a->r,a->z);arz0=G4TwoVector(a0->r,a0->z);
01107 arz1=G4TwoVector(a1->r,a1->z);brz=G4TwoVector(b->r,b->z);
01108
01109
01110 if(LeftOn(arz,arz1,arz0))
01111 {
01112 return Left(arz,brz,arz0)&&Left(brz,arz,arz1);
01113 }
01114 else
01115 {
01116 return !( LeftOn(arz,brz,arz1)&&LeftOn(brz,arz,arz0));
01117 }
01118 }
01119
01120
01121
01122
01123
01124 G4bool G4PolyPhiFace::Diagonal( G4PolyPhiFaceVertex *a, G4PolyPhiFaceVertex *b )
01125 {
01126 return InCone(a,b) && InCone(b,a) && Diagonalie(a,b);
01127 }
01128
01129
01130
01131
01132
01133 void G4PolyPhiFace::EarInit()
01134 {
01135 G4PolyPhiFaceVertex *corner = triangles;
01136 G4PolyPhiFaceVertex *c_prev,*c_next;
01137
01138 do
01139 {
01140
01141
01142 c_next=corner->next;
01143 c_prev=corner->prev;
01144
01145
01146
01147 corner->ear=Diagonal(c_prev,c_next);
01148 corner=corner->next;
01149
01150 } while( corner!=triangles );
01151 }
01152
01153
01154
01155
01156
01157 void G4PolyPhiFace::Triangulate()
01158 {
01159
01160
01161
01162 G4PolyPhiFaceVertex *tri_help = new G4PolyPhiFaceVertex[numEdges];
01163 triangles = tri_help;
01164 G4PolyPhiFaceVertex *triang = triangles;
01165
01166 std::vector<G4double> areas;
01167 std::vector<G4ThreeVector> points;
01168 G4double area=0.;
01169 G4PolyPhiFaceVertex *v0,*v1,*v2,*v3,*v4;
01170 v2=triangles;
01171
01172
01173
01174 G4PolyPhiFaceVertex *helper = corners;
01175 G4PolyPhiFaceVertex *helper2 = corners;
01176 do
01177 {
01178 triang->r = helper->r;
01179 triang->z = helper->z;
01180 triang->x = helper->x;
01181 triang->y= helper->y;
01182
01183
01184
01185 if( helper==corners )
01186 { triang->prev=triangles+numEdges-1; }
01187 else
01188 { triang->prev=helper2; }
01189
01190
01191
01192 if( helper<corners+numEdges-1 )
01193 { triang->next=triang+1; }
01194 else
01195 { triang->next=triangles; }
01196 helper2=triang;
01197 helper=helper->next;
01198 triang=triang->next;
01199
01200 } while( helper!=corners );
01201
01202 EarInit();
01203
01204 G4int n=numEdges;
01205 G4int i=0;
01206 G4ThreeVector p1,p2,p3,p4;
01207 const G4int max_n_loops=numEdges*10000;
01208
01209
01210
01211 while(n>3)
01212 {
01213 v2=triangles;
01214 do
01215 {
01216 if(v2->ear)
01217 {
01218
01219
01220 v3=v2->next; v4=v3->next;
01221 v1=v2->prev; v0=v1->prev;
01222
01223
01224
01225 p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
01226 p2=G4ThreeVector((v1)->x,(v1)->y,(v1)->z);
01227 p3=G4ThreeVector((v3)->x,(v3)->y,(v3)->z);
01228
01229 G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
01230 points.push_back(p4);
01231 areas.push_back(result1);
01232 area=area+result1;
01233
01234
01235
01236 v1->ear=Diagonal(v0,v3);
01237 v3->ear=Diagonal(v1,v4);
01238
01239
01240
01241
01242 v1->next=v3;
01243 v3->prev=v1;
01244 triangles=v3;
01245 n--;
01246
01247 break;
01248 }
01249
01250 v2=v2->next;
01251
01252 } while( v2!=triangles );
01253
01254 i++;
01255 if(i>=max_n_loops)
01256 {
01257 G4Exception( "G4PolyPhiFace::Triangulation()",
01258 "GeomSolids0003", FatalException,
01259 "Maximum number of steps is reached for triangulation!" );
01260 }
01261 }
01262
01263 if(v2->next)
01264 {
01265
01266
01267 v2=v2->next;
01268 p1=G4ThreeVector((v2)->x,(v2)->y,(v2)->z);
01269 p2=G4ThreeVector((v2->next)->x,(v2->next)->y,(v2->next)->z);
01270 p3=G4ThreeVector((v2->prev)->x,(v2->prev)->y,(v2->prev)->z);
01271 G4double result1 = SurfaceTriangle(p1,p2,p3,&p4 );
01272 points.push_back(p4);
01273 areas.push_back(result1);
01274 area=area+result1;
01275 }
01276
01277
01278
01279 fSurfaceArea = area;
01280
01281
01282
01283 G4double chose = area*G4UniformRand();
01284
01285
01286
01287 G4double Achose1, Achose2;
01288 Achose1=0; Achose2=0.;
01289 i=0;
01290 do
01291 {
01292 Achose2+=areas[i];
01293 if(chose>=Achose1 && chose<Achose2)
01294 {
01295 G4ThreeVector point;
01296 point=points[i] ;
01297 surface_point=point;
01298 break;
01299 }
01300 i++; Achose1=Achose2;
01301 } while( i<numEdges-2 );
01302
01303 delete [] tri_help;
01304 tri_help = 0;
01305 }