source: trunk/source/geometry/solids/specific/src/G4TwistTrapFlatSide.cc @ 1202

Last change on this file since 1202 was 1058, checked in by garnier, 15 years ago

file release beta
File size: 17.4 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: G4TwistTrapFlatSide.cc,v 1.6 2007/05/23 09:31:02 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30//
31// --------------------------------------------------------------------
32// GEANT 4 class source file
33//
34//
35// G4TwistTrapFlatSide.cc
36//
37// Author:
38//   30-Aug-2002 - O.Link (Oliver.Link@cern.ch)
39//
40// --------------------------------------------------------------------
41
42#include "G4TwistTrapFlatSide.hh"
43
44//=====================================================================
45//* constructors ------------------------------------------------------
46
47G4TwistTrapFlatSide::G4TwistTrapFlatSide( const G4String        &name,
48                              G4double      PhiTwist,
49                              G4double      pDx1,
50                              G4double      pDx2,
51                              G4double      pDy,
52                              G4double      pDz,
53                              G4double      pAlpha,
54                              G4double      pPhi,
55                              G4double      pTheta,
56                              G4int         handedness) 
57
58  : G4VTwistSurface(name)
59{
60   fHandedness = handedness;   // +z = +ve, -z = -ve
61
62   fDx1 = pDx1 ;
63   fDx2 = pDx2 ;
64   fDy = pDy ;
65   fDz = pDz ;
66   fAlpha = pAlpha ;
67   fTAlph = std::tan(fAlpha) ;
68   fPhi  = pPhi ;
69   fTheta = pTheta ;
70
71   fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi)  ;
72     // dx in surface equation
73   fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi)  ;
74     // dy in surface equation
75
76   fPhiTwist = PhiTwist ;
77
78   fCurrentNormal.normal.set( 0, 0, (fHandedness < 0 ? -1 : 1)); 
79         // Unit vector, in local coordinate system
80   fRot.rotateZ( fHandedness > 0 
81                 ? 0.5 * fPhiTwist
82                 : -0.5 * fPhiTwist );
83
84   fTrans.set(
85              fHandedness > 0 ? 0.5*fdeltaX : -0.5*fdeltaX , 
86              fHandedness > 0 ? 0.5*fdeltaY : -0.5*fdeltaY ,
87              fHandedness > 0 ? fDz : -fDz ) ;
88
89   fIsValidNorm = true;
90
91
92   fAxis[0] = kXAxis ;
93   fAxis[1] = kYAxis ;
94   fAxisMin[0] = kInfinity ;  // x-Axis cannot be fixed, because it
95   fAxisMax[0] =  kInfinity ; // depends on y
96   fAxisMin[1] = -fDy ;  // y - axis
97   fAxisMax[1] =  fDy ;
98
99   SetCorners();
100   SetBoundaries();
101}
102
103
104//=====================================================================
105//* Fake default constructor ------------------------------------------
106
107G4TwistTrapFlatSide::G4TwistTrapFlatSide( __void__& a )
108  : G4VTwistSurface(a)
109{
110}
111
112
113//=====================================================================
114//* destructor --------------------------------------------------------
115
116G4TwistTrapFlatSide::~G4TwistTrapFlatSide()
117{
118}
119
120//=====================================================================
121//* GetNormal ---------------------------------------------------------
122
123G4ThreeVector G4TwistTrapFlatSide::GetNormal(const G4ThreeVector & /* xx */ , 
124                                             G4bool isGlobal)
125{
126   if (isGlobal) {
127      return ComputeGlobalDirection(fCurrentNormal.normal);
128   } else {
129      return fCurrentNormal.normal;
130   }
131}
132
133//=====================================================================
134//* DistanceToSurface(p, v) -------------------------------------------
135
136G4int G4TwistTrapFlatSide::DistanceToSurface(const G4ThreeVector &gp,
137                                       const G4ThreeVector &gv,
138                                             G4ThreeVector  gxx[],
139                                             G4double       distance[],
140                                             G4int          areacode[],
141                                             G4bool         isvalid[],
142                                             EValidate      validate) 
143{
144   fCurStatWithV.ResetfDone(validate, &gp, &gv);
145   
146   if (fCurStatWithV.IsDone()) {
147      G4int i;
148      for (i=0; i<fCurStatWithV.GetNXX(); i++) {
149         gxx[i] = fCurStatWithV.GetXX(i);
150         distance[i] = fCurStatWithV.GetDistance(i);
151         areacode[i] = fCurStatWithV.GetAreacode(i);
152         isvalid[i]  = fCurStatWithV.IsValid(i);
153      }
154      return fCurStatWithV.GetNXX();
155   } else {
156      // initialize
157      G4int i;
158      for (i=0; i<2; i++) {
159         distance[i] = kInfinity;
160         areacode[i] = sOutside;
161         isvalid[i]  = false;
162         gxx[i].set(kInfinity, kInfinity, kInfinity);
163      }
164   }
165
166   G4ThreeVector p = ComputeLocalPoint(gp);
167   G4ThreeVector v = ComputeLocalDirection(gv);
168
169   //
170   // special case!
171   // if p is on surface, distance = 0.
172   //
173
174   if (std::fabs(p.z()) == 0.) {   // if p is on the plane
175      distance[0] = 0;
176      G4ThreeVector xx = p;
177      gxx[0] = ComputeGlobalPoint(xx);
178     
179      if (validate == kValidateWithTol) {
180         areacode[0] = GetAreaCode(xx);
181         if (!IsOutside(areacode[0])) {
182            isvalid[0] = true;
183         }
184      } else if (validate == kValidateWithoutTol) {
185         areacode[0] = GetAreaCode(xx, false);
186         if (IsInside(areacode[0])) {
187            isvalid[0] = true;
188         }
189      } else { // kDontValidate
190         areacode[0] = sInside;
191         isvalid[0] = true;
192      }
193
194      return 1;
195   }
196   //
197   // special case end
198   //
199   
200   if (v.z() == 0) { 
201
202      fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], 
203                                     isvalid[0], 0, validate, &gp, &gv);
204      return 0;
205   }
206   
207   distance[0] = - (p.z() / v.z());
208     
209   G4ThreeVector xx = p + distance[0]*v;
210   gxx[0] = ComputeGlobalPoint(xx);
211
212   if (validate == kValidateWithTol) {
213      areacode[0] = GetAreaCode(xx);
214      if (!IsOutside(areacode[0])) {
215         if (distance[0] >= 0) isvalid[0] = true;
216      }
217   } else if (validate == kValidateWithoutTol) {
218      areacode[0] = GetAreaCode(xx, false);
219      if (IsInside(areacode[0])) {
220         if (distance[0] >= 0) isvalid[0] = true;
221      }
222   } else { // kDontValidate
223      areacode[0] = sInside;
224         if (distance[0] >= 0) isvalid[0] = true;
225   }
226
227
228   fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
229                                  isvalid[0], 1, validate, &gp, &gv);
230
231#ifdef G4TWISTDEBUG
232   G4cerr << "ERROR - G4TwistTrapFlatSide::DistanceToSurface(p,v)" << G4endl;
233   G4cerr << "        Name        : " << GetName() << G4endl;
234   G4cerr << "        xx          : " << xx << G4endl;
235   G4cerr << "        gxx[0]      : " << gxx[0] << G4endl;
236   G4cerr << "        dist[0]     : " << distance[0] << G4endl;
237   G4cerr << "        areacode[0] : " << areacode[0] << G4endl;
238   G4cerr << "        isvalid[0]  : " << isvalid[0]  << G4endl;
239#endif
240   return 1;
241}
242
243//=====================================================================
244//* DistanceToSurface(p) ----------------------------------------------
245
246G4int G4TwistTrapFlatSide::DistanceToSurface(const G4ThreeVector &gp,
247                                             G4ThreeVector  gxx[],
248                                             G4double       distance[],
249                                             G4int          areacode[])
250{
251   // Calculate distance to plane in local coordinate,
252   // then return distance and global intersection points.
253   // 
254
255   fCurStat.ResetfDone(kDontValidate, &gp);
256
257   if (fCurStat.IsDone()) {
258      G4int i;
259      for (i=0; i<fCurStat.GetNXX(); i++) {
260         gxx[i] = fCurStat.GetXX(i);
261         distance[i] = fCurStat.GetDistance(i);
262         areacode[i] = fCurStat.GetAreacode(i);
263      }
264      return fCurStat.GetNXX();
265   } else {
266      // initialize
267      G4int i;
268      for (i=0; i<2; i++) {
269         distance[i] = kInfinity;
270         areacode[i] = sOutside;
271         gxx[i].set(kInfinity, kInfinity, kInfinity);
272      }
273   }
274   
275   G4ThreeVector p = ComputeLocalPoint(gp);
276   G4ThreeVector xx;
277
278   // The plane is placed on origin with making its normal
279   // parallel to z-axis.
280   if (std::fabs(p.z()) <= 0.5 * kCarTolerance)
281   {   // if p is on the plane, return 1
282      distance[0] = 0;
283      xx = p;
284   } else {
285      distance[0] = std::fabs(p.z());
286      xx.set(p.x(), p.y(), 0); 
287   }
288
289   gxx[0] = ComputeGlobalPoint(xx);
290   areacode[0] = sInside;
291   G4bool isvalid = true;
292   fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
293                             isvalid, 1, kDontValidate, &gp);
294   return 1;
295
296}
297
298G4int G4TwistTrapFlatSide::GetAreaCode(const G4ThreeVector &xx, 
299                                       G4bool withTol)
300{
301
302  static const G4double ctol = 0.5 * kCarTolerance;
303  G4int areacode = sInside;
304 
305  if (fAxis[0] == kXAxis && fAxis[1] == kYAxis) {
306
307    G4int yaxis = 1;
308   
309    G4double wmax = xAxisMax(xx.y(), fTAlph ) ;
310    G4double wmin = -xAxisMax(xx.y(), -fTAlph ) ;
311
312    if (withTol) {
313     
314      G4bool isoutside   = false;
315     
316      // test boundary of x-axis
317     
318      if (xx.x() < wmin + ctol) {
319        areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary; 
320        if (xx.x() <= wmin - ctol) isoutside = true;
321       
322      } else if (xx.x() > wmax - ctol) {
323        areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
324        if (xx.x() >= wmax + ctol)  isoutside = true;
325      }
326     
327      // test boundary of y-axis
328     
329      if (xx.y() < fAxisMin[yaxis] + ctol) {
330        areacode |= (sAxis1 & (sAxisY | sAxisMin)); 
331       
332        if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
333        else                        areacode |= sBoundary;
334        if (xx.y() <= fAxisMin[yaxis] - ctol) isoutside = true;
335       
336      } else if (xx.y() > fAxisMax[yaxis] - ctol) {
337        areacode |= (sAxis1 & (sAxisY | sAxisMax));
338       
339        if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
340        else                        areacode |= sBoundary; 
341        if (xx.y() >= fAxisMax[yaxis] + ctol) isoutside = true;
342      }
343     
344      // if isoutside = true, clear inside bit.             
345      // if not on boundary, add axis information.             
346     
347      if (isoutside) {
348        G4int tmpareacode = areacode & (~sInside);
349        areacode = tmpareacode;
350      } else if ((areacode & sBoundary) != sBoundary) {
351        areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisY);
352      }           
353     
354    } else {
355     
356      // boundary of x-axis
357     
358      if (xx.x() < wmin ) {
359        areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
360      } else if (xx.x() > wmax) {
361        areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
362      }
363     
364      // boundary of y-axis
365     
366      if (xx.y() < fAxisMin[yaxis]) {
367        areacode |= (sAxis1 & (sAxisY | sAxisMin));
368        if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
369        else                        areacode |= sBoundary; 
370       
371      } else if (xx.y() > fAxisMax[yaxis]) {
372        areacode |= (sAxis1 & (sAxisY | sAxisMax)) ;
373        if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
374        else                        areacode |= sBoundary; 
375      }
376     
377      if ((areacode & sBoundary) != sBoundary) {
378        areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisY);
379      }           
380    }
381    return areacode;
382  } else {
383    G4Exception("G4TwistTrapFlatSide::GetAreaCode()",
384                "NotImplemented", FatalException,
385                "Feature NOT implemented !");
386  }
387 
388  return areacode;
389}
390
391
392//=====================================================================
393//* SetCorners --------------------------------------------------------
394
395void G4TwistTrapFlatSide::SetCorners()
396{
397   // Set Corner points in local coodinate.
398
399   if (fAxis[0] == kXAxis && fAxis[1] == kYAxis) {
400   
401     G4double x, y, z;
402     
403     // corner of Axis0min and Axis1min
404     x = -fDx1 + fDy * fTAlph ;
405     y = -fDy ;
406     z = 0 ;
407     SetCorner(sC0Min1Min, x, y, z);
408     
409     // corner of Axis0max and Axis1min
410     x = fDx1 + fDy * fTAlph ;
411     y = -fDy ;
412     z = 0 ;
413     SetCorner(sC0Max1Min, x, y, z);
414     
415     // corner of Axis0max and Axis1max
416     x = fDx2 - fDy * fTAlph ;
417     y = fDy ;
418     z = 0 ;
419     SetCorner(sC0Max1Max, x, y, z);
420     
421     // corner of Axis0min and Axis1max
422     x = -fDx2 - fDy * fTAlph ;
423     y = fDy ;
424     z = 0 ;
425     SetCorner(sC0Min1Max, x, y, z);
426     
427   } else {
428     G4cerr << "ERROR - G4TwistTrapFlatSide::SetCorners()" << G4endl
429            << "        fAxis[0] = " << fAxis[0] << G4endl
430            << "        fAxis[1] = " << fAxis[1] << G4endl;
431     G4Exception("G4TwistTrapFlatSide::SetCorners()",
432                 "NotImplemented", FatalException,
433                 "Feature NOT implemented !");
434   }
435}
436
437//=====================================================================
438//* SetBoundaries() ---------------------------------------------------
439
440void G4TwistTrapFlatSide::SetBoundaries()
441{
442   // Set direction-unit vector of phi-boundary-lines in local coodinate.
443   // Don't call the function twice.
444   
445  G4ThreeVector direction ;
446
447  if (fAxis[0] == kXAxis && fAxis[1] == kYAxis) {
448   
449    // sAxis0 & sAxisMin
450    direction = - ( GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min) ) ;
451    direction = direction.unit();
452    SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction, 
453                GetCorner(sC0Min1Max), sAxisY) ;
454   
455    // sAxis0 & sAxisMax
456    direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min)  ; // inverse
457    direction = direction.unit();
458    SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction, 
459                GetCorner(sC0Max1Min), sAxisY);
460   
461    // sAxis1 & sAxisMin
462    direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
463    direction = direction.unit();
464    SetBoundary(sAxis1 & (sAxisY | sAxisMin), direction, 
465                GetCorner(sC0Min1Min), sAxisX);
466   
467    // sAxis1 & sAxisMax
468    direction = - ( GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max) ) ;
469    direction = direction.unit();
470    SetBoundary(sAxis1 & (sAxisY | sAxisMax), direction, 
471                GetCorner(sC0Max1Max), sAxisX);
472   
473  } else {
474    G4cerr << "ERROR - G4TwistTrapFlatSide::SetBoundaries()" << G4endl
475           << "        fAxis[0] = " << fAxis[0] << G4endl
476           << "        fAxis[1] = " << fAxis[1] << G4endl;
477    G4Exception("G4TwistTrapFlatSide::SetCorners()",
478                "NotImplemented", FatalException,
479                "Feature NOT implemented !");
480   }
481}
482
483//=====================================================================
484//* GetFacets() -------------------------------------------------------
485
486void G4TwistTrapFlatSide::GetFacets( G4int m, G4int n, G4double xyz[][3],
487                                     G4int faces[][4], G4int iside ) 
488{
489
490  G4double x,y    ;     // the two parameters for the surface equation
491  G4ThreeVector p ;  // a point on the surface, given by (z,u)
492
493  G4int nnode ;
494  G4int nface ;
495
496  G4double xmin,xmax ;
497
498  // calculate the (n-1)*(m-1) vertices
499
500  G4int i,j ;
501
502  for ( i = 0 ; i<n ; i++ ) {
503
504    y = -fDy + i*(2*fDy)/(n-1) ;
505
506    for ( j = 0 ; j<m ; j++ ) {
507
508      xmin = GetBoundaryMin(y) ;
509      xmax = GetBoundaryMax(y) ;
510      x = xmin + j*(xmax-xmin)/(m-1) ;
511
512      nnode = GetNode(i,j,m,n,iside) ;
513      p = SurfacePoint(x,y,true) ;  // surface point in global coordinate system
514
515      xyz[nnode][0] = p.x() ;
516      xyz[nnode][1] = p.y() ;
517      xyz[nnode][2] = p.z() ;
518
519      if ( i<n-1 && j<m-1 ) {   
520
521        nface = GetFace(i,j,m,n,iside) ;
522
523        if (fHandedness < 0) {  // lower side
524          faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,1) * ( GetNode(i  ,j  ,m,n,iside)+1) ; 
525          faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,1) * ( GetNode(i+1,j  ,m,n,iside)+1) ;
526          faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
527          faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,1) * ( GetNode(i  ,j+1,m,n,iside)+1) ;
528        } else {                // upper side
529          faces[nface][0] = GetEdgeVisibility(i,j,m,n,0,-1) * ( GetNode(i  ,j  ,m,n,iside)+1) ; 
530          faces[nface][1] = GetEdgeVisibility(i,j,m,n,1,-1) * ( GetNode(i  ,j+1,m,n,iside)+1) ;
531          faces[nface][2] = GetEdgeVisibility(i,j,m,n,2,-1) * ( GetNode(i+1,j+1,m,n,iside)+1) ;
532          faces[nface][3] = GetEdgeVisibility(i,j,m,n,3,-1) * ( GetNode(i+1,j  ,m,n,iside)+1) ;
533        }
534
535      }
536    }
537  }
538}
Note: See TracBrowser for help on using the repository browser.