Context Navigation

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

Last change on this file since 1298 was 1228, checked in by garnier, 15 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
59	G4UnionSolid:: G4UnionSolid( const G4String& pName,
60	G4VSolid* pSolidA ,
61	G4VSolid* pSolidB )
62	: G4BooleanSolid(pName,pSolidA,pSolidB)
63	{
64	}
65
66	/////////////////////////////////////////////////////////////////////
67	//
68	// Constructor
69
70	G4UnionSolid::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
84	G4UnionSolid::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
97	G4UnionSolid::G4UnionSolid( __void__& a )
98	: G4BooleanSolid(a)
99	{
100	}
101
102	///////////////////////////////////////////////////////////
103	//
104	// Destructor
105
106	G4UnionSolid::~G4UnionSolid()
107	{
108	}
109
110	///////////////////////////////////////////////////////////////
111	//
112	//
113
114	G4bool
115	G4UnionSolid::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
146	EInside 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
174	G4ThreeVector
175	G4UnionSolid::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
226	G4double
227	G4UnionSolid::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
255	G4double
256	G4UnionSolid::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
282	G4double
283	G4UnionSolid::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
372	G4double
373	G4UnionSolid::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
423	G4GeometryType G4UnionSolid::GetEntityType() const
424	{
425	return G4String("G4UnionSolid");
426	}
427
428	//////////////////////////////////////////////////////////////
429	//
430	//
431
432	void
433	G4UnionSolid::ComputeDimensions( G4VPVParameterisation*,
434	const G4int,
435	const G4VPhysicalVolume* )
436	{
437	}
438
439	/////////////////////////////////////////////////
440	//
441	//
442
443	void
444	G4UnionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
445	{
446	scene.AddSolid (*this);
447	}
448
449	////////////////////////////////////////////////////
450	//
451	//
452
453	G4Polyhedron*
454	G4UnionSolid::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
475	G4NURBS*
476	G4UnionSolid::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.

Download in other formats: