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