source: trunk/source/geometry/solids/Boolean/src/G4SubtractionSolid.cc@ 1347

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

tag geant4.9.4 beta 1 + modifs locales
File size: 14.5 KB
RevLine 
[831]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//
[1315]27// $Id: G4SubtractionSolid.cc,v 1.33 2010/05/11 15:03:45 gcosmo Exp $
[1337]28// GEANT4 tag $Name: geant4-09-04-beta-01 $
[831]29//
30// Implementation of methods for the class G4IntersectionSolid
31//
32// History:
33//
34// 14.10.98 V.Grichine: implementation of the first version
35// 19.10.98 V.Grichine: new algorithm of DistanceToIn(p,v)
36// 02.08.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
37// while -> do-while & surfaceA limitations
38// 13.09.00 V.Grichine: bug fixed in SurfaceNormal(p), p can be inside
39//
40// --------------------------------------------------------------------
41
42#include "G4SubtractionSolid.hh"
43
44#include "G4VoxelLimits.hh"
45#include "G4VPVParameterisation.hh"
46#include "G4GeometryTolerance.hh"
47
48#include "G4VGraphicsScene.hh"
49#include "G4Polyhedron.hh"
[1315]50#include "HepPolyhedronProcessor.h"
[831]51#include "G4NURBS.hh"
52// #include "G4NURBSbox.hh"
53
54#include <sstream>
55
56///////////////////////////////////////////////////////////////////
57//
58// Transfer all data members to G4BooleanSolid which is responsible
59// for them. pName will be in turn sent to G4VSolid
60
61G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
62 G4VSolid* pSolidA ,
63 G4VSolid* pSolidB )
64 : G4BooleanSolid(pName,pSolidA,pSolidB)
65{
66}
67
68///////////////////////////////////////////////////////////////
69//
70// Constructor
71
72G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
73 G4VSolid* pSolidA ,
74 G4VSolid* pSolidB ,
75 G4RotationMatrix* rotMatrix,
76 const G4ThreeVector& transVector )
77 : G4BooleanSolid(pName,pSolidA,pSolidB,rotMatrix,transVector)
78{
79}
80
81///////////////////////////////////////////////////////////////
82//
83// Constructor
84
85G4SubtractionSolid::G4SubtractionSolid( 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
98G4SubtractionSolid::G4SubtractionSolid( __void__& a )
99 : G4BooleanSolid(a)
100{
101}
102
103///////////////////////////////////////////////////////////////
104//
105// Destructor
106
107G4SubtractionSolid::~G4SubtractionSolid()
108{
109}
110
111///////////////////////////////////////////////////////////////
112//
113// CalculateExtent
114
115G4bool
116G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
117 const G4VoxelLimits& pVoxelLimit,
118 const G4AffineTransform& pTransform,
119 G4double& pMin,
120 G4double& pMax ) const
121{
122 // Since we cannot be sure how much the second solid subtracts
123 // from the first, we must use the first solid's extent!
124
125 return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
126 pTransform, pMin, pMax );
127}
128
129/////////////////////////////////////////////////////
130//
131// Touching ? Empty subtraction ?
132
133EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const
134{
135 EInside positionA = fPtrSolidA->Inside(p);
136 if (positionA == kOutside) return kOutside;
137
138 EInside positionB = fPtrSolidB->Inside(p);
139
140 if(positionA == kInside && positionB == kOutside)
141 {
142 return kInside ;
143 }
144 else
145 {
146 if(( positionA == kInside && positionB == kSurface) ||
147 ( positionB == kOutside && positionA == kSurface) ||
148 ( positionA == kSurface && positionB == kSurface &&
149 ( fPtrSolidA->SurfaceNormal(p) -
150 fPtrSolidB->SurfaceNormal(p) ).mag2() >
151 1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
152 {
153 return kSurface;
154 }
155 else
156 {
157 return kOutside;
158 }
159 }
160}
161
162//////////////////////////////////////////////////////////////
163//
164// SurfaceNormal
165
166G4ThreeVector
167G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const
168{
169 G4ThreeVector normal;
170 if( Inside(p) == kOutside )
171 {
172#ifdef G4BOOLDEBUG
173 G4cout << "WARNING - Invalid call [1] in "
174 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
175 << " Point p is inside !" << G4endl;
176 G4cout << " p = " << p << G4endl;
177 G4cerr << "WARNING - Invalid call [1] in "
178 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
179 << " Point p is inside !" << G4endl;
180 G4cerr << " p = " << p << G4endl;
181#endif
182 }
183 else
184 {
185 if( fPtrSolidA->Inside(p) == kSurface &&
186 fPtrSolidB->Inside(p) != kInside )
187 {
188 normal = fPtrSolidA->SurfaceNormal(p) ;
189 }
190 else if( fPtrSolidA->Inside(p) == kInside &&
191 fPtrSolidB->Inside(p) != kOutside )
192 {
193 normal = -fPtrSolidB->SurfaceNormal(p) ;
194 }
195 else
196 {
197 if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) )
198 {
199 normal = fPtrSolidA->SurfaceNormal(p) ;
200 }
201 else
202 {
203 normal = -fPtrSolidB->SurfaceNormal(p) ;
204 }
205#ifdef G4BOOLDEBUG
206 if(Inside(p) == kInside)
207 {
208 G4cout << "WARNING - Invalid call [2] in "
209 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
210 << " Point p is inside !" << G4endl;
211 G4cout << " p = " << p << G4endl;
212 G4cerr << "WARNING - Invalid call [2] in "
213 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
214 << " Point p is inside !" << G4endl;
215 G4cerr << " p = " << p << G4endl;
216 }
217#endif
218 }
219 }
220 return normal;
221}
222
223/////////////////////////////////////////////////////////////
224//
225// The same algorithm as in DistanceToIn(p)
226
227G4double
228G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p,
229 const G4ThreeVector& v ) const
230{
231 G4double dist = 0.0,disTmp = 0.0 ;
232
233#ifdef G4BOOLDEBUG
234 if( Inside(p) == kInside )
235 {
236 G4cout << "WARNING - Invalid call in "
237 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
238 << " Point p is inside !" << G4endl;
239 G4cout << " p = " << p << G4endl;
240 G4cout << " v = " << v << G4endl;
241 G4cerr << "WARNING - Invalid call in "
242 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
243 << " Point p is inside !" << G4endl;
244 G4cerr << " p = " << p << G4endl;
245 G4cerr << " v = " << v << G4endl;
246 }
247#endif
248
249 // if( // ( fPtrSolidA->Inside(p) != kOutside) && // case1:p in both A&B
250 if ( fPtrSolidB->Inside(p) != kOutside ) // start: out of B
251 {
252 dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ;
253
254 if( fPtrSolidA->Inside(p+dist*v) != kInside )
255 {
256 do
257 {
258 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
259
260 if(disTmp == kInfinity)
261 {
262 return kInfinity ;
263 }
264 dist += disTmp ;
265
266 if( Inside(p+dist*v) == kOutside )
267 {
268 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
269 dist += disTmp ;
270 }
271 }
272 while( Inside(p+dist*v) == kOutside ) ;
273 }
274 }
275 else // p outside A, start in A
276 {
277 dist = fPtrSolidA->DistanceToIn(p,v) ;
278
279 if( dist == kInfinity ) // past A, hence past A\B
280 {
281 return kInfinity ;
282 }
283 else
284 {
285 while( Inside(p+dist*v) == kOutside ) // pushing loop
286 {
287 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
288 dist += disTmp ;
289
290 if( Inside(p+dist*v) == kOutside )
291 {
292 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
293
294 if(disTmp == kInfinity) // past A, hence past A\B
295 {
296 return kInfinity ;
297 }
298 dist += disTmp ;
299 }
300 }
301 }
302 }
303
304 return dist ;
305}
306
307////////////////////////////////////////////////////////
308//
309// Approximate nearest distance from the point p to the intersection of
310// two solids. It is usually underestimated from the point of view of
311// isotropic safety
312
313G4double
314G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p ) const
315{
316 G4double dist=0.0;
317
318#ifdef G4BOOLDEBUG
319 if( Inside(p) == kInside )
320 {
321 G4cout << "WARNING - Invalid call in "
322 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
323 << " Point p is inside !" << G4endl;
324 G4cout << " p = " << p << G4endl;
325 G4cerr << "WARNING - Invalid call in "
326 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
327 << " Point p is inside !" << G4endl;
328 G4cerr << " p = " << p << G4endl;
329 }
330#endif
331
332 if( ( fPtrSolidA->Inside(p) != kOutside) && // case 1
333 ( fPtrSolidB->Inside(p) != kOutside) )
334 {
335 dist= fPtrSolidB->DistanceToOut(p) ;
336 }
337 else
338 {
339 dist= fPtrSolidA->DistanceToIn(p) ;
340 }
341
342 return dist;
343}
344
345//////////////////////////////////////////////////////////
346//
347// The same algorithm as DistanceToOut(p)
348
349G4double
350G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p,
351 const G4ThreeVector& v,
352 const G4bool calcNorm,
353 G4bool *validNorm,
354 G4ThreeVector *n ) const
355{
356#ifdef G4BOOLDEBUG
357 if( Inside(p) == kOutside )
358 {
359 G4cout << "Position:" << G4endl << G4endl;
360 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
361 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
362 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
363 G4cout << "Direction:" << G4endl << G4endl;
364 G4cout << "v.x() = " << v.x() << G4endl;
365 G4cout << "v.y() = " << v.y() << G4endl;
366 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
367 G4cout << "WARNING - Invalid call in "
368 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
369 << " Point p is outside !" << G4endl;
370 G4cout << " p = " << p << G4endl;
371 G4cout << " v = " << v << G4endl;
372 G4cerr << "WARNING - Invalid call in "
373 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
374 << " Point p is outside !" << G4endl;
375 G4cerr << " p = " << p << G4endl;
376 G4cerr << " v = " << v << G4endl;
377 }
378#endif
379
380 G4double distout;
381 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
382 G4double distB = fPtrSolidB->DistanceToIn(p,v) ;
383 if(distB < distA)
384 {
385 if(calcNorm)
386 {
387 *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
388 *validNorm = false ;
389 }
390 distout= distB ;
391 }
392 else
393 {
394 distout= distA ;
395 }
396 return distout;
397}
398
399//////////////////////////////////////////////////////////////
400//
401// Inverted algorithm of DistanceToIn(p)
402
403G4double
404G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const
405{
406 G4double dist=0.0;
407
408 if( Inside(p) == kOutside )
409 {
410#ifdef G4BOOLDEBUG
411 G4cout << "WARNING - Invalid call in "
412 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
413 << " Point p is outside" << G4endl;
414 G4cout << " p = " << p << G4endl;
415 G4cerr << "WARNING - Invalid call in "
416 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
417 << " Point p is outside" << G4endl;
418 G4cerr << " p = " << p << G4endl;
419#endif
420 }
421 else
422 {
423 dist= std::min(fPtrSolidA->DistanceToOut(p),
424 fPtrSolidB->DistanceToIn(p) ) ;
425 }
426 return dist;
427}
428
429//////////////////////////////////////////////////////////////
430//
431//
432
433G4GeometryType G4SubtractionSolid::GetEntityType() const
434{
435 return G4String("G4SubtractionSolid");
436}
437
438//////////////////////////////////////////////////////////////
439//
440//
441
442void
443G4SubtractionSolid::ComputeDimensions( G4VPVParameterisation*,
444 const G4int,
445 const G4VPhysicalVolume* )
446{
447}
448
449/////////////////////////////////////////////////
450//
451//
452
453void
454G4SubtractionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
455{
456 scene.AddSolid (*this);
457}
458
459////////////////////////////////////////////////////
460//
461//
462
463G4Polyhedron*
464G4SubtractionSolid::CreatePolyhedron () const
465{
[1315]466 HepPolyhedronProcessor processor;
467 // Stack components and components of components recursively
468 // See G4BooleanSolid::StackPolyhedron
469 G4Polyhedron* top = StackPolyhedron(processor, this);
470 G4Polyhedron* result = new G4Polyhedron(*top);
471 if (processor.execute(*result)) { return result; }
472 else { return 0; }
[831]473}
474
475/////////////////////////////////////////////////////////
476//
477//
478
479G4NURBS*
480G4SubtractionSolid::CreateNURBS () const
481{
482 // Take into account boolean operation - see CreatePolyhedron.
483 // return new G4NURBSbox (1.0, 1.0, 1.0);
484 return 0;
485}
Note: See TracBrowser for help on using the repository browser.