Context Navigation

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

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

Download in other formats: