| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration and of QinetiQ Ltd, *
20 // * subject to DEFCON 705 IPR conditions. 20 // * subject to DEFCON 705 IPR conditions. *
21 // * By using, copying, modifying or distri 21 // * By using, copying, modifying or distributing the software (or *
22 // * any work based on the software) you ag 22 // * any work based on the software) you agree to acknowledge its *
23 // * use in resulting scientific publicati 23 // * use in resulting scientific publications, and indicate your *
24 // * acceptance of all terms of the Geant4 Sof 24 // * acceptance of all terms of the Geant4 Software license. *
25 // ******************************************* 25 // ********************************************************************
26 // 26 //
27 // G4TessellatedSolid implementation << 27 // $Id: G4TessellatedSolid.cc,v 1.19 2009/04/27 08:06:27 gcosmo Exp $
>> 28 // GEANT4 tag $Name: geant4-09-03-patch-01 $
28 // 29 //
29 // 31.10.2004, P R Truscott, QinetiQ Ltd, UK - << 30 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 // 17.09.2007, P R Truscott, QinetiQ Ltd & Ric << 31 //
>> 32 // MODULE: G4TessellatedSolid.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 // 14 November 2007 P R Truscott, QinetiQ & Stan Seibert, U Texas
>> 46 // Bug fixes to CalculateExtent
>> 47 //
>> 48 // 17 September 2007, P R Truscott, QinetiQ Ltd & Richard Holmberg
31 // Updated extensively prio 49 // Updated extensively prior to this date to deal with
32 // concaved tessellated sur 50 // concaved tessellated surfaces, based on the algorithm
33 // of Richard Holmberg. Th 51 // of Richard Holmberg. This had been slightly modified
34 // to determine with inside 52 // to determine with inside the geometry by projecting
35 // random rays from the poi 53 // random rays from the point provided. Now random rays
36 // are predefined rather th 54 // are predefined rather than making use of random
37 // number generator at run- 55 // number generator at run-time.
38 // 12.10.2012, M Gayer, CERN, complete rewrite << 56 //
39 // requirements more than 5 << 57 // 22 November 2005, F Lei
40 // tens or more depending o << 58 // - Changed ::DescribeYourselfTo(), line 464
41 // to voxelization of surfa << 59 // - added GetPolyHedron()
42 // Speedup factor of thousa << 60 //
43 // facets in hundreds of th << 61 // 31 October 2004, P R Truscott, QinetiQ Ltd, UK
44 // 23.10.2016, E Tcherniaev, reimplemented Cal << 62 // - Created.
45 // use of G4BoundingEnvelop << 63 //
46 // ------------------------------------------- << 64 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47 65
48 #include "G4TessellatedSolid.hh" 66 #include "G4TessellatedSolid.hh"
49 << 67 #include "G4PolyhedronArbitrary.hh"
50 #if !defined(G4GEOM_USE_UTESSELLATEDSOLID) << 68 #include "globals.hh"
51 <<
52 #include <algorithm> <<
53 #include <fstream> <<
54 #include <iomanip> <<
55 #include <iostream> <<
56 #include <list> <<
57 #include <random> <<
58 #include <stack> <<
59 <<
60 #include "geomdefs.hh" <<
61 #include "Randomize.hh" 69 #include "Randomize.hh"
62 #include "G4SystemOfUnits.hh" <<
63 #include "G4PhysicalConstants.hh" <<
64 #include "G4GeometryTolerance.hh" <<
65 #include "G4VoxelLimits.hh" <<
66 #include "G4AffineTransform.hh" <<
67 #include "G4BoundingEnvelope.hh" <<
68 70
69 #include "G4VGraphicsScene.hh" << 71 #include <iostream>
70 #include "G4VisExtent.hh" <<
71 <<
72 #include "G4AutoLock.hh" <<
73 <<
74 namespace <<
75 { <<
76 G4Mutex polyhedronMutex = G4MUTEX_INITIALIZE <<
77 } <<
78 <<
79 using namespace std; <<
80 72
81 ////////////////////////////////////////////// 73 ///////////////////////////////////////////////////////////////////////////////
82 // 74 //
83 // Standard contructor has blank name and defi << 75 // Standard contructor has blank name and defines no facets.
84 // 76 //
85 G4TessellatedSolid::G4TessellatedSolid () : G4 << 77 G4TessellatedSolid::G4TessellatedSolid ()
>> 78 : G4VSolid("dummy"), fpPolyhedron(0), cubicVolume(0.), surfaceArea(0.)
86 { 79 {
87 Initialize(); << 80 dirTolerance = 1.0E-14;
>> 81
>> 82 geometryType = "G4TessellatedSolid";
>> 83 facets.clear();
>> 84 solidClosed = false;
>> 85
>> 86 xMinExtent = kInfinity;
>> 87 xMaxExtent = -kInfinity;
>> 88 yMinExtent = kInfinity;
>> 89 yMaxExtent = -kInfinity;
>> 90 zMinExtent = kInfinity;
>> 91 zMaxExtent = -kInfinity;
>> 92
>> 93 SetRandomVectorSet();
88 } 94 }
89 95
90 ////////////////////////////////////////////// 96 ///////////////////////////////////////////////////////////////////////////////
91 // 97 //
92 // Alternative constructor. Simple define name << 98 // Alternative constructor. Simple define name and geometry type - no facets
93 // to detine. 99 // to detine.
94 // 100 //
95 G4TessellatedSolid::G4TessellatedSolid (const << 101 G4TessellatedSolid::G4TessellatedSolid (const G4String &name)
96 : G4VSolid(name) << 102 : G4VSolid(name), fpPolyhedron(0), cubicVolume(0.), surfaceArea(0.)
97 { 103 {
98 Initialize(); << 104 dirTolerance = 1.0E-14;
>> 105
>> 106 geometryType = "G4TessellatedSolid";
>> 107 facets.clear();
>> 108 solidClosed = false;
>> 109
>> 110 xMinExtent = kInfinity;
>> 111 xMaxExtent = -kInfinity;
>> 112 yMinExtent = kInfinity;
>> 113 yMaxExtent = -kInfinity;
>> 114 zMinExtent = kInfinity;
>> 115 zMaxExtent = -kInfinity;
>> 116
>> 117 SetRandomVectorSet();
99 } 118 }
100 119
101 ////////////////////////////////////////////// 120 ///////////////////////////////////////////////////////////////////////////////
102 // 121 //
103 // Fake default constructor - sets only member 122 // Fake default constructor - sets only member data and allocates memory
104 // for usage restri 123 // for usage restricted to object persistency.
105 // 124 //
106 G4TessellatedSolid::G4TessellatedSolid( __void << 125 G4TessellatedSolid::G4TessellatedSolid( __void__& a )
>> 126 : G4VSolid(a), fpPolyhedron(0), facets(0),
>> 127 geometryType("G4TessellatedSolid"), cubicVolume(0.), surfaceArea(0.),
>> 128 vertexList(), xMinExtent(0.), xMaxExtent(0.),
>> 129 yMinExtent(0.), yMaxExtent(0.), zMinExtent(0.), zMaxExtent(0.),
>> 130 solidClosed(false)
107 { 131 {
108 Initialize(); << 132 SetRandomVectorSet();
109 fMinExtent.set(0,0,0); <<
110 fMaxExtent.set(0,0,0); <<
111 } <<
112 <<
113 <<
114 ////////////////////////////////////////////// <<
115 G4TessellatedSolid::~G4TessellatedSolid() <<
116 { <<
117 DeleteObjects(); <<
118 } 133 }
119 134
120 ////////////////////////////////////////////// 135 ///////////////////////////////////////////////////////////////////////////////
121 // 136 //
122 // Copy constructor. << 137 // Destructor.
123 // 138 //
124 G4TessellatedSolid::G4TessellatedSolid (const << 139 G4TessellatedSolid::~G4TessellatedSolid ()
125 : G4VSolid(ts) <<
126 { 140 {
127 Initialize(); << 141 DeleteObjects ();
128 <<
129 CopyObjects(ts); <<
130 } 142 }
131 143
132 ////////////////////////////////////////////// 144 ///////////////////////////////////////////////////////////////////////////////
133 // 145 //
134 // Assignment operator. << 146 // Define copy constructor.
135 // 147 //
136 G4TessellatedSolid& << 148 G4TessellatedSolid::G4TessellatedSolid (const G4TessellatedSolid &s)
137 G4TessellatedSolid::operator= (const G4Tessell << 149 : G4VSolid(s)
138 { 150 {
139 if (&ts == this) return *this; << 151 if (&s == this) { return; }
140 <<
141 // Copy base class data <<
142 G4VSolid::operator=(ts); <<
143 152
144 DeleteObjects (); << 153 dirTolerance = 1.0E-14;
>> 154
>> 155 geometryType = "G4TessellatedSolid";
>> 156 facets.clear();
>> 157 solidClosed = false;
>> 158
>> 159 xMinExtent = kInfinity;
>> 160 xMaxExtent = -kInfinity;
>> 161 yMinExtent = kInfinity;
>> 162 yMaxExtent = -kInfinity;
>> 163 zMinExtent = kInfinity;
>> 164 zMaxExtent = -kInfinity;
145 165
146 Initialize(); << 166 SetRandomVectorSet();
147 167
148 CopyObjects (ts); << 168 CopyObjects (s);
149 <<
150 return *this; <<
151 } 169 }
152 170
153 ////////////////////////////////////////////// 171 ///////////////////////////////////////////////////////////////////////////////
154 // 172 //
155 void G4TessellatedSolid::Initialize() << 173 // Define assignment operator.
>> 174 //
>> 175 const G4TessellatedSolid &
>> 176 G4TessellatedSolid::operator= (const G4TessellatedSolid &s)
156 { 177 {
157 kCarToleranceHalf = 0.5*kCarTolerance; << 178 if (&s == this) { return *this; }
158 << 179
159 fRebuildPolyhedron = false; fpPolyhedron = n << 180 DeleteObjects ();
160 fCubicVolume = 0.; fSurfaceArea = 0.; << 181 CopyObjects (s);
161 << 182
162 fGeometryType = "G4TessellatedSolid"; << 183 return *this;
163 fSolidClosed = false; <<
164 <<
165 fMinExtent.set(kInfinity,kInfinity,kInfinity <<
166 fMaxExtent.set(-kInfinity,-kInfinity,-kInfin <<
167 <<
168 SetRandomVectors(); <<
169 } 184 }
170 185
171 ////////////////////////////////////////////// 186 ///////////////////////////////////////////////////////////////////////////////
172 // 187 //
173 void G4TessellatedSolid::DeleteObjects() << 188 void G4TessellatedSolid::DeleteObjects ()
174 { 189 {
175 std::size_t size = fFacets.size(); << 190 for (std::vector<G4VFacet *>::iterator f=facets.begin(); f!=facets.end(); f++)
176 for (std::size_t i = 0; i < size; ++i) { de << 191 {
177 fFacets.clear(); << 192 delete *f;
178 delete fpPolyhedron; fpPolyhedron = nullptr; << 193 }
>> 194 facets.clear();
179 } 195 }
180 196
181 ////////////////////////////////////////////// 197 ///////////////////////////////////////////////////////////////////////////////
182 // 198 //
183 void G4TessellatedSolid::CopyObjects (const G4 << 199 void G4TessellatedSolid::CopyObjects (const G4TessellatedSolid &s)
184 { 200 {
185 G4ThreeVector reductionRatio; << 201 size_t n = s.GetNumberOfFacets();
186 G4int fmaxVoxels = fVoxels.GetMaxVoxels(redu << 202 for (size_t i=0; i<n; i++)
187 if (fmaxVoxels < 0) <<
188 fVoxels.SetMaxVoxels(reductionRatio); <<
189 else <<
190 fVoxels.SetMaxVoxels(fmaxVoxels); <<
191 <<
192 G4int n = ts.GetNumberOfFacets(); <<
193 for (G4int i = 0; i < n; ++i) <<
194 { 203 {
195 G4VFacet *facetClone = (ts.GetFacet(i))->G << 204 G4VFacet *facetClone = (s.GetFacet(i))->GetClone();
196 AddFacet(facetClone); 205 AddFacet(facetClone);
197 } 206 }
198 if (ts.GetSolidClosed()) SetSolidClosed(true << 207
>> 208 if ( s.GetSolidClosed() ) { SetSolidClosed(true); }
>> 209
>> 210 // cubicVolume = s.GetCubicVolume();
199 } 211 }
200 212
201 ////////////////////////////////////////////// 213 ///////////////////////////////////////////////////////////////////////////////
202 // 214 //
203 // Add a facet to the facet list. << 215 // Add a facet to the facet list. Note that you can add, but you cannot
204 // Note that you can add, but you cannot delet << 216 // delete.
205 // 217 //
206 G4bool G4TessellatedSolid::AddFacet (G4VFacet* << 218 G4bool G4TessellatedSolid::AddFacet (G4VFacet *aFacet)
207 { 219 {
208 // Add the facet to the vector. 220 // Add the facet to the vector.
209 // << 221
210 if (fSolidClosed) << 222 if (solidClosed)
211 { 223 {
212 G4Exception("G4TessellatedSolid::AddFacet( << 224 G4Exception("G4TessellatedSolid::AddFacet()", "InvalidSetup",
213 JustWarning, "Attempt to add f 225 JustWarning, "Attempt to add facets when solid is closed.");
214 return false; 226 return false;
215 } 227 }
216 else if (aFacet->IsDefined()) 228 else if (aFacet->IsDefined())
217 { 229 {
218 set<G4VertexInfo,G4VertexComparator>::iter << 230 if (facets.size() == 0)
219 = fFacetList.begin(), end = fFacetList.e <<
220 G4ThreeVector p = aFacet->GetCircumcentre( <<
221 G4VertexInfo value; <<
222 value.id = (G4int)fFacetList.size(); <<
223 value.mag2 = p.x() + p.y() + p.z(); <<
224 <<
225 G4bool found = false; <<
226 if (!OutsideOfExtent(p, kCarTolerance)) <<
227 { 231 {
228 G4double kCarTolerance3 = 3 * kCarTolera << 232 facets.push_back(aFacet);
229 pos = fFacetList.lower_bound(value); << 233 }
230 << 234 else
231 it = pos; << 235 {
232 while (!found && it != end) // Loop c << 236 G4bool found = false;
>> 237 FacetI it = facets.begin();
>> 238 do
>> 239 {
>> 240 found = (**it == *aFacet);
>> 241 } while (!found && ++it!=facets.end());
>> 242
>> 243 if (found)
233 { 244 {
234 G4int id = (*it).id; << 245 delete *it;
235 G4VFacet *facet = fFacets[id]; << 246 facets.erase(it);
236 G4ThreeVector q = facet->GetCircumcent <<
237 if ((found = (facet == aFacet))) break <<
238 G4double dif = q.x() + q.y() + q.z() - <<
239 if (dif > kCarTolerance3) break; <<
240 it++; <<
241 } 247 }
242 << 248 else
243 if (fFacets.size() > 1) <<
244 { 249 {
245 it = pos; << 250 facets.push_back(aFacet);
246 while (!found && it != begin) // Lo <<
247 { <<
248 --it; <<
249 G4int id = (*it).id; <<
250 G4VFacet *facet = fFacets[id]; <<
251 G4ThreeVector q = facet->GetCircumce <<
252 found = (facet == aFacet); <<
253 if (found) break; <<
254 G4double dif = value.mag2 - (q.x() + <<
255 if (dif > kCarTolerance3) break; <<
256 } <<
257 } 251 }
258 } 252 }
259 << 253
260 if (!found) <<
261 { <<
262 fFacets.push_back(aFacet); <<
263 fFacetList.insert(value); <<
264 } <<
265 return true; 254 return true;
266 } 255 }
267 else 256 else
268 { 257 {
269 G4Exception("G4TessellatedSolid::AddFacet( << 258 G4Exception("G4TessellatedSolid::AddFacet()", "InvalidSetup",
270 JustWarning, "Attempt to add f 259 JustWarning, "Attempt to add facet not properly defined.");
271 aFacet->StreamInfo(G4cout); << 260 G4cerr << "Facet attributes:" << G4endl;
>> 261 aFacet->StreamInfo(G4cerr);
>> 262 G4cerr << G4endl;
>> 263
272 return false; 264 return false;
273 } 265 }
274 } 266 }
275 267
276 ////////////////////////////////////////////// 268 ///////////////////////////////////////////////////////////////////////////////
277 // 269 //
278 G4int G4TessellatedSolid::SetAllUsingStack(con << 270 void G4TessellatedSolid::SetSolidClosed (const G4bool t)
279 con << 271 {
280 G4b << 272 if (t)
281 { << 273 {
282 vector<G4int> xyz = voxel; << 274 vertexList.clear();
283 stack<vector<G4int> > pos; << 275 for (FacetCI it=facets.begin(); it!=facets.end(); it++)
284 pos.push(xyz); <<
285 G4int filled = 0; <<
286 <<
287 vector<G4int> candidates; <<
288 <<
289 while (!pos.empty()) // Loop checking, 13 <<
290 { <<
291 xyz = pos.top(); <<
292 pos.pop(); <<
293 G4int index = fVoxels.GetVoxelsIndex(xyz); <<
294 if (!checked[index]) <<
295 { 276 {
296 checked.SetBitNumber(index, true); << 277 size_t m = vertexList.size();
297 << 278 G4ThreeVector p(0.0,0.0,0.0);
298 if (fVoxels.IsEmpty(index)) << 279 for (size_t i=0; i<(*it)->GetNumberOfVertices(); i++)
299 { << 280 {
300 ++filled; << 281 p = (*it)->GetVertex(i);
301 << 282 G4bool found = false;
302 fInsides.SetBitNumber(index, status); << 283 size_t j = 0;
303 << 284 while (j < m && !found)
304 for (auto i = 0; i <= 2; ++i) <<
305 { 285 {
306 if (xyz[i] < max[i] - 1) << 286 G4ThreeVector q = vertexList[j];
307 { << 287 found = (q-p).mag() < 0.5*kCarTolerance;
308 xyz[i]++; << 288 if (!found) j++;
309 pos.push(xyz); << 289 }
310 xyz[i]--; << 290
311 } << 291 if (!found)
312 << 292 {
313 if (xyz[i] > 0) << 293 vertexList.push_back(p);
314 { << 294 (*it)->SetVertexIndex(i,vertexList.size()-1);
315 xyz[i]--; << 295 }
316 pos.push(xyz); << 296 else
317 xyz[i]++; << 297 {
318 } << 298 (*it)->SetVertexIndex(i,j);
319 } 299 }
320 } 300 }
321 } 301 }
322 } << 302 //
323 return filled; << 303 // Now update the maximum x, y and z limits of the volume.
324 } << 304 //
325 << 305 for (size_t i=0; i<vertexList.size(); i++)
326 ////////////////////////////////////////////// <<
327 // <<
328 void G4TessellatedSolid::PrecalculateInsides() <<
329 { <<
330 vector<G4int> voxel(3), maxVoxels(3); <<
331 for (auto i = 0; i <= 2; ++i) <<
332 maxVoxels[i] = (G4int)fVoxels.GetBoundary( <<
333 G4int size = maxVoxels[0] * maxVoxels[1] * m <<
334 <<
335 G4SurfBits checked(size-1); <<
336 fInsides.Clear(); <<
337 fInsides.ResetBitNumber(size-1); <<
338 <<
339 G4ThreeVector point; <<
340 for (voxel[2] = 0; voxel[2] < maxVoxels[2] - <<
341 { <<
342 for (voxel[1] = 0; voxel[1] < maxVoxels[1] <<
343 { 306 {
344 for (voxel[0] = 0; voxel[0] < maxVoxels[ << 307 G4ThreeVector p = vertexList[i];
>> 308 G4double x = p.x();
>> 309 G4double y = p.y();
>> 310 G4double z = p.z();
>> 311
>> 312 if (i > 0)
>> 313 {
>> 314 if (x > xMaxExtent) xMaxExtent = x;
>> 315 if (x < xMinExtent) xMinExtent = x;
>> 316 if (y > yMaxExtent) yMaxExtent = y;
>> 317 if (y < yMinExtent) yMinExtent = y;
>> 318 if (z > zMaxExtent) zMaxExtent = z;
>> 319 if (z < zMinExtent) zMinExtent = z;
>> 320 }
>> 321 else
345 { 322 {
346 G4int index = fVoxels.GetVoxelsIndex(v << 323 xMaxExtent = x;
347 if (!checked[index] && fVoxels.IsEmpty << 324 xMinExtent = x;
348 { << 325 yMaxExtent = y;
349 for (auto i = 0; i <= 2; ++i) << 326 yMinExtent = y;
350 { << 327 zMaxExtent = z;
351 point[i] = fVoxels.GetBoundary(i)[ << 328 zMinExtent = z;
352 } <<
353 auto inside = (G4bool) (InsideNoVoxe <<
354 SetAllUsingStack(voxel, maxVoxels, i <<
355 } <<
356 else checked.SetBitNumber(index); <<
357 } 329 }
358 } 330 }
359 } <<
360 } <<
361 <<
362 ////////////////////////////////////////////// <<
363 // 331 //
364 void G4TessellatedSolid::Voxelize () <<
365 { <<
366 #ifdef G4SPECSDEBUG <<
367 G4cout << "Voxelizing..." << G4endl; <<
368 #endif <<
369 fVoxels.Voxelize(fFacets); <<
370 <<
371 if (fVoxels.Empty().GetNbits() != 0u) <<
372 { <<
373 #ifdef G4SPECSDEBUG <<
374 G4cout << "Precalculating Insides..." << G <<
375 #endif <<
376 PrecalculateInsides(); <<
377 } <<
378 } <<
379 <<
380 ////////////////////////////////////////////// <<
381 // 332 //
382 // Compute extremeFacets, i.e. find those face 333 // Compute extremeFacets, i.e. find those facets that have surface
383 // planes that bound the volume. 334 // planes that bound the volume.
384 // Note that this is going to reject concaved << 335 // Note that this is going to reject concaved surfaces as being extreme. Also
385 // note that if the vertex is on the facet, di 336 // note that if the vertex is on the facet, displacement is zero, so IsInside
386 // returns true. So will this work?? Need non << 337 // returns true. So will this work?? Need non-equality
387 // "G4bool inside = displacement < 0.0;" 338 // "G4bool inside = displacement < 0.0;"
388 // or 339 // or
389 // "G4bool inside = displacement <= -0.5*kCarT << 340 // "G4bool inside = displacement <= -0.5*kCarTolerance"
390 // (Notes from PT 13/08/2007). 341 // (Notes from PT 13/08/2007).
391 // 342 //
392 void G4TessellatedSolid::SetExtremeFacets() << 343 for (FacetCI it=facets.begin(); it!=facets.end(); it++)
393 { <<
394 // Copy vertices to local array <<
395 std::size_t vsize = fVertexList.size(); <<
396 std::vector<G4ThreeVector> vertices(vsize); <<
397 for (std::size_t i = 0; i < vsize; ++i) { ve <<
398 <<
399 // Shuffle vertices <<
400 std::mt19937 gen(12345678); <<
401 std::shuffle(vertices.begin(), vertices.end( <<
402 <<
403 // Select six extreme vertices in different <<
404 G4ThreeVector points[6]; <<
405 for (auto & point : points) { point = vertic <<
406 for (std::size_t i=1; i < vsize; ++i) <<
407 { <<
408 if (vertices[i].x() < points[0].x()) point <<
409 if (vertices[i].x() > points[1].x()) point <<
410 if (vertices[i].y() < points[2].y()) point <<
411 if (vertices[i].y() > points[3].y()) point <<
412 if (vertices[i].z() < points[4].z()) point <<
413 if (vertices[i].z() > points[5].z()) point <<
414 } <<
415 <<
416 // Find extreme facets <<
417 std::size_t size = fFacets.size(); <<
418 for (std::size_t j = 0; j < size; ++j) <<
419 { <<
420 G4VFacet &facet = *fFacets[j]; <<
421 <<
422 // Check extreme vertices <<
423 if (!facet.IsInside(points[0])) continue; <<
424 if (!facet.IsInside(points[1])) continue; <<
425 if (!facet.IsInside(points[2])) continue; <<
426 if (!facet.IsInside(points[3])) continue; <<
427 if (!facet.IsInside(points[4])) continue; <<
428 if (!facet.IsInside(points[5])) continue; <<
429 <<
430 // Check vertices <<
431 G4bool isExtreme = true; <<
432 for (std::size_t i=0; i < vsize; ++i) <<
433 { 344 {
434 if (!facet.IsInside(vertices[i])) << 345 G4bool isExtreme = true;
>> 346 for (size_t i=0; i<vertexList.size(); i++)
435 { 347 {
436 isExtreme = false; << 348 if (!(*it)->IsInside(vertexList[i]))
437 break; <<
438 } <<
439 } <<
440 if (isExtreme) fExtremeFacets.insert(&face <<
441 } <<
442 } <<
443 <<
444 ////////////////////////////////////////////// <<
445 // <<
446 void G4TessellatedSolid::CreateVertexList() <<
447 { <<
448 // The algorithm: <<
449 // we will have additional vertexListSorted, <<
450 // sorted by magnitude of vertice vector. <<
451 // New candidate for fVertexList - we will d <<
452 // item which would be within its magnitude <<
453 // We will go trough until we will reach > + <<
454 // Comparison (q-p).mag() < 0.5*kCarToleranc <<
455 // They can be just stored in std::vector, w <<
456 // on binary search. <<
457 <<
458 set<G4VertexInfo,G4VertexComparator> vertexL <<
459 set<G4VertexInfo,G4VertexComparator>::iterat <<
460 = vertexListSorted.begin(), end = vertexL <<
461 G4ThreeVector p; <<
462 G4VertexInfo value; <<
463 <<
464 fVertexList.clear(); <<
465 std::size_t size = fFacets.size(); <<
466 <<
467 G4double kCarTolerance24 = kCarTolerance * k <<
468 G4double kCarTolerance3 = 3 * kCarTolerance; <<
469 vector<G4int> newIndex(100); <<
470 <<
471 for (std::size_t k = 0; k < size; ++k) <<
472 { <<
473 G4VFacet &facet = *fFacets[k]; <<
474 G4int max = facet.GetNumberOfVertices(); <<
475 <<
476 for (G4int i = 0; i < max; ++i) <<
477 { <<
478 p = facet.GetVertex(i); <<
479 value.id = (G4int)fVertexList.size(); <<
480 value.mag2 = p.x() + p.y() + p.z(); <<
481 <<
482 G4bool found = false; <<
483 G4int id = 0; <<
484 if (!OutsideOfExtent(p, kCarTolerance)) <<
485 { <<
486 pos = vertexListSorted.lower_bound(val <<
487 it = pos; <<
488 while (it != end) // Loop checking, <<
489 { 349 {
490 id = (*it).id; << 350 isExtreme = false;
491 G4ThreeVector q = fVertexList[id]; << 351 break;
492 G4double dif = (q-p).mag2(); <<
493 found = (dif < kCarTolerance24); <<
494 if (found) break; <<
495 dif = q.x() + q.y() + q.z() - value. <<
496 if (dif > kCarTolerance3) break; <<
497 ++it; <<
498 } <<
499 <<
500 if (!found && (fVertexList.size() > 1) <<
501 { <<
502 it = pos; <<
503 while (it != begin) // Loop check <<
504 { <<
505 --it; <<
506 id = (*it).id; <<
507 G4ThreeVector q = fVertexList[id]; <<
508 G4double dif = (q-p).mag2(); <<
509 found = (dif < kCarTolerance24); <<
510 if (found) break; <<
511 dif = value.mag2 - (q.x() + q.y() <<
512 if (dif > kCarTolerance3) break; <<
513 } <<
514 } 352 }
515 } 353 }
516 << 354 if (isExtreme)
517 if (!found) << 355 extremeFacets.insert(*it);
518 { <<
519 #ifdef G4SPECSDEBUG <<
520 G4cout << p.x() << ":" << p.y() << ":" <<
521 G4cout << "Adding new vertex #" << i < <<
522 << " id " << value.id << G4endl <<
523 G4cout << "===" << G4endl; <<
524 #endif <<
525 fVertexList.push_back(p); <<
526 vertexListSorted.insert(value); <<
527 begin = vertexListSorted.begin(); <<
528 end = vertexListSorted.end(); <<
529 newIndex[i] = value.id; <<
530 // <<
531 // Now update the maximum x, y and z l <<
532 // <<
533 if (value.id == 0) fMinExtent = fMaxEx <<
534 else <<
535 { <<
536 if (p.x() > fMaxExtent.x()) fMaxExte <<
537 else if (p.x() < fMinExtent.x()) fMi <<
538 if (p.y() > fMaxExtent.y()) fMaxExte <<
539 else if (p.y() < fMinExtent.y()) fMi <<
540 if (p.z() > fMaxExtent.z()) fMaxExte <<
541 else if (p.z() < fMinExtent.z()) fMi <<
542 } <<
543 } <<
544 else <<
545 { <<
546 #ifdef G4SPECSDEBUG <<
547 G4cout << p.x() << ":" << p.y() << ":" <<
548 G4cout << "Vertex #" << i << " of face <<
549 << " found, redirecting to " << <<
550 G4cout << "===" << G4endl; <<
551 #endif <<
552 newIndex[i] = id; <<
553 } <<
554 } 356 }
555 // only now it is possible to change verti << 357 solidClosed = true;
556 // <<
557 facet.SetVertices(&fVertexList); <<
558 for (G4int i = 0; i < max; ++i) <<
559 facet.SetVertexIndex(i,newIndex[i]); <<
560 } <<
561 vector<G4ThreeVector>(fVertexList).swap(fVer <<
562 <<
563 #ifdef G4SPECSDEBUG <<
564 G4double previousValue = 0.; <<
565 for (auto res=vertexListSorted.cbegin(); res <<
566 { <<
567 G4int id = (*res).id; <<
568 G4ThreeVector vec = fVertexList[id]; <<
569 G4double mvalue = vec.x() + vec.y() + vec. <<
570 if (previousValue && (previousValue - 1e-9 <<
571 G4cout << "Error in CreateVertexList: pr <<
572 << " is smaller than mvalue " << <<
573 previousValue = mvalue; <<
574 } 358 }
575 #endif << 359 else
576 } <<
577 <<
578 ////////////////////////////////////////////// <<
579 // <<
580 void G4TessellatedSolid::DisplayAllocatedMemor <<
581 { <<
582 G4int without = AllocatedMemoryWithoutVoxels <<
583 G4int with = AllocatedMemory(); <<
584 G4double ratio = (G4double) with / without; <<
585 G4cout << "G4TessellatedSolid - Allocated me <<
586 << without << "; with " << with << "; <<
587 } <<
588 <<
589 ////////////////////////////////////////////// <<
590 // <<
591 void G4TessellatedSolid::SetSolidClosed (const <<
592 { <<
593 if (t) <<
594 { 360 {
595 #ifdef G4SPECSDEBUG << 361 solidClosed = false;
596 G4cout << "Creating vertex list..." << G4e <<
597 #endif <<
598 CreateVertexList(); <<
599 <<
600 #ifdef G4SPECSDEBUG <<
601 G4cout << "Setting extreme facets..." << G <<
602 #endif <<
603 SetExtremeFacets(); <<
604 <<
605 #ifdef G4SPECSDEBUG <<
606 G4cout << "Voxelizing..." << G4endl; <<
607 #endif <<
608 Voxelize(); <<
609 <<
610 #ifdef G4SPECSDEBUG <<
611 DisplayAllocatedMemory(); <<
612 #endif <<
613 <<
614 #ifdef G4SPECSDEBUG <<
615 G4cout << "Checking Structure..." << G4end <<
616 #endif <<
617 G4int irep = CheckStructure(); <<
618 if (irep != 0) <<
619 { <<
620 if ((irep & 1) != 0) <<
621 { <<
622 std::ostringstream message; <<
623 message << "Defects in solid: " << Ge <<
624 << " - negative cubic volume, <<
625 G4Exception("G4TessellatedSolid::SetS <<
626 "GeomSolids1001", JustWar <<
627 } <<
628 if ((irep & 2) != 0) <<
629 { <<
630 std::ostringstream message; <<
631 message << "Defects in solid: " << Ge <<
632 << " - some facets have wrong <<
633 G4Exception("G4TessellatedSolid::SetS <<
634 "GeomSolids1001", JustWar <<
635 } <<
636 if ((irep & 4) != 0) <<
637 { <<
638 std::ostringstream message; <<
639 message << "Defects in solid: " << Ge <<
640 << " - there are holes in the <<
641 G4Exception("G4TessellatedSolid::SetS <<
642 "GeomSolids1001", JustWar <<
643 } <<
644 } <<
645 } 362 }
646 fSolidClosed = t; <<
647 } 363 }
648 364
649 ////////////////////////////////////////////// 365 ///////////////////////////////////////////////////////////////////////////////
650 // 366 //
651 // GetSolidClosed 367 // GetSolidClosed
652 // 368 //
653 // Used to determine whether the solid is clos << 369 // Used to determine whether the solid is closed to adding further facets.
654 // 370 //
655 G4bool G4TessellatedSolid::GetSolidClosed () c 371 G4bool G4TessellatedSolid::GetSolidClosed () const
656 { << 372 {return solidClosed;}
657 return fSolidClosed; <<
658 } <<
659 <<
660 ////////////////////////////////////////////// <<
661 // <<
662 // CheckStructure <<
663 // <<
664 // Checks structure of the solid. Return value <<
665 // defect indicators, if any (0 means no defec <<
666 // 1 - cubic volume is negative, wrong orien <<
667 // 2 - some facets have wrong orientation <<
668 // 4 - holes in the surface <<
669 // <<
670 G4int G4TessellatedSolid::CheckStructure() con <<
671 { <<
672 G4int nedge = 0; <<
673 std::size_t nface = fFacets.size(); <<
674 <<
675 // Calculate volume <<
676 // <<
677 G4double volume = 0.; <<
678 for (std::size_t i = 0; i < nface; ++i) <<
679 { <<
680 G4VFacet& facet = *fFacets[i]; <<
681 nedge += facet.GetNumberOfVertices(); <<
682 volume += facet.GetArea()*(facet.GetVertex <<
683 } <<
684 auto ivolume = static_cast<G4int>(volume <= <<
685 <<
686 // Create sorted vector of edges <<
687 // <<
688 std::vector<int64_t> iedge(nedge); <<
689 G4int kk = 0; <<
690 for (std::size_t i = 0; i < nface; ++i) <<
691 { <<
692 G4VFacet& facet = *fFacets[i]; <<
693 G4int nnode = facet.GetNumberOfVertices(); <<
694 for (G4int k = 0; k < nnode; ++k) <<
695 { <<
696 int64_t i1 = facet.GetVertexIndex((k == <<
697 int64_t i2 = facet.GetVertexIndex(k); <<
698 auto inverse = static_cast<int64_t>(i2 <<
699 if (inverse != 0) std::swap(i1, i2); <<
700 iedge[kk++] = i1*1000000000 + i2*2 + inv <<
701 } <<
702 } <<
703 std::sort(iedge.begin(), iedge.end()); <<
704 <<
705 // Check edges, correct structure should con <<
706 // with different orientation <<
707 // <<
708 G4int iorder = 0; <<
709 G4int ihole = 0; <<
710 G4int i = 0; <<
711 while (i < nedge - 1) <<
712 { <<
713 if (iedge[i + 1] - iedge[i] == 1) // paire <<
714 { <<
715 i += 2; <<
716 } <<
717 else if (iedge[i + 1] == iedge[i]) // pair <<
718 { <<
719 iorder = 2; <<
720 i += 2; <<
721 } <<
722 else // unpaired edge <<
723 { <<
724 ihole = 4; <<
725 i++; <<
726 } <<
727 } <<
728 return ivolume + iorder + ihole; <<
729 } <<
730 373
731 ////////////////////////////////////////////// 374 ///////////////////////////////////////////////////////////////////////////////
732 // 375 //
733 // operator+= 376 // operator+=
734 // 377 //
735 // This operator allows the user to add two te 378 // This operator allows the user to add two tessellated solids together, so
736 // that the solid on the left then includes al 379 // that the solid on the left then includes all of the facets in the solid
737 // on the right. Note that copies of the face 380 // on the right. Note that copies of the facets are generated, rather than
738 // using the original facet set of the solid o 381 // using the original facet set of the solid on the right.
739 // 382 //
740 G4TessellatedSolid& << 383 const G4TessellatedSolid &G4TessellatedSolid::operator+=
741 G4TessellatedSolid::operator+=(const G4Tessell << 384 (const G4TessellatedSolid &right)
742 { 385 {
743 G4int size = right.GetNumberOfFacets(); << 386 for (size_t i=0; i<right.GetNumberOfFacets(); i++)
744 for (G4int i = 0; i < size; ++i) <<
745 AddFacet(right.GetFacet(i)->GetClone()); 387 AddFacet(right.GetFacet(i)->GetClone());
746 <<
747 return *this; 388 return *this;
748 } 389 }
749 390
750 ////////////////////////////////////////////// 391 ///////////////////////////////////////////////////////////////////////////////
751 // 392 //
752 // GetNumberOfFacets << 393 // GetFacet
753 // 394 //
754 G4int G4TessellatedSolid::GetNumberOfFacets() << 395 // Access pointer to facet in solid, indexed by integer i.
>> 396 //
>> 397 G4VFacet *G4TessellatedSolid::GetFacet (size_t i) const
755 { 398 {
756 return (G4int)fFacets.size(); << 399 return facets[i];
757 } 400 }
758 401
759 ////////////////////////////////////////////// 402 ///////////////////////////////////////////////////////////////////////////////
760 // 403 //
761 EInside G4TessellatedSolid::InsideVoxels(const << 404 // GetNumberOfFacets
>> 405 //
>> 406 size_t G4TessellatedSolid::GetNumberOfFacets () const
762 { 407 {
763 // << 408 return facets.size();
764 // First the simple test - check if we're ou <<
765 // of the tessellated solid. <<
766 // <<
767 if (OutsideOfExtent(p, kCarTolerance)) <<
768 return kOutside; <<
769 <<
770 vector<G4int> startingVoxel(3); <<
771 fVoxels.GetVoxel(startingVoxel, p); <<
772 <<
773 const G4double dirTolerance = 1.0E-14; <<
774 <<
775 const vector<G4int> &startingCandidates = <<
776 fVoxels.GetCandidates(startingVoxel); <<
777 std::size_t limit = startingCandidates.size( <<
778 if (limit == 0 && (fInsides.GetNbits() != 0u <<
779 { <<
780 G4int index = fVoxels.GetPointIndex(p); <<
781 EInside location = fInsides[index] ? kInsi <<
782 return location; <<
783 } <<
784 <<
785 G4double minDist = kInfinity; <<
786 <<
787 for(std::size_t i = 0; i < limit; ++i) <<
788 { <<
789 G4int candidate = startingCandidates[i]; <<
790 G4VFacet &facet = *fFacets[candidate]; <<
791 G4double dist = facet.Distance(p,minDist); <<
792 if (dist < minDist) minDist = dist; <<
793 if (dist <= kCarToleranceHalf) <<
794 return kSurface; <<
795 } <<
796 <<
797 // The following is something of an adaptati <<
798 // Rickard Holmberg augmented with informati <<
799 // "Geometric Tools for Computer Graphics," <<
800 // we're trying to determine whether we're i <<
801 // a few rays and determining if the first s <<
802 // vector between 0 to pi/2 (out-going) or p <<
803 // We should also avoid rays which are nearl <<
804 // tessellated surface, and therefore produc <<
805 // For the moment, this is a bit over-engine <<
806 // <<
807 G4double distOut = kInfinity; <<
808 G4double distIn = kInfinity; <<
809 G4double distO = 0.0; <<
810 G4double distI = 0.0; <<
811 G4double distFromSurfaceO = 0.0; <<
812 G4double distFromSurfaceI = 0.0; <<
813 G4ThreeVector normalO, normalI; <<
814 G4bool crossingO = false; <<
815 G4bool crossingI = false; <<
816 EInside location = kOutside; <<
817 G4int sm = 0; <<
818 <<
819 G4bool nearParallel = false; <<
820 do // Loop checking, 13.08.2015, G.Cosmo <<
821 { <<
822 // We loop until we find direction where t <<
823 // to the surface of any facet since this <<
824 // case is that the angles should be suffi <<
825 // are 20 random directions to select from <<
826 // <<
827 distOut = distIn = kInfinity; <<
828 const G4ThreeVector& v = fRandir[sm]; <<
829 ++sm; <<
830 // <<
831 // This code could be voxelized by the sam <<
832 // DistanceToOut(). We will traverse throu <<
833 // intersect only for those, which would b <<
834 // checked before. <<
835 // <<
836 G4ThreeVector currentPoint = p; <<
837 G4ThreeVector direction = v.unit(); <<
838 // G4SurfBits exclusion(fVoxels.GetBitsPer <<
839 vector<G4int> curVoxel(3); <<
840 curVoxel = startingVoxel; <<
841 G4double shiftBonus = kCarTolerance; <<
842 <<
843 G4bool crossed = false; <<
844 G4bool started = true; <<
845 <<
846 do // Loop checking, 13.08.2015, G.Cosm <<
847 { <<
848 const vector<G4int> &candidates = <<
849 started ? startingCandidates : fVoxels <<
850 started = false; <<
851 if (auto candidatesCount = (G4int)candid <<
852 { <<
853 for (G4int i = 0 ; i < candidatesCount <<
854 { <<
855 G4int candidate = candidates[i]; <<
856 // bits.SetBitNumber(candidate); <<
857 G4VFacet& facet = *fFacets[candidate <<
858 <<
859 crossingO = facet.Intersect(p,v,true <<
860 crossingI = facet.Intersect(p,v,fals <<
861 <<
862 if (crossingO || crossingI) <<
863 { <<
864 crossed = true; <<
865 <<
866 nearParallel = (crossingO <<
867 && std::fabs(normalO.dot( <<
868 || (crossingI && std::fab <<
869 if (!nearParallel) <<
870 { <<
871 if (crossingO && distO > 0.0 && <<
872 distOut = distO; <<
873 if (crossingI && distI > 0.0 && <<
874 distIn = distI; <<
875 } <<
876 else break; <<
877 } <<
878 } <<
879 if (nearParallel) break; <<
880 } <<
881 else <<
882 { <<
883 if (!crossed) <<
884 { <<
885 G4int index = fVoxels.GetVoxelsIndex <<
886 G4bool inside = fInsides[index]; <<
887 location = inside ? kInside : kOutsi <<
888 return location; <<
889 } <<
890 } <<
891 <<
892 G4double shift=fVoxels.DistanceToNext(cu <<
893 if (shift == kInfinity) break; <<
894 <<
895 currentPoint += direction * (shift + shi <<
896 } <<
897 while (fVoxels.UpdateCurrentVoxel(currentP <<
898 <<
899 } <<
900 while (nearParallel && sm != fMaxTries); <<
901 // <<
902 // Here we loop through the facets to find o <<
903 // between the ray and that facet. The test <<
904 // the ray is entering the facet or exiting. <<
905 // <<
906 #ifdef G4VERBOSE <<
907 if (sm == fMaxTries) <<
908 { <<
909 // <<
910 // We've run out of random vector directio <<
911 // low (nTries <= 0.5*maxTries) then this <<
912 // something wrong with geometry. <<
913 // <<
914 std::ostringstream message; <<
915 G4long oldprc = message.precision(16); <<
916 message << "Cannot determine whether point <<
917 << G4endl <<
918 << "Solid name = " << GetName() < <<
919 << "Geometry Type = " << fGeometryTyp <<
920 << "Number of facets = " << fFacets.size <<
921 << "Position:" << G4endl << G4endl <<
922 << "p.x() = " << p.x()/mm << " mm" << <<
923 << "p.y() = " << p.y()/mm << " mm" << <<
924 << "p.z() = " << p.z()/mm << " mm"; <<
925 message.precision(oldprc); <<
926 G4Exception("G4TessellatedSolid::Inside()" <<
927 "GeomSolids1002", JustWarning, <<
928 } <<
929 #endif <<
930 <<
931 // In the next if-then-elseif G4String the l <<
932 // (1) You don't hit anything so cannot be i <<
933 // constructed correctly! <<
934 // (2) Distance to inside (ie. nearest facet <<
935 // shorter than distance to outside (nea <<
936 // facet) - on condition of safety dista <<
937 // (3) Distance to outside is shorter than d <<
938 // we're inside. <<
939 // <<
940 if (distIn == kInfinity && distOut == kInfin <<
941 location = kOutside; <<
942 else if (distIn <= distOut - kCarToleranceHa <<
943 location = kOutside; <<
944 else if (distOut <= distIn - kCarToleranceHa <<
945 location = kInside; <<
946 <<
947 return location; <<
948 } 409 }
949 410
950 ////////////////////////////////////////////// 411 ///////////////////////////////////////////////////////////////////////////////
951 // 412 //
952 EInside G4TessellatedSolid::InsideNoVoxels (co << 413 // EInside G4TessellatedSolid::Inside (const G4ThreeVector &p) const
>> 414 //
>> 415 // This method must return:
>> 416 // * kOutside if the point at offset p is outside the shape
>> 417 // boundaries plus kCarTolerance/2,
>> 418 // * kSurface if the point is <= kCarTolerance/2 from a surface, or
>> 419 // * kInside otherwise.
>> 420 //
>> 421 EInside G4TessellatedSolid::Inside (const G4ThreeVector &p) const
953 { 422 {
954 // << 423 //
955 // First the simple test - check if we're ou << 424 // First the simple test - check if we're outside of the X-Y-Z extremes
956 // of the tessellated solid. << 425 // of the tessellated solid.
957 // << 426 //
958 if (OutsideOfExtent(p, kCarTolerance)) << 427 if ( p.x() < xMinExtent - kCarTolerance ||
>> 428 p.x() > xMaxExtent + kCarTolerance ||
>> 429 p.y() < yMinExtent - kCarTolerance ||
>> 430 p.y() > yMaxExtent + kCarTolerance ||
>> 431 p.z() < zMinExtent - kCarTolerance ||
>> 432 p.z() > zMaxExtent + kCarTolerance )
>> 433 {
959 return kOutside; 434 return kOutside;
960 << 435 }
961 const G4double dirTolerance = 1.0E-14; <<
962 436
963 G4double minDist = kInfinity; 437 G4double minDist = kInfinity;
964 // << 438 //
965 // Check if we are close to a surface << 439 //
966 // << 440 // Check if we are close to a surface
967 std::size_t size = fFacets.size(); << 441 //
968 for (std::size_t i = 0; i < size; ++i) << 442 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
969 { 443 {
970 G4VFacet& facet = *fFacets[i]; << 444 G4double dist = (*f)->Distance(p,minDist);
971 G4double dist = facet.Distance(p,minDist); <<
972 if (dist < minDist) minDist = dist; 445 if (dist < minDist) minDist = dist;
973 if (dist <= kCarToleranceHalf) << 446 if (dist <= 0.5*kCarTolerance)
974 { 447 {
975 return kSurface; 448 return kSurface;
976 } 449 }
977 } 450 }
978 // << 451 //
979 // The following is something of an adaptati << 452 //
980 // Rickard Holmberg augmented with informati << 453 // The following is something of an adaptation of the method implemented by
981 // "Geometric Tools for Computer Graphics," << 454 // Rickard Holmberg augmented with information from Schneider & Eberly,
982 // trying to determine whether we're inside << 455 // "Geometric Tools for Computer Graphics," pp700-701, 2003. In essence, we're
983 // rays and determining if the first surface << 456 // trying to determine whether we're inside the volume by projecting a few rays
984 // between 0 to pi/2 (out-going) or pi/2 to << 457 // and determining if the first surface crossed is has a normal vector between
985 // avoid rays which are nearly within the pl << 458 // 0 to pi/2 (out-going) or pi/2 to pi (in-going). We should also avoid rays
986 // and therefore produce rays randomly. For << 459 // which are nearly within the plane of the tessellated surface, and therefore
987 // over-engineered (belt-braces-and-ducttape << 460 // produce rays randomly. For the moment, this is a bit over-engineered
988 // << 461 // (belt-braces-and-ducttape).
>> 462 //
989 #if G4SPECSDEBUG 463 #if G4SPECSDEBUG
990 G4int nTry = 7; 464 G4int nTry = 7;
991 #else 465 #else
992 G4int nTry = 3; 466 G4int nTry = 3;
993 #endif 467 #endif
994 G4double distOut = kInfinity; 468 G4double distOut = kInfinity;
995 G4double distIn = kInfinity; 469 G4double distIn = kInfinity;
996 G4double distO = 0.0; 470 G4double distO = 0.0;
997 G4double distI = 0.0; 471 G4double distI = 0.0;
998 G4double distFromSurfaceO = 0.0; 472 G4double distFromSurfaceO = 0.0;
999 G4double distFromSurfaceI = 0.0; 473 G4double distFromSurfaceI = 0.0;
1000 G4ThreeVector normalO(0.0,0.0,0.0); 474 G4ThreeVector normalO(0.0,0.0,0.0);
1001 G4ThreeVector normalI(0.0,0.0,0.0); 475 G4ThreeVector normalI(0.0,0.0,0.0);
1002 G4bool crossingO = false; 476 G4bool crossingO = false;
1003 G4bool crossingI = false; 477 G4bool crossingI = false;
1004 EInside location = kOutside; 478 EInside location = kOutside;
1005 EInside locationprime = kOutside; 479 EInside locationprime = kOutside;
1006 G4int sm = 0; << 480 G4int m = 0;
1007 481
1008 for (G4int i=0; i<nTry; ++i) << 482 for (G4int i=0; i<nTry; i++)
1009 { 483 {
1010 G4bool nearParallel = false; 484 G4bool nearParallel = false;
1011 do // Loop checking, 13.08.2015, G.Cos << 485 do
1012 { 486 {
1013 // << 487 //
1014 // We loop until we find direction wher << 488 //
1015 // to the surface of any facet since th << 489 // We loop until we find direction where the vector is not nearly parallel
1016 // case is that the angles should be su << 490 // to the surface of any facet since this causes ambiguities. The usual
1017 // are 20 random directions to select f << 491 // case is that the angles should be sufficiently different, but there are 20
1018 // << 492 // random directions to select from - hopefully sufficient.
1019 distOut = distIn = kInfinity; << 493 //
1020 G4ThreeVector v = fRandir[sm]; << 494 distOut = kInfinity;
1021 sm++; << 495 distIn = kInfinity;
1022 auto f = fFacets.cbegin(); << 496 G4ThreeVector v = randir[m];
1023 << 497 m++;
1024 do // Loop checking, 13.08.2015, G.C << 498 FacetCI f = facets.begin();
>> 499 do
1025 { 500 {
1026 // << 501 //
1027 // Here we loop through the facets to << 502 //
1028 // intersection between the ray and t << 503 // Here we loop through the facets to find out if there is an intersection
1029 // separately whether the ray is ente << 504 // between the ray and that facet. The test if performed separately whether
1030 // << 505 // the ray is entering the facet or exiting.
1031 crossingO = ((*f)->Intersect(p,v,true << 506 //
1032 crossingI = ((*f)->Intersect(p,v,fals << 507 crossingO = ((*f)->Intersect(p,v,true,distO,distFromSurfaceO,normalO));
>> 508 crossingI = ((*f)->Intersect(p,v,false,distI,distFromSurfaceI,normalI));
1033 if (crossingO || crossingI) 509 if (crossingO || crossingI)
1034 { 510 {
1035 nearParallel = (crossingO && std::f << 511 nearParallel = (crossingO && std::abs(normalO.dot(v))<dirTolerance) ||
1036 || (crossingI && std::f << 512 (crossingI && std::abs(normalI.dot(v))<dirTolerance);
1037 if (!nearParallel) 513 if (!nearParallel)
1038 { 514 {
1039 if (crossingO && distO > 0.0 && d 515 if (crossingO && distO > 0.0 && distO < distOut) distOut = distO;
1040 if (crossingI && distI > 0.0 && d 516 if (crossingI && distI > 0.0 && distI < distIn) distIn = distI;
1041 } 517 }
1042 } 518 }
1043 } while (!nearParallel && ++f != fFacet << 519 } while (!nearParallel && ++f!=facets.end());
1044 } while (nearParallel && sm != fMaxTries) << 520 } while (nearParallel && m!=maxTries);
1045 521
1046 #ifdef G4VERBOSE 522 #ifdef G4VERBOSE
1047 if (sm == fMaxTries) << 523 if (m == maxTries)
1048 { 524 {
1049 // << 525 //
1050 // We've run out of random vector direc << 526 //
1051 // sufficiently low (nTries <= 0.5*maxT << 527 // We've run out of random vector directions. If nTries is set sufficiently
1052 // that there is something wrong with g << 528 // low (nTries <= 0.5*maxTries) then this would indicate that there is
1053 // << 529 // something wrong with geometry.
1054 std::ostringstream message; << 530 //
1055 G4long oldprc = message.precision(16); << 531 G4cout.precision(16) ;
1056 message << "Cannot determine whether po << 532 G4cout << G4endl ;
1057 << G4endl << 533 G4cout << "Solid name = " << GetName() << G4endl;
1058 << "Solid name = " << GetName() << 534 G4cout << "Geometry Type = " << geometryType << G4endl;
1059 << "Geometry Type = " << fGeometry << 535 G4cout << "Number of facets = " << facets.size() << G4endl;
1060 << "Number of facets = " << fFacets.s << 536 G4cout << "Position:" << G4endl << G4endl ;
1061 << "Position:" << G4endl << G4endl << 537 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1062 << "p.x() = " << p.x()/mm << " mm" << 538 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1063 << "p.y() = " << p.y()/mm << " mm" << 539 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1064 << "p.z() = " << p.z()/mm << " mm"; <<
1065 message.precision(oldprc); <<
1066 G4Exception("G4TessellatedSolid::Inside 540 G4Exception("G4TessellatedSolid::Inside()",
1067 "GeomSolids1002", JustWarning, messag << 541 "UnknownInsideOutside-MaxTries", JustWarning,
>> 542 "Cannot determine whether point is inside or outside volume!");
1068 } 543 }
1069 #endif 544 #endif
1070 // << 545 //
1071 // In the next if-then-elseif G4String th << 546 //
1072 // (1) You don't hit anything so cannot b << 547 // In the next if-then-elseif string the logic is as follows:
1073 // constructed correctly! << 548 // (1) You don't hit anything so cannot be inside volume, provided volume
1074 // (2) Distance to inside (ie. nearest fa << 549 // constructed correctly!
1075 // shorter than distance to outside ( << 550 // (2) Distance to inside (ie. nearest facet such that you enter facet) is
1076 // facet) - on condition of safety di << 551 // shorter than distance to outside (nearest facet such that you exit
1077 // (3) Distance to outside is shorter tha << 552 // facet) - on condition of safety distance - therefore we're outside.
1078 // we're inside. << 553 // (3) Distance to outside is shorter than distance to inside therefore we're
1079 // << 554 // inside.
>> 555 //
1080 if (distIn == kInfinity && distOut == kIn 556 if (distIn == kInfinity && distOut == kInfinity)
1081 locationprime = kOutside; 557 locationprime = kOutside;
1082 else if (distIn <= distOut - kCarToleranc << 558 else if (distIn <= distOut - kCarTolerance*0.5)
1083 locationprime = kOutside; 559 locationprime = kOutside;
1084 else if (distOut <= distIn - kCarToleranc << 560 else if (distOut <= distIn - kCarTolerance*0.5)
1085 locationprime = kInside; 561 locationprime = kInside;
1086 562
1087 if (i == 0) location = locationprime; << 563 if (i == 0) { location = locationprime; }
1088 } << 564 #ifdef G4VERBOSE
1089 << 565 else if (locationprime != location)
1090 return location; << 566 {
1091 } <<
1092 <<
1093 ///////////////////////////////////////////// <<
1094 // 567 //
1095 // Return index of the facet closest to the p <<
1096 // be located on the surface. Return -1 if no <<
1097 // 568 //
1098 G4int G4TessellatedSolid::GetFacetIndex (cons << 569 // Different ray directions result in different answer. Seems like the
1099 { << 570 // geometry is not constructed correctly.
1100 G4int index = -1; << 571 //
1101 << 572 G4cout.precision(16) ;
1102 if (fVoxels.GetCountOfVoxels() > 1) << 573 G4cout << G4endl ;
1103 { << 574 G4cout << "Solid name = " << GetName() << G4endl;
1104 vector<G4int> curVoxel(3); << 575 G4cout << "Geometry Type = " << geometryType << G4endl;
1105 fVoxels.GetVoxel(curVoxel, p); << 576 G4cout << "Number of facets = " << facets.size() << G4endl;
1106 const vector<G4int> &candidates = fVoxels << 577 G4cout << "Position:" << G4endl << G4endl ;
1107 if (auto limit = (G4int)candidates.size() << 578 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1108 { << 579 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1109 G4double minDist = kInfinity; << 580 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1110 for(G4int i = 0 ; i < limit ; ++i) << 581 G4Exception("G4TessellatedSolid::Inside()",
1111 { << 582 "UnknownInsideOutside", JustWarning,
1112 G4int candidate = candidates[i]; << 583 "Cannot determine whether point is inside or outside volume!");
1113 G4VFacet& facet = *fFacets[candidate] <<
1114 G4double dist = facet.Distance(p, min <<
1115 if (dist <= kCarToleranceHalf) return <<
1116 if (dist < minDist) <<
1117 { <<
1118 minDist = dist; <<
1119 index = candidate; <<
1120 } <<
1121 } <<
1122 } <<
1123 } <<
1124 else <<
1125 { <<
1126 G4double minDist = kInfinity; <<
1127 std::size_t size = fFacets.size(); <<
1128 for (std::size_t i = 0; i < size; ++i) <<
1129 { <<
1130 G4VFacet& facet = *fFacets[i]; <<
1131 G4double dist = facet.Distance(p, minDi <<
1132 if (dist < minDist) <<
1133 { <<
1134 minDist = dist; <<
1135 index = (G4int)i; <<
1136 } <<
1137 } 584 }
>> 585 #endif
1138 } 586 }
1139 return index; << 587
>> 588 return location;
1140 } 589 }
1141 590
1142 ///////////////////////////////////////////// 591 ///////////////////////////////////////////////////////////////////////////////
1143 // 592 //
>> 593 // G4ThreeVector G4TessellatedSolid::SurfaceNormal (const G4ThreeVector &p) const
>> 594 //
1144 // Return the outwards pointing unit normal o 595 // Return the outwards pointing unit normal of the shape for the
1145 // surface closest to the point at offset p. 596 // surface closest to the point at offset p.
1146 // <<
1147 G4bool G4TessellatedSolid::Normal (const G4Th <<
1148 G4Th <<
1149 { <<
1150 G4double minDist; <<
1151 G4VFacet* facet = nullptr; <<
1152 <<
1153 if (fVoxels.GetCountOfVoxels() > 1) <<
1154 { <<
1155 vector<G4int> curVoxel(3); <<
1156 fVoxels.GetVoxel(curVoxel, p); <<
1157 const vector<G4int> &candidates = fVoxels <<
1158 // fVoxels.GetCandidatesVoxelArray(p, can <<
1159 597
1160 if (auto limit = (G4int)candidates.size() << 598 G4ThreeVector G4TessellatedSolid::SurfaceNormal (const G4ThreeVector &p) const
1161 { << 599 {
1162 minDist = kInfinity; << 600 FacetCI minFacet;
1163 for(G4int i = 0 ; i < limit ; ++i) << 601 G4double minDist = kInfinity;
1164 { << 602 G4double dist = 0.0;
1165 G4int candidate = candidates[i]; << 603 G4ThreeVector normal;
1166 G4VFacet &fct = *fFacets[candidate]; << 604
1167 G4double dist = fct.Distance(p,minDis << 605 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
1168 if (dist < minDist) minDist = dist; <<
1169 if (dist <= kCarToleranceHalf) <<
1170 { <<
1171 aNormal = fct.GetSurfaceNormal(); <<
1172 return true; <<
1173 } <<
1174 } <<
1175 } <<
1176 minDist = MinDistanceFacet(p, true, facet <<
1177 } <<
1178 else <<
1179 { 606 {
1180 minDist = kInfinity; << 607 dist = (*f)->Distance(p,minDist);
1181 std::size_t size = fFacets.size(); << 608 if (dist < minDist)
1182 for (std::size_t i = 0; i < size; ++i) <<
1183 { 609 {
1184 G4VFacet& f = *fFacets[i]; << 610 minDist = dist;
1185 G4double dist = f.Distance(p, minDist); << 611 minFacet = f;
1186 if (dist < minDist) <<
1187 { <<
1188 minDist = dist; <<
1189 facet = &f; <<
1190 } <<
1191 } 612 }
1192 } 613 }
1193 << 614
1194 if (minDist != kInfinity) 615 if (minDist != kInfinity)
1195 { 616 {
1196 if (facet != nullptr) { aNormal = facet- << 617 normal = (*minFacet)->GetSurfaceNormal();
1197 return minDist <= kCarToleranceHalf; <<
1198 } 618 }
1199 else 619 else
1200 { 620 {
1201 #ifdef G4VERBOSE 621 #ifdef G4VERBOSE
1202 std::ostringstream message; << 622 G4cout << "WARNING - G4TessellatedSolid::SurfaceNormal(p)" << G4endl
1203 message << "Point p is not on surface !?" << 623 << " No facets found for point: " << p << " !" << G4endl
1204 << " No facets found for point << 624 << " Returning approximated value for normal." << G4endl;
1205 << " Returning approximated va << 625 G4Exception("G4TessellatedSolid::SurfaceNormal(p)", "Notification",
1206 << 626 JustWarning, "Point p is not on surface !?" );
1207 G4Exception("G4TessellatedSolid::SurfaceN <<
1208 "GeomSolids1002", JustWarning <<
1209 #endif 627 #endif
1210 aNormal = (p.z() > 0 ? G4ThreeVector(0,0, << 628 normal = (p.z()>0 ? G4ThreeVector(0,0,1) : G4ThreeVector(0,0,-1));
1211 return false; <<
1212 } 629 }
>> 630
>> 631 return normal;
1213 } 632 }
1214 633
1215 ///////////////////////////////////////////// 634 ///////////////////////////////////////////////////////////////////////////////
1216 // 635 //
1217 // G4double DistanceToIn(const G4ThreeVector& 636 // G4double DistanceToIn(const G4ThreeVector& p, const G4ThreeVector& v)
1218 // 637 //
1219 // Return the distance along the normalised v 638 // Return the distance along the normalised vector v to the shape,
1220 // from the point at offset p. If there is no 639 // from the point at offset p. If there is no intersection, return
1221 // kInfinity. The first intersection resultin << 640 // kInfinity. The first intersection resulting from ‘leaving’ a
1222 // surface/volume is discarded. Hence, this i 641 // surface/volume is discarded. Hence, this is tolerant of points on
1223 // surface of shape. 642 // surface of shape.
1224 // << 643
1225 G4double << 644 G4double G4TessellatedSolid::DistanceToIn (const G4ThreeVector &p,
1226 G4TessellatedSolid::DistanceToInNoVoxels (con << 645 const G4ThreeVector &v) const
1227 con <<
1228 <<
1229 { 646 {
1230 G4double minDist = kInfinity; 647 G4double minDist = kInfinity;
1231 G4double dist = 0.0; 648 G4double dist = 0.0;
1232 G4double distFromSurface = 0.0; 649 G4double distFromSurface = 0.0;
1233 G4ThreeVector normal; << 650 G4ThreeVector normal(0.0,0.0,0.0);
1234 << 651
1235 #if G4SPECSDEBUG 652 #if G4SPECSDEBUG
1236 if (Inside(p) == kInside ) << 653 if ( Inside(p) == kInside )
1237 { 654 {
1238 std::ostringstream message; << 655 G4cout.precision(16) ;
1239 G4int oldprc = message.precision(16) ; << 656 G4cout << G4endl ;
1240 message << "Point p is already inside!?" << 657 // DumpInfo();
1241 << "Position:" << G4endl << G4endl << 658 G4cout << "Position:" << G4endl << G4endl ;
1242 << " p.x() = " << p.x()/mm << " mm" << 659 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1243 << " p.y() = " << p.y()/mm << " mm" << 660 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1244 << " p.z() = " << p.z()/mm << " mm" << 661 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1245 << "DistanceToOut(p) == " << DistanceTo << 662 G4cout << "DistanceToOut(p) == " << DistanceToOut(p) << G4endl;
1246 message.precision(oldprc) ; << 663 G4Exception("G4TriangularFacet::DistanceToIn(p,v)", "Notification", JustWarning,
1247 G4Exception("G4TriangularFacet::DistanceT << 664 "Point p is already inside!?" );
1248 "GeomSolids1002", JustWarning <<
1249 } 665 }
1250 #endif 666 #endif
1251 667
1252 std::size_t size = fFacets.size(); << 668 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
1253 for (std::size_t i = 0; i < size; ++i) <<
1254 { 669 {
1255 G4VFacet& facet = *fFacets[i]; << 670 if ((*f)->Intersect(p,v,false,dist,distFromSurface,normal))
1256 if (facet.Intersect(p,v,false,dist,distFr <<
1257 { 671 {
1258 // << 672 if (distFromSurface > 0.5*kCarTolerance && dist >= 0.0 && dist < minDist)
1259 // set minDist to the new distance to c <<
1260 // is in positive direction and point i <<
1261 // within 0.5*kCarTolerance of the surf <<
1262 // zero and leave member function immed <<
1263 // proposed by & credit to Akira Okumur <<
1264 // <<
1265 if (distFromSurface > kCarToleranceHalf <<
1266 { 673 {
1267 minDist = dist; 674 minDist = dist;
1268 } 675 }
1269 else <<
1270 { <<
1271 if (-kCarToleranceHalf <= dist && dis <<
1272 { <<
1273 return 0.0; <<
1274 } <<
1275 else <<
1276 { <<
1277 if (distFromSurface > -kCarToleran <<
1278 && distFromSurface < kCarToleran <<
1279 { <<
1280 minDist = dist; <<
1281 } <<
1282 } <<
1283 } <<
1284 } 676 }
1285 } 677 }
>> 678
1286 return minDist; 679 return minDist;
1287 } 680 }
1288 681
1289 ///////////////////////////////////////////// 682 ///////////////////////////////////////////////////////////////////////////////
1290 // 683 //
1291 G4double << 684 // G4double DistanceToIn(const G4ThreeVector& p)
1292 G4TessellatedSolid::DistanceToOutNoVoxels (co << 685 //
1293 co << 686 // Calculate distance to nearest surface of shape from an outside point p. The
1294 << 687 // distance can be an underestimate.
1295 <<
1296 <<
1297 { <<
1298 G4double minDist = kInfinity; <<
1299 G4double dist = 0.0; <<
1300 G4double distFromSurface = 0.0; <<
1301 G4ThreeVector normal, minNormal; <<
1302 688
>> 689 G4double G4TessellatedSolid::DistanceToIn (const G4ThreeVector &p) const
>> 690 {
>> 691 G4double minDist = kInfinity;
>> 692 G4double dist = 0.0;
>> 693
1303 #if G4SPECSDEBUG 694 #if G4SPECSDEBUG
1304 if ( Inside(p) == kOutside ) << 695 if ( Inside(p) == kInside )
1305 { 696 {
1306 std::ostringstream message; << 697 G4cout.precision(16) ;
1307 G4int oldprc = message.precision(16) ; << 698 G4cout << G4endl ;
1308 message << "Point p is already outside!?" << 699 // DumpInfo();
1309 << "Position:" << G4endl << G4endl << 700 G4cout << "Position:" << G4endl << G4endl ;
1310 << " p.x() = " << p.x()/mm << " mm" << 701 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1311 << " p.y() = " << p.y()/mm << " mm" << 702 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1312 << " p.z() = " << p.z()/mm << " mm" << 703 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1313 << "DistanceToIn(p) == " << DistanceToI << 704 G4cout << "DistanceToOut(p) == " << DistanceToOut(p) << G4endl;
1314 message.precision(oldprc) ; << 705 G4Exception("G4TriangularFacet::DistanceToIn(p)", "Notification", JustWarning,
1315 G4Exception("G4TriangularFacet::DistanceT << 706 "Point p is already inside!?" );
1316 "GeomSolids1002", JustWarning <<
1317 } 707 }
1318 #endif 708 #endif
1319 709
1320 G4bool isExtreme = false; << 710 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
1321 std::size_t size = fFacets.size(); <<
1322 for (std::size_t i = 0; i < size; ++i) <<
1323 { 711 {
1324 G4VFacet& facet = *fFacets[i]; << 712 dist = (*f)->Distance(p,minDist,false);
1325 if (facet.Intersect(p,v,true,dist,distFro << 713 if (dist < minDist) { minDist = dist; }
1326 { <<
1327 if (distFromSurface > 0.0 && distFromSu <<
1328 && facet.Distance(p,kCarTolerance) <= <<
1329 { <<
1330 // We are on a surface. Return zero. <<
1331 aConvex = (fExtremeFacets.find(&facet <<
1332 // Normal(p, aNormalVector); <<
1333 // aNormalVector = facet.GetSurfaceNo <<
1334 aNormalVector = normal; <<
1335 return 0.0; <<
1336 } <<
1337 if (dist >= 0.0 && dist < minDist) <<
1338 { <<
1339 minDist = dist; <<
1340 minNormal = normal; <<
1341 isExtreme = (fExtremeFacets.find(&fac <<
1342 } <<
1343 } <<
1344 } <<
1345 if (minDist < kInfinity) <<
1346 { <<
1347 aNormalVector = minNormal; <<
1348 aConvex = isExtreme; <<
1349 return minDist; <<
1350 } <<
1351 else <<
1352 { <<
1353 // No intersection found <<
1354 aConvex = false; <<
1355 Normal(p, aNormalVector); <<
1356 return 0.0; <<
1357 } 714 }
>> 715
>> 716 return minDist;
1358 } 717 }
1359 718
1360 ///////////////////////////////////////////// 719 ///////////////////////////////////////////////////////////////////////////////
1361 // 720 //
1362 void G4TessellatedSolid:: << 721 // G4double DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
1363 DistanceToOutCandidates(const std::vector<G4i << 722 // const G4bool calcNorm=false,
1364 const G4ThreeVector& << 723 // G4bool *validNorm=0, G4ThreeVector *n=0);
1365 const G4ThreeVector& << 724 //
1366 G4double& minDi << 725 // Return distance along the normalised vector v to the shape, from a
1367 G4int& minCandi << 726 // point at an offset p inside or on the surface of the
>> 727 // shape. Intersections with surfaces, when the point is not greater
>> 728 // than kCarTolerance/2 from a surface, must be ignored.
>> 729 // If calcNorm is true, then it must also set validNorm to either
>> 730 // * true, if the solid lies entirely behind or on the exiting
>> 731 // surface. Then it must set n to the outwards normal vector
>> 732 // (the Magnitude of the vector is not defined).
>> 733 // * false, if the solid does not lie entirely behind or on the
>> 734 // exiting surface.
>> 735 // If calcNorm is false, then validNorm and n are unused.
>> 736
>> 737 G4double G4TessellatedSolid::DistanceToOut (const G4ThreeVector &p,
>> 738 const G4ThreeVector &v, const G4bool calcNorm,
>> 739 G4bool *validNorm, G4ThreeVector *n) const
1368 { 740 {
1369 auto candidatesCount = (G4int)candidates.si << 741 G4double minDist = kInfinity;
1370 G4double dist = 0.0; 742 G4double dist = 0.0;
1371 G4double distFromSurface = 0.0; 743 G4double distFromSurface = 0.0;
1372 G4ThreeVector normal; << 744 G4ThreeVector normal(0.0,0.0,0.0);
>> 745 G4ThreeVector minNormal(0.0,0.0,0.0);
>> 746
>> 747 #if G4SPECSDEBUG
>> 748 if ( Inside(p) == kOutside )
>> 749 {
>> 750 G4cout.precision(16) ;
>> 751 G4cout << G4endl ;
>> 752 // DumpInfo();
>> 753 G4cout << "Position:" << G4endl << G4endl ;
>> 754 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
>> 755 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
>> 756 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
>> 757 G4cout << "DistanceToIn(p) == " << DistanceToIn(p) << G4endl;
>> 758 G4Exception("G4TriangularFacet::DistanceToOut(p)", "Notification", JustWarning,
>> 759 "Point p is already outside !?" );
>> 760 }
>> 761 #endif
1373 762
1374 for (G4int i = 0 ; i < candidatesCount; ++i << 763 G4bool isExtreme = false;
>> 764 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
1375 { 765 {
1376 G4int candidate = candidates[i]; << 766 if ((*f)->Intersect(p,v,true,dist,distFromSurface,normal))
1377 G4VFacet& facet = *fFacets[candidate]; << 767 {
1378 if (facet.Intersect(aPoint,direction,true << 768 if (distFromSurface > 0.0 && distFromSurface <= 0.5*kCarTolerance &&
1379 { << 769 (*f)->Distance(p,kCarTolerance) <= 0.5*kCarTolerance)
1380 if (distFromSurface > 0.0 && distFromSu <<
1381 && facet.Distance(aPoint,kCarTolerance <<
1382 { 770 {
1383 // We are on a surface << 771 // We are on a surface. Return zero.
1384 // << 772 if (calcNorm) {
1385 minDist = 0.0; << 773 *validNorm = extremeFacets.count(*f);
1386 minNormal = normal; << 774 *n = SurfaceNormal(p);
1387 minCandidate = candidate; << 775 }
1388 break; << 776 return 0.0;
1389 } 777 }
1390 if (dist >= 0.0 && dist < minDist) 778 if (dist >= 0.0 && dist < minDist)
1391 { 779 {
1392 minDist = dist; << 780 minDist = dist;
1393 minNormal = normal; 781 minNormal = normal;
1394 minCandidate = candidate; << 782 isExtreme = extremeFacets.count(*f);
1395 } 783 }
1396 } 784 }
1397 } 785 }
1398 } << 786
1399 << 787 if (minDist < kInfinity)
1400 ///////////////////////////////////////////// <<
1401 // <<
1402 G4double <<
1403 G4TessellatedSolid::DistanceToOutCore(const G <<
1404 const G <<
1405 G <<
1406 G <<
1407 G <<
1408 { <<
1409 G4double minDistance; <<
1410 <<
1411 if (fVoxels.GetCountOfVoxels() > 1) <<
1412 { 788 {
1413 minDistance = kInfinity; << 789 if (calcNorm)
1414 <<
1415 G4ThreeVector currentPoint = aPoint; <<
1416 G4ThreeVector direction = aDirection.unit <<
1417 G4double totalShift = 0.; <<
1418 vector<G4int> curVoxel(3); <<
1419 if (!fVoxels.Contains(aPoint)) return 0.; <<
1420 <<
1421 fVoxels.GetVoxel(curVoxel, currentPoint); <<
1422 <<
1423 G4double shiftBonus = kCarTolerance; <<
1424 <<
1425 const vector<G4int>* old = nullptr; <<
1426 <<
1427 G4int minCandidate = -1; <<
1428 do // Loop checking, 13.08.2015, G.Cos <<
1429 { <<
1430 const vector<G4int>& candidates = fVoxe <<
1431 if (old == &candidates) <<
1432 ++old; <<
1433 if (old != &candidates && !candidates.e <<
1434 { <<
1435 DistanceToOutCandidates(candidates, a <<
1436 aNormalVector <<
1437 if (minDistance <= totalShift) break; <<
1438 } <<
1439 <<
1440 G4double shift=fVoxels.DistanceToNext(c <<
1441 if (shift == kInfinity) break; <<
1442 <<
1443 totalShift += shift; <<
1444 if (minDistance <= totalShift) break; <<
1445 <<
1446 currentPoint += direction * (shift + sh <<
1447 <<
1448 old = &candidates; <<
1449 } <<
1450 while (fVoxels.UpdateCurrentVoxel(current <<
1451 <<
1452 if (minCandidate < 0) <<
1453 { <<
1454 // No intersection found <<
1455 minDistance = 0.; <<
1456 aConvex = false; <<
1457 Normal(aPoint, aNormalVector); <<
1458 } <<
1459 else <<
1460 { 790 {
1461 aConvex = (fExtremeFacets.find(fFacets[ << 791 *validNorm = isExtreme;
1462 != fExtremeFacets.end()); << 792 *n = minNormal;
1463 } 793 }
>> 794 return minDist;
1464 } 795 }
1465 else 796 else
1466 { 797 {
1467 minDistance = DistanceToOutNoVoxels(aPoin << 798 // No intersection found
1468 aConv << 799 if (calcNorm)
1469 } <<
1470 return minDistance; <<
1471 } <<
1472 <<
1473 ///////////////////////////////////////////// <<
1474 // <<
1475 G4double G4TessellatedSolid:: <<
1476 DistanceToInCandidates(const std::vector<G4in <<
1477 const G4ThreeVector& a <<
1478 const G4ThreeVector& d <<
1479 { <<
1480 auto candidatesCount = (G4int)candidates.si <<
1481 G4double dist = 0.0; <<
1482 G4double distFromSurface = 0.0; <<
1483 G4ThreeVector normal; <<
1484 <<
1485 G4double minDistance = kInfinity; <<
1486 for (G4int i = 0 ; i < candidatesCount; ++i <<
1487 { <<
1488 G4int candidate = candidates[i]; <<
1489 G4VFacet& facet = *fFacets[candidate]; <<
1490 if (facet.Intersect(aPoint,direction,fals <<
1491 { 800 {
1492 // << 801 *validNorm = false;
1493 // Set minDist to the new distance to c << 802 *n = SurfaceNormal(p);
1494 // in positive direction and point is n <<
1495 // within 0.5*kCarTolerance of the surf <<
1496 // zero and leave member function immed <<
1497 // proposed by & credit to Akira Okumur <<
1498 // <<
1499 if ( (distFromSurface > kCarToleranceHa <<
1500 && (dist >= 0.0) && (dist < minDistan <<
1501 { <<
1502 minDistance = dist; <<
1503 } <<
1504 else <<
1505 { <<
1506 if (-kCarToleranceHalf <= dist && dis <<
1507 { <<
1508 return 0.0; <<
1509 } <<
1510 else if (distFromSurface > -kCarTole <<
1511 && distFromSurface < kCarTole <<
1512 { <<
1513 minDistance = dist; <<
1514 } <<
1515 } <<
1516 } 803 }
>> 804 return 0.0;
1517 } 805 }
1518 return minDistance; <<
1519 } 806 }
1520 807
1521 ///////////////////////////////////////////// 808 ///////////////////////////////////////////////////////////////////////////////
1522 // 809 //
1523 G4double << 810 // G4double DistanceToOut(const G4ThreeVector& p)
1524 G4TessellatedSolid::DistanceToInCore(const G4 <<
1525 const G4 <<
1526 G4 <<
1527 { <<
1528 G4double minDistance; <<
1529 <<
1530 if (fVoxels.GetCountOfVoxels() > 1) <<
1531 { <<
1532 minDistance = kInfinity; <<
1533 G4ThreeVector currentPoint = aPoint; <<
1534 G4ThreeVector direction = aDirection.unit <<
1535 G4double shift = fVoxels.DistanceToFirst( <<
1536 if (shift == kInfinity) return shift; <<
1537 G4double shiftBonus = kCarTolerance; <<
1538 if (shift != 0.0) <<
1539 currentPoint += direction * (shift + sh <<
1540 // if (!fVoxels.Contains(currentPoint)) <<
1541 G4double totalShift = shift; <<
1542 <<
1543 // G4SurfBits exclusion; // (1/*fVoxels.G <<
1544 vector<G4int> curVoxel(3); <<
1545 <<
1546 fVoxels.GetVoxel(curVoxel, currentPoint); <<
1547 do // Loop checking, 13.08.2015, G.Cos <<
1548 { <<
1549 const vector<G4int>& candidates = fVoxe <<
1550 if (!candidates.empty()) <<
1551 { <<
1552 G4double distance=DistanceToInCandida <<
1553 if (minDistance > distance) minDistan <<
1554 if (distance < totalShift) break; <<
1555 } <<
1556 <<
1557 shift = fVoxels.DistanceToNext(currentP <<
1558 if (shift == kInfinity /*|| shift == 0* <<
1559 <<
1560 totalShift += shift; <<
1561 if (minDistance < totalShift) break; <<
1562 <<
1563 currentPoint += direction * (shift + sh <<
1564 } <<
1565 while (fVoxels.UpdateCurrentVoxel(current <<
1566 } <<
1567 else <<
1568 { <<
1569 minDistance = DistanceToInNoVoxels(aPoint <<
1570 } <<
1571 <<
1572 return minDistance; <<
1573 } <<
1574 <<
1575 ///////////////////////////////////////////// <<
1576 // 811 //
1577 G4bool << 812 // Calculate distance to nearest surface of shape from an inside
1578 G4TessellatedSolid::CompareSortedVoxel(const << 813 // point. The distance can be an underestimate.
1579 const <<
1580 { <<
1581 return l.second < r.second; <<
1582 } <<
1583 814
1584 ///////////////////////////////////////////// << 815 G4double G4TessellatedSolid::DistanceToOut (const G4ThreeVector &p) const
1585 // <<
1586 G4double <<
1587 G4TessellatedSolid::MinDistanceFacet(const G4 <<
1588 G4 <<
1589 G4 <<
1590 { 816 {
1591 G4double minDist = kInfinity; 817 G4double minDist = kInfinity;
1592 << 818 G4double dist = 0.0;
1593 G4int size = fVoxels.GetVoxelBoxesSize(); << 819
1594 vector<pair<G4int, G4double> > voxelsSorted <<
1595 <<
1596 pair<G4int, G4double> info; <<
1597 <<
1598 for (G4int i = 0; i < size; ++i) <<
1599 { <<
1600 const G4VoxelBox& voxelBox = fVoxels.GetV <<
1601 <<
1602 G4ThreeVector pointShifted = p - voxelBox <<
1603 G4double safety = fVoxels.MinDistanceToBo <<
1604 info.first = i; <<
1605 info.second = safety; <<
1606 voxelsSorted[i] = info; <<
1607 } <<
1608 <<
1609 std::sort(voxelsSorted.begin(), voxelsSorte <<
1610 &G4TessellatedSolid::CompareSorte <<
1611 <<
1612 for (G4int i = 0; i < size; ++i) <<
1613 { <<
1614 const pair<G4int,G4double>& inf = voxelsS <<
1615 G4double dist = inf.second; <<
1616 if (dist > minDist) break; <<
1617 <<
1618 const vector<G4int>& candidates = fVoxels <<
1619 auto csize = (G4int)candidates.size(); <<
1620 for (G4int j = 0; j < csize; ++j) <<
1621 { <<
1622 G4int candidate = candidates[j]; <<
1623 G4VFacet& facet = *fFacets[candidate]; <<
1624 dist = simple ? facet.Distance(p,minDis <<
1625 : facet.Distance(p,minDis <<
1626 if (dist < minDist) <<
1627 { <<
1628 minDist = dist; <<
1629 minFacet = &facet; <<
1630 } <<
1631 } <<
1632 } <<
1633 return minDist; <<
1634 } <<
1635 <<
1636 ///////////////////////////////////////////// <<
1637 // <<
1638 G4double G4TessellatedSolid::SafetyFromOutsid <<
1639 <<
1640 { <<
1641 #if G4SPECSDEBUG <<
1642 if ( Inside(p) == kInside ) <<
1643 { <<
1644 std::ostringstream message; <<
1645 G4int oldprc = message.precision(16) ; <<
1646 message << "Point p is already inside!?" <<
1647 << "Position:" << G4endl << G4endl <<
1648 << "p.x() = " << p.x()/mm << " mm" << <<
1649 << "p.y() = " << p.y()/mm << " mm" << <<
1650 << "p.z() = " << p.z()/mm << " mm" << <<
1651 << "DistanceToOut(p) == " << DistanceTo <<
1652 message.precision(oldprc) ; <<
1653 G4Exception("G4TriangularFacet::DistanceT <<
1654 "GeomSolids1002", JustWarning <<
1655 } <<
1656 #endif <<
1657 <<
1658 G4double minDist; <<
1659 <<
1660 if (fVoxels.GetCountOfVoxels() > 1) <<
1661 { <<
1662 if (!aAccurate) <<
1663 return fVoxels.DistanceToBoundingBox(p) <<
1664 <<
1665 if (!OutsideOfExtent(p, kCarTolerance)) <<
1666 { <<
1667 vector<G4int> startingVoxel(3); <<
1668 fVoxels.GetVoxel(startingVoxel, p); <<
1669 const vector<G4int> &candidates = fVoxe <<
1670 if (candidates.empty() && (fInsides.Get <<
1671 { <<
1672 G4int index = fVoxels.GetPointIndex(p <<
1673 if (fInsides[index]) return 0.; <<
1674 } <<
1675 } <<
1676 <<
1677 G4VFacet* facet; <<
1678 minDist = MinDistanceFacet(p, true, facet <<
1679 } <<
1680 else <<
1681 { <<
1682 minDist = kInfinity; <<
1683 std::size_t size = fFacets.size(); <<
1684 for (std::size_t i = 0; i < size; ++i) <<
1685 { <<
1686 G4VFacet& facet = *fFacets[i]; <<
1687 G4double dist = facet.Distance(p,minDis <<
1688 if (dist < minDist) minDist = dist; <<
1689 } <<
1690 } <<
1691 return minDist; <<
1692 } <<
1693 <<
1694 ///////////////////////////////////////////// <<
1695 // <<
1696 G4double <<
1697 G4TessellatedSolid::SafetyFromInside (const G <<
1698 { <<
1699 #if G4SPECSDEBUG 820 #if G4SPECSDEBUG
1700 if ( Inside(p) == kOutside ) 821 if ( Inside(p) == kOutside )
1701 { 822 {
1702 std::ostringstream message; << 823 G4cout.precision(16) ;
1703 G4int oldprc = message.precision(16) ; << 824 G4cout << G4endl ;
1704 message << "Point p is already outside!?" << 825 // DumpInfo();
1705 << "Position:" << G4endl << G4endl << 826 G4cout << "Position:" << G4endl << G4endl ;
1706 << "p.x() = " << p.x()/mm << " mm" << << 827 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
1707 << "p.y() = " << p.y()/mm << " mm" << << 828 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
1708 << "p.z() = " << p.z()/mm << " mm" << << 829 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
1709 << "DistanceToIn(p) == " << DistanceToI << 830 G4cout << "DistanceToIn(p) == " << DistanceToIn(p) << G4endl;
1710 message.precision(oldprc) ; << 831 G4Exception("G4TriangularFacet::DistanceToOut(p)", "Notification", JustWarning,
1711 G4Exception("G4TriangularFacet::DistanceT << 832 "Point p is already outside !?" );
1712 "GeomSolids1002", JustWarning <<
1713 } 833 }
1714 #endif 834 #endif
1715 835
1716 G4double minDist; << 836 for (FacetCI f=facets.begin(); f!=facets.end(); f++)
1717 <<
1718 if (OutsideOfExtent(p, kCarTolerance)) retu <<
1719 <<
1720 if (fVoxels.GetCountOfVoxels() > 1) <<
1721 { <<
1722 G4VFacet* facet; <<
1723 minDist = MinDistanceFacet(p, true, facet <<
1724 } <<
1725 else <<
1726 { 837 {
1727 minDist = kInfinity; << 838 dist = (*f)->Distance(p,minDist,true);
1728 G4double dist = 0.0; << 839 if (dist < minDist) minDist = dist;
1729 std::size_t size = fFacets.size(); <<
1730 for (std::size_t i = 0; i < size; ++i) <<
1731 { <<
1732 G4VFacet& facet = *fFacets[i]; <<
1733 dist = facet.Distance(p,minDist); <<
1734 if (dist < minDist) minDist = dist; <<
1735 } <<
1736 } 840 }
>> 841
1737 return minDist; 842 return minDist;
1738 } 843 }
1739 844
1740 ///////////////////////////////////////////// 845 ///////////////////////////////////////////////////////////////////////////////
1741 // 846 //
1742 // G4GeometryType GetEntityType() const; 847 // G4GeometryType GetEntityType() const;
1743 // 848 //
1744 // Provide identification of the class of an << 849 // Provide identification of the class of an object (required for
>> 850 // persistency and STEP interface).
1745 // 851 //
1746 G4GeometryType G4TessellatedSolid::GetEntityT 852 G4GeometryType G4TessellatedSolid::GetEntityType () const
1747 { 853 {
1748 return fGeometryType; << 854 return geometryType;
1749 } 855 }
1750 856
1751 ///////////////////////////////////////////// 857 ///////////////////////////////////////////////////////////////////////////////
1752 // 858 //
1753 // IsFaceted << 859 void G4TessellatedSolid::DescribeYourselfTo (G4VGraphicsScene& scene) const
1754 // <<
1755 G4bool G4TessellatedSolid::IsFaceted () const <<
1756 { <<
1757 return true; <<
1758 } <<
1759 <<
1760 ///////////////////////////////////////////// <<
1761 // <<
1762 std::ostream &G4TessellatedSolid::StreamInfo( <<
1763 { 860 {
1764 os << G4endl; << 861 scene.AddSolid (*this);
1765 os << "Solid name = " << GetName() <<
1766 os << "Geometry Type = " << fGeometryTyp <<
1767 os << "Number of facets = " << fFacets.size <<
1768 <<
1769 std::size_t size = fFacets.size(); <<
1770 for (std::size_t i = 0; i < size; ++i) <<
1771 { <<
1772 os << "FACET # = " << i + 1 << G <<
1773 G4VFacet &facet = *fFacets[i]; <<
1774 facet.StreamInfo(os); <<
1775 } <<
1776 os << G4endl; <<
1777 <<
1778 return os; <<
1779 } 862 }
1780 863
1781 ///////////////////////////////////////////// 864 ///////////////////////////////////////////////////////////////////////////////
1782 // 865 //
1783 // Make a clone of the object << 866 // Dispatch to parameterisation for replication mechanism dimension
1784 // << 867 // computation & modification.
1785 G4VSolid* G4TessellatedSolid::Clone() const << 868 //
1786 { << 869 //void G4TessellatedSolid::ComputeDimensions (G4VPVParameterisation* p,
1787 return new G4TessellatedSolid(*this); << 870 // const G4int n, const G4VPhysicalVolume* pRep) const
1788 } << 871 //{
>> 872 // G4VSolid *ptr = 0;
>> 873 // ptr = *this;
>> 874 // p->ComputeDimensions(ptr,n,pRep);
>> 875 //}
1789 876
1790 ///////////////////////////////////////////// 877 ///////////////////////////////////////////////////////////////////////////////
1791 // 878 //
1792 // EInside G4TessellatedSolid::Inside (const << 879 std::ostream &G4TessellatedSolid::StreamInfo(std::ostream &os) const
1793 // <<
1794 // This method must return: <<
1795 // * kOutside if the point at offset p is <<
1796 // boundaries plus kCarTolerance/2, <<
1797 // * kSurface if the point is <= kCarToler <<
1798 // * kInside otherwise. <<
1799 // <<
1800 EInside G4TessellatedSolid::Inside (const G4T <<
1801 { 880 {
1802 EInside location; << 881 os << G4endl;
1803 << 882 os << "Geometry Type = " << geometryType << G4endl;
1804 if (fVoxels.GetCountOfVoxels() > 1) << 883 os << "Number of facets = " << facets.size() << G4endl;
1805 { << 884
1806 location = InsideVoxels(aPoint); << 885 for (FacetCI f = facets.begin(); f != facets.end(); f++)
1807 } <<
1808 else <<
1809 { 886 {
1810 location = InsideNoVoxels(aPoint); << 887 os << "FACET # = " << f-facets.begin()+1 << G4endl;
>> 888 (*f)->StreamInfo(os);
1811 } 889 }
1812 return location; << 890 os <<G4endl;
1813 } << 891
1814 << 892 return os;
1815 ///////////////////////////////////////////// <<
1816 // <<
1817 G4ThreeVector G4TessellatedSolid::SurfaceNorm <<
1818 { <<
1819 G4ThreeVector n; <<
1820 Normal(p, n); <<
1821 return n; <<
1822 } <<
1823 <<
1824 ///////////////////////////////////////////// <<
1825 // <<
1826 // G4double DistanceToIn(const G4ThreeVector& <<
1827 // <<
1828 // Calculate distance to nearest surface of s <<
1829 // distance can be an underestimate. <<
1830 // <<
1831 G4double G4TessellatedSolid::DistanceToIn(con <<
1832 { <<
1833 return SafetyFromOutside(p, false); <<
1834 } 893 }
1835 894
1836 ///////////////////////////////////////////// 895 ///////////////////////////////////////////////////////////////////////////////
1837 // 896 //
1838 G4double G4TessellatedSolid::DistanceToIn(con << 897 G4Polyhedron *G4TessellatedSolid::CreatePolyhedron () const
1839 con <<
1840 { 898 {
1841 G4double dist = DistanceToInCore(p,v,kInfin << 899 size_t nVertices = vertexList.size();
1842 #ifdef G4SPECSDEBUG << 900 size_t nFacets = facets.size();
1843 if (dist < kInfinity) << 901 G4PolyhedronArbitrary *polyhedron =
>> 902 new G4PolyhedronArbitrary (nVertices, nFacets);
>> 903 for (G4ThreeVectorList::const_iterator v = vertexList.begin();
>> 904 v!=vertexList.end(); v++) polyhedron->AddVertex(*v);
>> 905
>> 906 for (FacetCI f=facets.begin(); f != facets.end(); f++)
1844 { 907 {
1845 if (Inside(p + dist*v) != kSurface) << 908 size_t v[4];
>> 909 for (size_t j=0; j<4; j++)
1846 { 910 {
1847 std::ostringstream message; << 911 size_t i = (*f)->GetVertexIndex(j);
1848 message << "Invalid response from facet << 912 if (i == 999999999) v[j] = 0;
1849 << G4endl << 913 else v[j] = i+1;
1850 << "at point: " << p << "and d <<
1851 G4Exception("G4TessellatedSolid::Distan <<
1852 "GeomSolids1002", JustWarni <<
1853 } 914 }
1854 } << 915 if ((*f)->GetEntityType() == "G4RectangularFacet")
1855 #endif <<
1856 return dist; <<
1857 } <<
1858 <<
1859 ///////////////////////////////////////////// <<
1860 // <<
1861 // G4double DistanceToOut(const G4ThreeVector <<
1862 // <<
1863 // Calculate distance to nearest surface of s <<
1864 // point. The distance can be an underestimat <<
1865 // <<
1866 G4double G4TessellatedSolid::DistanceToOut(co <<
1867 { <<
1868 return SafetyFromInside(p, false); <<
1869 } <<
1870 <<
1871 ///////////////////////////////////////////// <<
1872 // <<
1873 // G4double DistanceToOut(const G4ThreeVector <<
1874 // const G4bool calcNo <<
1875 // G4bool *validNorm=0 <<
1876 // <<
1877 // Return distance along the normalised vecto <<
1878 // point at an offset p inside or on the surf <<
1879 // shape. Intersections with surfaces, when t <<
1880 // than kCarTolerance/2 from a surface, must <<
1881 // If calcNorm is true, then it must also <<
1882 // * true, if the solid lies entirely beh <<
1883 // surface. Then it must set n to the <<
1884 // (the Magnitude of the vector is not <<
1885 // * false, if the solid does not lie ent <<
1886 // exiting surface. <<
1887 // If calcNorm is false, then validNorm and n <<
1888 // <<
1889 G4double G4TessellatedSolid::DistanceToOut(co <<
1890 co <<
1891 co <<
1892 <<
1893 <<
1894 { <<
1895 G4ThreeVector n; <<
1896 G4bool valid; <<
1897 <<
1898 G4double dist = DistanceToOutCore(p, v, n, <<
1899 if (calcNorm) <<
1900 { <<
1901 *norm = n; <<
1902 *validNorm = valid; <<
1903 } <<
1904 #ifdef G4SPECSDEBUG <<
1905 if (dist < kInfinity) <<
1906 { <<
1907 if (Inside(p + dist*v) != kSurface) <<
1908 { 916 {
1909 std::ostringstream message; << 917 size_t i = v[3];
1910 message << "Invalid response from facet << 918 v[3] = v[2];
1911 << G4endl << 919 v[2] = i;
1912 << "at point: " << p << "and d <<
1913 G4Exception("G4TessellatedSolid::Distan <<
1914 "GeomSolids1002", JustWarni <<
1915 } 920 }
>> 921 polyhedron->AddFacet(v[0],v[1],v[2],v[3]);
1916 } 922 }
1917 #endif << 923
1918 return dist; << 924 return (G4Polyhedron*) polyhedron;
1919 } <<
1920 <<
1921 ///////////////////////////////////////////// <<
1922 // <<
1923 void G4TessellatedSolid::DescribeYourselfTo ( <<
1924 { <<
1925 scene.AddSolid (*this); <<
1926 } 925 }
1927 926
1928 ///////////////////////////////////////////// 927 ///////////////////////////////////////////////////////////////////////////////
1929 // 928 //
1930 G4Polyhedron* G4TessellatedSolid::CreatePolyh << 929 G4NURBS *G4TessellatedSolid::CreateNURBS () const
1931 { 930 {
1932 auto nVertices = (G4int)fVertexList.size(); << 931 return 0;
1933 auto nFacets = (G4int)fFacets.size(); <<
1934 auto polyhedron = new G4Polyhedron(nVertice <<
1935 for (auto i = 0; i < nVertices; ++i) <<
1936 { <<
1937 polyhedron->SetVertex(i+1, fVertexList[i] <<
1938 } <<
1939 <<
1940 for (auto i = 0; i < nFacets; ++i) <<
1941 { <<
1942 G4VFacet* facet = fFacets[i]; <<
1943 G4int v[4] = {0}; <<
1944 G4int n = facet->GetNumberOfVertices(); <<
1945 if (n > 4) n = 4; <<
1946 for (auto j = 0; j < n; ++j) <<
1947 { <<
1948 v[j] = facet->GetVertexIndex(j) + 1; <<
1949 } <<
1950 polyhedron->SetFacet(i+1, v[0], v[1], v[2 <<
1951 } <<
1952 polyhedron->SetReferences(); <<
1953 <<
1954 return polyhedron; <<
1955 } 932 }
1956 933
1957 ///////////////////////////////////////////// 934 ///////////////////////////////////////////////////////////////////////////////
1958 // 935 //
1959 // GetPolyhedron 936 // GetPolyhedron
1960 // 937 //
1961 G4Polyhedron* G4TessellatedSolid::GetPolyhedr << 938 G4Polyhedron* G4TessellatedSolid::GetPolyhedron () const
1962 { 939 {
1963 if (fpPolyhedron == nullptr || << 940 if (!fpPolyhedron ||
1964 fRebuildPolyhedron || <<
1965 fpPolyhedron->GetNumberOfRotationStepsA 941 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
1966 fpPolyhedron->GetNumberOfRotationSteps( 942 fpPolyhedron->GetNumberOfRotationSteps())
1967 { << 943 {
1968 G4AutoLock l(&polyhedronMutex); << 944 delete fpPolyhedron;
1969 delete fpPolyhedron; << 945 fpPolyhedron = CreatePolyhedron();
1970 fpPolyhedron = CreatePolyhedron(); << 946 }
1971 fRebuildPolyhedron = false; <<
1972 l.unlock(); <<
1973 } <<
1974 return fpPolyhedron; 947 return fpPolyhedron;
1975 } 948 }
1976 949
1977 ///////////////////////////////////////////// 950 ///////////////////////////////////////////////////////////////////////////////
1978 // 951 //
1979 // Get bounding box << 952 // CalculateExtent
1980 // 953 //
1981 void G4TessellatedSolid::BoundingLimits(G4Thr << 954 // Based on correction provided by Stan Seibert, University of Texas.
1982 G4Thr <<
1983 { <<
1984 pMin = fMinExtent; <<
1985 pMax = fMaxExtent; <<
1986 <<
1987 // Check correctness of the bounding box <<
1988 // <<
1989 if (pMin.x() >= pMax.x() || pMin.y() >= pMa <<
1990 { <<
1991 std::ostringstream message; <<
1992 message << "Bad bounding box (min >= max) <<
1993 << GetName() << " !" <<
1994 << "\npMin = " << pMin <<
1995 << "\npMax = " << pMax; <<
1996 G4Exception("G4TessellatedSolid::Bounding <<
1997 "GeomMgt0001", JustWarning, m <<
1998 DumpInfo(); <<
1999 } <<
2000 } <<
2001 <<
2002 ///////////////////////////////////////////// <<
2003 // <<
2004 // Calculate extent under transform and speci <<
2005 // 955 //
2006 G4bool 956 G4bool
2007 G4TessellatedSolid::CalculateExtent(const EAx 957 G4TessellatedSolid::CalculateExtent(const EAxis pAxis,
2008 const G4V 958 const G4VoxelLimits& pVoxelLimit,
2009 const G4A 959 const G4AffineTransform& pTransform,
2010 G4d 960 G4double& pMin, G4double& pMax) const
2011 { 961 {
2012 G4ThreeVector bmin, bmax; << 962 G4ThreeVectorList transVertexList(vertexList);
2013 963
2014 // Check bounding box (bbox) << 964 // Put solid into transformed frame
2015 // << 965 for (size_t i=0; i<vertexList.size(); i++)
2016 BoundingLimits(bmin,bmax); << 966 { pTransform.ApplyPointTransform(transVertexList[i]); }
2017 G4BoundingEnvelope bbox(bmin,bmax); << 967
2018 << 968 // Find min and max extent in each dimension
2019 // Use simple bounding-box to help in the c << 969 G4ThreeVector minExtent(kInfinity, kInfinity, kInfinity);
2020 // << 970 G4ThreeVector maxExtent(-kInfinity, -kInfinity, -kInfinity);
2021 return bbox.CalculateExtent(pAxis,pVoxelLim << 971 for (size_t i=0; i<transVertexList.size(); i++)
2022 << 972 {
2023 #if 0 << 973 for (G4int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; axis++)
2024 // Precise extent computation (disabled by << 974 {
2025 // << 975 G4double coordinate = transVertexList[i][axis];
2026 if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVo << 976 if (coordinate < minExtent[axis])
2027 { << 977 { minExtent[axis] = coordinate; }
2028 return (pMin < pMax) ? true : false; << 978 if (coordinate > maxExtent[axis])
2029 } << 979 { maxExtent[axis] = coordinate; }
2030 << 980 }
2031 // The extent is calculated as cumulative e << 981 }
2032 // formed by facets and the center of the b << 982
2033 // << 983 // Check for containment and clamp to voxel boundaries
2034 G4double eminlim = pVoxelLimit.GetMinExtent << 984 for (G4int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; axis++)
2035 G4double emaxlim = pVoxelLimit.GetMaxExtent << 985 {
2036 << 986 EAxis geomAxis = kXAxis; // G4 geom classes use different index type
2037 G4ThreeVectorList base; << 987 switch(axis)
2038 G4ThreeVectorList apex(1); << 988 {
2039 std::vector<const G4ThreeVectorList *> pyra << 989 case G4ThreeVector::X: geomAxis = kXAxis; break;
2040 pyramid[0] = &base; << 990 case G4ThreeVector::Y: geomAxis = kYAxis; break;
2041 pyramid[1] = &apex; << 991 case G4ThreeVector::Z: geomAxis = kZAxis; break;
2042 apex[0] = (bmin+bmax)*0.5; << 992 }
2043 << 993 G4bool isLimited = pVoxelLimit.IsLimited(geomAxis);
2044 // main loop along facets << 994 G4double voxelMinExtent = pVoxelLimit.GetMinExtent(geomAxis);
2045 pMin = kInfinity; << 995 G4double voxelMaxExtent = pVoxelLimit.GetMaxExtent(geomAxis);
2046 pMax = -kInfinity; << 996
2047 for (G4int i=0; i<GetNumberOfFacets(); ++i) << 997 if (isLimited)
2048 { << 998 {
2049 G4VFacet* facet = GetFacet(i); << 999 if ( minExtent[axis] > voxelMaxExtent+kCarTolerance ||
2050 if (std::abs((facet->GetSurfaceNormal()). << 1000 maxExtent[axis] < voxelMinExtent-kCarTolerance )
2051 < kCarToleranceHalf) continue; << 1001 {
2052 << 1002 return false ;
2053 G4int nv = facet->GetNumberOfVertices(); << 1003 }
2054 base.resize(nv); << 1004 else
2055 for (G4int k=0; k<nv; ++k) { base[k] = fa << 1005 {
2056 << 1006 if (minExtent[axis] < voxelMinExtent)
2057 G4double emin,emax; << 1007 {
2058 G4BoundingEnvelope benv(pyramid); << 1008 minExtent[axis] = voxelMinExtent ;
2059 if (!benv.CalculateExtent(pAxis,pVoxelLim << 1009 }
2060 if (emin < pMin) pMin = emin; << 1010 if (maxExtent[axis] > voxelMaxExtent)
2061 if (emax > pMax) pMax = emax; << 1011 {
2062 if (eminlim > pMin && emaxlim < pMax) bre << 1012 maxExtent[axis] = voxelMaxExtent;
2063 } << 1013 }
2064 return (pMin < pMax); << 1014 }
2065 #endif << 1015 }
>> 1016 }
>> 1017
>> 1018 // Convert pAxis into G4ThreeVector index
>> 1019 G4int vecAxis=0;
>> 1020 switch(pAxis)
>> 1021 {
>> 1022 case kXAxis: vecAxis = G4ThreeVector::X; break;
>> 1023 case kYAxis: vecAxis = G4ThreeVector::Y; break;
>> 1024 case kZAxis: vecAxis = G4ThreeVector::Z; break;
>> 1025 default: break;
>> 1026 }
>> 1027
>> 1028 pMin = minExtent[vecAxis] - kCarTolerance;
>> 1029 pMax = maxExtent[vecAxis] + kCarTolerance;
>> 1030
>> 1031 return true;
2066 } 1032 }
2067 1033
2068 ///////////////////////////////////////////// 1034 ///////////////////////////////////////////////////////////////////////////////
2069 // 1035 //
2070 G4double G4TessellatedSolid::GetMinXExtent () 1036 G4double G4TessellatedSolid::GetMinXExtent () const
2071 { << 1037 {return xMinExtent;}
2072 return fMinExtent.x(); <<
2073 } <<
2074 1038
2075 ///////////////////////////////////////////// 1039 ///////////////////////////////////////////////////////////////////////////////
2076 // 1040 //
2077 G4double G4TessellatedSolid::GetMaxXExtent () 1041 G4double G4TessellatedSolid::GetMaxXExtent () const
2078 { << 1042 {return xMaxExtent;}
2079 return fMaxExtent.x(); <<
2080 } <<
2081 1043
2082 ///////////////////////////////////////////// 1044 ///////////////////////////////////////////////////////////////////////////////
2083 // 1045 //
2084 G4double G4TessellatedSolid::GetMinYExtent () 1046 G4double G4TessellatedSolid::GetMinYExtent () const
2085 { << 1047 {return yMinExtent;}
2086 return fMinExtent.y(); <<
2087 } <<
2088 1048
2089 ///////////////////////////////////////////// 1049 ///////////////////////////////////////////////////////////////////////////////
2090 // 1050 //
2091 G4double G4TessellatedSolid::GetMaxYExtent () 1051 G4double G4TessellatedSolid::GetMaxYExtent () const
2092 { << 1052 {return yMaxExtent;}
2093 return fMaxExtent.y(); <<
2094 } <<
2095 1053
2096 ///////////////////////////////////////////// 1054 ///////////////////////////////////////////////////////////////////////////////
2097 // 1055 //
2098 G4double G4TessellatedSolid::GetMinZExtent () 1056 G4double G4TessellatedSolid::GetMinZExtent () const
2099 { << 1057 {return zMinExtent;}
2100 return fMinExtent.z(); <<
2101 } <<
2102 1058
2103 ///////////////////////////////////////////// 1059 ///////////////////////////////////////////////////////////////////////////////
2104 // 1060 //
2105 G4double G4TessellatedSolid::GetMaxZExtent () 1061 G4double G4TessellatedSolid::GetMaxZExtent () const
2106 { << 1062 {return zMaxExtent;}
2107 return fMaxExtent.z(); <<
2108 } <<
2109 1063
2110 ///////////////////////////////////////////// 1064 ///////////////////////////////////////////////////////////////////////////////
2111 // 1065 //
2112 G4VisExtent G4TessellatedSolid::GetExtent () 1066 G4VisExtent G4TessellatedSolid::GetExtent () const
2113 { 1067 {
2114 return { fMinExtent.x(), fMaxExtent.x(), << 1068 return G4VisExtent (xMinExtent, xMaxExtent, yMinExtent, yMaxExtent,
2115 fMinExtent.y(), fMaxExtent.y(), << 1069 zMinExtent, zMaxExtent);
2116 fMinExtent.z(), fMaxExtent.z() }; <<
2117 } 1070 }
2118 1071
2119 ///////////////////////////////////////////// 1072 ///////////////////////////////////////////////////////////////////////////////
2120 // 1073 //
2121 G4double G4TessellatedSolid::GetCubicVolume ( 1074 G4double G4TessellatedSolid::GetCubicVolume ()
2122 { 1075 {
2123 if (fCubicVolume != 0.) return fCubicVolume << 1076 if(cubicVolume != 0.) {;}
2124 << 1077 else { cubicVolume = G4VSolid::GetCubicVolume(); }
2125 // For explanation of the following algorit << 1078 return cubicVolume;
2126 // https://en.wikipedia.org/wiki/Polyhedron <<
2127 // http://wwwf.imperial.ac.uk/~rn/centroid. <<
2128 <<
2129 std::size_t size = fFacets.size(); <<
2130 for (std::size_t i = 0; i < size; ++i) <<
2131 { <<
2132 G4VFacet &facet = *fFacets[i]; <<
2133 G4double area = facet.GetArea(); <<
2134 G4ThreeVector unit_normal = facet.GetSurf <<
2135 fCubicVolume += area * (facet.GetVertex(0 <<
2136 } <<
2137 fCubicVolume /= 3.; <<
2138 return fCubicVolume; <<
2139 } 1079 }
2140 1080
2141 ///////////////////////////////////////////// 1081 ///////////////////////////////////////////////////////////////////////////////
2142 // 1082 //
2143 G4double G4TessellatedSolid::GetSurfaceArea ( 1083 G4double G4TessellatedSolid::GetSurfaceArea ()
2144 { 1084 {
2145 if (fSurfaceArea != 0.) return fSurfaceArea << 1085 if(surfaceArea != 0.) { return surfaceArea; }
2146 1086
2147 std::size_t size = fFacets.size(); << 1087 for (FacetI f=facets.begin(); f!=facets.end(); f++)
2148 for (std::size_t i = 0; i < size; ++i) <<
2149 { 1088 {
2150 G4VFacet &facet = *fFacets[i]; << 1089 surfaceArea += (*f)->GetArea();
2151 fSurfaceArea += facet.GetArea(); <<
2152 } 1090 }
2153 return fSurfaceArea; << 1091 return surfaceArea;
2154 } 1092 }
2155 1093
2156 ///////////////////////////////////////////// 1094 ///////////////////////////////////////////////////////////////////////////////
2157 // 1095 //
2158 G4ThreeVector G4TessellatedSolid::GetPointOnS 1096 G4ThreeVector G4TessellatedSolid::GetPointOnSurface() const
2159 { 1097 {
2160 // Select randomly a facet and return a ran 1098 // Select randomly a facet and return a random point on it
2161 1099
2162 auto i = (G4int) G4RandFlat::shoot(0., fFac << 1100 G4int i = CLHEP::RandFlat::shootInt(facets.size());
2163 return fFacets[i]->GetPointOnFace(); << 1101 return facets[i]->GetPointOnFace();
2164 } 1102 }
2165 <<
2166 ///////////////////////////////////////////// 1103 ///////////////////////////////////////////////////////////////////////////////
2167 // 1104 //
2168 // SetRandomVectorSet 1105 // SetRandomVectorSet
2169 // 1106 //
2170 // This is a set of predefined random vectors 1107 // This is a set of predefined random vectors (if that isn't a contradition
2171 // in terms!) used to generate rays from a us 1108 // in terms!) used to generate rays from a user-defined point. The member
2172 // function Inside uses these to determine wh 1109 // function Inside uses these to determine whether the point is inside or
2173 // outside of the tessellated solid. All vec 1110 // outside of the tessellated solid. All vectors should be unit vectors.
2174 // 1111 //
2175 void G4TessellatedSolid::SetRandomVectors () << 1112 void G4TessellatedSolid::SetRandomVectorSet()
2176 { 1113 {
2177 fRandir.resize(20); << 1114 randir[0] = G4ThreeVector(-0.9577428892113370, 0.2732676269591740, 0.0897405271949221);
2178 fRandir[0] = << 1115 randir[1] = G4ThreeVector(-0.8331264504940770,-0.5162067214954600,-0.1985722492445700);
2179 G4ThreeVector(-0.9577428892113370, 0.2732 << 1116 randir[2] = G4ThreeVector(-0.1516671651108820, 0.9666292616127460, 0.2064580868390110);
2180 fRandir[1] = << 1117 randir[3] = G4ThreeVector( 0.6570250350323190,-0.6944539025883300, 0.2933460081893360);
2181 G4ThreeVector(-0.8331264504940770,-0.5162 << 1118 randir[4] = G4ThreeVector(-0.4820456281280320,-0.6331060000098690,-0.6056474264406270);
2182 fRandir[2] = << 1119 randir[5] = G4ThreeVector( 0.7629032554236800, 0.1016854697539910,-0.6384658864065180);
2183 G4ThreeVector(-0.1516671651108820, 0.9666 << 1120 randir[6] = G4ThreeVector( 0.7689540409061150, 0.5034929891988220, 0.3939600142169160);
2184 fRandir[3] = << 1121 randir[7] = G4ThreeVector( 0.5765188359255740, 0.5997271636278330,-0.5549354566343150);
2185 G4ThreeVector( 0.6570250350323190,-0.6944 << 1122 randir[8] = G4ThreeVector( 0.6660632777862070,-0.6362809868288380, 0.3892379937580790);
2186 fRandir[4] = << 1123 randir[9] = G4ThreeVector( 0.3824415020414780, 0.6541792713761380,-0.6525243125110690);
2187 G4ThreeVector(-0.4820456281280320,-0.6331 << 1124 randir[10] = G4ThreeVector(-0.5107726564526760, 0.6020905056811610, 0.6136760679616570);
2188 fRandir[5] = << 1125 randir[11] = G4ThreeVector( 0.7459135439578050, 0.6618796061649330, 0.0743530220183488);
2189 G4ThreeVector( 0.7629032554236800 , 0.101 << 1126 randir[12] = G4ThreeVector( 0.1536405855311580, 0.8117477913978260,-0.5634359711967240);
2190 fRandir[6] = << 1127 randir[13] = G4ThreeVector( 0.0744395301705579,-0.8707110101772920,-0.4861286795736560);
2191 G4ThreeVector( 0.7689540409061150, 0.5034 << 1128 randir[14] = G4ThreeVector(-0.1665874645185400, 0.6018553940549240,-0.7810369397872780);
2192 fRandir[7] = << 1129 randir[15] = G4ThreeVector( 0.7766902003633100, 0.6014617505959970,-0.1870724331097450);
2193 G4ThreeVector( 0.5765188359255740, 0.5997 << 1130 randir[16] = G4ThreeVector(-0.8710128685847430,-0.1434320216603030,-0.4698551243971010);
2194 fRandir[8] = << 1131 randir[17] = G4ThreeVector( 0.8901082092766820,-0.4388411398893870, 0.1229871120030100);
2195 G4ThreeVector( 0.6660632777862070,-0.6362 << 1132 randir[18] = G4ThreeVector(-0.6430417431544370,-0.3295938228697690, 0.6912779675984150);
2196 fRandir[9] = << 1133 randir[19] = G4ThreeVector( 0.6331124368380410, 0.6306211461665000, 0.4488714875425340);
2197 G4ThreeVector( 0.3824415020414780, 0.6541 <<
2198 fRandir[10] = <<
2199 G4ThreeVector(-0.5107726564526760, 0.6020 <<
2200 fRandir[11] = <<
2201 G4ThreeVector( 0.7459135439578050, 0.6618 <<
2202 fRandir[12] = <<
2203 G4ThreeVector( 0.1536405855311580, 0.8117 <<
2204 fRandir[13] = <<
2205 G4ThreeVector( 0.0744395301705579,-0.8707 <<
2206 fRandir[14] = <<
2207 G4ThreeVector(-0.1665874645185400, 0.6018 <<
2208 fRandir[15] = <<
2209 G4ThreeVector( 0.7766902003633100, 0.6014 <<
2210 fRandir[16] = <<
2211 G4ThreeVector(-0.8710128685847430,-0.1434 <<
2212 fRandir[17] = <<
2213 G4ThreeVector( 0.8901082092766820,-0.4388 <<
2214 fRandir[18] = <<
2215 G4ThreeVector(-0.6430417431544370,-0.3295 <<
2216 fRandir[19] = <<
2217 G4ThreeVector( 0.6331124368380410, 0.6306 <<
2218 <<
2219 fMaxTries = 20; <<
2220 } <<
2221 <<
2222 ///////////////////////////////////////////// <<
2223 // <<
2224 G4int G4TessellatedSolid::AllocatedMemoryWith <<
2225 { <<
2226 G4int base = sizeof(*this); <<
2227 base += fVertexList.capacity() * sizeof(G4T <<
2228 base += fRandir.capacity() * sizeof(G4Three <<
2229 <<
2230 std::size_t limit = fFacets.size(); <<
2231 for (std::size_t i = 0; i < limit; ++i) <<
2232 { <<
2233 G4VFacet& facet = *fFacets[i]; <<
2234 base += facet.AllocatedMemory(); <<
2235 } <<
2236 <<
2237 for (const auto & fExtremeFacet : fExtremeF <<
2238 { <<
2239 G4VFacet &facet = *fExtremeFacet; <<
2240 base += facet.AllocatedMemory(); <<
2241 } <<
2242 return base; <<
2243 } <<
2244 <<
2245 ///////////////////////////////////////////// <<
2246 // <<
2247 G4int G4TessellatedSolid::AllocatedMemory() <<
2248 { <<
2249 G4int size = AllocatedMemoryWithoutVoxels() <<
2250 G4int sizeInsides = fInsides.GetNbytes(); <<
2251 G4int sizeVoxels = fVoxels.AllocatedMemory( <<
2252 size += sizeInsides + sizeVoxels; <<
2253 return size; <<
2254 } <<
2255 1134
2256 #endif << 1135 maxTries = 20;
>> 1136 }
2257 1137