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

Last change on this file since 1042 was 850, checked in by garnier, 17 years ago

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