Context Navigation

G4IntersectionSolid.cc @ 1337

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

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

Context Navigation

source: trunk/source/geometry/solids/Boolean/src/G4IntersectionSolid.cc @ 1337

Download in other formats: