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