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