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

Last change on this file since 1038 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
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: G4SubtractionSolid.cc,v 1.31 2007/10/23 14:42:31 grichine Exp $
28// GEANT4 tag $Name: HEAD $
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"
50#include "G4NURBS.hh"
51// #include "G4NURBSbox.hh"
52
53#include <sstream>
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
60G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
61 G4VSolid* pSolidA ,
62 G4VSolid* pSolidB )
63 : G4BooleanSolid(pName,pSolidA,pSolidB)
64{
65}
66
67///////////////////////////////////////////////////////////////
68//
69// Constructor
70
71G4SubtractionSolid::G4SubtractionSolid( 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// Constructor
83
84G4SubtractionSolid::G4SubtractionSolid( const G4String& pName,
85 G4VSolid* pSolidA ,
86 G4VSolid* pSolidB ,
87 const G4Transform3D& transform )
88 : G4BooleanSolid(pName,pSolidA,pSolidB,transform)
89{
90}
91
92//////////////////////////////////////////////////////////////////
93//
94// Fake default constructor - sets only member data and allocates memory
95// for usage restricted to object persistency.
96
97G4SubtractionSolid::G4SubtractionSolid( __void__& a )
98 : G4BooleanSolid(a)
99{
100}
101
102///////////////////////////////////////////////////////////////
103//
104// Destructor
105
106G4SubtractionSolid::~G4SubtractionSolid()
107{
108}
109
110///////////////////////////////////////////////////////////////
111//
112// CalculateExtent
113
114G4bool
115G4SubtractionSolid::CalculateExtent( const EAxis pAxis,
116 const G4VoxelLimits& pVoxelLimit,
117 const G4AffineTransform& pTransform,
118 G4double& pMin,
119 G4double& pMax ) const
120{
121 // Since we cannot be sure how much the second solid subtracts
122 // from the first, we must use the first solid's extent!
123
124 return fPtrSolidA->CalculateExtent( pAxis, pVoxelLimit,
125 pTransform, pMin, pMax );
126}
127
128/////////////////////////////////////////////////////
129//
130// Touching ? Empty subtraction ?
131
132EInside G4SubtractionSolid::Inside( const G4ThreeVector& p ) const
133{
134 EInside positionA = fPtrSolidA->Inside(p);
135 if (positionA == kOutside) return kOutside;
136
137 EInside positionB = fPtrSolidB->Inside(p);
138
139 if(positionA == kInside && positionB == kOutside)
140 {
141 return kInside ;
142 }
143 else
144 {
145 if(( positionA == kInside && positionB == kSurface) ||
146 ( positionB == kOutside && positionA == kSurface) ||
147 ( positionA == kSurface && positionB == kSurface &&
148 ( fPtrSolidA->SurfaceNormal(p) -
149 fPtrSolidB->SurfaceNormal(p) ).mag2() >
150 1000*G4GeometryTolerance::GetInstance()->GetRadialTolerance() ) )
151 {
152 return kSurface;
153 }
154 else
155 {
156 return kOutside;
157 }
158 }
159}
160
161//////////////////////////////////////////////////////////////
162//
163// SurfaceNormal
164
165G4ThreeVector
166G4SubtractionSolid::SurfaceNormal( const G4ThreeVector& p ) const
167{
168 G4ThreeVector normal;
169 if( Inside(p) == kOutside )
170 {
171#ifdef G4BOOLDEBUG
172 G4cout << "WARNING - Invalid call [1] in "
173 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
174 << " Point p is inside !" << G4endl;
175 G4cout << " p = " << p << G4endl;
176 G4cerr << "WARNING - Invalid call [1] in "
177 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
178 << " Point p is inside !" << G4endl;
179 G4cerr << " p = " << p << G4endl;
180#endif
181 }
182 else
183 {
184 if( fPtrSolidA->Inside(p) == kSurface &&
185 fPtrSolidB->Inside(p) != kInside )
186 {
187 normal = fPtrSolidA->SurfaceNormal(p) ;
188 }
189 else if( fPtrSolidA->Inside(p) == kInside &&
190 fPtrSolidB->Inside(p) != kOutside )
191 {
192 normal = -fPtrSolidB->SurfaceNormal(p) ;
193 }
194 else
195 {
196 if ( fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToIn(p) )
197 {
198 normal = fPtrSolidA->SurfaceNormal(p) ;
199 }
200 else
201 {
202 normal = -fPtrSolidB->SurfaceNormal(p) ;
203 }
204#ifdef G4BOOLDEBUG
205 if(Inside(p) == kInside)
206 {
207 G4cout << "WARNING - Invalid call [2] in "
208 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
209 << " Point p is inside !" << G4endl;
210 G4cout << " p = " << p << G4endl;
211 G4cerr << "WARNING - Invalid call [2] in "
212 << "G4SubtractionSolid::SurfaceNormal(p)" << G4endl
213 << " Point p is inside !" << G4endl;
214 G4cerr << " p = " << p << G4endl;
215 }
216#endif
217 }
218 }
219 return normal;
220}
221
222/////////////////////////////////////////////////////////////
223//
224// The same algorithm as in DistanceToIn(p)
225
226G4double
227G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p,
228 const G4ThreeVector& v ) const
229{
230 G4double dist = 0.0,disTmp = 0.0 ;
231
232#ifdef G4BOOLDEBUG
233 if( Inside(p) == kInside )
234 {
235 G4cout << "WARNING - Invalid call in "
236 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
237 << " Point p is inside !" << G4endl;
238 G4cout << " p = " << p << G4endl;
239 G4cout << " v = " << v << G4endl;
240 G4cerr << "WARNING - Invalid call in "
241 << "G4SubtractionSolid::DistanceToIn(p,v)" << G4endl
242 << " Point p is inside !" << G4endl;
243 G4cerr << " p = " << p << G4endl;
244 G4cerr << " v = " << v << G4endl;
245 }
246#endif
247
248 // if( // ( fPtrSolidA->Inside(p) != kOutside) && // case1:p in both A&B
249 if ( fPtrSolidB->Inside(p) != kOutside ) // start: out of B
250 {
251 dist = fPtrSolidB->DistanceToOut(p,v) ; // ,calcNorm,validNorm,n) ;
252
253 if( fPtrSolidA->Inside(p+dist*v) != kInside )
254 {
255 do
256 {
257 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
258
259 if(disTmp == kInfinity)
260 {
261 return kInfinity ;
262 }
263 dist += disTmp ;
264
265 if( Inside(p+dist*v) == kOutside )
266 {
267 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
268 dist += disTmp ;
269 }
270 }
271 while( Inside(p+dist*v) == kOutside ) ;
272 }
273 }
274 else // p outside A, start in A
275 {
276 dist = fPtrSolidA->DistanceToIn(p,v) ;
277
278 if( dist == kInfinity ) // past A, hence past A\B
279 {
280 return kInfinity ;
281 }
282 else
283 {
284 while( Inside(p+dist*v) == kOutside ) // pushing loop
285 {
286 disTmp = fPtrSolidB->DistanceToOut(p+dist*v,v) ;
287 dist += disTmp ;
288
289 if( Inside(p+dist*v) == kOutside )
290 {
291 disTmp = fPtrSolidA->DistanceToIn(p+dist*v,v) ;
292
293 if(disTmp == kInfinity) // past A, hence past A\B
294 {
295 return kInfinity ;
296 }
297 dist += disTmp ;
298 }
299 }
300 }
301 }
302
303 return dist ;
304}
305
306////////////////////////////////////////////////////////
307//
308// Approximate nearest distance from the point p to the intersection of
309// two solids. It is usually underestimated from the point of view of
310// isotropic safety
311
312G4double
313G4SubtractionSolid::DistanceToIn( const G4ThreeVector& p ) const
314{
315 G4double dist=0.0;
316
317#ifdef G4BOOLDEBUG
318 if( Inside(p) == kInside )
319 {
320 G4cout << "WARNING - Invalid call in "
321 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
322 << " Point p is inside !" << G4endl;
323 G4cout << " p = " << p << G4endl;
324 G4cerr << "WARNING - Invalid call in "
325 << "G4SubtractionSolid::DistanceToIn(p)" << G4endl
326 << " Point p is inside !" << G4endl;
327 G4cerr << " p = " << p << G4endl;
328 }
329#endif
330
331 if( ( fPtrSolidA->Inside(p) != kOutside) && // case 1
332 ( fPtrSolidB->Inside(p) != kOutside) )
333 {
334 dist= fPtrSolidB->DistanceToOut(p) ;
335 }
336 else
337 {
338 dist= fPtrSolidA->DistanceToIn(p) ;
339 }
340
341 return dist;
342}
343
344//////////////////////////////////////////////////////////
345//
346// The same algorithm as DistanceToOut(p)
347
348G4double
349G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p,
350 const G4ThreeVector& v,
351 const G4bool calcNorm,
352 G4bool *validNorm,
353 G4ThreeVector *n ) const
354{
355#ifdef G4BOOLDEBUG
356 if( Inside(p) == kOutside )
357 {
358 G4cout << "Position:" << G4endl << G4endl;
359 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
360 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
361 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
362 G4cout << "Direction:" << G4endl << G4endl;
363 G4cout << "v.x() = " << v.x() << G4endl;
364 G4cout << "v.y() = " << v.y() << G4endl;
365 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
366 G4cout << "WARNING - Invalid call in "
367 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
368 << " Point p is outside !" << G4endl;
369 G4cout << " p = " << p << G4endl;
370 G4cout << " v = " << v << G4endl;
371 G4cerr << "WARNING - Invalid call in "
372 << "G4SubtractionSolid::DistanceToOut(p,v)" << G4endl
373 << " Point p is outside !" << G4endl;
374 G4cerr << " p = " << p << G4endl;
375 G4cerr << " v = " << v << G4endl;
376 }
377#endif
378
379 G4double distout;
380 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,validNorm,n) ;
381 G4double distB = fPtrSolidB->DistanceToIn(p,v) ;
382 if(distB < distA)
383 {
384 if(calcNorm)
385 {
386 *n = -(fPtrSolidB->SurfaceNormal(p+distB*v)) ;
387 *validNorm = false ;
388 }
389 distout= distB ;
390 }
391 else
392 {
393 distout= distA ;
394 }
395 return distout;
396}
397
398//////////////////////////////////////////////////////////////
399//
400// Inverted algorithm of DistanceToIn(p)
401
402G4double
403G4SubtractionSolid::DistanceToOut( const G4ThreeVector& p ) const
404{
405 G4double dist=0.0;
406
407 if( Inside(p) == kOutside )
408 {
409#ifdef G4BOOLDEBUG
410 G4cout << "WARNING - Invalid call in "
411 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
412 << " Point p is outside" << G4endl;
413 G4cout << " p = " << p << G4endl;
414 G4cerr << "WARNING - Invalid call in "
415 << "G4SubtractionSolid::DistanceToOut(p)" << G4endl
416 << " Point p is outside" << G4endl;
417 G4cerr << " p = " << p << G4endl;
418#endif
419 }
420 else
421 {
422 dist= std::min(fPtrSolidA->DistanceToOut(p),
423 fPtrSolidB->DistanceToIn(p) ) ;
424 }
425 return dist;
426}
427
428//////////////////////////////////////////////////////////////
429//
430//
431
432G4GeometryType G4SubtractionSolid::GetEntityType() const
433{
434 return G4String("G4SubtractionSolid");
435}
436
437//////////////////////////////////////////////////////////////
438//
439//
440
441void
442G4SubtractionSolid::ComputeDimensions( G4VPVParameterisation*,
443 const G4int,
444 const G4VPhysicalVolume* )
445{
446}
447
448/////////////////////////////////////////////////
449//
450//
451
452void
453G4SubtractionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
454{
455 scene.AddSolid (*this);
456}
457
458////////////////////////////////////////////////////
459//
460//
461
462G4Polyhedron*
463G4SubtractionSolid::CreatePolyhedron () const
464{
465 G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
466 G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
467 if (pA && pB)
468 {
469 G4Polyhedron* resultant = new G4Polyhedron (pA->subtract(*pB));
470 return resultant;
471 }
472 else
473 {
474 std::ostringstream oss;
475 oss << "Solid - " << GetName()
476 << " - one of the Boolean components has no" << G4endl
477 << " corresponding polyhedron. Returning NULL !";
478 G4Exception("G4SubtractionSolid::CreatePolyhedron()", "InvalidSetup",
479 JustWarning, oss.str().c_str());
480 return 0;
481 }
482}
483
484/////////////////////////////////////////////////////////
485//
486//
487
488G4NURBS*
489G4SubtractionSolid::CreateNURBS () const
490{
491 // Take into account boolean operation - see CreatePolyhedron.
492 // return new G4NURBSbox (1.0, 1.0, 1.0);
493 return 0;
494}
Note: See TracBrowser for help on using the repository browser.