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