Context Navigation

source: trunk/source/geometry/solids/Boolean/src/G4IntersectionSolid.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: 15.3 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: G4IntersectionSolid.cc,v 1.30 2006/11/08 09:37:41 gcosmo Exp $
28	// GEANT4 tag $Name: geant4-09-03 $
29	//
30	// Implementation of methods for the class G4IntersectionSolid
31	//
32	// History:
33	//
34	// 17.02.05 V.Grichine: bug was fixed in DistanceToIn(p,v) based on algorithm
35	// proposed by Dino Bazzacco <dino.bazzacco@pd.infn.it>
36	// 29.05.01 V.Grichine: bug was fixed in DistanceToIn(p,v)
37	// 16.03.01 V.Grichine: modifications in CalculateExtent() and Inside()
38	// 29.07.99 V.Grichine: modifications in DistanceToIn(p,v)
39	// 12.09.98 V.Grichine: first implementation
40	//
41	// --------------------------------------------------------------------
42
43	#include "G4IntersectionSolid.hh"
44
45	#include <sstream>
46
47	#include "G4VoxelLimits.hh"
48	#include "G4VPVParameterisation.hh"
49
50	#include "G4VGraphicsScene.hh"
51	#include "G4Polyhedron.hh"
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
61	G4IntersectionSolid::G4IntersectionSolid( const G4String& pName,
62	G4VSolid* pSolidA ,
63	G4VSolid* pSolidB )
64	: G4BooleanSolid(pName,pSolidA,pSolidB)
65	{
66	}
67
68	///////////////////////////////////////////////////////////////////
69	//
70
71	G4IntersectionSolid::G4IntersectionSolid( 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
84	G4IntersectionSolid::G4IntersectionSolid( 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	G4IntersectionSolid::G4IntersectionSolid( __void__& a )
98	: G4BooleanSolid(a)
99	{
100	}
101
102	///////////////////////////////////////////////////////////////
103	//
104	//
105
106	G4IntersectionSolid::~G4IntersectionSolid()
107	{
108	}
109
110	///////////////////////////////////////////////////////////////
111	//
112	//
113
114	G4bool
115	G4IntersectionSolid::CalculateExtent(const EAxis pAxis,
116	const G4VoxelLimits& pVoxelLimit,
117	const G4AffineTransform& pTransform,
118	G4double& pMin,
119	G4double& pMax) const
120	{
121	G4bool retA, retB, out;
122	G4double minA, minB, maxA, maxB;
123
124	retA = fPtrSolidA
125	->CalculateExtent( pAxis, pVoxelLimit, pTransform, minA, maxA);
126	retB = fPtrSolidB
127	->CalculateExtent( pAxis, pVoxelLimit, pTransform, minB, maxB);
128
129	if( retA && retB )
130	{
131	pMin = std::max( minA, minB );
132	pMax = std::min( maxA, maxB );
133	out = (pMax > pMin); // true;
134	#ifdef G4BOOLDEBUG
135	// G4cout.precision(16);
136	// G4cout<<"pMin = "<<pMin<<"; pMax = "<<pMax<<G4endl;
137	#endif
138	}
139	else out = false;
140
141	return out; // It exists in this slice only if both exist in it.
142	}
143
144	/////////////////////////////////////////////////////
145	//
146	// Touching ? Empty intersection ?
147
148	EInside G4IntersectionSolid::Inside(const G4ThreeVector& p) const
149	{
150	EInside positionA = fPtrSolidA->Inside(p) ;
151
152	if( positionA == kOutside ) return kOutside ;
153
154	EInside positionB = fPtrSolidB->Inside(p) ;
155
156	if(positionA == kInside && positionB == kInside)
157	{
158	return kInside ;
159	}
160	else
161	{
162	if((positionA == kInside && positionB == kSurface) \|\|
163	(positionB == kInside && positionA == kSurface) \|\|
164	(positionA == kSurface && positionB == kSurface) )
165	{
166	return kSurface ;
167	}
168	else
169	{
170	return kOutside ;
171	}
172	}
173	}
174
175	//////////////////////////////////////////////////////////////
176	//
177
178	G4ThreeVector
179	G4IntersectionSolid::SurfaceNormal( const G4ThreeVector& p ) const
180	{
181	G4ThreeVector normal;
182	EInside insideA, insideB;
183
184	insideA= fPtrSolidA->Inside(p);
185	insideB= fPtrSolidB->Inside(p);
186
187	#ifdef G4BOOLDEBUG
188	if( (insideA == kOutside) \|\| (insideB == kOutside) )
189	{
190	G4cout << "WARNING - Invalid call in "
191	<< "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
192	<< " Point p is outside !" << G4endl;
193	G4cout << " p = " << p << G4endl;
194	G4cerr << "WARNING - Invalid call in "
195	<< "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
196	<< " Point p is outside !" << G4endl;
197	G4cerr << " p = " << p << G4endl;
198	}
199	#endif
200
201	// OLD: if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) )
202
203	// On the surface of both is difficult ... treat it like on A now!
204	//
205	// if( (insideA == kSurface) && (insideB == kSurface) )
206	// normal= fPtrSolidA->SurfaceNormal(p) ;
207	// else
208	if( insideA == kSurface )
209	{
210	normal= fPtrSolidA->SurfaceNormal(p) ;
211	}
212	else if( insideB == kSurface )
213	{
214	normal= fPtrSolidB->SurfaceNormal(p) ;
215	}
216	// We are on neither surface, so we should generate an exception
217	else
218	{
219	if(fPtrSolidA->DistanceToOut(p) <= fPtrSolidB->DistanceToOut(p) )
220	normal= fPtrSolidA->SurfaceNormal(p) ;
221	else
222	normal= fPtrSolidB->SurfaceNormal(p) ;
223	#ifdef G4BOOLDEBUG
224	G4cout << "WARNING - Invalid call in "
225	<< "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
226	<< " Point p is out of surface !" << G4endl;
227	G4cout << " p = " << p << G4endl;
228	G4cerr << "WARNING - Invalid call in "
229	<< "G4IntersectionSolid::SurfaceNormal(p)" << G4endl
230	<< " Point p is out of surface !" << G4endl;
231	G4cerr << " p = " << p << G4endl;
232	#endif
233	}
234
235	return normal;
236	}
237
238	/////////////////////////////////////////////////////////////
239	//
240	// The same algorithm as in DistanceToIn(p)
241
242	G4double
243	G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p,
244	const G4ThreeVector& v ) const
245	{
246	G4double dist = 0.0;
247	if( Inside(p) == kInside )
248	{
249	#ifdef G4BOOLDEBUG
250	G4cout << "WARNING - Invalid call in "
251	<< "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
252	<< " Point p is inside !" << G4endl;
253	G4cout << " p = " << p << G4endl;
254	G4cout << " v = " << v << G4endl;
255	G4cerr << "WARNING - Invalid call in "
256	<< "G4IntersectionSolid::DistanceToIn(p,v)" << G4endl
257	<< " Point p is inside !" << G4endl;
258	G4cerr << " p = " << p << G4endl;
259	G4cerr << " v = " << v << G4endl;
260	#endif
261	}
262	else // if( Inside(p) == kSurface )
263	{
264	EInside wA = fPtrSolidA->Inside(p);
265	EInside wB = fPtrSolidB->Inside(p);
266
267	G4ThreeVector pA = p, pB = p;
268	G4double dA = 0., dA1=0., dA2=0.;
269	G4double dB = 0., dB1=0., dB2=0.;
270	G4bool doA = true, doB = true;
271
272	while(true)
273	{
274	if(doA)
275	{
276	// find next valid range for A
277
278	dA1 = 0.;
279
280	if( wA != kInside )
281	{
282	dA1 = fPtrSolidA->DistanceToIn(pA, v);
283
284	if( dA1 == kInfinity ) return kInfinity;
285
286	pA += dA1*v;
287	}
288	dA2 = dA1 + fPtrSolidA->DistanceToOut(pA, v);
289	}
290	dA1 += dA;
291	dA2 += dA;
292
293	if(doB)
294	{
295	// find next valid range for B
296
297	dB1 = 0.;
298	if(wB != kInside)
299	{
300	dB1 = fPtrSolidB->DistanceToIn(pB, v);
301
302	if(dB1 == kInfinity) return kInfinity;
303
304	pB += dB1*v;
305	}
306	dB2 = dB1 + fPtrSolidB->DistanceToOut(pB, v);
307	}
308	dB1 += dB;
309	dB2 += dB;
310
311	// check if they overlap
312
313	if( dA1 < dB1 )
314	{
315	if( dB1 < dA2 ) return dB1;
316
317	dA = dA2;
318	pA = p + dA*v; // continue from here
319	wA = kSurface;
320	doA = true;
321	doB = false;
322	}
323	else
324	{
325	if( dA1 < dB2 ) return dA1;
326
327	dB = dB2;
328	pB = p + dB*v; // continue from here
329	wB = kSurface;
330	doB = true;
331	doA = false;
332	}
333	}
334	}
335	return dist ;
336	}
337
338	////////////////////////////////////////////////////////
339	//
340	// Approximate nearest distance from the point p to the intersection of
341	// two solids
342
343	G4double
344	G4IntersectionSolid::DistanceToIn( const G4ThreeVector& p) const
345	{
346	#ifdef G4BOOLDEBUG
347	if( Inside(p) == kInside )
348	{
349	G4cout << "WARNING - Invalid call in "
350	<< "G4IntersectionSolid::DistanceToIn(p)" << G4endl
351	<< " Point p is inside !" << G4endl;
352	G4cout << " p = " << p << G4endl;
353	G4cerr << "WARNING - Invalid call in "
354	<< "G4IntersectionSolid::DistanceToIn(p)" << G4endl
355	<< " Point p is inside !" << G4endl;
356	G4cerr << " p = " << p << G4endl;
357	}
358	#endif
359	EInside sideA = fPtrSolidA->Inside(p) ;
360	EInside sideB = fPtrSolidB->Inside(p) ;
361	G4double dist=0.0 ;
362
363	if( sideA != kInside && sideB != kOutside )
364	{
365	dist = fPtrSolidA->DistanceToIn(p) ;
366	}
367	else
368	{
369	if( sideB != kInside && sideA != kOutside )
370	{
371	dist = fPtrSolidB->DistanceToIn(p) ;
372	}
373	else
374	{
375	dist = std::min(fPtrSolidA->DistanceToIn(p),
376	fPtrSolidB->DistanceToIn(p) ) ;
377	}
378	}
379	return dist ;
380	}
381
382	//////////////////////////////////////////////////////////
383	//
384	// The same algorithm as DistanceToOut(p)
385
386	G4double
387	G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p,
388	const G4ThreeVector& v,
389	const G4bool calcNorm,
390	G4bool *validNorm,
391	G4ThreeVector *n ) const
392	{
393	G4bool validNormA, validNormB;
394	G4ThreeVector nA, nB;
395
396	#ifdef G4BOOLDEBUG
397	if( Inside(p) == kOutside )
398	{
399	G4cout << "Position:" << G4endl << G4endl;
400	G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
401	G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
402	G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
403	G4cout << "Direction:" << G4endl << G4endl;
404	G4cout << "v.x() = " << v.x() << G4endl;
405	G4cout << "v.y() = " << v.y() << G4endl;
406	G4cout << "v.z() = " << v.z() << G4endl << G4endl;
407	G4cout << "WARNING - Invalid call in "
408	<< "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
409	<< " Point p is outside !" << G4endl;
410	G4cout << " p = " << p << G4endl;
411	G4cout << " v = " << v << G4endl;
412	G4cerr << "WARNING - Invalid call in "
413	<< "G4IntersectionSolid::DistanceToOut(p,v)" << G4endl
414	<< " Point p is outside !" << G4endl;
415	G4cerr << " p = " << p << G4endl;
416	G4cerr << " v = " << v << G4endl;
417	}
418	#endif
419	G4double distA = fPtrSolidA->DistanceToOut(p,v,calcNorm,&validNormA,&nA) ;
420	G4double distB = fPtrSolidB->DistanceToOut(p,v,calcNorm,&validNormB,&nB) ;
421
422	G4double dist = std::min(distA,distB) ;
423
424	if( calcNorm )
425	{
426	if ( distA < distB )
427	{
428	*validNorm = validNormA;
429	*n = nA;
430	}
431	else
432	{
433	*validNorm = validNormB;
434	*n = nB;
435	}
436	}
437
438	return dist ;
439	}
440
441	//////////////////////////////////////////////////////////////
442	//
443	// Inverted algorithm of DistanceToIn(p)
444
445	G4double
446	G4IntersectionSolid::DistanceToOut( const G4ThreeVector& p ) const
447	{
448	#ifdef G4BOOLDEBUG
449	if( Inside(p) == kOutside )
450	{
451	G4cout << "WARNING - Invalid call in "
452	<< "G4IntersectionSolid::DistanceToOut(p)" << G4endl
453	<< " Point p is outside !" << G4endl;
454	G4cout << " p = " << p << G4endl;
455	G4cerr << "WARNING - Invalid call in "
456	<< "G4IntersectionSolid::DistanceToOut(p)" << G4endl
457	<< " Point p is outside !" << G4endl;
458	G4cerr << " p = " << p << G4endl;
459	}
460	#endif
461
462	return std::min(fPtrSolidA->DistanceToOut(p),
463	fPtrSolidB->DistanceToOut(p) ) ;
464
465	}
466
467	//////////////////////////////////////////////////////////////
468	//
469	//
470
471	void
472	G4IntersectionSolid::ComputeDimensions( G4VPVParameterisation*,
473	const G4int,
474	const G4VPhysicalVolume* )
475	{
476	}
477
478	/////////////////////////////////////////////////
479	//
480	//
481
482	G4GeometryType G4IntersectionSolid::GetEntityType() const
483	{
484	return G4String("G4IntersectionSolid");
485	}
486
487	/////////////////////////////////////////////////
488	//
489	//
490
491	void
492	G4IntersectionSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
493	{
494	scene.AddSolid (*this);
495	}
496
497	////////////////////////////////////////////////////
498	//
499	//
500
501	G4Polyhedron*
502	G4IntersectionSolid::CreatePolyhedron () const
503	{
504	G4Polyhedron* pA = fPtrSolidA->GetPolyhedron();
505	G4Polyhedron* pB = fPtrSolidB->GetPolyhedron();
506	if (pA && pB)
507	{
508	G4Polyhedron* resultant = new G4Polyhedron (pA->intersect(*pB));
509	return resultant;
510	}
511	else
512	{
513	std::ostringstream oss;
514	oss << "Solid - " << GetName()
515	<< " - one of the Boolean components has no" << G4endl
516	<< " corresponding polyhedron. Returning NULL !";
517	G4Exception("G4IntersectionSolid::CreatePolyhedron()", "InvalidSetup",
518	JustWarning, oss.str().c_str());
519	return 0;
520	}
521	}
522
523	/////////////////////////////////////////////////////////
524	//
525	//
526
527	G4NURBS*
528	G4IntersectionSolid::CreateNURBS () const
529	{
530	// Take into account boolean operation - see CreatePolyhedron.
531	// return new G4NURBSbox (1.0, 1.0, 1.0);
532	return 0;
533	}

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

Download in other formats: