| 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. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Implementation of methods for the class G4U <<
27 // 26 //
28 // 23.04.18 E.Tcherniaev: added extended BBox, << 27 // $Id: G4UnionSolid.cc,v 1.35 2007/10/23 14:42:31 grichine Exp $
29 // 17.03.17 E.Tcherniaev: revision of SurfaceN << 28 // GEANT4 tag $Name: geant4-09-01-patch-03 $
>> 29 //
>> 30 // Implementation of methods for the class G4IntersectionSolid
>> 31 //
>> 32 // History:
>> 33 //
30 // 12.09.98 V.Grichine: first implementation 34 // 12.09.98 V.Grichine: first implementation
>> 35 // 28.11.98 V.Grichine: fix while loops in DistToIn/Out
>> 36 // 27.07.99 V.Grichine: modifications in DistToOut(p,v,...), while -> do-while
>> 37 // 16.03.01 V.Grichine: modifications in CalculateExtent()
>> 38 //
31 // ------------------------------------------- 39 // --------------------------------------------------------------------
32 40
33 #include <sstream> <<
34 <<
35 #include "G4UnionSolid.hh" 41 #include "G4UnionSolid.hh"
36 42
37 #include "G4SystemOfUnits.hh" << 43 #include <sstream>
>> 44
38 #include "G4VoxelLimits.hh" 45 #include "G4VoxelLimits.hh"
39 #include "G4VPVParameterisation.hh" 46 #include "G4VPVParameterisation.hh"
40 #include "G4GeometryTolerance.hh" 47 #include "G4GeometryTolerance.hh"
41 48
42 #include "G4VGraphicsScene.hh" 49 #include "G4VGraphicsScene.hh"
43 #include "G4Polyhedron.hh" 50 #include "G4Polyhedron.hh"
44 #include "G4PolyhedronArbitrary.hh" << 51 #include "G4NURBS.hh"
45 #include "HepPolyhedronProcessor.h" << 52 // #include "G4NURBSbox.hh"
46 <<
47 #include "G4IntersectionSolid.hh" <<
48 53
49 ////////////////////////////////////////////// << 54 ///////////////////////////////////////////////////////////////////
50 // 55 //
51 // Transfer all data members to G4BooleanSolid 56 // Transfer all data members to G4BooleanSolid which is responsible
52 // for them. pName will be in turn sent to G4V 57 // for them. pName will be in turn sent to G4VSolid
53 58
54 G4UnionSolid:: G4UnionSolid( const G4String& p 59 G4UnionSolid:: G4UnionSolid( const G4String& pName,
55 G4VSolid* p 60 G4VSolid* pSolidA ,
56 G4VSolid* p 61 G4VSolid* pSolidB )
57 : G4BooleanSolid(pName,pSolidA,pSolidB) 62 : G4BooleanSolid(pName,pSolidA,pSolidB)
58 { 63 {
59 Init(); <<
60 } 64 }
61 65
62 ////////////////////////////////////////////// << 66 /////////////////////////////////////////////////////////////////////
63 // 67 //
64 // Constructor 68 // Constructor
65 69
66 G4UnionSolid::G4UnionSolid( const G4String& pN 70 G4UnionSolid::G4UnionSolid( const G4String& pName,
67 G4VSolid* pS 71 G4VSolid* pSolidA ,
68 G4VSolid* pS 72 G4VSolid* pSolidB ,
69 G4RotationMa 73 G4RotationMatrix* rotMatrix,
70 const G4ThreeVecto 74 const G4ThreeVector& transVector )
71 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMa 75 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
72 76
73 { 77 {
74 Init(); <<
75 } 78 }
76 79
77 ////////////////////////////////////////////// << 80 ///////////////////////////////////////////////////////////
78 // 81 //
79 // Constructor 82 // Constructor
80 83
81 G4UnionSolid::G4UnionSolid( const G4String& pN 84 G4UnionSolid::G4UnionSolid( const G4String& pName,
82 G4VSolid* pS 85 G4VSolid* pSolidA ,
83 G4VSolid* pS 86 G4VSolid* pSolidB ,
84 const G4Transform3 87 const G4Transform3D& transform )
85 : G4BooleanSolid(pName,pSolidA,pSolidB,trans 88 : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
86 { 89 {
87 Init(); <<
88 } 90 }
89 91
90 ////////////////////////////////////////////// << 92 //////////////////////////////////////////////////////////////////
91 // 93 //
92 // Fake default constructor - sets only member 94 // Fake default constructor - sets only member data and allocates memory
93 // for usage restri 95 // for usage restricted to object persistency.
94 96
95 G4UnionSolid::G4UnionSolid( __void__& a ) 97 G4UnionSolid::G4UnionSolid( __void__& a )
96 : G4BooleanSolid(a) 98 : G4BooleanSolid(a)
97 { 99 {
98 } 100 }
99 101
100 ////////////////////////////////////////////// << 102 ///////////////////////////////////////////////////////////
101 // 103 //
102 // Destructor 104 // Destructor
103 105
104 G4UnionSolid::~G4UnionSolid() 106 G4UnionSolid::~G4UnionSolid()
105 = default; <<
106 <<
107 ////////////////////////////////////////////// <<
108 // <<
109 // Copy constructor <<
110 <<
111 G4UnionSolid::G4UnionSolid(const G4UnionSolid& <<
112 : G4BooleanSolid (rhs) <<
113 { 107 {
114 fPMin = rhs.fPMin; <<
115 fPMax = rhs.fPMax; <<
116 halfCarTolerance=0.5*kCarTolerance; <<
117 } 108 }
118 109
119 ////////////////////////////////////////////// << 110 ///////////////////////////////////////////////////////////////
120 // 111 //
121 // Assignment operator <<
122 <<
123 G4UnionSolid& G4UnionSolid::operator = (const <<
124 { <<
125 // Check assignment to self <<
126 // <<
127 if (this == &rhs) { return *this; } <<
128 <<
129 // Copy base class data <<
130 // <<
131 G4BooleanSolid::operator=(rhs); <<
132 <<
133 fPMin = rhs.fPMin; <<
134 fPMax = rhs.fPMax; <<
135 halfCarTolerance = rhs.halfCarTolerance; <<
136 <<
137 return *this; <<
138 } <<
139 <<
140 ////////////////////////////////////////////// <<
141 // <<
142 // Initialisation <<
143 <<
144 void G4UnionSolid::Init() <<
145 { <<
146 G4ThreeVector pdelta(kCarTolerance,kCarToler <<
147 G4ThreeVector pmin, pmax; <<
148 BoundingLimits(pmin, pmax); <<
149 fPMin = pmin - pdelta; <<
150 fPMax = pmax + pdelta; <<
151 halfCarTolerance=0.5*kCarTolerance; <<
152 } <<
153 <<
154 ////////////////////////////////////////////// <<
155 // 112 //
156 // Get bounding box <<
157 <<
158 void G4UnionSolid::BoundingLimits(G4ThreeVecto <<
159 G4ThreeVecto <<
160 { <<
161 G4ThreeVector minA,maxA, minB,maxB; <<
162 fPtrSolidA->BoundingLimits(minA,maxA); <<
163 fPtrSolidB->BoundingLimits(minB,maxB); <<
164 <<
165 pMin.set(std::min(minA.x(),minB.x()), <<
166 std::min(minA.y(),minB.y()), <<
167 std::min(minA.z(),minB.z())); <<
168 <<
169 pMax.set(std::max(maxA.x(),maxB.x()), <<
170 std::max(maxA.y(),maxB.y()), <<
171 std::max(maxA.z(),maxB.z())); <<
172 <<
173 // Check correctness of the bounding box <<
174 // <<
175 if (pMin.x() >= pMax.x() || pMin.y() >= pMax <<
176 { <<
177 std::ostringstream message; <<
178 message << "Bad bounding box (min >= max) <<
179 << GetName() << " !" <<
180 << "\npMin = " << pMin <<
181 << "\npMax = " << pMax; <<
182 G4Exception("G4UnionSolid::BoundingLimits( <<
183 JustWarning, message); <<
184 DumpInfo(); <<
185 } <<
186 } <<
187 <<
188 ////////////////////////////////////////////// <<
189 // <<
190 // Calculate extent under transform and specif <<
191 113
192 G4bool 114 G4bool
193 G4UnionSolid::CalculateExtent( const EAxis pAx 115 G4UnionSolid::CalculateExtent( const EAxis pAxis,
194 const G4VoxelLi 116 const G4VoxelLimits& pVoxelLimit,
195 const G4AffineT 117 const G4AffineTransform& pTransform,
196 G4double& 118 G4double& pMin,
197 G4double& 119 G4double& pMax ) const
198 { 120 {
199 G4bool touchesA, touchesB, out ; 121 G4bool touchesA, touchesB, out ;
200 G4double minA = kInfinity, minB = kInfinit 122 G4double minA = kInfinity, minB = kInfinity,
201 maxA = -kInfinity, maxB = -kInfinit 123 maxA = -kInfinity, maxB = -kInfinity;
202 124
203 touchesA = fPtrSolidA->CalculateExtent( pAxi 125 touchesA = fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
204 pTra 126 pTransform, minA, maxA);
205 touchesB = fPtrSolidB->CalculateExtent( pAxi << 127 touchesB= fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit,
206 pTra << 128 pTransform, minB, maxB);
207 if( touchesA || touchesB ) 129 if( touchesA || touchesB )
208 { 130 {
209 pMin = std::min( minA, minB ); 131 pMin = std::min( minA, minB );
210 pMax = std::max( maxA, maxB ); 132 pMax = std::max( maxA, maxB );
211 out = true ; 133 out = true ;
212 } 134 }
213 else << 135 else out = false ;
214 { <<
215 out = false ; <<
216 } <<
217 136
218 return out ; // It exists in this slice if 137 return out ; // It exists in this slice if either one does.
219 } 138 }
220 139
221 ////////////////////////////////////////////// << 140 /////////////////////////////////////////////////////
222 // 141 //
223 // Important comment: When solids A and B touc 142 // Important comment: When solids A and B touch together along flat
224 // surface the surface points will be consider 143 // surface the surface points will be considered as kSurface, while points
225 // located around will correspond to kInside 144 // located around will correspond to kInside
226 145
227 EInside G4UnionSolid::Inside( const G4ThreeVec 146 EInside G4UnionSolid::Inside( const G4ThreeVector& p ) const
228 { 147 {
229 if (std::max(p.z()-fPMax.z(), fPMin.z()-p.z( <<
230 <<
231 EInside positionA = fPtrSolidA->Inside(p); 148 EInside positionA = fPtrSolidA->Inside(p);
232 if (positionA == kInside) { return position << 149 if (positionA == kInside) return kInside;
233 EInside positionB = fPtrSolidB->Inside(p); <<
234 if (positionA == kOutside) { return position <<
235 <<
236 if (positionB == kInside) { return position <<
237 if (positionB == kOutside) { return position <<
238 150
239 // Both points are on surface << 151 EInside positionB = fPtrSolidB->Inside(p);
240 // <<
241 static const G4double rtol <<
242 = 1000*G4GeometryTolerance::GetInstance()- <<
243 152
244 return ((fPtrSolidA->SurfaceNormal(p) + << 153 if( positionA == kInside || positionB == kInside ||
245 fPtrSolidB->SurfaceNormal(p)).mag2( << 154 ( positionA == kSurface && positionB == kSurface &&
>> 155 ( fPtrSolidA->SurfaceNormal(p) +
>> 156 fPtrSolidB->SurfaceNormal(p) ).mag2() <
>> 157 1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
>> 158 {
>> 159 return kInside;
>> 160 }
>> 161 else
>> 162 {
>> 163 if( ( positionA != kInside && positionB == kSurface ) ||
>> 164 ( positionB != kInside && positionA == kSurface ) ||
>> 165 ( positionA == kSurface && positionB == kSurface ) ) return kSurface;
>> 166 else return kOutside;
>> 167 }
246 } 168 }
247 169
248 ////////////////////////////////////////////// << 170 //////////////////////////////////////////////////////////////
>> 171 //
249 // 172 //
250 // Get surface normal <<
251 173
252 G4ThreeVector 174 G4ThreeVector
253 G4UnionSolid::SurfaceNormal( const G4ThreeVect 175 G4UnionSolid::SurfaceNormal( const G4ThreeVector& p ) const
254 { 176 {
255 EInside positionA = fPtrSolidA->Inside(p); << 177 G4ThreeVector normal;
256 EInside positionB = fPtrSolidB->Inside(p); <<
257 <<
258 if (positionA == kSurface && <<
259 positionB == kOutside) return fPtrSolidA <<
260 178
261 if (positionA == kOutside && << 179 #ifdef G4BOOLDEBUG
262 positionB == kSurface) return fPtrSolidB << 180 if( Inside(p) == kOutside )
>> 181 {
>> 182 G4cout << "WARNING - Invalid call in "
>> 183 << "G4UnionSolid::SurfaceNormal(p)" << G4endl
>> 184 << " Point p is outside !" << G4endl;
>> 185 G4cout << " p = " << p << G4endl;
>> 186 G4cerr << "WARNING - Invalid call in "
>> 187 << "G4UnionSolid::SurfaceNormal(p)" << G4endl
>> 188 << " Point p is outside !" << G4endl;
>> 189 G4cerr << " p = " << p << G4endl;
>> 190 }
>> 191 #endif
263 192
264 if (positionA == kSurface && << 193 if(fPtrSolidA->Inside(p) == kSurface && fPtrSolidB->Inside(p) != kInside)
265 positionB == kSurface) <<
266 { <<
267 if (Inside(p) == kSurface) <<
268 { 194 {
269 G4ThreeVector normalA = fPtrSolidA->Surf << 195 normal= fPtrSolidA->SurfaceNormal(p) ;
270 G4ThreeVector normalB = fPtrSolidB->Surf <<
271 return (normalA + normalB).unit(); <<
272 } 196 }
273 } << 197 else if(fPtrSolidB->Inside(p) == kSurface &&
>> 198 fPtrSolidA->Inside(p) != kInside)
>> 199 {
>> 200 normal= fPtrSolidB->SurfaceNormal(p) ;
>> 201 }
>> 202 else
>> 203 {
>> 204 normal= fPtrSolidA->SurfaceNormal(p) ;
274 #ifdef G4BOOLDEBUG 205 #ifdef G4BOOLDEBUG
275 G4String surf[3] = { "OUTSIDE", "SURFACE", " << 206 if(Inside(p)==kInside)
276 std::ostringstream message; << 207 {
277 G4int oldprc = message.precision(16); << 208 G4cout << "WARNING - Invalid call in "
278 message << "Invalid call of SurfaceNormal(p) << 209 << "G4UnionSolid::SurfaceNormal(p)" << G4endl
279 << GetName() << " !" << 210 << " Point p is inside !" << G4endl;
280 << "\nPoint p" << p << " is " << sur << 211 G4cout << " p = " << p << G4endl;
281 message.precision(oldprc); << 212 G4cerr << "WARNING - Invalid call in "
282 G4Exception("G4UnionSolid::SurfaceNormal()", << 213 << "G4UnionSolid::SurfaceNormal(p)" << G4endl
283 JustWarning, message); << 214 << " Point p is inside !" << G4endl;
>> 215 G4cerr << " p = " << p << G4endl;
>> 216 }
284 #endif 217 #endif
285 return fPtrSolidA->SurfaceNormal(p); << 218 }
>> 219 return normal;
286 } 220 }
287 221
288 ////////////////////////////////////////////// << 222 /////////////////////////////////////////////////////////////
289 // 223 //
290 // The same algorithm as in DistanceToIn(p) 224 // The same algorithm as in DistanceToIn(p)
291 225
292 G4double 226 G4double
293 G4UnionSolid::DistanceToIn( const G4ThreeVecto 227 G4UnionSolid::DistanceToIn( const G4ThreeVector& p,
294 const G4ThreeVecto << 228 const G4ThreeVector& v ) const
295 { 229 {
296 #ifdef G4BOOLDEBUG 230 #ifdef G4BOOLDEBUG
297 if( Inside(p) == kInside ) 231 if( Inside(p) == kInside )
298 { 232 {
299 G4cout << "WARNING - Invalid call in " 233 G4cout << "WARNING - Invalid call in "
300 << "G4UnionSolid::DistanceToIn(p,v) 234 << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
301 << " Point p is inside !" << G4end 235 << " Point p is inside !" << G4endl;
302 G4cout << " p = " << p << G4endl; 236 G4cout << " p = " << p << G4endl;
303 G4cout << " v = " << v << G4endl; 237 G4cout << " v = " << v << G4endl;
304 G4cerr << "WARNING - Invalid call in " 238 G4cerr << "WARNING - Invalid call in "
305 << "G4UnionSolid::DistanceToIn(p,v) 239 << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
306 << " Point p is inside !" << G4end 240 << " Point p is inside !" << G4endl;
307 G4cerr << " p = " << p << G4endl; 241 G4cerr << " p = " << p << G4endl;
308 G4cerr << " v = " << v << G4endl; 242 G4cerr << " v = " << v << G4endl;
309 } 243 }
310 #endif 244 #endif
311 245
312 return std::min(fPtrSolidA->DistanceToIn(p,v 246 return std::min(fPtrSolidA->DistanceToIn(p,v),
313 fPtrSolidB->DistanceToIn(p,v << 247 fPtrSolidB->DistanceToIn(p,v) ) ;
314 } 248 }
315 249
316 ////////////////////////////////////////////// << 250 ////////////////////////////////////////////////////////
317 // 251 //
318 // Approximate nearest distance from the point 252 // Approximate nearest distance from the point p to the union of
319 // two solids 253 // two solids
320 254
321 G4double 255 G4double
322 G4UnionSolid::DistanceToIn( const G4ThreeVecto << 256 G4UnionSolid::DistanceToIn( const G4ThreeVector& p) const
323 { 257 {
324 #ifdef G4BOOLDEBUG 258 #ifdef G4BOOLDEBUG
325 if( Inside(p) == kInside ) 259 if( Inside(p) == kInside )
326 { 260 {
327 G4cout << "WARNING - Invalid call in " 261 G4cout << "WARNING - Invalid call in "
328 << "G4UnionSolid::DistanceToIn(p)" 262 << "G4UnionSolid::DistanceToIn(p)" << G4endl
329 << " Point p is inside !" << G4end 263 << " Point p is inside !" << G4endl;
330 G4cout << " p = " << p << G4endl; 264 G4cout << " p = " << p << G4endl;
331 G4cerr << "WARNING - Invalid call in " 265 G4cerr << "WARNING - Invalid call in "
332 << "G4UnionSolid::DistanceToIn(p)" 266 << "G4UnionSolid::DistanceToIn(p)" << G4endl
333 << " Point p is inside !" << G4end 267 << " Point p is inside !" << G4endl;
334 G4cerr << " p = " << p << G4endl; 268 G4cerr << " p = " << p << G4endl;
335 } 269 }
336 #endif 270 #endif
337 G4double distA = fPtrSolidA->DistanceToIn(p) 271 G4double distA = fPtrSolidA->DistanceToIn(p) ;
338 G4double distB = fPtrSolidB->DistanceToIn(p) 272 G4double distB = fPtrSolidB->DistanceToIn(p) ;
339 G4double safety = std::min(distA,distB) ; 273 G4double safety = std::min(distA,distB) ;
340 if(safety < 0.0) safety = 0.0 ; 274 if(safety < 0.0) safety = 0.0 ;
341 return safety ; 275 return safety ;
342 } 276 }
343 277
344 ////////////////////////////////////////////// << 278 //////////////////////////////////////////////////////////
345 // 279 //
346 // The same algorithm as DistanceToOut(p) 280 // The same algorithm as DistanceToOut(p)
347 281
348 G4double 282 G4double
349 G4UnionSolid::DistanceToOut( const G4ThreeVect 283 G4UnionSolid::DistanceToOut( const G4ThreeVector& p,
350 const G4ThreeVect << 284 const G4ThreeVector& v,
351 const G4bool calc << 285 const G4bool calcNorm,
352 G4bool* val << 286 G4bool *validNorm,
353 G4ThreeVect << 287 G4ThreeVector *n ) const
354 { 288 {
355 G4double dist = 0.0, disTmp = 0.0 ; 289 G4double dist = 0.0, disTmp = 0.0 ;
356 G4ThreeVector normTmp; 290 G4ThreeVector normTmp;
357 G4ThreeVector* nTmp = &normTmp; << 291 G4ThreeVector* nTmp= &normTmp;
358 292
359 if( Inside(p) == kOutside ) 293 if( Inside(p) == kOutside )
360 { 294 {
361 #ifdef G4BOOLDEBUG 295 #ifdef G4BOOLDEBUG
362 G4cout << "Position:" << G4endl << G4en 296 G4cout << "Position:" << G4endl << G4endl;
363 G4cout << "p.x() = " << p.x()/mm << " 297 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
364 G4cout << "p.y() = " << p.y()/mm << " 298 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
365 G4cout << "p.z() = " << p.z()/mm << " 299 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
366 G4cout << "Direction:" << G4endl << G4en 300 G4cout << "Direction:" << G4endl << G4endl;
367 G4cout << "v.x() = " << v.x() << G4end 301 G4cout << "v.x() = " << v.x() << G4endl;
368 G4cout << "v.y() = " << v.y() << G4end 302 G4cout << "v.y() = " << v.y() << G4endl;
369 G4cout << "v.z() = " << v.z() << G4end 303 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
370 G4cout << "WARNING - Invalid call in " 304 G4cout << "WARNING - Invalid call in "
371 << "G4UnionSolid::DistanceToOut(p 305 << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
372 << " Point p is outside !" << G4 306 << " Point p is outside !" << G4endl;
373 G4cout << " p = " << p << G4end 307 G4cout << " p = " << p << G4endl;
374 G4cout << " v = " << v << G4end 308 G4cout << " v = " << v << G4endl;
375 G4cerr << "WARNING - Invalid call in " 309 G4cerr << "WARNING - Invalid call in "
376 << "G4UnionSolid::DistanceToOut(p 310 << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
377 << " Point p is outside !" << G4 311 << " Point p is outside !" << G4endl;
378 G4cerr << " p = " << p << G4end 312 G4cerr << " p = " << p << G4endl;
379 G4cerr << " v = " << v << G4end 313 G4cerr << " v = " << v << G4endl;
380 #endif 314 #endif
381 } 315 }
382 else 316 else
383 { 317 {
384 EInside positionA = fPtrSolidA->Inside(p) 318 EInside positionA = fPtrSolidA->Inside(p) ;
>> 319 // EInside positionB = fPtrSolidB->Inside(p) ;
385 320
386 if( positionA != kOutside ) 321 if( positionA != kOutside )
387 { 322 {
388 do // Loop checking, 13.08.2015, G.Cosm << 323 do
389 { 324 {
390 disTmp = fPtrSolidA->DistanceToOut(p+d 325 disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
391 val << 326 validNorm,nTmp) ;
392 dist += disTmp ; 327 dist += disTmp ;
393 328
394 if(fPtrSolidB->Inside(p+dist*v) != kOu 329 if(fPtrSolidB->Inside(p+dist*v) != kOutside)
395 { 330 {
396 disTmp = fPtrSolidB->DistanceToOut(p 331 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
397 v << 332 validNorm,nTmp) ;
398 dist += disTmp ; 333 dist += disTmp ;
399 } 334 }
400 } 335 }
401 while( (fPtrSolidA->Inside(p+dist*v) != << 336 // while( Inside(p+dist*v) == kInside ) ;
402 && (disTmp > halfCarTolerance) ); << 337 while( fPtrSolidA->Inside(p+dist*v) != kOutside &&
>> 338 disTmp > 0.5*kCarTolerance ) ;
403 } 339 }
404 else // if( positionB != kOutside ) 340 else // if( positionB != kOutside )
405 { 341 {
406 do // Loop checking, 13.08.2015, G.Cosm << 342 do
407 { 343 {
408 disTmp = fPtrSolidB->DistanceToOut(p+d 344 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
409 val << 345 validNorm,nTmp) ;
410 dist += disTmp ; 346 dist += disTmp ;
411 347
412 if(fPtrSolidA->Inside(p+dist*v) != kOu 348 if(fPtrSolidA->Inside(p+dist*v) != kOutside)
413 { 349 {
414 disTmp = fPtrSolidA->DistanceToOut(p 350 disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
415 v << 351 validNorm,nTmp) ;
416 dist += disTmp ; 352 dist += disTmp ;
417 } 353 }
418 } 354 }
419 while( (fPtrSolidB->Inside(p+dist*v) != << 355 // while( Inside(p+dist*v) == kInside ) ;
420 && (disTmp > halfCarTolerance) ); << 356 while( (fPtrSolidB->Inside(p+dist*v) != kOutside)
>> 357 && (disTmp > 0.5*kCarTolerance) ) ;
421 } 358 }
422 } 359 }
423 if( calcNorm ) 360 if( calcNorm )
424 { 361 {
425 *validNorm = false ; 362 *validNorm = false ;
426 *n = *nTmp ; 363 *n = *nTmp ;
427 } 364 }
428 return dist ; 365 return dist ;
429 } 366 }
430 367
431 ////////////////////////////////////////////// << 368 //////////////////////////////////////////////////////////////
432 // 369 //
433 // Inverted algorithm of DistanceToIn(p) 370 // Inverted algorithm of DistanceToIn(p)
434 371
435 G4double 372 G4double
436 G4UnionSolid::DistanceToOut( const G4ThreeVect 373 G4UnionSolid::DistanceToOut( const G4ThreeVector& p ) const
437 { 374 {
438 G4double distout = 0.0; 375 G4double distout = 0.0;
439 if( Inside(p) == kOutside ) 376 if( Inside(p) == kOutside )
440 { 377 {
441 #ifdef G4BOOLDEBUG 378 #ifdef G4BOOLDEBUG
442 G4cout << "WARNING - Invalid call in " 379 G4cout << "WARNING - Invalid call in "
443 << "G4UnionSolid::DistanceToOut(p)" 380 << "G4UnionSolid::DistanceToOut(p)" << G4endl
444 << " Point p is outside !" << G4en 381 << " Point p is outside !" << G4endl;
445 G4cout << " p = " << p << G4endl; 382 G4cout << " p = " << p << G4endl;
446 G4cerr << "WARNING - Invalid call in " 383 G4cerr << "WARNING - Invalid call in "
447 << "G4UnionSolid::DistanceToOut(p)" 384 << "G4UnionSolid::DistanceToOut(p)" << G4endl
448 << " Point p is outside !" << G4en 385 << " Point p is outside !" << G4endl;
449 G4cerr << " p = " << p << G4endl; 386 G4cerr << " p = " << p << G4endl;
450 #endif 387 #endif
451 } 388 }
452 else 389 else
453 { 390 {
454 EInside positionA = fPtrSolidA->Inside(p) 391 EInside positionA = fPtrSolidA->Inside(p) ;
455 EInside positionB = fPtrSolidB->Inside(p) 392 EInside positionB = fPtrSolidB->Inside(p) ;
456 393
457 // Is this equivalent ?? 394 // Is this equivalent ??
458 // if( ! ( (positionA == kOutside)) && 395 // if( ! ( (positionA == kOutside)) &&
459 // (positionB == kOutside)) ) 396 // (positionB == kOutside)) )
460 if((positionA == kInside && positionB == 397 if((positionA == kInside && positionB == kInside ) ||
461 (positionA == kInside && positionB == 398 (positionA == kInside && positionB == kSurface ) ||
462 (positionA == kSurface && positionB == 399 (positionA == kSurface && positionB == kInside ) )
463 { 400 {
464 distout= std::max(fPtrSolidA->DistanceTo 401 distout= std::max(fPtrSolidA->DistanceToOut(p),
465 fPtrSolidB->DistanceTo << 402 fPtrSolidB->DistanceToOut(p) ) ;
466 } 403 }
467 else 404 else
468 { 405 {
469 if(positionA == kOutside) 406 if(positionA == kOutside)
470 { 407 {
471 distout= fPtrSolidB->DistanceToOut(p) 408 distout= fPtrSolidB->DistanceToOut(p) ;
472 } 409 }
473 else 410 else
474 { 411 {
475 distout= fPtrSolidA->DistanceToOut(p) 412 distout= fPtrSolidA->DistanceToOut(p) ;
476 } 413 }
477 } 414 }
478 } 415 }
479 return distout; 416 return distout;
480 } 417 }
481 418
482 ////////////////////////////////////////////// << 419 //////////////////////////////////////////////////////////////
>> 420 //
483 // 421 //
484 // GetEntityType <<
485 422
486 G4GeometryType G4UnionSolid::GetEntityType() c 423 G4GeometryType G4UnionSolid::GetEntityType() const
487 { 424 {
488 return {"G4UnionSolid"}; << 425 return G4String("G4UnionSolid");
489 } 426 }
490 427
491 ////////////////////////////////////////////// << 428 //////////////////////////////////////////////////////////////
492 // 429 //
493 // Make a clone of the object <<
494 <<
495 G4VSolid* G4UnionSolid::Clone() const <<
496 { <<
497 return new G4UnionSolid(*this); <<
498 } <<
499 <<
500 ////////////////////////////////////////////// <<
501 // 430 //
502 // ComputeDimensions <<
503 431
504 void 432 void
505 G4UnionSolid::ComputeDimensions( G4VPVPa 433 G4UnionSolid::ComputeDimensions( G4VPVParameterisation*,
506 const G4int, 434 const G4int,
507 const G4VPhys 435 const G4VPhysicalVolume* )
508 { 436 {
509 } 437 }
510 438
511 ////////////////////////////////////////////// << 439 /////////////////////////////////////////////////
512 // 440 //
513 // DescribeYourselfTo << 441 //
514 442
515 void 443 void
516 G4UnionSolid::DescribeYourselfTo ( G4VGraphics 444 G4UnionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
517 { 445 {
518 scene.AddSolid (*this); 446 scene.AddSolid (*this);
519 } 447 }
520 448
521 ////////////////////////////////////////////// << 449 ////////////////////////////////////////////////////
>> 450 //
522 // 451 //
523 // CreatePolyhedron <<
524 452
525 G4Polyhedron* 453 G4Polyhedron*
526 G4UnionSolid::CreatePolyhedron () const << 454 G4UnionSolid::CreatePolyhedron () const
527 { 455 {
528 if (fExternalBoolProcessor == nullptr) << 456 G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
529 { << 457 G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
530 HepPolyhedronProcessor processor; << 458 if (pA && pB) {
531 // Stack components and components of comp << 459 G4Polyhedron* resultant = new G4Polyhedron (pA->add(*pB));
532 // See G4BooleanSolid::StackPolyhedron << 460 return resultant;
533 G4Polyhedron* top = StackPolyhedron(proces << 461 } else {
534 auto result = new G4Polyhedron(*top); << 462 std::ostringstream oss;
535 if (processor.execute(*result)) << 463 oss << GetName() <<
536 { << 464 ": one of the Boolean components has no corresponding polyhedron.";
537 return result; << 465 G4Exception("G4UnionSolid::CreatePolyhedron",
538 } << 466 "", JustWarning, oss.str().c_str());
539 else << 467 return 0;
540 { <<
541 return nullptr; <<
542 } <<
543 } <<
544 else <<
545 { <<
546 return fExternalBoolProcessor->Process(thi <<
547 } 468 }
548 } 469 }
549 470
550 ////////////////////////////////////////////// << 471 /////////////////////////////////////////////////////////
>> 472 //
551 // 473 //
552 // GetCubicVolume <<
553 474
554 G4double G4UnionSolid::GetCubicVolume() << 475 G4NURBS*
>> 476 G4UnionSolid::CreateNURBS () const
555 { 477 {
556 if( fCubicVolume >= 0. ) << 478 // Take into account boolean operation - see CreatePolyhedron.
557 { << 479 // return new G4NURBSbox (1.0, 1.0, 1.0);
558 return fCubicVolume; << 480 return 0;
559 } <<
560 G4ThreeVector bminA, bmaxA, bminB, bmaxB; <<
561 fPtrSolidA->BoundingLimits(bminA, bmaxA); <<
562 fPtrSolidB->BoundingLimits(bminB, bmaxB); <<
563 G4bool noIntersection = <<
564 bminA.x() >= bmaxB.x() || bminA.y() >= bm <<
565 bminB.x() >= bmaxA.x() || bminB.y() >= bm <<
566 <<
567 if (noIntersection) <<
568 { <<
569 fCubicVolume = fPtrSolidA->GetCubicVolume( <<
570 } <<
571 else <<
572 { <<
573 if (GetNumOfConstituents() > 10) <<
574 { <<
575 fCubicVolume = G4BooleanSolid::GetCubicV <<
576 } <<
577 else <<
578 { <<
579 G4IntersectionSolid intersectVol("Tempor <<
580 fPtrSol <<
581 intersectVol.SetCubVolStatistics(GetCubV <<
582 intersectVol.SetCubVolEpsilon(GetCubVolE <<
583 <<
584 fCubicVolume = fPtrSolidA->GetCubicVolum <<
585 - intersectVol.GetCubicVolume(); <<
586 } <<
587 } <<
588 return fCubicVolume; <<
589 } 481 }
590 482