source: trunk/source/geometry/solids/Boolean/src/G4UnionSolid.cc@ 1341

Last change on this file since 1341 was 1337, checked in by garnier, 15 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 14.7 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4UnionSolid.cc,v 1.37 2010/05/11 15:03:45 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30// Implementation of methods for the class G4IntersectionSolid
31//
32// History:
33//
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//
39// --------------------------------------------------------------------
40
41#include "G4UnionSolid.hh"
42
43#include <sstream>
44
45#include "G4VoxelLimits.hh"
46#include "G4VPVParameterisation.hh"
47#include "G4GeometryTolerance.hh"
48
49#include "G4VGraphicsScene.hh"
50#include "G4Polyhedron.hh"
51#include "HepPolyhedronProcessor.h"
52#include "G4NURBS.hh"
53// #include "G4NURBSbox.hh"
54
55///////////////////////////////////////////////////////////////////
56//
57// Transfer all data members to G4BooleanSolid which is responsible
58// for them. pName will be in turn sent to G4VSolid
59
60G4UnionSolid:: G4UnionSolid( const G4String& pName,
61 G4VSolid* pSolidA ,
62 G4VSolid* pSolidB )
63 : G4BooleanSolid(pName,pSolidA,pSolidB)
64{
65}
66
67/////////////////////////////////////////////////////////////////////
68//
69// Constructor
70
71G4UnionSolid::G4UnionSolid( const G4String& pName,
72 G4VSolid* pSolidA ,
73 G4VSolid* pSolidB ,
74 G4RotationMatrix* rotMatrix,
75 const G4ThreeVector& transVector )
76 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
77
78{
79}
80
81///////////////////////////////////////////////////////////
82//
83// Constructor
84
85G4UnionSolid::G4UnionSolid( const G4String& pName,
86 G4VSolid* pSolidA ,
87 G4VSolid* pSolidB ,
88 const G4Transform3D& transform )
89 : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
90{
91}
92
93//////////////////////////////////////////////////////////////////
94//
95// Fake default constructor - sets only member data and allocates memory
96// for usage restricted to object persistency.
97
98G4UnionSolid::G4UnionSolid( __void__& a )
99 : G4BooleanSolid(a)
100{
101}
102
103///////////////////////////////////////////////////////////
104//
105// Destructor
106
107G4UnionSolid::~G4UnionSolid()
108{
109}
110
111///////////////////////////////////////////////////////////////
112//
113//
114
115G4bool
116G4UnionSolid::CalculateExtent( const EAxis pAxis,
117 const G4VoxelLimits& pVoxelLimit,
118 const G4AffineTransform& pTransform,
119 G4double& pMin,
120 G4double& pMax ) const
121{
122 G4bool touchesA, touchesB, out ;
123 G4double minA = kInfinity, minB = kInfinity,
124 maxA = -kInfinity, maxB = -kInfinity;
125
126 touchesA = fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
127 pTransform, minA, maxA);
128 touchesB= fPtrSolidB->CalculateExtent( pAxis, pVoxelLimit,
129 pTransform, minB, maxB);
130 if( touchesA || touchesB )
131 {
132 pMin = std::min( minA, minB );
133 pMax = std::max( maxA, maxB );
134 out = true ;
135 }
136 else out = false ;
137
138 return out ; // It exists in this slice if either one does.
139}
140
141/////////////////////////////////////////////////////
142//
143// Important comment: When solids A and B touch together along flat
144// surface the surface points will be considered as kSurface, while points
145// located around will correspond to kInside
146
147EInside G4UnionSolid::Inside( const G4ThreeVector& p ) const
148{
149 EInside positionA = fPtrSolidA->Inside(p);
150 if (positionA == kInside) return kInside;
151
152 EInside positionB = fPtrSolidB->Inside(p);
153
154 if( positionA == kInside || positionB == kInside ||
155 ( positionA == kSurface && positionB == kSurface &&
156 ( fPtrSolidA->SurfaceNormal(p) +
157 fPtrSolidB->SurfaceNormal(p) ).mag2() <
158 1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
159 {
160 return kInside;
161 }
162 else
163 {
164 if( ( positionA != kInside && positionB == kSurface ) ||
165 ( positionB != kInside && positionA == kSurface ) ||
166 ( positionA == kSurface && positionB == kSurface ) ) return kSurface;
167 else return kOutside;
168 }
169}
170
171//////////////////////////////////////////////////////////////
172//
173//
174
175G4ThreeVector
176G4UnionSolid::SurfaceNormal( const G4ThreeVector& p ) const
177{
178 G4ThreeVector normal;
179
180#ifdef G4BOOLDEBUG
181 if( Inside(p) == kOutside )
182 {
183 G4cout << "WARNING - Invalid call in "
184 << "G4UnionSolid::SurfaceNormal(p)" << G4endl
185 << " Point p is outside !" << G4endl;
186 G4cout << " p = " << p << G4endl;
187 G4cerr << "WARNING - Invalid call in "
188 << "G4UnionSolid::SurfaceNormal(p)" << G4endl
189 << " Point p is outside !" << G4endl;
190 G4cerr << " p = " << p << G4endl;
191 }
192#endif
193
194 if(fPtrSolidA->Inside(p) == kSurface && fPtrSolidB->Inside(p) != kInside)
195 {
196 normal= fPtrSolidA->SurfaceNormal(p) ;
197 }
198 else if(fPtrSolidB->Inside(p) == kSurface &&
199 fPtrSolidA->Inside(p) != kInside)
200 {
201 normal= fPtrSolidB->SurfaceNormal(p) ;
202 }
203 else
204 {
205 normal= fPtrSolidA->SurfaceNormal(p) ;
206#ifdef G4BOOLDEBUG
207 if(Inside(p)==kInside)
208 {
209 G4cout << "WARNING - Invalid call in "
210 << "G4UnionSolid::SurfaceNormal(p)" << G4endl
211 << " Point p is inside !" << G4endl;
212 G4cout << " p = " << p << G4endl;
213 G4cerr << "WARNING - Invalid call in "
214 << "G4UnionSolid::SurfaceNormal(p)" << G4endl
215 << " Point p is inside !" << G4endl;
216 G4cerr << " p = " << p << G4endl;
217 }
218#endif
219 }
220 return normal;
221}
222
223/////////////////////////////////////////////////////////////
224//
225// The same algorithm as in DistanceToIn(p)
226
227G4double
228G4UnionSolid::DistanceToIn( const G4ThreeVector& p,
229 const G4ThreeVector& v ) const
230{
231#ifdef G4BOOLDEBUG
232 if( Inside(p) == kInside )
233 {
234 G4cout << "WARNING - Invalid call in "
235 << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
236 << " Point p is inside !" << G4endl;
237 G4cout << " p = " << p << G4endl;
238 G4cout << " v = " << v << G4endl;
239 G4cerr << "WARNING - Invalid call in "
240 << "G4UnionSolid::DistanceToIn(p,v)" << G4endl
241 << " Point p is inside !" << G4endl;
242 G4cerr << " p = " << p << G4endl;
243 G4cerr << " v = " << v << G4endl;
244 }
245#endif
246
247 return std::min(fPtrSolidA->DistanceToIn(p,v),
248 fPtrSolidB->DistanceToIn(p,v) ) ;
249}
250
251////////////////////////////////////////////////////////
252//
253// Approximate nearest distance from the point p to the union of
254// two solids
255
256G4double
257G4UnionSolid::DistanceToIn( const G4ThreeVector& p) const
258{
259#ifdef G4BOOLDEBUG
260 if( Inside(p) == kInside )
261 {
262 G4cout << "WARNING - Invalid call in "
263 << "G4UnionSolid::DistanceToIn(p)" << G4endl
264 << " Point p is inside !" << G4endl;
265 G4cout << " p = " << p << G4endl;
266 G4cerr << "WARNING - Invalid call in "
267 << "G4UnionSolid::DistanceToIn(p)" << G4endl
268 << " Point p is inside !" << G4endl;
269 G4cerr << " p = " << p << G4endl;
270 }
271#endif
272 G4double distA = fPtrSolidA->DistanceToIn(p) ;
273 G4double distB = fPtrSolidB->DistanceToIn(p) ;
274 G4double safety = std::min(distA,distB) ;
275 if(safety < 0.0) safety = 0.0 ;
276 return safety ;
277}
278
279//////////////////////////////////////////////////////////
280//
281// The same algorithm as DistanceToOut(p)
282
283G4double
284G4UnionSolid::DistanceToOut( const G4ThreeVector& p,
285 const G4ThreeVector& v,
286 const G4bool calcNorm,
287 G4bool *validNorm,
288 G4ThreeVector *n ) const
289{
290 G4double dist = 0.0, disTmp = 0.0 ;
291 G4ThreeVector normTmp;
292 G4ThreeVector* nTmp= &normTmp;
293
294 if( Inside(p) == kOutside )
295 {
296#ifdef G4BOOLDEBUG
297 G4cout << "Position:" << G4endl << G4endl;
298 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
299 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
300 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
301 G4cout << "Direction:" << G4endl << G4endl;
302 G4cout << "v.x() = " << v.x() << G4endl;
303 G4cout << "v.y() = " << v.y() << G4endl;
304 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
305 G4cout << "WARNING - Invalid call in "
306 << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
307 << " Point p is outside !" << G4endl;
308 G4cout << " p = " << p << G4endl;
309 G4cout << " v = " << v << G4endl;
310 G4cerr << "WARNING - Invalid call in "
311 << "G4UnionSolid::DistanceToOut(p,v)" << G4endl
312 << " Point p is outside !" << G4endl;
313 G4cerr << " p = " << p << G4endl;
314 G4cerr << " v = " << v << G4endl;
315#endif
316 }
317 else
318 {
319 EInside positionA = fPtrSolidA->Inside(p) ;
320 // EInside positionB = fPtrSolidB->Inside(p) ;
321
322 if( positionA != kOutside )
323 {
324 do
325 {
326 disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
327 validNorm,nTmp) ;
328 dist += disTmp ;
329
330 if(fPtrSolidB->Inside(p+dist*v) != kOutside)
331 {
332 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
333 validNorm,nTmp) ;
334 dist += disTmp ;
335 }
336 }
337 // while( Inside(p+dist*v) == kInside ) ;
338 while( fPtrSolidA->Inside(p+dist*v) != kOutside &&
339 disTmp > 0.5*kCarTolerance ) ;
340 }
341 else // if( positionB != kOutside )
342 {
343 do
344 {
345 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v,calcNorm,
346 validNorm,nTmp) ;
347 dist += disTmp ;
348
349 if(fPtrSolidA->Inside(p+dist*v) != kOutside)
350 {
351 disTmp = fPtrSolidA->DistanceToOut(p+dist*v,v,calcNorm,
352 validNorm,nTmp) ;
353 dist += disTmp ;
354 }
355 }
356 // while( Inside(p+dist*v) == kInside ) ;
357 while( (fPtrSolidB->Inside(p+dist*v) != kOutside)
358 && (disTmp > 0.5*kCarTolerance) ) ;
359 }
360 }
361 if( calcNorm )
362 {
363 *validNorm = false ;
364 *n = *nTmp ;
365 }
366 return dist ;
367}
368
369//////////////////////////////////////////////////////////////
370//
371// Inverted algorithm of DistanceToIn(p)
372
373G4double
374G4UnionSolid::DistanceToOut( const G4ThreeVector& p ) const
375{
376 G4double distout = 0.0;
377 if( Inside(p) == kOutside )
378 {
379#ifdef G4BOOLDEBUG
380 G4cout << "WARNING - Invalid call in "
381 << "G4UnionSolid::DistanceToOut(p)" << G4endl
382 << " Point p is outside !" << G4endl;
383 G4cout << " p = " << p << G4endl;
384 G4cerr << "WARNING - Invalid call in "
385 << "G4UnionSolid::DistanceToOut(p)" << G4endl
386 << " Point p is outside !" << G4endl;
387 G4cerr << " p = " << p << G4endl;
388#endif
389 }
390 else
391 {
392 EInside positionA = fPtrSolidA->Inside(p) ;
393 EInside positionB = fPtrSolidB->Inside(p) ;
394
395 // Is this equivalent ??
396 // if( ! ( (positionA == kOutside)) &&
397 // (positionB == kOutside)) )
398 if((positionA == kInside && positionB == kInside ) ||
399 (positionA == kInside && positionB == kSurface ) ||
400 (positionA == kSurface && positionB == kInside ) )
401 {
402 distout= std::max(fPtrSolidA->DistanceToOut(p),
403 fPtrSolidB->DistanceToOut(p) ) ;
404 }
405 else
406 {
407 if(positionA == kOutside)
408 {
409 distout= fPtrSolidB->DistanceToOut(p) ;
410 }
411 else
412 {
413 distout= fPtrSolidA->DistanceToOut(p) ;
414 }
415 }
416 }
417 return distout;
418}
419
420//////////////////////////////////////////////////////////////
421//
422//
423
424G4GeometryType G4UnionSolid::GetEntityType() const
425{
426 return G4String("G4UnionSolid");
427}
428
429//////////////////////////////////////////////////////////////
430//
431//
432
433void
434G4UnionSolid::ComputeDimensions( G4VPVParameterisation*,
435 const G4int,
436 const G4VPhysicalVolume* )
437{
438}
439
440/////////////////////////////////////////////////
441//
442//
443
444void
445G4UnionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
446{
447 scene.AddSolid (*this);
448}
449
450////////////////////////////////////////////////////
451//
452//
453
454G4Polyhedron*
455G4UnionSolid::CreatePolyhedron () const
456{
457 HepPolyhedronProcessor processor;
458 // Stack components and components of components recursively
459 // See G4BooleanSolid::StackPolyhedron
460 G4Polyhedron* top = StackPolyhedron(processor, this);
461 G4Polyhedron* result = new G4Polyhedron(*top);
462 if (processor.execute(*result)) { return result; }
463 else { return 0; }
464}
465
466/////////////////////////////////////////////////////////
467//
468//
469
470G4NURBS*
471G4UnionSolid::CreateNURBS () const
472{
473 // Take into account boolean operation - see CreatePolyhedron.
474 // return new G4NURBSbox (1.0, 1.0, 1.0);
475 return 0;
476}
Note: See TracBrowser for help on using the repository browser.