| Geant4 Cross Reference |
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 // G4TwistTubsFlatSide implementation
27 //
28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp), created.
29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
30 // from original version in Jupiter-2.5.02 application.
31 // --------------------------------------------------------------------
32
33 #include "G4TwistTubsFlatSide.hh"
34 #include "G4GeometryTolerance.hh"
35
36 //=====================================================================
37 //* constructors ------------------------------------------------------
38
39 G4TwistTubsFlatSide::G4TwistTubsFlatSide(const G4String& name,
40 const G4RotationMatrix& rot,
41 const G4ThreeVector& tlate,
42 const G4ThreeVector& n,
43 const EAxis axis0 ,
44 const EAxis axis1 ,
45 G4double axis0min,
46 G4double axis1min,
47 G4double axis0max,
48 G4double axis1max )
49 : G4VTwistSurface(name, rot, tlate, 0, axis0, axis1,
50 axis0min, axis1min, axis0max, axis1max)
51 {
52 if (axis0 == kPhi && axis1 == kRho)
53 {
54 G4Exception("G4TwistTubsFlatSide::G4TwistTubsFlatSide()",
55 "GeomSolids0002", FatalErrorInArgument,
56 "Should swap axis0 and axis1!");
57 }
58
59 G4ThreeVector normal = rot.inverse()*n;
60 fCurrentNormal.normal = normal.unit(); // in local coordinate system
61 fIsValidNorm = true;
62
63 SetCorners();
64 SetBoundaries();
65
66 fSurfaceArea = 1. ; // NOTE: not yet implemented!
67 }
68
69 G4TwistTubsFlatSide::G4TwistTubsFlatSide( const G4String& name,
70 G4double EndInnerRadius[2],
71 G4double EndOuterRadius[2],
72 G4double DPhi,
73 G4double EndPhi[2],
74 G4double EndZ[2],
75 G4int handedness )
76 : G4VTwistSurface(name)
77 {
78 fHandedness = handedness; // +z = +ve, -z = -ve
79 fAxis[0] = kRho; // in local coordinate system
80 fAxis[1] = kPhi;
81 G4int i = (handedness < 0 ? 0 : 1);
82 fAxisMin[0] = EndInnerRadius[i]; // Inner-hype radius at z=0
83 fAxisMax[0] = EndOuterRadius[i]; // Outer-hype radius at z=0
84 fAxisMin[1] = -0.5*DPhi;
85 fAxisMax[1] = -fAxisMin[1];
86 fCurrentNormal.normal.set(0, 0, (fHandedness < 0 ? -1 : 1));
87 // Unit vector, in local coordinate system
88 fRot.rotateZ(EndPhi[i]);
89 fTrans.set(0, 0, EndZ[i]);
90 fIsValidNorm = true;
91
92 SetCorners();
93 SetBoundaries();
94
95 fSurfaceArea = 0.5*DPhi * (EndOuterRadius[i]*EndOuterRadius[i]
96 - EndInnerRadius[i]*EndInnerRadius[i] ) ;
97 }
98
99 //=====================================================================
100 //* Fake default constructor ------------------------------------------
101
102 G4TwistTubsFlatSide::G4TwistTubsFlatSide( __void__& a )
103 : G4VTwistSurface(a)
104 {
105 }
106
107 //=====================================================================
108 //* destructor --------------------------------------------------------
109
110 G4TwistTubsFlatSide::~G4TwistTubsFlatSide() = default;
111
112 //=====================================================================
113 //* GetNormal ---------------------------------------------------------
114
115 G4ThreeVector G4TwistTubsFlatSide::GetNormal(const G4ThreeVector& /* xx */ ,
116 G4bool isGlobal)
117 {
118 if (isGlobal)
119 {
120 return ComputeGlobalDirection(fCurrentNormal.normal);
121 }
122 else
123 {
124 return fCurrentNormal.normal;
125 }
126 }
127
128 //=====================================================================
129 //* DistanceToSurface(p, v) -------------------------------------------
130
131 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector& gp,
132 const G4ThreeVector& gv,
133 G4ThreeVector gxx[],
134 G4double distance[],
135 G4int areacode[],
136 G4bool isvalid[],
137 EValidate validate)
138 {
139 fCurStatWithV.ResetfDone(validate, &gp, &gv);
140
141 if (fCurStatWithV.IsDone())
142 {
143 for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
144 {
145 gxx[i] = fCurStatWithV.GetXX(i);
146 distance[i] = fCurStatWithV.GetDistance(i);
147 areacode[i] = fCurStatWithV.GetAreacode(i);
148 isvalid[i] = fCurStatWithV.IsValid(i);
149 }
150 return fCurStatWithV.GetNXX();
151 }
152 else // initialize
153 {
154 for (G4int i=0; i<2; ++i)
155 {
156 distance[i] = kInfinity;
157 areacode[i] = sOutside;
158 isvalid[i] = false;
159 gxx[i].set(kInfinity, kInfinity, kInfinity);
160 }
161 }
162
163 G4ThreeVector p = ComputeLocalPoint(gp);
164 G4ThreeVector v = ComputeLocalDirection(gv);
165
166 //
167 // special case!
168 // if p is on surface, distance = 0.
169 //
170
171 if (std::fabs(p.z()) == 0.) // if p is on the plane
172 {
173 distance[0] = 0;
174 G4ThreeVector xx = p;
175 gxx[0] = ComputeGlobalPoint(xx);
176
177 if (validate == kValidateWithTol)
178 {
179 areacode[0] = GetAreaCode(xx);
180 if (!IsOutside(areacode[0]))
181 {
182 isvalid[0] = true;
183 }
184 }
185 else if (validate == kValidateWithoutTol)
186 {
187 areacode[0] = GetAreaCode(xx, false);
188 if (IsInside(areacode[0]))
189 {
190 isvalid[0] = true;
191 }
192 }
193 else // kDontValidate
194 {
195 areacode[0] = sInside;
196 isvalid[0] = true;
197 }
198 return 1;
199 }
200 //
201 // special case end
202 //
203
204 if (v.z() == 0)
205 {
206 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
207 isvalid[0], 0, validate, &gp, &gv);
208 return 0;
209 }
210
211 distance[0] = - (p.z() / v.z());
212
213 G4ThreeVector xx = p + distance[0]*v;
214 gxx[0] = ComputeGlobalPoint(xx);
215
216 if (validate == kValidateWithTol)
217 {
218 areacode[0] = GetAreaCode(xx);
219 if (!IsOutside(areacode[0]))
220 {
221 if (distance[0] >= 0) isvalid[0] = true;
222 }
223 }
224 else if (validate == kValidateWithoutTol)
225 {
226 areacode[0] = GetAreaCode(xx, false);
227 if (IsInside(areacode[0]))
228 {
229 if (distance[0] >= 0) isvalid[0] = true;
230 }
231 }
232 else // kDontValidate
233 {
234 areacode[0] = sInside;
235 if (distance[0] >= 0) isvalid[0] = true;
236 }
237
238 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
239 isvalid[0], 1, validate, &gp, &gv);
240
241 #ifdef G4TWISTDEBUG
242 G4cerr << "ERROR - G4TwistTubsFlatSide::DistanceToSurface(p,v)" << G4endl;
243 G4cerr << " Name : " << GetName() << G4endl;
244 G4cerr << " xx : " << xx << G4endl;
245 G4cerr << " gxx[0] : " << gxx[0] << G4endl;
246 G4cerr << " dist[0] : " << distance[0] << G4endl;
247 G4cerr << " areacode[0] : " << areacode[0] << G4endl;
248 G4cerr << " isvalid[0] : " << isvalid[0] << G4endl;
249 #endif
250 return 1;
251 }
252
253 //=====================================================================
254 //* DistanceToSurface(p) ----------------------------------------------
255
256 G4int G4TwistTubsFlatSide::DistanceToSurface(const G4ThreeVector& gp,
257 G4ThreeVector gxx[],
258 G4double distance[],
259 G4int areacode[])
260 {
261 // Calculate distance to plane in local coordinate,
262 // then return distance and global intersection points.
263 //
264
265 fCurStat.ResetfDone(kDontValidate, &gp);
266
267 if (fCurStat.IsDone())
268 {
269 for (G4int i=0; i<fCurStat.GetNXX(); ++i)
270 {
271 gxx[i] = fCurStat.GetXX(i);
272 distance[i] = fCurStat.GetDistance(i);
273 areacode[i] = fCurStat.GetAreacode(i);
274 }
275 return fCurStat.GetNXX();
276 }
277 else // initialize
278 {
279 for (auto i=0; i<2; ++i)
280 {
281 distance[i] = kInfinity;
282 areacode[i] = sOutside;
283 gxx[i].set(kInfinity, kInfinity, kInfinity);
284 }
285 }
286
287 G4ThreeVector p = ComputeLocalPoint(gp);
288 G4ThreeVector xx;
289
290 // The plane is placed on origin with making its normal
291 // parallel to z-axis.
292 if (std::fabs(p.z()) <= 0.5 * kCarTolerance) // if p is on plane, return 1
293 {
294 distance[0] = 0;
295 xx = p;
296 }
297 else
298 {
299 distance[0] = std::fabs(p.z());
300 xx.set(p.x(), p.y(), 0);
301 }
302
303 gxx[0] = ComputeGlobalPoint(xx);
304 areacode[0] = sInside;
305 G4bool isvalid = true;
306 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
307 isvalid, 1, kDontValidate, &gp);
308 return 1;
309 }
310
311 //=====================================================================
312 //* GetAreaCode -------------------------------------------------------
313
314 G4int G4TwistTubsFlatSide::GetAreaCode(const G4ThreeVector &xx,
315 G4bool withTol)
316 {
317 const G4double rtol
318 = 0.5*G4GeometryTolerance::GetInstance()->GetRadialTolerance();
319
320 G4int areacode = sInside;
321
322 if (fAxis[0] == kRho && fAxis[1] == kPhi)
323 {
324 G4int rhoaxis = 0;
325
326 G4ThreeVector dphimin; // direction of phi-minimum boundary
327 G4ThreeVector dphimax; // direction of phi-maximum boundary
328 dphimin = GetCorner(sC0Max1Min);
329 dphimax = GetCorner(sC0Max1Max);
330
331 if (withTol)
332 {
333 G4bool isoutside = false;
334
335 // test boundary of rho-axis
336
337 if (xx.getRho() <= fAxisMin[rhoaxis] + rtol)
338 {
339 areacode |= (sAxis0 & (sAxisRho|sAxisMin)) | sBoundary; // rho-min
340 if (xx.getRho() < fAxisMin[rhoaxis] - rtol) isoutside = true;
341
342 }
343 else if (xx.getRho() >= fAxisMax[rhoaxis] - rtol)
344 {
345 areacode |= (sAxis0 & (sAxisRho|sAxisMax)) | sBoundary; // rho-max
346 if (xx.getRho() > fAxisMax[rhoaxis] + rtol) isoutside = true;
347 }
348
349 // test boundary of phi-axis
350
351 if (AmIOnLeftSide(xx, dphimin) >= 0) // xx is on dphimin
352 {
353 areacode |= (sAxis1 & (sAxisPhi | sAxisMin));
354 if ((areacode & sBoundary) != 0) areacode |= sCorner; // xx is on corner.
355 else areacode |= sBoundary;
356
357 if (AmIOnLeftSide(xx, dphimin) > 0) isoutside = true;
358
359 }
360 else if (AmIOnLeftSide(xx, dphimax) <= 0) // xx is on dphimax
361 {
362 areacode |= (sAxis1 & (sAxisPhi | sAxisMax));
363 if ((areacode & sBoundary) != 0) areacode |= sCorner; // xx is on corner.
364 else areacode |= sBoundary;
365
366 if (AmIOnLeftSide(xx, dphimax) < 0) isoutside = true;
367 }
368
369 // if isoutside = true, clear inside bit.
370 // if not on boundary, add axis information.
371
372 if (isoutside)
373 {
374 G4int tmpareacode = areacode & (~sInside);
375 areacode = tmpareacode;
376 }
377 else if ((areacode & sBoundary) != sBoundary)
378 {
379 areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
380 }
381
382 }
383 else
384 {
385 // out of boundary of rho-axis
386
387 if (xx.getRho() < fAxisMin[rhoaxis])
388 {
389 areacode |= (sAxis0 & (sAxisRho | sAxisMin)) | sBoundary;
390 }
391 else if (xx.getRho() > fAxisMax[rhoaxis])
392 {
393 areacode |= (sAxis0 & (sAxisRho | sAxisMax)) | sBoundary;
394 }
395
396 // out of boundary of phi-axis
397
398 if (AmIOnLeftSide(xx, dphimin, false) >= 0) // xx is leftside or
399 {
400 areacode |= (sAxis1 & (sAxisPhi | sAxisMin)) ; // boundary of dphimin
401 if ((areacode & sBoundary) != 0) areacode |= sCorner; // xx is on corner
402 else areacode |= sBoundary;
403
404 }
405 else if (AmIOnLeftSide(xx, dphimax, false) <= 0) // xx is rightside or
406 {
407 areacode |= (sAxis1 & (sAxisPhi | sAxisMax)); // boundary of dphimax
408 if ((areacode & sBoundary) != 0) areacode |= sCorner; // xx is on corner
409 else areacode |= sBoundary;
410
411 }
412
413 if ((areacode & sBoundary) != sBoundary) {
414 areacode |= (sAxis0 & sAxisRho) | (sAxis1 & sAxisPhi);
415 }
416
417 }
418 return areacode;
419 }
420 else
421 {
422 std::ostringstream message;
423 message << "Feature NOT implemented !" << G4endl
424 << " fAxis[0] = " << fAxis[0] << G4endl
425 << " fAxis[1] = " << fAxis[1];
426 G4Exception("G4TwistTubsFlatSide::GetAreaCode()", "GeomSolids0001",
427 FatalException, message);
428 }
429 return areacode;
430 }
431
432 //=====================================================================
433 //* SetCorners --------------------------------------------------------
434
435 void G4TwistTubsFlatSide::SetCorners()
436 {
437 // Set Corner points in local coodinate.
438
439 if (fAxis[0] == kRho && fAxis[1] == kPhi)
440 {
441 G4int rhoaxis = 0; // kRho
442 G4int phiaxis = 1; // kPhi
443
444 G4double x, y, z;
445 // corner of Axis0min and Axis1min
446 x = fAxisMin[rhoaxis]*std::cos(fAxisMin[phiaxis]);
447 y = fAxisMin[rhoaxis]*std::sin(fAxisMin[phiaxis]);
448 z = 0;
449 SetCorner(sC0Min1Min, x, y, z);
450 // corner of Axis0max and Axis1min
451 x = fAxisMax[rhoaxis]*std::cos(fAxisMin[phiaxis]);
452 y = fAxisMax[rhoaxis]*std::sin(fAxisMin[phiaxis]);
453 z = 0;
454 SetCorner(sC0Max1Min, x, y, z);
455 // corner of Axis0max and Axis1max
456 x = fAxisMax[rhoaxis]*std::cos(fAxisMax[phiaxis]);
457 y = fAxisMax[rhoaxis]*std::sin(fAxisMax[phiaxis]);
458 z = 0;
459 SetCorner(sC0Max1Max, x, y, z);
460 // corner of Axis0min and Axis1max
461 x = fAxisMin[rhoaxis]*std::cos(fAxisMax[phiaxis]);
462 y = fAxisMin[rhoaxis]*std::sin(fAxisMax[phiaxis]);
463 z = 0;
464 SetCorner(sC0Min1Max, x, y, z);
465
466 }
467 else
468 {
469 std::ostringstream message;
470 message << "Feature NOT implemented !" << G4endl
471 << " fAxis[0] = " << fAxis[0] << G4endl
472 << " fAxis[1] = " << fAxis[1];
473 G4Exception("G4TwistTubsFlatSide::SetCorners()", "GeomSolids0001",
474 FatalException, message);
475 }
476 }
477
478 //=====================================================================
479 //* SetBoundaries() ---------------------------------------------------
480
481 void G4TwistTubsFlatSide::SetBoundaries()
482 {
483 // Set direction-unit vector of phi-boundary-lines in local coodinate.
484 // Don't call the function twice.
485
486 if (fAxis[0] == kRho && fAxis[1] == kPhi)
487 {
488 G4ThreeVector direction;
489 // sAxis0 & sAxisMin
490 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
491 direction = direction.unit();
492 SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction,
493 GetCorner(sC0Min1Min), sAxisPhi);
494
495 // sAxis0 & sAxisMax
496 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
497 direction = direction.unit();
498 SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction,
499 GetCorner(sC0Max1Min), sAxisPhi);
500
501 // sAxis1 & sAxisMin
502 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
503 direction = direction.unit();
504 SetBoundary(sAxis1 & (sAxisRho | sAxisMin), direction,
505 GetCorner(sC0Min1Min), sAxisRho);
506
507 // sAxis1 & sAxisMax
508 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
509 direction = direction.unit();
510 SetBoundary(sAxis1 & (sAxisRho | sAxisMax), direction,
511 GetCorner(sC0Min1Max), sAxisPhi);
512 }
513 else
514 {
515 std::ostringstream message;
516 message << "Feature NOT implemented !" << G4endl
517 << " fAxis[0] = " << fAxis[0] << G4endl
518 << " fAxis[1] = " << fAxis[1];
519 G4Exception("G4TwistTubsFlatSide::SetBoundaries()", "GeomSolids0001",
520 FatalException, message);
521 }
522 }
523
524 //=====================================================================
525 //* GetFacets() -------------------------------------------------------
526
527 void G4TwistTubsFlatSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
528 G4int faces[][4], G4int iside )
529 {
530 G4ThreeVector p ;
531
532 G4double rmin = fAxisMin[0] ;
533 G4double rmax = fAxisMax[0] ;
534 G4double phimin, phimax ;
535
536 G4double r,phi ;
537 G4int nnode,nface ;
538
539 for ( G4int i = 0 ; i<n ; ++i )
540 {
541 r = rmin + i*(rmax-rmin)/(n-1) ;
542
543 phimin = GetBoundaryMin(r) ;
544 phimax = GetBoundaryMax(r) ;
545
546 for ( G4int j = 0 ; j<k ; ++j )
547 {
548 phi = phimin + j*(phimax-phimin)/(k-1) ;
549
550 nnode = GetNode(i,j,k,n,iside) ;
551 p = SurfacePoint(phi,r,true) ; // surface point in global coord.system
552
553 xyz[nnode][0] = p.x() ;
554 xyz[nnode][1] = p.y() ;
555 xyz[nnode][2] = p.z() ;
556
557 if ( i<n-1 && j<k-1 ) // conterclock wise filling
558 {
559 nface = GetFace(i,j,k,n,iside) ;
560
561 if (fHandedness < 0) // lower side
562 {
563 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1)
564 * ( GetNode(i ,j ,k,n,iside)+1) ;
565 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1)
566 * ( GetNode(i ,j+1,k,n,iside)+1) ;
567 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1)
568 * ( GetNode(i+1,j+1,k,n,iside)+1) ;
569 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1)
570 * ( GetNode(i+1,j ,k,n,iside)+1) ;
571 }
572 else // upper side
573 {
574 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1)
575 * ( GetNode(i ,j ,k,n,iside)+1) ;
576 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1)
577 * ( GetNode(i+1,j ,k,n,iside)+1) ;
578 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1)
579 * ( GetNode(i+1,j+1,k,n,iside)+1) ;
580 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1)
581 * ( GetNode(i ,j+1,k,n,iside)+1) ;
582 }
583 }
584 }
585 }
586 }
587