source: trunk/source/geometry/solids/Boolean/src/G4IntersectionSolid.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: 15.2 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: G4IntersectionSolid.cc,v 1.32 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// 17.02.05 V.Grichine: bug was fixed in DistanceToIn(p,v) based on algorithm
35// proposed by Dino Bazzacco <dino.bazzacco@pd.infn.it>
36// 29.05.01 V.Grichine: bug was fixed in DistanceToIn(p,v)
37// 16.03.01 V.Grichine: modifications in CalculateExtent() and Inside()
38// 29.07.99 V.Grichine: modifications in DistanceToIn(p,v)
39// 12.09.98 V.Grichine: first implementation
40//
41// --------------------------------------------------------------------
42
43#include "G4IntersectionSolid.hh"
44
45#include <sstream>
46
47#include "G4VoxelLimits.hh"
48#include "G4VPVParameterisation.hh"
49
50#include "G4VGraphicsScene.hh"
51#include "G4Polyhedron.hh"
52#include "HepPolyhedronProcessor.h"
53#include "G4NURBS.hh"
54// #include "G4NURBSbox.hh"
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//
61
62G4IntersectionSolid::G4IntersectionSolid( const G4String& pName,
63 G4VSolid* pSolidA ,
64 G4VSolid* pSolidB )
65 : G4BooleanSolid(pName,pSolidA,pSolidB)
66{
67}
68
69///////////////////////////////////////////////////////////////////
70//
71
72G4IntersectionSolid::G4IntersectionSolid( 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//
84
85G4IntersectionSolid::G4IntersectionSolid( 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
98G4IntersectionSolid::G4IntersectionSolid( __void__& a )
99 : G4BooleanSolid(a)
100{
101}
102
103///////////////////////////////////////////////////////////////
104//
105//
106
107G4IntersectionSolid::~G4IntersectionSolid()
108{
109}
110
111///////////////////////////////////////////////////////////////
112//
113//
114
115G4bool
116G4IntersectionSolid::CalculateExtent(const EAxis pAxis,
117 const G4VoxelLimits& pVoxelLimit,
118 const G4AffineTransform& pTransform,
119 G4double& pMin,
120 G4double& pMax) const
121{
122 G4bool retA, retB, out;
123 G4double minA, minB, maxA, maxB;
124
125 retA = fPtrSolidA
126 ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minA, maxA);
127 retB = fPtrSolidB
128 ->CalculateExtent( pAxis, pVoxelLimit, pTransform, minB, maxB);
129
130 if( retA && retB )
131 {
132 pMin = std::max( minA, minB );
133 pMax = std::min( maxA, maxB );
134 out = (pMax > pMin); // true;
135#ifdef G4BOOLDEBUG
136 // G4cout.precision(16);
137 // G4cout<<"pMin = "<<pMin<<"; pMax = "<<pMax<<G4endl;
138#endif
139 }
140 else out = false;
141
142 return out; // It exists in this slice only if both exist in it.
143}
144
145/////////////////////////////////////////////////////
146//
147// Touching ? Empty intersection ?
148
149EInside G4IntersectionSolid::Inside(const G4ThreeVector& p) const
150{
151 EInside positionA = fPtrSolidA->Inside(p) ;
152
153 if( positionA == kOutside ) return kOutside ;
154
155 EInside positionB = fPtrSolidB->Inside(p) ;
156
157 if(positionA == kInside && positionB == kInside)
158 {
159 return kInside ;
160 }
161 else
162 {
163 if((positionA == kInside && positionB == kSurface) ||
164 (positionB == kInside && positionA == kSurface) ||
165 (positionA == kSurface && positionB == kSurface) )
166 {
167 return kSurface ;
168 }
169 else
170 {
171 return kOutside ;
172 }
173 }
174}
175
176//////////////////////////////////////////////////////////////
177//
178
179G4ThreeVector
180G4IntersectionSolid::SurfaceNormal( const G4ThreeVector& p ) const
181{
182 G4ThreeVector normal;
183 EInside insideA, insideB;
184
185 insideA= fPtrSolidA->Inside(p);
186 insideB= fPtrSolidB->Inside(p);
187
188#ifdef G4BOOLDEBUG
189 if( (insideA == kOutside) || (insideB == kOutside) )
190 {
191 G4cout << "WARNING - Invalid call in "
192 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
193 << " Point p is outside !" << G4endl;
194 G4cout << " p = " << p << G4endl;
195 G4cerr << "WARNING - Invalid call in "
196 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
197 << " Point p is outside !" << G4endl;
198 G4cerr << " p = " << p << G4endl;
199 }
200#endif
201
202 // OLD: if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) )
203
204 // On the surface of both is difficult ... treat it like on A now!
205 //
206 // if( (insideA == kSurface) && (insideB == kSurface) )
207 // normal= fPtrSolidA->SurfaceNormal(p) ;
208 // else
209 if( insideA == kSurface )
210 {
211 normal= fPtrSolidA->SurfaceNormal(p) ;
212 }
213 else if( insideB == kSurface )
214 {
215 normal= fPtrSolidB->SurfaceNormal(p) ;
216 }
217 // We are on neither surface, so we should generate an exception
218 else
219 {
220 if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) )
221 normal= fPtrSolidA->SurfaceNormal(p) ;
222 else
223 normal= fPtrSolidB->SurfaceNormal(p) ;
224#ifdef G4BOOLDEBUG
225 G4cout << "WARNING - Invalid call in "
226 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
227 << " Point p is out of surface !" << G4endl;
228 G4cout << " p = " << p << G4endl;
229 G4cerr << "WARNING - Invalid call in "
230 << "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
231 << " Point p is out of surface !" << G4endl;
232 G4cerr << " p = " << p << G4endl;
233#endif
234 }
235
236 return normal;
237}
238
239/////////////////////////////////////////////////////////////
240//
241// The same algorithm as in DistanceToIn(p)
242
243G4double
244G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p,
245 const G4ThreeVector& v ) const
246{
247 G4double dist = 0.0;
248 if( Inside(p) == kInside )
249 {
250#ifdef G4BOOLDEBUG
251 G4cout << "WARNING - Invalid call in "
252 << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
253 << " Point p is inside !" << G4endl;
254 G4cout << " p = " << p << G4endl;
255 G4cout << " v = " << v << G4endl;
256 G4cerr << "WARNING - Invalid call in "
257 << "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
258 << " Point p is inside !" << G4endl;
259 G4cerr << " p = " << p << G4endl;
260 G4cerr << " v = " << v << G4endl;
261#endif
262 }
263 else // if( Inside(p) == kSurface )
264 {
265 EInside wA = fPtrSolidA->Inside(p);
266 EInside wB = fPtrSolidB->Inside(p);
267
268 G4ThreeVector pA = p, pB = p;
269 G4double dA = 0., dA1=0., dA2=0.;
270 G4double dB = 0., dB1=0., dB2=0.;
271 G4bool doA = true, doB = true;
272
273 while(true)
274 {
275 if(doA)
276 {
277 // find next valid range for A
278
279 dA1 = 0.;
280
281 if( wA != kInside )
282 {
283 dA1 = fPtrSolidA->DistanceToIn(pA, v);
284
285 if( dA1 == kInfinity ) return kInfinity;
286
287 pA += dA1*v;
288 }
289 dA2 = dA1 + fPtrSolidA->DistanceToOut(pA, v);
290 }
291 dA1 += dA;
292 dA2 += dA;
293
294 if(doB)
295 {
296 // find next valid range for B
297
298 dB1 = 0.;
299 if(wB != kInside)
300 {
301 dB1 = fPtrSolidB->DistanceToIn(pB, v);
302
303 if(dB1 == kInfinity) return kInfinity;
304
305 pB += dB1*v;
306 }
307 dB2 = dB1 + fPtrSolidB->DistanceToOut(pB, v);
308 }
309 dB1 += dB;
310 dB2 += dB;
311
312 // check if they overlap
313
314 if( dA1 < dB1 )
315 {
316 if( dB1 < dA2 ) return dB1;
317
318 dA = dA2;
319 pA = p + dA*v; // continue from here
320 wA = kSurface;
321 doA = true;
322 doB = false;
323 }
324 else
325 {
326 if( dA1 < dB2 ) return dA1;
327
328 dB = dB2;
329 pB = p + dB*v; // continue from here
330 wB = kSurface;
331 doB = true;
332 doA = false;
333 }
334 }
335 }
336 return dist ;
337}
338
339////////////////////////////////////////////////////////
340//
341// Approximate nearest distance from the point p to the intersection of
342// two solids
343
344G4double
345G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p) const
346{
347#ifdef G4BOOLDEBUG
348 if( Inside(p) == kInside )
349 {
350 G4cout << "WARNING - Invalid call in "
351 << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
352 << " Point p is inside !" << G4endl;
353 G4cout << " p = " << p << G4endl;
354 G4cerr << "WARNING - Invalid call in "
355 << "G4IntersectionSolid::DistanceToIn(p)" << G4endl
356 << " Point p is inside !" << G4endl;
357 G4cerr << " p = " << p << G4endl;
358 }
359#endif
360 EInside sideA = fPtrSolidA->Inside(p) ;
361 EInside sideB = fPtrSolidB->Inside(p) ;
362 G4double dist=0.0 ;
363
364 if( sideA != kInside && sideB != kOutside )
365 {
366 dist = fPtrSolidA->DistanceToIn(p) ;
367 }
368 else
369 {
370 if( sideB != kInside && sideA != kOutside )
371 {
372 dist = fPtrSolidB->DistanceToIn(p) ;
373 }
374 else
375 {
376 dist = std::min(fPtrSolidA->DistanceToIn(p),
377 fPtrSolidB->DistanceToIn(p) ) ;
378 }
379 }
380 return dist ;
381}
382
383//////////////////////////////////////////////////////////
384//
385// The same algorithm as DistanceToOut(p)
386
387G4double
388G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p,
389 const G4ThreeVector& v,
390 const G4bool calcNorm,
391 G4bool *validNorm,
392 G4ThreeVector *n ) const
393{
394 G4bool validNormA, validNormB;
395 G4ThreeVector nA, nB;
396
397#ifdef G4BOOLDEBUG
398 if( Inside(p) == kOutside )
399 {
400 G4cout << "Position:" << G4endl << G4endl;
401 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
402 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
403 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
404 G4cout << "Direction:" << G4endl << G4endl;
405 G4cout << "v.x() = " << v.x() << G4endl;
406 G4cout << "v.y() = " << v.y() << G4endl;
407 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
408 G4cout << "WARNING - Invalid call in "
409 << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
410 << " Point p is outside !" << G4endl;
411 G4cout << " p = " << p << G4endl;
412 G4cout << " v = " << v << G4endl;
413 G4cerr << "WARNING - Invalid call in "
414 << "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
415 << " Point p is outside !" << G4endl;
416 G4cerr << " p = " << p << G4endl;
417 G4cerr << " v = " << v << G4endl;
418 }
419#endif
420 G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,&validNormA,&nA) ;
421 G4double distB = fPtrSolidB->DistanceToOut(p,v,calcNorm,&validNormB,&nB) ;
422
423 G4double dist = std::min(distA,distB) ;
424
425 if( calcNorm )
426 {
427 if ( distA < distB )
428 {
429 *validNorm = validNormA;
430 *n = nA;
431 }
432 else
433 {
434 *validNorm = validNormB;
435 *n = nB;
436 }
437 }
438
439 return dist ;
440}
441
442//////////////////////////////////////////////////////////////
443//
444// Inverted algorithm of DistanceToIn(p)
445
446G4double
447G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p ) const
448{
449#ifdef G4BOOLDEBUG
450 if( Inside(p) == kOutside )
451 {
452 G4cout << "WARNING - Invalid call in "
453 << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
454 << " Point p is outside !" << G4endl;
455 G4cout << " p = " << p << G4endl;
456 G4cerr << "WARNING - Invalid call in "
457 << "G4IntersectionSolid::DistanceToOut(p)" << G4endl
458 << " Point p is outside !" << G4endl;
459 G4cerr << " p = " << p << G4endl;
460 }
461#endif
462
463 return std::min(fPtrSolidA->DistanceToOut(p),
464 fPtrSolidB->DistanceToOut(p) ) ;
465
466}
467
468//////////////////////////////////////////////////////////////
469//
470//
471
472void
473G4IntersectionSolid::ComputeDimensions( G4VPVParameterisation*,
474 const G4int,
475 const G4VPhysicalVolume* )
476{
477}
478
479/////////////////////////////////////////////////
480//
481//
482
483G4GeometryType G4IntersectionSolid::GetEntityType() const
484{
485 return G4String("G4IntersectionSolid");
486}
487
488/////////////////////////////////////////////////
489//
490//
491
492void
493G4IntersectionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
494{
495 scene.AddSolid (*this);
496}
497
498////////////////////////////////////////////////////
499//
500//
501
502G4Polyhedron*
503G4IntersectionSolid::CreatePolyhedron () const
504{
505 HepPolyhedronProcessor processor;
506 // Stack components and components of components recursively
507 // See G4BooleanSolid::StackPolyhedron
508 G4Polyhedron* top = StackPolyhedron(processor, this);
509 G4Polyhedron* result = new G4Polyhedron(*top);
510 if (processor.execute(*result)) { return result; }
511 else { return 0; }
512}
513
514/////////////////////////////////////////////////////////
515//
516//
517
518G4NURBS*
519G4IntersectionSolid::CreateNURBS () const
520{
521 // Take into account boolean operation - see CreatePolyhedron.
522 // return new G4NURBSbox (1.0, 1.0, 1.0);
523 return 0;
524}
Note: See TracBrowser for help on using the repository browser.