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

Last change on this file since 1202 was 1058, checked in by garnier, 17 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.