source: trunk/source/geometry/solids/Boolean/src/G4IntersectionSolid.cc@ 883

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

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