source: trunk/source/geometry/solids/specific/src/G4QuadrangularFacet.cc@ 1200

Last change on this file since 1200 was 1058, checked in by garnier, 17 years ago

file release beta
File size: 9.0 KB
RevLine 
[831]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 and of QinetiQ Ltd, *
[1058]20// * subject to DEFCON 705 IPR conditions. *
[831]21// * By using, copying, modifying or distributing the software (or *
22// * any work based on the software) you agree to acknowledge its *
23// * use in resulting scientific publications, and indicate your *
24// * acceptance of all terms of the Geant4 Software license. *
25// ********************************************************************
26//
27//
[1058]28// $Id: G4QuadrangularFacet.cc,v 1.7 2008/12/18 12:57:34 gunter Exp $
29// GEANT4 tag $Name: geant4-09-02-ref-02 $
30//
[831]31// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32//
33// MODULE: G4QuadrangularFacet.cc
34//
35// Date: 15/06/2005
36// Author: P R Truscott
37// Organisation: QinetiQ Ltd, UK
38// Customer: UK Ministry of Defence : RAO CRP TD Electronic Systems
39// Contract: C/MAT/N03517
40//
41// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42//
43// CHANGE HISTORY
44// --------------
45//
46// 31 October 2004, P R Truscott, QinetiQ Ltd, UK - Created.
47//
48// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49
50#include "G4QuadrangularFacet.hh"
51#include "globals.hh"
52#include "Randomize.hh"
53
54///////////////////////////////////////////////////////////////////////////////
55//
56// !!!THIS IS A FUDGE!!! IT'S TWO ADJACENT G4TRIANGULARFACETS
57// --- NOT EFFICIENT BUT PRACTICAL.
58//
59G4QuadrangularFacet::G4QuadrangularFacet (const G4ThreeVector Pt0,
60 const G4ThreeVector vt1, const G4ThreeVector vt2,
61 const G4ThreeVector vt3, G4FacetVertexType vertexType)
62 : G4VFacet()
63{
64 P0 = Pt0;
65 nVertices = 4;
66 if (vertexType == ABSOLUTE)
67 {
68 P.push_back(vt1);
69 P.push_back(vt2);
70 P.push_back(vt3);
71
72 E.push_back(vt1 - P0);
73 E.push_back(vt2 - P0);
74 E.push_back(vt3 - P0);
75 }
76 else
77 {
78 P.push_back(P0 + vt1);
79 P.push_back(P0 + vt2);
80 P.push_back(P0 + vt3);
81
82 E.push_back(vt1);
83 E.push_back(vt2);
84 E.push_back(vt3);
85 }
86
87 G4double length1 = E[0].mag();
88 G4double length2 = (P[1]-P[0]).mag();
89 G4double length3 = (P[2]-P[1]).mag();
90 G4double length4 = E[2].mag();
91
92 G4ThreeVector normal1 = E[0].cross(E[1]).unit();
93 G4ThreeVector normal2 = E[1].cross(E[2]).unit();
94
95 if (length1 <= kCarTolerance || length2 <= kCarTolerance ||
96 length3 <= kCarTolerance || length4 <= kCarTolerance ||
97 normal1.dot(normal2) < 0.9999999999)
98 {
99 G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
100 "InvalidSetup", JustWarning,
101 "Length of sides of facet are too small or sides not planar.");
102 G4cerr << G4endl;
103 G4cerr << "P0 = " << P0 << G4endl;
104 G4cerr << "P1 = " << P[0] << G4endl;
105 G4cerr << "P2 = " << P[1] << G4endl;
106 G4cerr << "P3 = " << P[2] << G4endl;
107 G4cerr << "Side lengths = P0->P1" << length1 << G4endl;
108 G4cerr << "Side lengths = P1->P2" << length2 << G4endl;
109 G4cerr << "Side lengths = P2->P3" << length3 << G4endl;
110 G4cerr << "Side lengths = P3->P0" << length4 << G4endl;
111 G4cerr << G4endl;
112
113 isDefined = false;
114 geometryType = "G4QuadragularFacet";
115 surfaceNormal = G4ThreeVector(0.0,0.0,0.0);
116 }
117 else
118 {
119 isDefined = true;
120 geometryType = "G4QuadrangularFacet";
121
122 facet1 = new G4TriangularFacet(P0,P[0],P[1],ABSOLUTE);
123 facet2 = new G4TriangularFacet(P0,P[1],P[2],ABSOLUTE);
124 surfaceNormal = normal1;
125
126 G4ThreeVector vtmp = 0.5 * (E[0] + E[1]);
127 circumcentre = P0 + vtmp;
128 radiusSqr = vtmp.mag2();
129 radius = std::sqrt(radiusSqr);
130
131 for (size_t i=0; i<4; i++) I.push_back(0);
132 }
133}
134
135///////////////////////////////////////////////////////////////////////////////
136//
137G4QuadrangularFacet::~G4QuadrangularFacet ()
138{
139 if (facet1 != 0) delete facet1;
140 if (facet2 != 0) delete facet2;
141
142 P.clear();
143 E.clear();
144 I.clear();
145}
146
147///////////////////////////////////////////////////////////////////////////////
148//
149G4VFacet *G4QuadrangularFacet::GetClone ()
150{
151 G4QuadrangularFacet *c =
152 new G4QuadrangularFacet (P0, P[0], P[1], P[2], ABSOLUTE);
153 G4VFacet *cc = 0;
154 cc = c;
155 return cc;
156}
157
158///////////////////////////////////////////////////////////////////////////////
159//
160G4ThreeVector G4QuadrangularFacet::Distance (const G4ThreeVector &p)
161{
162 G4ThreeVector v1 = facet1->Distance(p);
163 G4ThreeVector v2 = facet2->Distance(p);
164
165 if (v1.mag2() < v2.mag2()) return v1;
166 else return v2;
167}
168
169///////////////////////////////////////////////////////////////////////////////
170//
171G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p,
172 const G4double)
173{
174 /*G4ThreeVector D = P0 - p;
175 G4double d = E[0].dot(D);
176 G4double e = E[1].dot(D);
177 G4double s = b*e - c*d;
178 G4double t = b*d - a*e;*/
179 G4double dist = kInfinity;
180
181 /*if (s+t > 1.0 || s < 0.0 || t < 0.0)
182 {
183 G4ThreeVector D0 = P0 - p;
184 G4ThreeVector D1 = P[0] - p;
185 G4ThreeVector D2 = P[1] - p;
186
187 G4double d0 = D0.mag();
188 G4double d1 = D1.mag();
189 G4double d2 = D2.mag();
190
191 dist = min(d0, min(d1, d2));
192 if (dist > minDist) return kInfinity;
193 }*/
194
195 dist = Distance(p).mag();
196
197 return dist;
198}
199
200///////////////////////////////////////////////////////////////////////////////
201//
202G4double G4QuadrangularFacet::Distance (const G4ThreeVector &p,
203 const G4double, const G4bool outgoing)
204{
205 /*G4ThreeVector D = P0 - p;
206 G4double d = E[0].dot(D);
207 G4double e = E[1].dot(D);
208 G4double s = b*e - c*d;
209 G4double t = b*d - a*e;*/
210 G4double dist = kInfinity;
211
212 /*if (s+t > 1.0 || s < 0.0 || t < 0.0)
213 {
214 G4ThreeVector D0 = P0 - p;
215 G4ThreeVector D1 = P[0] - p;
216 G4ThreeVector D2 = P[1] - p;
217
218 G4double d0 = D0.mag();
219 G4double d1 = D1.mag();
220 G4double d2 = D2.mag();
221
222 dist = min(d0, min(d1, d2));
223 if (dist > minDist ||
224 (D0.dot(surfaceNormal) > 0.0 && !outgoing) ||
225 (D0.dot(surfaceNormal) < 0.0 && outgoing)) return kInfinity;
226 }*/
227
228 G4ThreeVector v = Distance(p);
229 G4double dir = v.dot(surfaceNormal);
230 if ((dir > dirTolerance && !outgoing) ||
231 (dir <-dirTolerance && outgoing)) dist = kInfinity;
232 else dist = v.mag();
233
234 return dist;
235}
236
237///////////////////////////////////////////////////////////////////////////////
238//
239G4double G4QuadrangularFacet::Extent (const G4ThreeVector axis)
240{
241 G4double s = P0.dot(axis);
242 for (G4ThreeVectorList::iterator it=P.begin(); it!=P.end(); it++)
243 {
244 G4double sp = it->dot(axis);
245 if (sp > s) s = sp;
246 }
247
248 return s;
249}
250
251///////////////////////////////////////////////////////////////////////////////
252//
253G4bool G4QuadrangularFacet::Intersect (const G4ThreeVector &p,
254 const G4ThreeVector &v, G4bool outgoing, G4double &distance,
255 G4double &distFromSurface, G4ThreeVector &normal)
256{
257 G4bool intersect =
258 facet1->Intersect(p,v,outgoing,distance,distFromSurface,normal);
259 if (!intersect)
260 {
261 intersect = facet2->Intersect(p,v,outgoing,distance,distFromSurface,normal);
262 }
263
264 if (!intersect)
265 {
266 distance = kInfinity;
267 distFromSurface = kInfinity;
268 normal = G4ThreeVector(0.0,0.0,0.0);
269 }
270
271 return intersect;
272}
273
274////////////////////////////////////////////////////////////////////////
275//
276// GetPointOnFace
277//
278// Auxiliary method for get a random point on surface
279
280G4ThreeVector G4QuadrangularFacet::GetPointOnFace() const
281{
282 G4ThreeVector pr;
283
284 if ( G4UniformRand() < 0.5 )
285 {
286 pr = facet1->GetPointOnFace();
287 }
288 else
289 {
290 pr = facet2->GetPointOnFace();
291 }
292
293 return pr;
294}
295
296////////////////////////////////////////////////////////////////////////
297//
298// GetArea
299//
300// Auxiliary method for returning the surface area
301
302G4double G4QuadrangularFacet::GetArea()
303{
304 if (!area) { area = facet1->GetArea() + facet2->GetArea(); }
305
306 return area;
307}
Note: See TracBrowser for help on using the repository browser.