source: trunk/source/geometry/management/src/G4ReflectedSolid.cc@ 1347

Last change on this file since 1347 was 1340, checked in by garnier, 15 years ago

update ti head
File size: 19.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: G4ReflectedSolid.cc,v 1.13 2010/10/19 15:20:18 gcosmo Exp $
28//
29// GEANT4 tag $Name: geommng-V09-03-05 $
30//
31// Implementation for G4ReflectedSolid class for boolean
32// operations between other solids
33//
34// Author: Vladimir Grichine, 23.07.01 (Vladimir.Grichine@cern.ch)
35//
36// --------------------------------------------------------------------
37
38#include "G4ReflectedSolid.hh"
39
40#include <sstream>
41
42#include "G4Point3D.hh"
43#include "G4Normal3D.hh"
44
45#include "G4VoxelLimits.hh"
46
47#include "G4VPVParameterisation.hh"
48
49#include "G4VGraphicsScene.hh"
50#include "G4Polyhedron.hh"
51#include "G4NURBS.hh"
52// #include "G4NURBSbox.hh"
53
54
55/////////////////////////////////////////////////////////////////
56//
57// Constructor using HepTransform3D, in fact HepReflect3D
58
59G4ReflectedSolid::G4ReflectedSolid( const G4String& pName,
60 G4VSolid* pSolid ,
61 const G4Transform3D& transform )
62 : G4VSolid(pName), fpPolyhedron(0)
63{
64 fPtrSolid = pSolid ;
65 G4RotationMatrix rotMatrix ;
66
67 fDirectTransform =
68 new G4AffineTransform(rotMatrix, transform.getTranslation()) ;
69 fPtrTransform =
70 new G4AffineTransform(rotMatrix, transform.getTranslation()) ;
71 fPtrTransform->Invert() ;
72
73 fDirectTransform3D = new G4Transform3D(transform) ;
74 fPtrTransform3D = new G4Transform3D(transform.inverse()) ;
75}
76
77///////////////////////////////////////////////////////////////////
78//
79
80G4ReflectedSolid::~G4ReflectedSolid()
81{
82 if(fPtrTransform)
83 {
84 delete fPtrTransform; fPtrTransform=0;
85 delete fDirectTransform; fDirectTransform=0;
86 }
87 if(fPtrTransform3D)
88 {
89 delete fPtrTransform3D; fPtrTransform3D=0;
90 delete fDirectTransform3D; fDirectTransform3D=0;
91 }
92 delete fpPolyhedron;
93}
94
95///////////////////////////////////////////////////////////////////
96//
97
98G4ReflectedSolid::G4ReflectedSolid(const G4ReflectedSolid& rhs)
99 : G4VSolid(rhs), fPtrSolid(rhs.fPtrSolid), fpPolyhedron(0)
100{
101 fPtrTransform = new G4AffineTransform(*rhs.fPtrTransform);
102 fDirectTransform = new G4AffineTransform(*rhs.fDirectTransform);
103 fPtrTransform3D = new G4Transform3D(*rhs.fPtrTransform3D);
104 fDirectTransform3D = new G4Transform3D(*rhs.fDirectTransform3D);
105}
106
107///////////////////////////////////////////////////////////////////
108//
109
110G4ReflectedSolid& G4ReflectedSolid::operator=(const G4ReflectedSolid& rhs)
111{
112 // Check assignment to self
113 //
114 if (this == &rhs) { return *this; }
115
116 // Copy base class data
117 //
118 G4VSolid::operator=(rhs);
119
120 // Copy data
121 //
122 fPtrSolid= rhs.fPtrSolid; fpPolyhedron= 0;
123 delete fPtrTransform;
124 fPtrTransform= new G4AffineTransform(*rhs.fPtrTransform);
125 delete fDirectTransform;
126 fDirectTransform= new G4AffineTransform(*rhs.fDirectTransform);
127 delete fPtrTransform3D;
128 fPtrTransform3D= new G4Transform3D(*rhs.fPtrTransform3D);
129 delete fDirectTransform3D;
130 fDirectTransform3D= new G4Transform3D(*rhs.fDirectTransform3D);
131
132 return *this;
133}
134
135///////////////////////////////////////////////////////////////////
136//
137
138G4GeometryType G4ReflectedSolid::GetEntityType() const
139{
140 return G4String("G4ReflectedSolid");
141}
142
143const G4ReflectedSolid* G4ReflectedSolid::GetReflectedSolidPtr() const
144{
145 return this;
146}
147
148G4ReflectedSolid* G4ReflectedSolid::GetReflectedSolidPtr()
149{
150 return this;
151}
152
153G4VSolid* G4ReflectedSolid::GetConstituentMovedSolid() const
154{
155 return fPtrSolid;
156}
157
158/////////////////////////////////////////////////////////////////////////////
159
160G4AffineTransform G4ReflectedSolid::GetTransform() const
161{
162 G4AffineTransform aTransform = *fPtrTransform;
163 return aTransform;
164}
165
166void G4ReflectedSolid::SetTransform(G4AffineTransform& transform)
167{
168 fPtrTransform = &transform ;
169 fpPolyhedron = 0;
170}
171
172//////////////////////////////////////////////////////////////////////////////
173
174G4AffineTransform G4ReflectedSolid::GetDirectTransform() const
175{
176 G4AffineTransform aTransform= *fDirectTransform;
177 return aTransform;
178}
179
180void G4ReflectedSolid::SetDirectTransform(G4AffineTransform& transform)
181{
182 fDirectTransform = &transform ;
183 fpPolyhedron = 0;
184}
185
186/////////////////////////////////////////////////////////////////////////////
187
188G4Transform3D G4ReflectedSolid::GetTransform3D() const
189{
190 G4Transform3D aTransform = *fPtrTransform3D;
191 return aTransform;
192}
193
194void G4ReflectedSolid::SetTransform3D(G4Transform3D& transform)
195{
196 fPtrTransform3D = &transform ;
197 fpPolyhedron = 0;
198}
199
200//////////////////////////////////////////////////////////////////////////////
201
202G4Transform3D G4ReflectedSolid::GetDirectTransform3D() const
203{
204 G4Transform3D aTransform= *fDirectTransform3D;
205 return aTransform;
206}
207
208void G4ReflectedSolid::SetDirectTransform3D(G4Transform3D& transform)
209{
210 fDirectTransform3D = &transform ;
211 fpPolyhedron = 0;
212}
213
214/////////////////////////////////////////////////////////////////////////////
215
216G4RotationMatrix G4ReflectedSolid::GetFrameRotation() const
217{
218 G4RotationMatrix InvRotation= fDirectTransform->NetRotation();
219 return InvRotation;
220}
221
222void G4ReflectedSolid::SetFrameRotation(const G4RotationMatrix& matrix)
223{
224 fDirectTransform->SetNetRotation(matrix);
225}
226
227/////////////////////////////////////////////////////////////////////////////
228
229G4ThreeVector G4ReflectedSolid::GetFrameTranslation() const
230{
231 return fPtrTransform->NetTranslation();
232}
233
234void G4ReflectedSolid::SetFrameTranslation(const G4ThreeVector& vector)
235{
236 fPtrTransform->SetNetTranslation(vector);
237}
238
239///////////////////////////////////////////////////////////////
240
241G4RotationMatrix G4ReflectedSolid::GetObjectRotation() const
242{
243 G4RotationMatrix Rotation= fPtrTransform->NetRotation();
244 return Rotation;
245}
246
247void G4ReflectedSolid::SetObjectRotation(const G4RotationMatrix& matrix)
248{
249 fPtrTransform->SetNetRotation(matrix);
250}
251
252///////////////////////////////////////////////////////////////////////
253
254G4ThreeVector G4ReflectedSolid::GetObjectTranslation() const
255{
256 return fDirectTransform->NetTranslation();
257}
258
259void G4ReflectedSolid::SetObjectTranslation(const G4ThreeVector& vector)
260{
261 fDirectTransform->SetNetTranslation(vector);
262}
263
264///////////////////////////////////////////////////////////////
265//
266//
267
268G4bool
269G4ReflectedSolid::CalculateExtent( const EAxis pAxis,
270 const G4VoxelLimits& pVoxelLimit,
271 const G4AffineTransform& pTransform,
272 G4double& pMin,
273 G4double& pMax ) const
274{
275
276 G4VoxelLimits unLimit;
277 G4AffineTransform unTransform;
278
279 G4double x1 = -kInfinity, x2 = kInfinity,
280 y1 = -kInfinity, y2 = kInfinity,
281 z1 = -kInfinity, z2 = kInfinity;
282
283 G4bool existsAfterClip = false ;
284 existsAfterClip =
285 fPtrSolid->CalculateExtent(kXAxis,unLimit,unTransform,x1,x2);
286 existsAfterClip =
287 fPtrSolid->CalculateExtent(kYAxis,unLimit,unTransform,y1,y2);
288 existsAfterClip =
289 fPtrSolid->CalculateExtent(kZAxis,unLimit,unTransform,z1,z2);
290
291 existsAfterClip = false;
292 pMin = +kInfinity ;
293 pMax = -kInfinity ;
294
295 G4Transform3D pTransform3D = G4Transform3D(pTransform.NetRotation().inverse(),
296 pTransform.NetTranslation());
297
298 G4Transform3D transform3D = pTransform3D*(*fDirectTransform3D);
299
300 G4Point3D tmpPoint;
301
302 // Calculate rotated vertex coordinates
303
304 G4ThreeVectorList* vertices = new G4ThreeVectorList();
305
306 if (vertices)
307 {
308 vertices->reserve(8);
309
310 G4ThreeVector vertex0(x1,y1,z1) ;
311 tmpPoint = transform3D*G4Point3D(vertex0);
312 vertex0 = G4ThreeVector(tmpPoint.x(),tmpPoint.y(),tmpPoint.z());
313 vertices->push_back(vertex0);
314
315 G4ThreeVector vertex1(x2,y1,z1) ;
316 tmpPoint = transform3D*G4Point3D(vertex1);
317 vertex1 = G4ThreeVector(tmpPoint.x(),tmpPoint.y(),tmpPoint.z());
318 vertices->push_back(vertex1);
319
320 G4ThreeVector vertex2(x2,y2,z1) ;
321 tmpPoint = transform3D*G4Point3D(vertex2);
322 vertex2 = G4ThreeVector(tmpPoint.x(),tmpPoint.y(),tmpPoint.z());
323 vertices->push_back(vertex2);
324
325 G4ThreeVector vertex3(x1,y2,z1) ;
326 tmpPoint = transform3D*G4Point3D(vertex3);
327 vertex3 = G4ThreeVector(tmpPoint.x(),tmpPoint.y(),tmpPoint.z());
328 vertices->push_back(vertex3);
329
330 G4ThreeVector vertex4(x1,y1,z2) ;
331 tmpPoint = transform3D*G4Point3D(vertex4);
332 vertex4 = G4ThreeVector(tmpPoint.x(),tmpPoint.y(),tmpPoint.z());
333 vertices->push_back(vertex4);
334
335 G4ThreeVector vertex5(x2,y1,z2) ;
336 tmpPoint = transform3D*G4Point3D(vertex5);
337 vertex5 = G4ThreeVector(tmpPoint.x(),tmpPoint.y(),tmpPoint.z());
338 vertices->push_back(vertex5);
339
340 G4ThreeVector vertex6(x2,y2,z2) ;
341 tmpPoint = transform3D*G4Point3D(vertex6);
342 vertex6 = G4ThreeVector(tmpPoint.x(),tmpPoint.y(),tmpPoint.z());
343 vertices->push_back(vertex6);
344
345 G4ThreeVector vertex7(x1,y2,z2) ;
346 tmpPoint = transform3D*G4Point3D(vertex7);
347 vertex7 = G4ThreeVector(tmpPoint.x(),tmpPoint.y(),tmpPoint.z());
348 vertices->push_back(vertex7);
349 }
350 else
351 {
352 DumpInfo();
353 G4Exception("G4ReflectedSolid::CalculateExtent()",
354 "FatalError", FatalException,
355 "Error in allocation of vertices. Out of memory !");
356 }
357
358 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
359 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax) ;
360 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
361
362 if (pVoxelLimit.IsLimited(pAxis) == false)
363 {
364 if ( pMin != kInfinity || pMax != -kInfinity )
365 {
366 existsAfterClip = true ;
367
368 // Add 2*tolerance to avoid precision troubles
369
370 pMin -= kCarTolerance;
371 pMax += kCarTolerance;
372 }
373 }
374 else
375 {
376 G4ThreeVector clipCentre(
377 ( pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
378 ( pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
379 ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
380
381 if ( pMin != kInfinity || pMax != -kInfinity )
382 {
383 existsAfterClip = true ;
384
385
386 // Check to see if endpoints are in the solid
387
388 clipCentre(pAxis) = pVoxelLimit.GetMinExtent(pAxis);
389
390 if (Inside(transform3D.inverse()*G4Point3D(clipCentre)) != kOutside)
391 {
392 pMin = pVoxelLimit.GetMinExtent(pAxis);
393 }
394 else
395 {
396 pMin -= kCarTolerance;
397 }
398 clipCentre(pAxis) = pVoxelLimit.GetMaxExtent(pAxis);
399
400 if (Inside(transform3D.inverse()*G4Point3D(clipCentre)) != kOutside)
401 {
402 pMax = pVoxelLimit.GetMaxExtent(pAxis);
403 }
404 else
405 {
406 pMax += kCarTolerance;
407 }
408 }
409 // Check for case where completely enveloping clipping volume
410 // If point inside then we are confident that the solid completely
411 // envelopes the clipping volume. Hence set min/max extents according
412 // to clipping volume extents along the specified axis.
413
414 else if (Inside(transform3D.inverse()*G4Point3D(clipCentre)) != kOutside)
415 {
416 existsAfterClip = true ;
417 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
418 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
419 }
420 }
421 delete vertices;
422 return existsAfterClip;
423}
424
425/////////////////////////////////////////////////////
426//
427//
428
429EInside G4ReflectedSolid::Inside(const G4ThreeVector& p) const
430{
431
432 G4Point3D newPoint = (*fDirectTransform3D)*G4Point3D(p) ;
433 // G4Point3D newPoint = (*fPtrTransform3D)*G4Point3D(p) ;
434
435 return fPtrSolid->Inside(G4ThreeVector(newPoint.x(),
436 newPoint.y(),
437 newPoint.z())) ;
438}
439
440//////////////////////////////////////////////////////////////
441//
442//
443
444G4ThreeVector
445G4ReflectedSolid::SurfaceNormal( const G4ThreeVector& p ) const
446{
447 G4Point3D newPoint = (*fDirectTransform3D)*G4Point3D(p) ;
448 G4ThreeVector normal =
449 fPtrSolid->SurfaceNormal(G4ThreeVector(newPoint.x(),
450 newPoint.y(),
451 newPoint.z() ) ) ;
452 G4Point3D newN = (*fDirectTransform3D)*G4Point3D(normal) ;
453 newN.unit() ;
454
455 return G4ThreeVector(newN.x(),newN.y(),newN.z()) ;
456}
457
458/////////////////////////////////////////////////////////////
459//
460// The same algorithm as in DistanceToIn(p)
461
462G4double
463G4ReflectedSolid::DistanceToIn( const G4ThreeVector& p,
464 const G4ThreeVector& v ) const
465{
466 G4Point3D newPoint = (*fDirectTransform3D)*G4Point3D(p) ;
467 G4Point3D newDirection = (*fDirectTransform3D)*G4Point3D(v) ;
468 newDirection.unit() ;
469 return fPtrSolid->DistanceToIn(
470 G4ThreeVector(newPoint.x(),newPoint.y(),newPoint.z()),
471 G4ThreeVector(newDirection.x(),newDirection.y(),newDirection.z())) ;
472}
473
474////////////////////////////////////////////////////////
475//
476// Approximate nearest distance from the point p to the intersection of
477// two solids
478
479G4double
480G4ReflectedSolid::DistanceToIn( const G4ThreeVector& p) const
481{
482 G4Point3D newPoint = (*fDirectTransform3D)*G4Point3D(p) ;
483 return fPtrSolid->DistanceToIn(
484 G4ThreeVector(newPoint.x(),newPoint.y(),newPoint.z())) ;
485}
486
487//////////////////////////////////////////////////////////
488//
489// The same algorithm as DistanceToOut(p)
490
491G4double
492G4ReflectedSolid::DistanceToOut( const G4ThreeVector& p,
493 const G4ThreeVector& v,
494 const G4bool calcNorm,
495 G4bool *validNorm,
496 G4ThreeVector *n ) const
497{
498 G4ThreeVector solNorm ;
499
500 G4Point3D newPoint = (*fDirectTransform3D)*G4Point3D(p) ;
501 G4Point3D newDirection = (*fDirectTransform3D)*G4Point3D(v);
502 newDirection.unit() ;
503
504 G4double dist =
505 fPtrSolid->DistanceToOut(
506 G4ThreeVector(newPoint.x(),newPoint.y(),newPoint.z()),
507 G4ThreeVector(newDirection.x(),newDirection.y(),newDirection.z()),
508 calcNorm, validNorm, &solNorm) ;
509 if(calcNorm)
510 {
511 G4Point3D newN = (*fDirectTransform3D)*G4Point3D(solNorm);
512 newN.unit() ;
513 *n = G4ThreeVector(newN.x(),newN.y(),newN.z());
514 }
515 return dist ;
516}
517
518//////////////////////////////////////////////////////////////
519//
520// Inverted algorithm of DistanceToIn(p)
521
522G4double
523G4ReflectedSolid::DistanceToOut( const G4ThreeVector& p ) const
524{
525 G4Point3D newPoint = (*fDirectTransform3D)*G4Point3D(p);
526 return fPtrSolid->DistanceToOut(
527 G4ThreeVector(newPoint.x(),newPoint.y(),newPoint.z()));
528}
529
530//////////////////////////////////////////////////////////////
531//
532//
533
534void
535G4ReflectedSolid::ComputeDimensions( G4VPVParameterisation*,
536 const G4int,
537 const G4VPhysicalVolume* )
538{
539 DumpInfo();
540 G4Exception("G4ReflectedSolid::ComputeDimensions()",
541 "NotApplicable", FatalException,
542 "Method not applicable in this context!");
543}
544
545//////////////////////////////////////////////////////////////
546//
547// Return a point (G4ThreeVector) randomly and uniformly selected
548// on the solid surface
549
550G4ThreeVector G4ReflectedSolid::GetPointOnSurface() const
551{
552 G4ThreeVector p = fPtrSolid->GetPointOnSurface();
553 G4Point3D newPoint = (*fDirectTransform3D)*G4Point3D(p);
554
555 return G4ThreeVector(newPoint.x(),newPoint.y(),newPoint.z());
556}
557
558//////////////////////////////////////////////////////////////////////////
559//
560// Make a clone of this object
561
562G4VSolid* G4ReflectedSolid::Clone() const
563{
564 return new G4ReflectedSolid(*this);
565}
566
567
568//////////////////////////////////////////////////////////////////////////
569//
570// Stream object contents to an output stream
571
572std::ostream& G4ReflectedSolid::StreamInfo(std::ostream& os) const
573{
574 os << "-----------------------------------------------------------\n"
575 << " *** Dump for Reflected solid - " << GetName() << " ***\n"
576 << " ===================================================\n"
577 << " Solid type: " << GetEntityType() << "\n"
578 << " Parameters of constituent solid: \n"
579 << "===========================================================\n";
580 fPtrSolid->StreamInfo(os);
581 os << "===========================================================\n"
582 << " Transformations: \n"
583 << " Direct transformation - translation : \n"
584 << " " << fDirectTransform->NetTranslation() << "\n"
585 << " - rotation : \n"
586 << " ";
587 fDirectTransform->NetRotation().print(os);
588 os << "\n"
589 << "===========================================================\n";
590
591 return os;
592}
593
594/////////////////////////////////////////////////
595//
596//
597
598void
599G4ReflectedSolid::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
600{
601 scene.AddSolid (*this);
602}
603
604////////////////////////////////////////////////////
605//
606//
607
608G4Polyhedron*
609G4ReflectedSolid::CreatePolyhedron () const
610{
611 G4Polyhedron* polyhedron = fPtrSolid->CreatePolyhedron();
612 if (polyhedron)
613 {
614 polyhedron->Transform(*fDirectTransform3D);
615 return polyhedron;
616 }
617 else
618 {
619 std::ostringstream oss;
620 oss << "Solid - " << GetName()
621 << " - original solid has no" << G4endl
622 << " corresponding polyhedron. Returning NULL!";
623 G4Exception("G4ReflectedSolid::CreatePolyhedron()", "InvalidSetup",
624 JustWarning, oss.str().c_str());
625 return 0;
626 }
627}
628
629/////////////////////////////////////////////////////////
630//
631//
632
633G4NURBS*
634G4ReflectedSolid::CreateNURBS () const
635{
636 // Take into account local transformation - see CreatePolyhedron.
637 // return fPtrSolid->CreateNURBS() ;
638 return 0;
639}
640
641/////////////////////////////////////////////////////////
642//
643//
644
645G4Polyhedron*
646G4ReflectedSolid::GetPolyhedron () const
647{
648 if (!fpPolyhedron ||
649 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
650 fpPolyhedron->GetNumberOfRotationSteps())
651 {
652 delete fpPolyhedron;
653 fpPolyhedron = CreatePolyhedron ();
654 }
655 return fpPolyhedron;
656}
Note: See TracBrowser for help on using the repository browser.