Context Navigation

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

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

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

Download in other formats: