| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4TwistTrapFlatSide implementation
27 //
28 // Author: 30-Aug-2002 - O.Link (Oliver.Link@c
29 // -------------------------------------------
30
31 #include "G4TwistTrapFlatSide.hh"
32
33 //============================================
34 //* constructors -----------------------------
35
36 G4TwistTrapFlatSide::G4TwistTrapFlatSide( cons
37
38
39
40
41
42
43
44
45
46 : G4VTwistSurface(name)
47 {
48 fHandedness = handedness; // +z = +ve, -z
49
50 fDx1 = pDx1 ;
51 fDx2 = pDx2 ;
52 fDy = pDy ;
53 fDz = pDz ;
54 fAlpha = pAlpha ;
55 fTAlph = std::tan(fAlpha) ;
56 fPhi = pPhi ;
57 fTheta = pTheta ;
58
59 fdeltaX = 2 * fDz * std::tan(fTheta) * std:
60 // dx in surface equation
61 fdeltaY = 2 * fDz * std::tan(fTheta) * std:
62 // dy in surface equation
63
64 fPhiTwist = PhiTwist ;
65
66 fCurrentNormal.normal.set( 0, 0, (fHandedne
67 // Unit vector, in local coordinate s
68 fRot.rotateZ( fHandedness > 0
69 ? 0.5 * fPhiTwist
70 : -0.5 * fPhiTwist );
71
72 fTrans.set(
73 fHandedness > 0 ? 0.5*fdeltaX :
74 fHandedness > 0 ? 0.5*fdeltaY :
75 fHandedness > 0 ? fDz : -fDz ) ;
76
77 fIsValidNorm = true;
78
79
80 fAxis[0] = kXAxis ;
81 fAxis[1] = kYAxis ;
82 fAxisMin[0] = kInfinity ; // x-Axis cannot
83 fAxisMax[0] = kInfinity ; // depends on y
84 fAxisMin[1] = -fDy ; // y - axis
85 fAxisMax[1] = fDy ;
86
87 SetCorners();
88 SetBoundaries();
89 }
90
91
92 //============================================
93 //* Fake default constructor -----------------
94
95 G4TwistTrapFlatSide::G4TwistTrapFlatSide( __vo
96 : G4VTwistSurface(a)
97 {
98 }
99
100
101 //============================================
102 //* destructor -------------------------------
103
104 G4TwistTrapFlatSide::~G4TwistTrapFlatSide() =
105
106 //============================================
107 //* GetNormal --------------------------------
108
109 G4ThreeVector G4TwistTrapFlatSide::GetNormal(c
110 G
111 {
112 if (isGlobal)
113 {
114 return ComputeGlobalDirection(fCurrentNo
115 }
116 else
117 {
118 return fCurrentNormal.normal;
119 }
120 }
121
122 //============================================
123 //* DistanceToSurface(p, v) ------------------
124
125 G4int G4TwistTrapFlatSide::DistanceToSurface(c
126 c
127
128
129
130
131
132 {
133 fCurStatWithV.ResetfDone(validate, &gp, &gv
134
135 if (fCurStatWithV.IsDone())
136 {
137 for (G4int i=0; i<fCurStatWithV.GetNXX()
138 {
139 gxx[i] = fCurStatWithV.GetXX(i);
140 distance[i] = fCurStatWithV.GetDistan
141 areacode[i] = fCurStatWithV.GetAreaco
142 isvalid[i] = fCurStatWithV.IsValid(i
143 }
144 return fCurStatWithV.GetNXX();
145 }
146 else // initialize
147 {
148 for (auto i=0; i<2; ++i)
149 {
150 distance[i] = kInfinity;
151 areacode[i] = sOutside;
152 isvalid[i] = false;
153 gxx[i].set(kInfinity, kInfinity, kInf
154 }
155 }
156
157 G4ThreeVector p = ComputeLocalPoint(gp);
158 G4ThreeVector v = ComputeLocalDirection(gv)
159
160 //
161 // special case!
162 // if p is on surface, distance = 0.
163 //
164
165 if (std::fabs(p.z()) == 0.) // if p is
166 {
167 distance[0] = 0;
168 G4ThreeVector xx = p;
169 gxx[0] = ComputeGlobalPoint(xx);
170
171 if (validate == kValidateWithTol)
172 {
173 areacode[0] = GetAreaCode(xx);
174 if (!IsOutside(areacode[0]))
175 {
176 isvalid[0] = true;
177 }
178 }
179 else if (validate == kValidateWithoutTol
180 {
181 areacode[0] = GetAreaCode(xx, false);
182 if (IsInside(areacode[0]))
183 {
184 isvalid[0] = true;
185 }
186 }
187 else // kDontValidate
188 {
189 areacode[0] = sInside;
190 isvalid[0] = true;
191 }
192 return 1;
193 }
194 //
195 // special case end
196 //
197
198 if (v.z() == 0) {
199
200 fCurStatWithV.SetCurrentStatus(0, gxx[0]
201 isvalid[0
202 return 0;
203 }
204
205 distance[0] = - (p.z() / v.z());
206
207 G4ThreeVector xx = p + distance[0]*v;
208 gxx[0] = ComputeGlobalPoint(xx);
209
210 if (validate == kValidateWithTol)
211 {
212 areacode[0] = GetAreaCode(xx);
213 if (!IsOutside(areacode[0]))
214 {
215 if (distance[0] >= 0) isvalid[0] = tr
216 }
217 }
218 else if (validate == kValidateWithoutTol)
219 {
220 areacode[0] = GetAreaCode(xx, false);
221 if (IsInside(areacode[0]))
222 {
223 if (distance[0] >= 0) isvalid[0] = tr
224 }
225 }
226 else // kDontValidate
227 {
228 areacode[0] = sInside;
229 if (distance[0] >= 0) isvalid[0] = tr
230 }
231
232 fCurStatWithV.SetCurrentStatus(0, gxx[0], d
233 isvalid[0],
234
235 #ifdef G4TWISTDEBUG
236 G4cerr << "ERROR - G4TwistTrapFlatSide::Dis
237 G4cerr << " Name : " << GetNa
238 G4cerr << " xx : " << xx <<
239 G4cerr << " gxx[0] : " << gxx[0
240 G4cerr << " dist[0] : " << dista
241 G4cerr << " areacode[0] : " << areac
242 G4cerr << " isvalid[0] : " << isval
243 #endif
244 return 1;
245 }
246
247 //============================================
248 //* DistanceToSurface(p) ---------------------
249
250 G4int G4TwistTrapFlatSide::DistanceToSurface(c
251 G
252 G
253 G
254 {
255 // Calculate distance to plane in local coo
256 // then return distance and global intersec
257 //
258
259 fCurStat.ResetfDone(kDontValidate, &gp);
260
261 if (fCurStat.IsDone())
262 {
263 for (G4int i=0; i<fCurStat.GetNXX(); ++i
264 {
265 gxx[i] = fCurStat.GetXX(i);
266 distance[i] = fCurStat.GetDistance(i)
267 areacode[i] = fCurStat.GetAreacode(i)
268 }
269 return fCurStat.GetNXX();
270 }
271 else // initialize
272 {
273 for (auto i=0; i<2; ++i)
274 {
275 distance[i] = kInfinity;
276 areacode[i] = sOutside;
277 gxx[i].set(kInfinity, kInfinity, kInf
278 }
279 }
280
281 G4ThreeVector p = ComputeLocalPoint(gp);
282 G4ThreeVector xx;
283
284 // The plane is placed on origin with makin
285 // parallel to z-axis.
286 if (std::fabs(p.z()) <= 0.5 * kCarTolerance
287 { // if p is on the plane, return 1
288 distance[0] = 0;
289 xx = p;
290 }
291 else
292 {
293 distance[0] = std::fabs(p.z());
294 xx.set(p.x(), p.y(), 0);
295 }
296
297 gxx[0] = ComputeGlobalPoint(xx);
298 areacode[0] = sInside;
299 G4bool isvalid = true;
300 fCurStat.SetCurrentStatus(0, gxx[0], distan
301 isvalid, 1, kDont
302 return 1;
303
304 }
305
306 //============================================
307 //* GetAreaCode() ----------------------------
308
309 G4int G4TwistTrapFlatSide::GetAreaCode(const G
310 G
311 {
312
313 static const G4double ctol = 0.5 * kCarToler
314 G4int areacode = sInside;
315
316 if (fAxis[0] == kXAxis && fAxis[1] == kYAxis
317 {
318 G4int yaxis = 1;
319
320 G4double wmax = xAxisMax(xx.y(), fTAlph )
321 G4double wmin = -xAxisMax(xx.y(), -fTAlph
322
323 if (withTol)
324 {
325 G4bool isoutside = false;
326
327 // test boundary of x-axis
328
329 if (xx.x() < wmin + ctol)
330 {
331 areacode |= (sAxis0 & (sAxisX | sAxisM
332 if (xx.x() <= wmin - ctol) isoutside =
333
334 }
335 else if (xx.x() > wmax - ctol)
336 {
337 areacode |= (sAxis0 & (sAxisX | sAxisM
338 if (xx.x() >= wmax + ctol) isoutside
339 }
340
341 // test boundary of y-axis
342
343 if (xx.y() < fAxisMin[yaxis] + ctol)
344 {
345 areacode |= (sAxis1 & (sAxisY | sAxisM
346
347 if ((areacode & sBoundary) != 0) are
348 else areacode |
349 if (xx.y() <= fAxisMin[yaxis] - ctol)
350
351 }
352 else if (xx.y() > fAxisMax[yaxis] - ctol
353 {
354 areacode |= (sAxis1 & (sAxisY | sAxisM
355
356 if ((areacode & sBoundary) != 0) are
357 else areacode |
358 if (xx.y() >= fAxisMax[yaxis] + ctol)
359 }
360
361 // if isoutside = true, clear inside bit
362 // if not on boundary, add axis informat
363
364 if (isoutside)
365 {
366 G4int tmpareacode = areacode & (~sInsi
367 areacode = tmpareacode;
368 }
369 else if ((areacode & sBoundary) != sBoun
370 {
371 areacode |= (sAxis0 & sAxisX) | (sAxis
372 }
373 }
374 else
375 {
376 // boundary of x-axis
377
378 if (xx.x() < wmin )
379 {
380 areacode |= (sAxis0 & (sAxisX | sAxisM
381 }
382 else if (xx.x() > wmax)
383 {
384 areacode |= (sAxis0 & (sAxisX | sAxisM
385 }
386
387 // boundary of y-axis
388
389 if (xx.y() < fAxisMin[yaxis])
390 {
391 areacode |= (sAxis1 & (sAxisY | sAxisM
392 if ((areacode & sBoundary) != 0) are
393 else areacode |
394
395 }
396 else if (xx.y() > fAxisMax[yaxis])
397 {
398 areacode |= (sAxis1 & (sAxisY | sAxisM
399 if ((areacode & sBoundary) != 0) are
400 else areacode |
401 }
402
403 if ((areacode & sBoundary) != sBoundary)
404 {
405 areacode |= (sAxis0 & sAxisX) | (sAxis
406 }
407 }
408 return areacode;
409 }
410 else
411 {
412 G4Exception("G4TwistTrapFlatSide::GetAreaC
413 "GeomSolids0001", FatalExcepti
414 "Feature NOT implemented !");
415 }
416
417 return areacode;
418 }
419
420 //============================================
421 //* SetCorners -------------------------------
422
423 void G4TwistTrapFlatSide::SetCorners()
424 {
425 // Set Corner points in local coodinate.
426
427 if (fAxis[0] == kXAxis && fAxis[1] == kYAxi
428 {
429 G4double x, y, z;
430
431 // corner of Axis0min and Axis1min
432 x = -fDx1 + fDy * fTAlph ;
433 y = -fDy ;
434 z = 0 ;
435 SetCorner(sC0Min1Min, x, y, z);
436
437 // corner of Axis0max and Axis1min
438 x = fDx1 + fDy * fTAlph ;
439 y = -fDy ;
440 z = 0 ;
441 SetCorner(sC0Max1Min, x, y, z);
442
443 // corner of Axis0max and Axis1max
444 x = fDx2 - fDy * fTAlph ;
445 y = fDy ;
446 z = 0 ;
447 SetCorner(sC0Max1Max, x, y, z);
448
449 // corner of Axis0min and Axis1max
450 x = -fDx2 - fDy * fTAlph ;
451 y = fDy ;
452 z = 0 ;
453 SetCorner(sC0Min1Max, x, y, z);
454
455 }
456 else
457 {
458 std::ostringstream message;
459 message << "Feature NOT implemented !" <<
460 << " fAxis[0] = " << fAxis
461 << " fAxis[1] = " << fAxis
462 G4Exception("G4TwistTrapFlatSide::SetCorn
463 "GeomSolids0001", FatalExcept
464 }
465 }
466
467 //============================================
468 //* SetBoundaries() --------------------------
469
470 void G4TwistTrapFlatSide::SetBoundaries()
471 {
472 // Set direction-unit vector of phi-boundar
473 // Don't call the function twice.
474
475 G4ThreeVector direction ;
476
477 if (fAxis[0] == kXAxis && fAxis[1] == kYAxis
478 {
479 // sAxis0 & sAxisMin
480 direction = - ( GetCorner(sC0Min1Max) - Ge
481 direction = direction.unit();
482 SetBoundary(sAxis0 & (sAxisX | sAxisMin),
483 GetCorner(sC0Min1Max), sAxisY)
484
485 // sAxis0 & sAxisMax
486 direction = GetCorner(sC0Max1Max) - GetCor
487 direction = direction.unit();
488 SetBoundary(sAxis0 & (sAxisX | sAxisMax),
489 GetCorner(sC0Max1Min), sAxisY)
490
491 // sAxis1 & sAxisMin
492 direction = GetCorner(sC0Max1Min) - GetCor
493 direction = direction.unit();
494 SetBoundary(sAxis1 & (sAxisY | sAxisMin),
495 GetCorner(sC0Min1Min), sAxisX)
496
497 // sAxis1 & sAxisMax
498 direction = - ( GetCorner(sC0Max1Max) - Ge
499 direction = direction.unit();
500 SetBoundary(sAxis1 & (sAxisY | sAxisMax),
501 GetCorner(sC0Max1Max), sAxisX)
502
503 }
504 else
505 {
506 std::ostringstream message;
507 message << "Feature NOT implemented !" <<
508 << " fAxis[0] = " << fAxis[
509 << " fAxis[1] = " << fAxis[
510 G4Exception("G4TwistTrapFlatSide::SetCorne
511 "GeomSolids0001", FatalExcepti
512 }
513 }
514
515 //============================================
516 //* GetFacets() ------------------------------
517
518 void G4TwistTrapFlatSide::GetFacets( G4int k,
519 G4int fac
520 {
521 G4double x,y ; // the two parameters
522 G4ThreeVector p ; // a point on the surface
523
524 G4int nnode ;
525 G4int nface ;
526
527 G4double xmin,xmax ;
528
529 // calculate the (n-1)*(k-1) vertices
530
531 for ( G4int i = 0 ; i<n ; ++i )
532 {
533 y = -fDy + i*(2*fDy)/(n-1) ;
534
535 for ( G4int j = 0 ; j<k ; ++j )
536 {
537 xmin = GetBoundaryMin(y) ;
538 xmax = GetBoundaryMax(y) ;
539 x = xmin + j*(xmax-xmin)/(k-1) ;
540
541 nnode = GetNode(i,j,k,n,iside) ;
542 p = SurfacePoint(x,y,true) ; // surface
543
544 xyz[nnode][0] = p.x() ;
545 xyz[nnode][1] = p.y() ;
546 xyz[nnode][2] = p.z() ;
547
548 if ( i<n-1 && j<k-1 )
549 {
550 nface = GetFace(i,j,k,n,iside) ;
551
552 if (fHandedness < 0) // lower side
553 {
554 faces[nface][0] = GetEdgeVisibility(
555 * ( GetNode(i ,j ,
556 faces[nface][1] = GetEdgeVisibility(
557 * ( GetNode(i+1,j ,
558 faces[nface][2] = GetEdgeVisibility(
559 * ( GetNode(i+1,j+1,
560 faces[nface][3] = GetEdgeVisibility(
561 * ( GetNode(i ,j+1,
562 }
563 else // upper side
564 {
565 faces[nface][0] = GetEdgeVisibility(
566 * ( GetNode(i ,j ,
567 faces[nface][1] = GetEdgeVisibility(
568 * ( GetNode(i ,j+1,
569 faces[nface][2] = GetEdgeVisibility(
570 * ( GetNode(i+1,j+1,
571 faces[nface][3] = GetEdgeVisibility(
572 * ( GetNode(i+1,j ,
573 }
574 }
575 }
576 }
577 }
578