Geant4  10.01.p02
UQuadrangularFacet.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * This Software is part of the AIDA Unified Solids Library package *
4 // * See: https://aidasoft.web.cern.ch/USolids *
5 // ********************************************************************
6 //
7 // $Id:$
8 //
9 // --------------------------------------------------------------------
10 //
11 // UQuadrangularFacet
12 //
13 // 17.10.12 Marek Gayer
14 // Created from original implementation in Geant4
15 // --------------------------------------------------------------------
16 
17 #include "UUtils.hh"
18 #include "VUSolid.hh"
19 #include "UQuadrangularFacet.hh"
20 
21 using namespace std;
22 
24 //
25 // !!!THIS IS A FUDGE!!! IT'S TWO ADJACENT UTRIANGULARFACETS
26 // --- NOT EFFICIENT BUT PRACTICAL.
27 //
29  const UVector3& vt1, const UVector3& vt2,
30  const UVector3& vt3, UFacetVertexType vertexType)
31 {
32  UVector3 e1, e2, e3;
33 
34  SetVertex(0, vt0);
35  if (vertexType == UABSOLUTE)
36  {
37  SetVertex(1, vt1);
38  SetVertex(2, vt2);
39  SetVertex(3, vt3);
40 
41  e1 = vt1 - vt0;
42  e2 = vt2 - vt0;
43  e3 = vt3 - vt0;
44  }
45  else
46  {
47  SetVertex(1, vt0 + vt1);
48  SetVertex(2, vt0 + vt2);
49  SetVertex(3, vt0 + vt3);
50 
51  e1 = vt1;
52  e2 = vt2;
53  e3 = vt3;
54  }
55  double length1 = e1.Mag();
56  double length2 = (GetVertex(2) - GetVertex(1)).Mag();
57  double length3 = (GetVertex(3) - GetVertex(2)).Mag();
58  double length4 = e3.Mag();
59 
60  UVector3 normal1 = e1.Cross(e2).Unit();
61  UVector3 normal2 = e2.Cross(e3).Unit();
62 
63  bool isDefined = (length1 > VUSolid::Tolerance() && length2 > VUSolid::Tolerance() &&
64  length3 > VUSolid::Tolerance() && length4 > VUSolid::Tolerance() &&
65  normal1.Dot(normal2) >= 0.9999999999);
66 
67  if (isDefined)
68  {
69  fFacet1 = UTriangularFacet(GetVertex(0), GetVertex(1), GetVertex(2), UABSOLUTE);
70  fFacet2 = UTriangularFacet(GetVertex(0), GetVertex(2), GetVertex(3), UABSOLUTE);
71 
72  UTriangularFacet facet3(GetVertex(0), GetVertex(1), GetVertex(3), UABSOLUTE);
73  UTriangularFacet facet4(GetVertex(1), GetVertex(2), GetVertex(3), UABSOLUTE);
74 
75  UVector3 normal12 = fFacet1.GetSurfaceNormal() + fFacet2.GetSurfaceNormal();
76  UVector3 normal34 = facet3.GetSurfaceNormal() + facet4.GetSurfaceNormal();
77  UVector3 normal = 0.25 * (normal12 + normal34);
78 
79  fFacet1.SetSurfaceNormal(normal);
80  fFacet2.SetSurfaceNormal(normal);
81 
82  UVector3 vtmp = 0.5 * (e1 + e2);
83  fCircumcentre = GetVertex(0) + vtmp;
84  double radiusSqr = vtmp.Mag2();
85  fRadius = std::sqrt(radiusSqr);
86  }
87  else
88  {
89  UUtils::Exception("UQuadrangularFacet::UQuadrangularFacet()", "GeomSolids1002",
90  UWarning, 1, "Length of sides of facet are too small or sides not planar.");
91  cerr << endl;
92  cerr << "P0 = " << GetVertex(0) << endl;
93  cerr << "P1 = " << GetVertex(1) << endl;
94  cerr << "P2 = " << GetVertex(2) << endl;
95  cerr << "P3 = " << GetVertex(3) << endl;
96  cerr << "Side lengths = P0->P1" << length1 << endl;
97  cerr << "Side lengths = P1->P2" << length2 << endl;
98  cerr << "Side lengths = P2->P3" << length3 << endl;
99  cerr << "Side lengths = P3->P0" << length4 << endl;
100  cerr << endl;
101  fRadius = 0;
102  }
103 }
104 
106 //
108 {
109 }
110 
112 //
114 {
115  fFacet1 = rhs.fFacet1;
116  fFacet2 = rhs.fFacet2;
117  fRadius = 0;
118 }
119 
121 //
123 {
124  if (this == &rhs)
125  return *this;
126 
127  fFacet1 = rhs.fFacet1;
128  fFacet2 = rhs.fFacet2;
129 
130  fRadius = 0;
131 
132  return *this;
133 }
134 
136 //
138 {
140  return c;
141 }
142 
144 //
146 {
147  UVector3 v1 = fFacet1.Distance(p);
148  UVector3 v2 = fFacet2.Distance(p);
149 
150  if (v1.Mag2() < v2.Mag2()) return v1;
151  else return v2;
152 }
153 
155 //
157  const double)
158 {
159  double dist = Distance(p).Mag();
160  return dist;
161 }
162 
164 //
165 double UQuadrangularFacet::Distance(const UVector3& p, const double, const bool outgoing)
166 {
167  double dist;
168 
169  UVector3 v = Distance(p);
170  double dir = v.Dot(GetSurfaceNormal());
171  if ((dir > dirTolerance && !outgoing) || (dir < -dirTolerance && outgoing))
172  dist = UUtils::kInfinity;
173  else
174  dist = v.Mag();
175  return dist;
176 }
177 
179 {
180  double ss = 0;
181 
182  for (int i = 0; i <= 3; ++i)
183  {
184  double sp = GetVertex(i).Dot(axis);
185  if (sp > ss) ss = sp;
186  }
187  return ss;
188 }
189 
190 bool UQuadrangularFacet::Intersect(const UVector3& p, const UVector3& v, bool outgoing, double& distance, double& distFromSurface, UVector3& normal)
191 {
192  bool intersect = fFacet1.Intersect(p, v, outgoing, distance, distFromSurface, normal);
193  if (!intersect) intersect = fFacet2.Intersect(p, v, outgoing, distance, distFromSurface, normal);
194  if (!intersect)
195  {
196  distance = distFromSurface = UUtils::kInfinity;
197  normal.Set(0);
198  }
199  return intersect;
200 }
201 
202 // Auxiliary method for get a random point on surface
204 {
206  return pr;
207 }
208 
209 // Auxiliary method for returning the surface area
211 {
212  double area = fFacet1.GetArea() + fFacet2.GetArea();
213  return area;
214 }
215 
217 {
218  return "QuadrangularFacet";
219 }
220 
222 {
223  return fFacet1.GetSurfaceNormal();
224 }
double Mag2() const
Definition: UVector3.hh:287
UTriangularFacet fFacet1
virtual UGeometryType GetEntityType() const
static const double dirTolerance
Definition: VUFacet.hh:88
UVector3 GetPointOnFace() const
UVector3 Distance(const UVector3 &p)
UVector3 GetSurfaceNormal() const
UVector3 Cross(const UVector3 &) const
Definition: UVector3.hh:281
static const G4double e2
UQuadrangularFacet & operator=(const UQuadrangularFacet &right)
UFacetVertexType
Definition: VUFacet.hh:49
static double Tolerance()
Definition: VUSolid.hh:127
static double normal(HepRandomEngine *eptr)
Definition: RandPoisson.cc:77
UQuadrangularFacet(const UVector3 &Pt0, const UVector3 &vt1, const UVector3 &vt2, const UVector3 &vt3, UFacetVertexType)
static const double kInfinity
Definition: UUtils.hh:54
double Extent(const UVector3 axis)
double Dot(const UVector3 &) const
Definition: UVector3.hh:276
double Mag() const
Definition: UVector3.cc:48
static const G4double e1
void Exception(const char *originOfException, const char *exceptionCode, UExceptionSeverity severity, int level, const char *description)
Definition: UUtils.cc:122
void Set(double xx, double yy, double zz)
Definition: UVector3.hh:264
UVector3 Unit() const
Definition: UVector3.cc:80
UTriangularFacet fFacet2
UVector3 GetPointOnFace() const
UVector3 GetVertex(int i) const
bool Intersect(const UVector3 &p, const UVector3 &v, const bool outgoing, double &distance, double &distFromSurface, UVector3 &normal)
static const G4double e3
UVector3 GetSurfaceNormal() const
double Random(double min=0.0, double max=1.0)
Definition: UUtils.cc:69
bool Intersect(const UVector3 &p, const UVector3 &v, const bool outgoing, double &distance, double &distFromSurface, UVector3 &normal)
UVector3 Distance(const UVector3 &p)