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 #include "G4QuadrangularFacet.hh"
00044 #include "geomdefs.hh"
00045 #include "Randomize.hh"
00046
00047 using namespace std;
00048
00050
00051
00052
00053
00054 G4QuadrangularFacet::G4QuadrangularFacet (const G4ThreeVector &vt0,
00055 const G4ThreeVector &vt1,
00056 const G4ThreeVector &vt2,
00057 const G4ThreeVector &vt3,
00058 G4FacetVertexType vertexType)
00059 {
00060 G4ThreeVector e1, e2, e3;
00061
00062 SetVertex(0, vt0);
00063 if (vertexType == ABSOLUTE)
00064 {
00065 SetVertex(1, vt1);
00066 SetVertex(2, vt2);
00067 SetVertex(3, vt3);
00068
00069 e1 = vt1 - vt0;
00070 e2 = vt2 - vt0;
00071 e3 = vt3 - vt0;
00072 }
00073 else
00074 {
00075 SetVertex(1, vt0 + vt1);
00076 SetVertex(2, vt0 + vt2);
00077 SetVertex(3, vt0 + vt3);
00078
00079 e1 = vt1;
00080 e2 = vt2;
00081 e3 = vt3;
00082 }
00083 G4double length1 = e1.mag();
00084 G4double length2 = (GetVertex(2)-GetVertex(1)).mag();
00085 G4double length3 = (GetVertex(3)-GetVertex(2)).mag();
00086 G4double length4 = e3.mag();
00087
00088 G4ThreeVector normal1 = e1.cross(e2).unit();
00089 G4ThreeVector normal2 = e2.cross(e3).unit();
00090
00091 bool isDefined = (length1 > kCarTolerance && length2 > kCarTolerance &&
00092 length3 > kCarTolerance && length4 > kCarTolerance &&
00093 normal1.dot(normal2) >= 0.9999999999);
00094
00095 if (isDefined)
00096 {
00097 fFacet1 = G4TriangularFacet (GetVertex(0),GetVertex(1),
00098 GetVertex(2),ABSOLUTE);
00099 fFacet2 = G4TriangularFacet (GetVertex(0),GetVertex(2),
00100 GetVertex(3),ABSOLUTE);
00101
00102 G4TriangularFacet facet3 (GetVertex(0),GetVertex(1),GetVertex(3),ABSOLUTE);
00103 G4TriangularFacet facet4 (GetVertex(1),GetVertex(2),GetVertex(3),ABSOLUTE);
00104
00105 G4ThreeVector normal12 = fFacet1.GetSurfaceNormal()
00106 + fFacet2.GetSurfaceNormal();
00107 G4ThreeVector normal34 = facet3.GetSurfaceNormal()
00108 + facet4.GetSurfaceNormal();
00109 G4ThreeVector normal = 0.25 * (normal12 + normal34);
00110
00111 fFacet1.SetSurfaceNormal (normal);
00112 fFacet2.SetSurfaceNormal (normal);
00113
00114 G4ThreeVector vtmp = 0.5 * (e1 + e2);
00115 fCircumcentre = GetVertex(0) + vtmp;
00116 G4double radiusSqr = vtmp.mag2();
00117 fRadius = std::sqrt(radiusSqr);
00118 }
00119 else
00120 {
00121 G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
00122 "GeomSolids0002", JustWarning,
00123 "Length of sides of facet are too small or sides not planar.");
00124 G4cout << G4endl;
00125 G4cout << "P0 = " << GetVertex(0) << G4endl;
00126 G4cout << "P1 = " << GetVertex(1) << G4endl;
00127 G4cout << "P2 = " << GetVertex(2) << G4endl;
00128 G4cout << "P3 = " << GetVertex(3) << G4endl;
00129 G4cout << "Side lengths = P0->P1" << length1 << G4endl;
00130 G4cout << "Side lengths = P1->P2" << length2 << G4endl;
00131 G4cout << "Side lengths = P2->P3" << length3 << G4endl;
00132 G4cout << "Side lengths = P3->P0" << length4 << G4endl;
00133 G4cout << G4endl;
00134 fRadius = 0.0;
00135 }
00136 }
00137
00139
00140 G4QuadrangularFacet::~G4QuadrangularFacet ()
00141 {
00142 }
00143
00145
00146 G4QuadrangularFacet::G4QuadrangularFacet (const G4QuadrangularFacet &rhs)
00147 : G4VFacet(rhs)
00148 {
00149 fFacet1 = rhs.fFacet1;
00150 fFacet2 = rhs.fFacet2;
00151 fRadius = 0.0;
00152 }
00153
00155
00156 G4QuadrangularFacet &
00157 G4QuadrangularFacet::operator=(const G4QuadrangularFacet &rhs)
00158 {
00159 if (this == &rhs)
00160 return *this;
00161
00162 fFacet1 = rhs.fFacet1;
00163 fFacet2 = rhs.fFacet2;
00164 fRadius = 0.0;
00165
00166 return *this;
00167 }
00168
00170
00171 G4VFacet *G4QuadrangularFacet::GetClone ()
00172 {
00173 G4QuadrangularFacet *c = new G4QuadrangularFacet (GetVertex(0), GetVertex(1),
00174 GetVertex(2), GetVertex(3),
00175 ABSOLUTE);
00176 return c;
00177 }
00178
00180
00181 G4ThreeVector G4QuadrangularFacet::Distance (const G4ThreeVector &p)
00182 {
00183 G4ThreeVector v1 = fFacet1.Distance(p);
00184 G4ThreeVector v2 = fFacet2.Distance(p);
00185
00186 if (v1.mag2() < v2.mag2()) return v1;
00187 else return v2;
00188 }
00189
00191
00192 G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p,
00193 G4double)
00194 {
00195 G4double dist = Distance(p).mag();
00196 return dist;
00197 }
00198
00200
00201 G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p, G4double,
00202 const G4bool outgoing)
00203 {
00204 G4double dist;
00205
00206 G4ThreeVector v = Distance(p);
00207 G4double dir = v.dot(GetSurfaceNormal());
00208 if ( ((dir > dirTolerance) && (!outgoing))
00209 || ((dir < -dirTolerance) && outgoing))
00210 dist = kInfinity;
00211 else
00212 dist = v.mag();
00213 return dist;
00214 }
00215
00217
00218 G4double G4QuadrangularFacet::Extent (const G4ThreeVector axis)
00219 {
00220 G4double ss = 0;
00221
00222 for (G4int i = 0; i <= 3; ++i)
00223 {
00224 G4double sp = GetVertex(i).dot(axis);
00225 if (sp > ss) ss = sp;
00226 }
00227 return ss;
00228 }
00229
00231
00232 G4bool G4QuadrangularFacet::Intersect (const G4ThreeVector &p,
00233 const G4ThreeVector &v,
00234 G4bool outgoing,
00235 G4double &distance,
00236 G4double &distFromSurface,
00237 G4ThreeVector &normal)
00238 {
00239 G4bool intersect =
00240 fFacet1.Intersect(p,v,outgoing,distance,distFromSurface,normal);
00241 if (!intersect) intersect =
00242 fFacet2.Intersect(p,v,outgoing,distance,distFromSurface,normal);
00243 if (!intersect)
00244 {
00245 distance = distFromSurface = kInfinity;
00246 normal.set(0,0,0);
00247 }
00248 return intersect;
00249 }
00250
00252
00253
00254
00255 G4ThreeVector G4QuadrangularFacet::GetPointOnFace() const
00256 {
00257 G4ThreeVector pr = (G4RandFlat::shoot(0.,1.) < 0.5)
00258 ? fFacet1.GetPointOnFace() : fFacet2.GetPointOnFace();
00259 return pr;
00260 }
00261
00263
00264
00265
00266 G4double G4QuadrangularFacet::GetArea()
00267 {
00268 G4double area = fFacet1.GetArea() + fFacet2.GetArea();
00269 return area;
00270 }
00271
00273
00274 G4String G4QuadrangularFacet::GetEntityType () const
00275 {
00276 return "G4QuadrangularFacet";
00277 }
00278
00280
00281 G4ThreeVector G4QuadrangularFacet::GetSurfaceNormal () const
00282 {
00283 return fFacet1.GetSurfaceNormal();
00284 }