Geant4  10.02
G4QuadrangularFacet.cc
Go to the documentation of this file.
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 to 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 //
28 // $Id: G4QuadrangularFacet.cc 66819 2013-01-12 16:20:10Z gcosmo $
29 //
30 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31 //
32 // CHANGE HISTORY
33 // --------------
34 //
35 // 31 October 2004, P R Truscott, QinetiQ Ltd, UK - Created.
36 // 12 October 2012, M Gayer, CERN
37 // New implementation reducing memory requirements by 50%,
38 // and considerable CPU speedup together with the new
39 // implementation of G4TessellatedSolid.
40 //
41 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42 
43 #include "G4QuadrangularFacet.hh"
44 #include "geomdefs.hh"
45 #include "Randomize.hh"
46 
47 using namespace std;
48 
50 //
51 // !!!THIS IS A FUDGE!!! IT'S TWO ADJACENT G4TRIANGULARFACETS
52 // --- NOT EFFICIENT BUT PRACTICAL.
53 //
55  const G4ThreeVector &vt1,
56  const G4ThreeVector &vt2,
57  const G4ThreeVector &vt3,
58  G4FacetVertexType vertexType)
59 {
61 
62  SetVertex(0, vt0);
63  if (vertexType == ABSOLUTE)
64  {
65  SetVertex(1, vt1);
66  SetVertex(2, vt2);
67  SetVertex(3, vt3);
68 
69  e1 = vt1 - vt0;
70  e2 = vt2 - vt0;
71  e3 = vt3 - vt0;
72  }
73  else
74  {
75  SetVertex(1, vt0 + vt1);
76  SetVertex(2, vt0 + vt2);
77  SetVertex(3, vt0 + vt3);
78 
79  e1 = vt1;
80  e2 = vt2;
81  e3 = vt3;
82  }
83  G4double length1 = e1.mag();
84  G4double length2 = (GetVertex(2)-GetVertex(1)).mag();
85  G4double length3 = (GetVertex(3)-GetVertex(2)).mag();
86  G4double length4 = e3.mag();
87 
88  G4ThreeVector normal1 = e1.cross(e2).unit();
89  G4ThreeVector normal2 = e2.cross(e3).unit();
90 
91  bool isDefined = (length1 > kCarTolerance && length2 > kCarTolerance &&
92  length3 > kCarTolerance && length4 > kCarTolerance &&
93  normal1.dot(normal2) >= 0.9999999999);
94 
95  if (isDefined)
96  {
97  fFacet1 = G4TriangularFacet (GetVertex(0),GetVertex(1),
98  GetVertex(2),ABSOLUTE);
99  fFacet2 = G4TriangularFacet (GetVertex(0),GetVertex(2),
100  GetVertex(3),ABSOLUTE);
101 
102  G4TriangularFacet facet3 (GetVertex(0),GetVertex(1),GetVertex(3),ABSOLUTE);
103  G4TriangularFacet facet4 (GetVertex(1),GetVertex(2),GetVertex(3),ABSOLUTE);
104 
105  G4ThreeVector normal12 = fFacet1.GetSurfaceNormal()
106  + fFacet2.GetSurfaceNormal();
107  G4ThreeVector normal34 = facet3.GetSurfaceNormal()
108  + facet4.GetSurfaceNormal();
109  G4ThreeVector normal = 0.25 * (normal12 + normal34);
110 
111  fFacet1.SetSurfaceNormal (normal);
112  fFacet2.SetSurfaceNormal (normal);
113 
114  G4ThreeVector vtmp = 0.5 * (e1 + e2);
115  fCircumcentre = GetVertex(0) + vtmp;
116  G4double radiusSqr = vtmp.mag2();
117  fRadius = std::sqrt(radiusSqr);
118  }
119  else
120  {
121  G4Exception("G4QuadrangularFacet::G4QuadrangularFacet()",
122  "GeomSolids0002", JustWarning,
123  "Length of sides of facet are too small or sides not planar.");
124  G4cout << G4endl;
125  G4cout << "P0 = " << GetVertex(0) << G4endl;
126  G4cout << "P1 = " << GetVertex(1) << G4endl;
127  G4cout << "P2 = " << GetVertex(2) << G4endl;
128  G4cout << "P3 = " << GetVertex(3) << G4endl;
129  G4cout << "Side lengths = P0->P1" << length1 << G4endl;
130  G4cout << "Side lengths = P1->P2" << length2 << G4endl;
131  G4cout << "Side lengths = P2->P3" << length3 << G4endl;
132  G4cout << "Side lengths = P3->P0" << length4 << G4endl;
133  G4cout << G4endl;
134  fRadius = 0.0;
135  }
136 }
137 
139 //
141 {
142 }
143 
145 //
147  : G4VFacet(rhs)
148 {
149  fFacet1 = rhs.fFacet1;
150  fFacet2 = rhs.fFacet2;
151  fRadius = 0.0;
152 }
153 
155 //
158 {
159  if (this == &rhs)
160  return *this;
161 
162  fFacet1 = rhs.fFacet1;
163  fFacet2 = rhs.fFacet2;
164  fRadius = 0.0;
165 
166  return *this;
167 }
168 
170 //
172 {
174  GetVertex(2), GetVertex(3),
175  ABSOLUTE);
176  return c;
177 }
178 
180 //
182 {
185 
186  if (v1.mag2() < v2.mag2()) return v1;
187  else return v2;
188 }
189 
191 //
193  G4double)
194 {
195  G4double dist = Distance(p).mag();
196  return dist;
197 }
198 
200 //
202  const G4bool outgoing)
203 {
204  G4double dist;
205 
206  G4ThreeVector v = Distance(p);
207  G4double dir = v.dot(GetSurfaceNormal());
208  if ( ((dir > dirTolerance) && (!outgoing))
209  || ((dir < -dirTolerance) && outgoing))
210  dist = kInfinity;
211  else
212  dist = v.mag();
213  return dist;
214 }
215 
217 //
219 {
220  G4double ss = 0;
221 
222  for (G4int i = 0; i <= 3; ++i)
223  {
224  G4double sp = GetVertex(i).dot(axis);
225  if (sp > ss) ss = sp;
226  }
227  return ss;
228 }
229 
231 //
233  const G4ThreeVector &v,
234  G4bool outgoing,
235  G4double &distance,
236  G4double &distFromSurface,
238 {
239  G4bool intersect =
240  fFacet1.Intersect(p,v,outgoing,distance,distFromSurface,normal);
241  if (!intersect) intersect =
242  fFacet2.Intersect(p,v,outgoing,distance,distFromSurface,normal);
243  if (!intersect)
244  {
245  distance = distFromSurface = kInfinity;
246  normal.set(0,0,0);
247  }
248  return intersect;
249 }
250 
252 //
253 // Auxiliary method to get a random point on surface
254 //
256 {
257  G4ThreeVector pr = (G4RandFlat::shoot(0.,1.) < 0.5)
259  return pr;
260 }
261 
263 //
264 // Auxiliary method for returning the surface area
265 //
267 {
268  G4double area = fFacet1.GetArea() + fFacet2.GetArea();
269  return area;
270 }
271 
273 //
275 {
276  return "G4QuadrangularFacet";
277 }
278 
280 //
282 {
283  return fFacet1.GetSurfaceNormal();
284 }
ThreeVector shoot(const G4int Ap, const G4int Af)
G4TriangularFacet fFacet1
G4bool Intersect(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool outgoing, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal)
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
static const G4double e2
G4double Extent(const G4ThreeVector axis)
G4ThreeVector GetSurfaceNormal() const
int G4int
Definition: G4Types.hh:78
G4bool Intersect(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool outgoing, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal)
G4QuadrangularFacet(const G4ThreeVector &Pt0, const G4ThreeVector &vt1, const G4ThreeVector &vt2, const G4ThreeVector &vt3, G4FacetVertexType)
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
G4GLOB_DLL std::ostream G4cout
G4FacetVertexType
Definition: G4VFacet.hh:56
G4ThreeVector GetVertex(G4int i) const
static const G4double dirTolerance
Definition: G4VFacet.hh:101
bool G4bool
Definition: G4Types.hh:79
G4ThreeVector GetPointOnFace() const
G4QuadrangularFacet & operator=(const G4QuadrangularFacet &right)
static const G4double e1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4TriangularFacet fFacet2
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4ThreeVector GetPointOnFace() const
static const G4double e3
G4ThreeVector GetSurfaceNormal() const
G4ThreeVector Distance(const G4ThreeVector &p)
G4GeometryType GetEntityType() const
G4ThreeVector Distance(const G4ThreeVector &p)