Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: