source: trunk/source/geometry/solids/Boolean/src/G4UnionSolid.cc@ 1279

Last change on this file since 1279 was 1228, checked in by garnier, 16 years ago

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