Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
HepPolyhedron.cc
Go to the documentation of this file.
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// G4 Polyhedron library
27//
28// History:
29// 23.07.96 E.Chernyaev <[email protected]> - initial version
30//
31// 30.09.96 E.Chernyaev
32// - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
33// - added GetNextUnitNormal, GetNextEdgeIndices, GetNextEdge
34//
35// 15.12.96 E.Chernyaev
36// - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
37// - rewritten G4PolyhedronCons;
38// - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
39//
40// 01.06.97 E.Chernyaev
41// - modified RotateAroundZ, added SetSideFacets
42//
43// 19.03.00 E.Chernyaev
44// - implemented boolean operations (add, subtract, intersect) on polyhedra;
45//
46// 25.05.01 E.Chernyaev
47// - added GetSurfaceArea() and GetVolume()
48//
49// 05.11.02 E.Chernyaev
50// - added createTwistedTrap() and createPolyhedron()
51//
52// 20.06.05 G.Cosmo
53// - added HepPolyhedronEllipsoid
54//
55// 18.07.07 T.Nikitina
56// - added HepPolyhedronParaboloid
57//
58// 22.02.20 E.Chernyaev
59// - added HepPolyhedronTet, HepPolyhedronHyberbolicMirror
60//
61// 12.05.21 E.Chernyaev
62// - added TriangulatePolygon(), RotateContourAroundZ()
63// - added HepPolyhedronPgon, HepPolyhedronPcon given by rz-contour
64//
65// 26.03.22 E.Chernyaev
66// - added SetVertex(), SetFacet()
67// - added HepPolyhedronTetMesh
68//
69// 04.04.22 E.Chernyaev
70// - added JoinCoplanarFacets()
71//
72// 07.04.22 E.Chernyaev
73// - added HepPolyhedronBoxMesh
74
75#include "HepPolyhedron.h"
76#include "G4PhysicalConstants.hh"
77#include "G4Vector3D.hh"
78
79#include <cstdlib> // Required on some compilers for std::abs(int) ...
80#include <cmath>
81#include <algorithm>
82
83using CLHEP::perMillion;
84using CLHEP::deg;
85using CLHEP::pi;
86using CLHEP::twopi;
87using CLHEP::nm;
88const G4double spatialTolerance = 0.01*nm;
89
90/***********************************************************************
91 * *
92 * Name: HepPolyhedron operator << Date: 09.05.96 *
93 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
94 * *
95 * Function: Print contents of G4 polyhedron *
96 * *
97 ***********************************************************************/
98std::ostream & operator<<(std::ostream & ostr, const G4Facet & facet) {
99 for (const auto& edge : facet.edge) {
100 ostr << " " << edge.v << "/" << edge.f;
101 }
102 return ostr;
103}
104
105std::ostream & operator<<(std::ostream & ostr, const HepPolyhedron & ph) {
106 ostr << std::endl;
107 ostr << "Nvertices=" << ph.nvert << ", Nfacets=" << ph.nface << std::endl;
108 G4int i;
109 for (i=1; i<=ph.nvert; i++) {
110 ostr << "xyz(" << i << ")="
111 << ph.pV[i].x() << ' ' << ph.pV[i].y() << ' ' << ph.pV[i].z()
112 << std::endl;
113 }
114 for (i=1; i<=ph.nface; i++) {
115 ostr << "face(" << i << ")=" << ph.pF[i] << std::endl;
116 }
117 return ostr;
118}
119
120HepPolyhedron::HepPolyhedron(G4int Nvert, G4int Nface)
121/***********************************************************************
122 * *
123 * Name: HepPolyhedron constructor with Date: 26.03.2022 *
124 * allocation of memory Revised: *
125 * Author: E.Tcherniaev (E.Chernyaev) *
126 * *
127 ***********************************************************************/
128: nvert(0), nface(0), pV(nullptr), pF(nullptr)
129{
130 AllocateMemory(Nvert, Nface);
131}
132
133HepPolyhedron::HepPolyhedron(const HepPolyhedron &from)
134/***********************************************************************
135 * *
136 * Name: HepPolyhedron copy constructor Date: 23.07.96 *
137 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
138 * *
139 ***********************************************************************/
140: nvert(0), nface(0), pV(nullptr), pF(nullptr)
141{
142 AllocateMemory(from.nvert, from.nface);
143 for (G4int i=1; i<=nvert; i++) pV[i] = from.pV[i];
144 for (G4int k=1; k<=nface; k++) pF[k] = from.pF[k];
145}
146
147HepPolyhedron::HepPolyhedron(HepPolyhedron&& from)
148/***********************************************************************
149 * *
150 * Name: HepPolyhedron move constructor Date: 04.11.2019 *
151 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
152 * *
153 ***********************************************************************/
154: nvert(0), nface(0), pV(nullptr), pF(nullptr)
155{
156 nvert = from.nvert;
157 nface = from.nface;
158 pV = from.pV;
159 pF = from.pF;
160
161 // Release the data from the source object
162 from.nvert = 0;
163 from.nface = 0;
164 from.pV = nullptr;
165 from.pF = nullptr;
166}
167
168HepPolyhedron & HepPolyhedron::operator=(const HepPolyhedron &from)
169/***********************************************************************
170 * *
171 * Name: HepPolyhedron operator = Date: 23.07.96 *
172 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
173 * *
174 * Function: Copy contents of one polyhedron to another *
175 * *
176 ***********************************************************************/
177{
178 if (this != &from) {
179 AllocateMemory(from.nvert, from.nface);
180 for (G4int i=1; i<=nvert; i++) pV[i] = from.pV[i];
181 for (G4int k=1; k<=nface; k++) pF[k] = from.pF[k];
182 }
183 return *this;
184}
185
186HepPolyhedron & HepPolyhedron::operator=(HepPolyhedron&& from)
187/***********************************************************************
188 * *
189 * Name: HepPolyhedron move operator = Date: 04.11.2019 *
190 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
191 * *
192 * Function: Move contents of one polyhedron to another *
193 * *
194 ***********************************************************************/
195{
196 if (this != &from) {
197 delete [] pV;
198 delete [] pF;
199 nvert = from.nvert;
200 nface = from.nface;
201 pV = from.pV;
202 pF = from.pF;
203
204 // Release the data from the source object
205 from.nvert = 0;
206 from.nface = 0;
207 from.pV = nullptr;
208 from.pF = nullptr;
209 }
210 return *this;
211}
212
213G4int
214HepPolyhedron::FindNeighbour(G4int iFace, G4int iNode, G4int iOrder) const
215/***********************************************************************
216 * *
217 * Name: HepPolyhedron::FindNeighbour Date: 22.11.99 *
218 * Author: E.Chernyaev Revised: *
219 * *
220 * Function: Find neighbouring face *
221 * *
222 ***********************************************************************/
223{
224 G4int i;
225 for (i=0; i<4; i++) {
226 if (iNode == std::abs(pF[iFace].edge[i].v)) break;
227 }
228 if (i == 4) {
229 std::cerr
230 << "HepPolyhedron::FindNeighbour: face " << iFace
231 << " has no node " << iNode
232 << std::endl;
233 return 0;
234 }
235 if (iOrder < 0) {
236 if ( --i < 0) i = 3;
237 if (pF[iFace].edge[i].v == 0) i = 2;
238 }
239 return (pF[iFace].edge[i].v > 0) ? 0 : pF[iFace].edge[i].f;
240}
241
242G4Normal3D HepPolyhedron::FindNodeNormal(G4int iFace, G4int iNode) const
243/***********************************************************************
244 * *
245 * Name: HepPolyhedron::FindNodeNormal Date: 22.11.99 *
246 * Author: E.Chernyaev Revised: *
247 * *
248 * Function: Find normal at given node *
249 * *
250 ***********************************************************************/
251{
252 G4Normal3D normal = GetUnitNormal(iFace);
253 G4int k = iFace, iOrder = 1;
254
255 for(;;) {
256 k = FindNeighbour(k, iNode, iOrder);
257 if (k == iFace) break;
258 if (k > 0) {
259 normal += GetUnitNormal(k);
260 }else{
261 if (iOrder < 0) break;
262 k = iFace;
263 iOrder = -iOrder;
264 }
265 }
266 return normal.unit();
267}
268
269G4int HepPolyhedron::GetNumberOfRotationSteps()
270/***********************************************************************
271 * *
272 * Name: HepPolyhedron::GetNumberOfRotationSteps Date: 24.06.97 *
273 * Author: J.Allison (Manchester University) Revised: *
274 * *
275 * Function: Get number of steps for whole circle *
276 * *
277 ***********************************************************************/
278{
279 return fNumberOfRotationSteps;
280}
281
282void HepPolyhedron::SetVertex(G4int index, const G4Point3D& v)
283/***********************************************************************
284 * *
285 * Name: HepPolyhedron::SetVertex Date: 26.03.22 *
286 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
287 * *
288 * Function: Set vertex *
289 * *
290 ***********************************************************************/
291{
292 if (index < 1 || index > nvert)
293 {
294 std::cerr
295 << "HepPolyhedron::SetVertex: vertex index = " << index
296 << " is out of range\n"
297 << " N. of vertices = " << nvert << "\n"
298 << " N. of facets = " << nface << std::endl;
299 return;
300 }
301 pV[index] = v;
302}
303
304void
305HepPolyhedron::SetFacet(G4int index, G4int iv1, G4int iv2, G4int iv3, G4int iv4)
306/***********************************************************************
307 * *
308 * Name: HepPolyhedron::SetFacet Date: 26.03.22 *
309 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
310 * *
311 * Function: Set facet *
312 * *
313 ***********************************************************************/
314{
315 if (index < 1 || index > nface)
316 {
317 std::cerr
318 << "HepPolyhedron::SetFacet: facet index = " << index
319 << " is out of range\n"
320 << " N. of vertices = " << nvert << "\n"
321 << " N. of facets = " << nface << std::endl;
322 return;
323 }
324 if (iv1 < 1 || iv1 > nvert ||
325 iv2 < 1 || iv2 > nvert ||
326 iv3 < 1 || iv3 > nvert ||
327 iv4 < 0 || iv4 > nvert)
328 {
329 std::cerr
330 << "HepPolyhedron::SetFacet: incorrectly specified facet"
331 << " (" << iv1 << ", " << iv2 << ", " << iv3 << ", " << iv4 << ")\n"
332 << " N. of vertices = " << nvert << "\n"
333 << " N. of facets = " << nface << std::endl;
334 return;
335 }
336 pF[index] = G4Facet(iv1, 0, iv2, 0, iv3, 0, iv4, 0);
337}
338
339void HepPolyhedron::SetNumberOfRotationSteps(G4int n)
340/***********************************************************************
341 * *
342 * Name: HepPolyhedron::SetNumberOfRotationSteps Date: 24.06.97 *
343 * Author: J.Allison (Manchester University) Revised: *
344 * *
345 * Function: Set number of steps for whole circle *
346 * *
347 ***********************************************************************/
348{
349 const G4int nMin = 3;
350 if (n < nMin) {
351 std::cerr
352 << "HepPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
353 << "number of steps per circle < " << nMin << "; forced to " << nMin
354 << std::endl;
355 fNumberOfRotationSteps = nMin;
356 }else{
357 fNumberOfRotationSteps = n;
358 }
359}
360
361void HepPolyhedron::ResetNumberOfRotationSteps()
362/***********************************************************************
363 * *
364 * Name: HepPolyhedron::GetNumberOfRotationSteps Date: 24.06.97 *
365 * Author: J.Allison (Manchester University) Revised: *
366 * *
367 * Function: Reset number of steps for whole circle to default value *
368 * *
369 ***********************************************************************/
370{
371 fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
372}
373
374void HepPolyhedron::AllocateMemory(G4int Nvert, G4int Nface)
375/***********************************************************************
376 * *
377 * Name: HepPolyhedron::AllocateMemory Date: 19.06.96 *
378 * Author: E.Chernyaev (IHEP/Protvino) Revised: 05.11.02 *
379 * *
380 * Function: Allocate memory for GEANT4 polyhedron *
381 * *
382 * Input: Nvert - number of nodes *
383 * Nface - number of faces *
384 * *
385 ***********************************************************************/
386{
387 if (nvert == Nvert && nface == Nface) return;
388 delete [] pV;
389 delete [] pF;
390 if (Nvert > 0 && Nface > 0) {
391 nvert = Nvert;
392 nface = Nface;
393 pV = new G4Point3D[nvert+1];
394 pF = new G4Facet[nface+1];
395 }else{
396 nvert = 0; nface = 0; pV = nullptr; pF = nullptr;
397 }
398}
399
400void HepPolyhedron::CreatePrism()
401/***********************************************************************
402 * *
403 * Name: HepPolyhedron::CreatePrism Date: 15.07.96 *
404 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
405 * *
406 * Function: Set facets for a prism *
407 * *
408 ***********************************************************************/
409{
410 enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
411
412 pF[1] = G4Facet(1,LEFT, 4,BACK, 3,RIGHT, 2,FRONT);
413 pF[2] = G4Facet(5,TOP, 8,BACK, 4,BOTTOM, 1,FRONT);
414 pF[3] = G4Facet(8,TOP, 7,RIGHT, 3,BOTTOM, 4,LEFT);
415 pF[4] = G4Facet(7,TOP, 6,FRONT, 2,BOTTOM, 3,BACK);
416 pF[5] = G4Facet(6,TOP, 5,LEFT, 1,BOTTOM, 2,RIGHT);
417 pF[6] = G4Facet(5,FRONT, 6,RIGHT, 7,BACK, 8,LEFT);
418}
419
420void HepPolyhedron::RotateEdge(G4int k1, G4int k2, G4double r1, G4double r2,
421 G4int v1, G4int v2, G4int vEdge,
422 G4bool ifWholeCircle, G4int nds, G4int &kface)
423/***********************************************************************
424 * *
425 * Name: HepPolyhedron::RotateEdge Date: 05.12.96 *
426 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
427 * *
428 * Function: Create set of facets by rotation of an edge around Z-axis *
429 * *
430 * Input: k1, k2 - end vertices of the edge *
431 * r1, r2 - radiuses of the end vertices *
432 * v1, v2 - visibility of edges produced by rotation of the end *
433 * vertices *
434 * vEdge - visibility of the edge *
435 * ifWholeCircle - is true in case of whole circle rotation *
436 * nds - number of discrete steps *
437 * r[] - r-coordinates *
438 * kface - current free cell in the pF array *
439 * *
440 ***********************************************************************/
441{
442 if (r1 == 0. && r2 == 0.) return;
443
444 G4int i;
445 G4int i1 = k1;
446 G4int i2 = k2;
447 G4int ii1 = ifWholeCircle ? i1 : i1+nds;
448 G4int ii2 = ifWholeCircle ? i2 : i2+nds;
449 G4int vv = ifWholeCircle ? vEdge : 1;
450
451 if (nds == 1) {
452 if (r1 == 0.) {
453 pF[kface++] = G4Facet(i1,0, v2*i2,0, (i2+1),0);
454 }else if (r2 == 0.) {
455 pF[kface++] = G4Facet(i1,0, i2,0, v1*(i1+1),0);
456 }else{
457 pF[kface++] = G4Facet(i1,0, v2*i2,0, (i2+1),0, v1*(i1+1),0);
458 }
459 }else{
460 if (r1 == 0.) {
461 pF[kface++] = G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0);
462 for (i2++,i=1; i<nds-1; i2++,i++) {
463 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
464 }
465 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
466 }else if (r2 == 0.) {
467 pF[kface++] = G4Facet(vv*i1,0, vEdge*i2,0, v1*(i1+1),0);
468 for (i1++,i=1; i<nds-1; i1++,i++) {
469 pF[kface++] = G4Facet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
470 }
471 pF[kface++] = G4Facet(vEdge*i1,0, vv*i2,0, v1*ii1,0);
472 }else{
473 pF[kface++] = G4Facet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
474 for (i1++,i2++,i=1; i<nds-1; i1++,i2++,i++) {
475 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
476 }
477 pF[kface++] = G4Facet(vEdge*i1,0, v2*i2,0, vv*ii2,0, v1*ii1,0);
478 }
479 }
480}
481
482void HepPolyhedron::SetSideFacets(G4int ii[4], G4int vv[4],
483 G4int *kk, G4double *r,
484 G4double dphi, G4int nds, G4int &kface)
485/***********************************************************************
486 * *
487 * Name: HepPolyhedron::SetSideFacets Date: 20.05.97 *
488 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
489 * *
490 * Function: Set side facets for the case of incomplete rotation *
491 * *
492 * Input: ii[4] - indices of original vertices *
493 * vv[4] - visibility of edges *
494 * kk[] - indices of nodes *
495 * r[] - radiuses *
496 * dphi - delta phi *
497 * nds - number of discrete steps *
498 * kface - current free cell in the pF array *
499 * *
500 ***********************************************************************/
501{
502 G4int k1, k2, k3, k4;
503
504 if (std::abs(dphi-pi) < perMillion) { // half a circle
505 for (G4int i=0; i<4; i++) {
506 k1 = ii[i];
507 k2 = ii[(i+1)%4];
508 if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;
509 }
510 }
511
512 if (ii[1] == ii[2]) {
513 k1 = kk[ii[0]];
514 k2 = kk[ii[2]];
515 k3 = kk[ii[3]];
516 pF[kface++] = G4Facet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
517 if (r[ii[0]] != 0.) k1 += nds;
518 if (r[ii[2]] != 0.) k2 += nds;
519 if (r[ii[3]] != 0.) k3 += nds;
520 pF[kface++] = G4Facet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
521 }else if (kk[ii[0]] == kk[ii[1]]) {
522 k1 = kk[ii[0]];
523 k2 = kk[ii[2]];
524 k3 = kk[ii[3]];
525 pF[kface++] = G4Facet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
526 if (r[ii[0]] != 0.) k1 += nds;
527 if (r[ii[2]] != 0.) k2 += nds;
528 if (r[ii[3]] != 0.) k3 += nds;
529 pF[kface++] = G4Facet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
530 }else if (kk[ii[2]] == kk[ii[3]]) {
531 k1 = kk[ii[0]];
532 k2 = kk[ii[1]];
533 k3 = kk[ii[2]];
534 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
535 if (r[ii[0]] != 0.) k1 += nds;
536 if (r[ii[1]] != 0.) k2 += nds;
537 if (r[ii[2]] != 0.) k3 += nds;
538 pF[kface++] = G4Facet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
539 }else{
540 k1 = kk[ii[0]];
541 k2 = kk[ii[1]];
542 k3 = kk[ii[2]];
543 k4 = kk[ii[3]];
544 pF[kface++] = G4Facet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
545 if (r[ii[0]] != 0.) k1 += nds;
546 if (r[ii[1]] != 0.) k2 += nds;
547 if (r[ii[2]] != 0.) k3 += nds;
548 if (r[ii[3]] != 0.) k4 += nds;
549 pF[kface++] = G4Facet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
550 }
551}
552
553void HepPolyhedron::RotateAroundZ(G4int nstep, G4double phi, G4double dphi,
554 G4int np1, G4int np2,
555 const G4double *z, G4double *r,
556 G4int nodeVis, G4int edgeVis)
557/***********************************************************************
558 * *
559 * Name: HepPolyhedron::RotateAroundZ Date: 27.11.96 *
560 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
561 * *
562 * Function: Create HepPolyhedron for a solid produced by rotation of *
563 * two polylines around Z-axis *
564 * *
565 * Input: nstep - number of discrete steps, if 0 then default *
566 * phi - starting phi angle *
567 * dphi - delta phi *
568 * np1 - number of points in external polyline *
569 * (must be negative in case of closed polyline) *
570 * np2 - number of points in internal polyline (may be 1) *
571 * z[] - z-coordinates (+z >>> -z for both polylines) *
572 * r[] - r-coordinates *
573 * nodeVis - how to Draw edges joing consecutive positions of *
574 * node during rotation *
575 * edgeVis - how to Draw edges *
576 * *
577 ***********************************************************************/
578{
579 static const G4double wholeCircle = twopi;
580
581 // S E T R O T A T I O N P A R A M E T E R S
582
583 G4bool ifWholeCircle = std::abs(dphi-wholeCircle) < perMillion;
584 G4double delPhi = ifWholeCircle ? wholeCircle : dphi;
585 G4int nSphi = nstep;
586 if (nSphi <= 0) nSphi = GetNumberOfRotationSteps()*delPhi/wholeCircle + 0.5;
587 if (nSphi == 0) nSphi = 1;
588 G4int nVphi = ifWholeCircle ? nSphi : nSphi + 1;
589 G4bool ifClosed = np1 <= 0; // true if external polyline is closed
590
591 // C O U N T V E R T I C E S
592
593 G4int absNp1 = std::abs(np1);
594 G4int absNp2 = std::abs(np2);
595 G4int i1beg = 0;
596 G4int i1end = absNp1-1;
597 G4int i2beg = absNp1;
598 G4int i2end = absNp1+absNp2-1;
599 G4int i, j, k;
600
601 for(i=i1beg; i<=i2end; i++) {
602 if (std::abs(r[i]) < spatialTolerance) r[i] = 0.;
603 }
604
605 // external polyline - check position of nodes relative to Z
606 //
607 G4int Nverts = 0;
608 for (i=i1beg; i<=i1end; i++) {
609 Nverts += (r[i] == 0.) ? 1 : nVphi;
610 }
611
612 // internal polyline
613 //
614 G4bool ifSide1 = false; // whether to create bottom faces
615 G4bool ifSide2 = false; // whether to create top faces
616
617 if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) { // first node
618 Nverts += (r[i2beg] == 0.) ? 1 : nVphi;
619 ifSide1 = true;
620 }
621
622 for(i=i2beg+1; i<i2end; i++) { // intermediate nodes
623 Nverts += (r[i] == 0.) ? 1 : nVphi;
624 }
625
626 if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) { // last node
627 if (absNp2 > 1) Nverts += (r[i2end] == 0.) ? 1 : nVphi;
628 ifSide2 = true;
629 }
630
631 // C O U N T F A C E S
632
633 // external lateral faces
634 //
635 G4int Nfaces = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi;
636
637 // internal lateral faces
638 //
639 if (absNp2 > 1) {
640 for(i=i2beg; i<i2end; i++) {
641 if (r[i] > 0. || r[i+1] > 0.) Nfaces += nSphi;
642 }
643
644 if (ifClosed) {
645 if (r[i2end] > 0. || r[i2beg] > 0.) Nfaces += nSphi;
646 }
647 }
648
649 // bottom and top faces
650 //
651 if (!ifClosed) {
652 if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) Nfaces += nSphi;
653 if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) Nfaces += nSphi;
654 }
655
656 // phi_wedge faces
657 //
658 if (!ifWholeCircle) {
659 Nfaces += ifClosed ? 2*absNp1 : 2*(absNp1-1);
660 }
661
662 // A L L O C A T E M E M O R Y
663
664 AllocateMemory(Nverts, Nfaces);
665 if (pV == nullptr || pF == nullptr) return;
666
667 // G E N E R A T E V E R T I C E S
668
669 G4int *kk; // array of start indices along polylines
670 kk = new G4int[absNp1+absNp2];
671
672 // external polyline
673 //
674 k = 1; // free position in array of vertices pV
675 for(i=i1beg; i<=i1end; i++) {
676 kk[i] = k;
677 if (r[i] == 0.)
678 { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
679 }
680
681 // first point of internal polyline
682 //
683 i = i2beg;
684 if (ifSide1) {
685 kk[i] = k;
686 if (r[i] == 0.)
687 { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
688 }else{
689 kk[i] = kk[i1beg];
690 }
691
692 // intermediate points of internal polyline
693 //
694 for(i=i2beg+1; i<i2end; i++) {
695 kk[i] = k;
696 if (r[i] == 0.)
697 { pV[k++] = G4Point3D(0, 0, z[i]); } else { k += nVphi; }
698 }
699
700 // last point of internal polyline
701 //
702 if (absNp2 > 1) {
703 i = i2end;
704 if (ifSide2) {
705 kk[i] = k;
706 if (r[i] == 0.) pV[k] = G4Point3D(0, 0, z[i]);
707 }else{
708 kk[i] = kk[i1end];
709 }
710 }
711
712 // set vertices
713 //
714 G4double cosPhi, sinPhi;
715
716 for(j=0; j<nVphi; j++) {
717 cosPhi = std::cos(phi+j*delPhi/nSphi);
718 sinPhi = std::sin(phi+j*delPhi/nSphi);
719 for(i=i1beg; i<=i2end; i++) {
720 if (r[i] != 0.)
721 pV[kk[i]+j] = G4Point3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
722 }
723 }
724
725 // G E N E R A T E F A C E S
726
727 // external faces
728 //
729 G4int v1,v2;
730
731 k = 1; // free position in array of faces pF
732 v2 = ifClosed ? nodeVis : 1;
733 for(i=i1beg; i<i1end; i++) {
734 v1 = v2;
735 if (!ifClosed && i == i1end-1) {
736 v2 = 1;
737 }else{
738 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
739 }
740 RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
741 edgeVis, ifWholeCircle, nSphi, k);
742 }
743 if (ifClosed) {
744 RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
745 edgeVis, ifWholeCircle, nSphi, k);
746 }
747
748 // internal faces
749 //
750 if (absNp2 > 1) {
751 v2 = ifClosed ? nodeVis : 1;
752 for(i=i2beg; i<i2end; i++) {
753 v1 = v2;
754 if (!ifClosed && i==i2end-1) {
755 v2 = 1;
756 }else{
757 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
758 }
759 RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
760 edgeVis, ifWholeCircle, nSphi, k);
761 }
762 if (ifClosed) {
763 RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
764 edgeVis, ifWholeCircle, nSphi, k);
765 }
766 }
767
768 // bottom and top faces
769 //
770 if (!ifClosed) {
771 if (ifSide1) {
772 RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
773 -1, ifWholeCircle, nSphi, k);
774 }
775 if (ifSide2) {
776 RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
777 -1, ifWholeCircle, nSphi, k);
778 }
779 }
780
781 // phi_wedge faces in case of incomplete circle
782 //
783 if (!ifWholeCircle) {
784
785 G4int ii[4], vv[4];
786
787 if (ifClosed) {
788 for (i=i1beg; i<=i1end; i++) {
789 ii[0] = i;
790 ii[3] = (i == i1end) ? i1beg : i+1;
791 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
792 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
793 vv[0] = -1;
794 vv[1] = 1;
795 vv[2] = -1;
796 vv[3] = 1;
797 SetSideFacets(ii, vv, kk, r, delPhi, nSphi, k);
798 }
799 }else{
800 for (i=i1beg; i<i1end; i++) {
801 ii[0] = i;
802 ii[3] = i+1;
803 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
804 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
805 vv[0] = (i == i1beg) ? 1 : -1;
806 vv[1] = 1;
807 vv[2] = (i == i1end-1) ? 1 : -1;
808 vv[3] = 1;
809 SetSideFacets(ii, vv, kk, r, delPhi, nSphi, k);
810 }
811 }
812 }
813
814 delete [] kk; // free memory
815
816 // final check
817 //
818 if (k-1 != nface) {
819 std::cerr
820 << "HepPolyhedron::RotateAroundZ: number of generated faces ("
821 << k-1 << ") is not equal to the number of allocated faces ("
822 << nface << ")"
823 << std::endl;
824 }
825}
826
827void
828HepPolyhedron::RotateContourAroundZ(G4int nstep,
829 G4double phi,
830 G4double dphi,
831 const std::vector<G4TwoVector> &rz,
832 G4int nodeVis,
833 G4int edgeVis)
834/***********************************************************************
835 * *
836 * Name: HepPolyhedron::RotateContourAroundZ Date: 12.05.21 *
837 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
838 * *
839 * Function: Create HepPolyhedron for a solid produced by rotation of *
840 * a closed polyline (rz-contour) around Z-axis *
841 * *
842 * Input: nstep - number of discrete steps, if 0 then default *
843 * phi - starting phi angle *
844 * dphi - delta phi *
845 * rz - rz-contour *
846 * nodeVis - how to Draw edges joing consecutive positions of *
847 * node during rotation *
848 * edgeVis - how to Draw edges *
849 * *
850 ***********************************************************************/
851{
852 // S E T R O T A T I O N P A R A M E T E R S
853
854 G4bool ifWholeCircle = std::abs(dphi - twopi) < perMillion;
855 G4double delPhi = (ifWholeCircle) ? twopi : dphi;
856 G4int nSphi = nstep;
857 if (nSphi <= 0) nSphi = GetNumberOfRotationSteps()*delPhi/twopi + 0.5;
858 if (nSphi == 0) nSphi = 1;
859 G4int nVphi = (ifWholeCircle) ? nSphi : nSphi + 1;
860
861 // C A L C U L A T E A R E A
862
863 G4int Nrz = (G4int)rz.size();
864 G4double area = 0;
865 for (G4int i = 0; i < Nrz; ++i)
866 {
867 G4int k = (i == 0) ? Nrz - 1 : i - 1;
868 area += rz[k].x()*rz[i].y() - rz[i].x()*rz[k].y();
869 }
870
871 // P R E P A R E P O L Y L I N E
872
873 auto r = new G4double[Nrz];
874 auto z = new G4double[Nrz];
875 for (G4int i = 0; i < Nrz; ++i)
876 {
877 r[i] = rz[i].x();
878 z[i] = rz[i].y();
879 if (std::abs(r[i]) < spatialTolerance) r[i] = 0.;
880 }
881
882 // C O U N T V E R T I C E S A N D F A C E S
883
884 G4int Nverts = 0;
885 for(G4int i = 0; i < Nrz; ++i) Nverts += (r[i] == 0.) ? 1 : nVphi;
886
887 G4int Nedges = Nrz;
888 for (G4int i = 0; i < Nrz; ++i)
889 {
890 G4int k = (i == 0) ? Nrz - 1 : i - 1;
891 Nedges -= static_cast<int>(r[k] == 0 && r[i] == 0);
892 }
893
894 G4int Nfaces = Nedges*nSphi; // lateral faces
895 if (!ifWholeCircle) Nfaces += 2*(Nrz - 2); // phi_wedge faces
896
897 // A L L O C A T E M E M O R Y
898
899 AllocateMemory(Nverts, Nfaces);
900 if (pV == nullptr || pF == nullptr)
901 {
902 delete [] r;
903 delete [] z;
904 return;
905 }
906
907 // S E T V E R T I C E S
908
909 auto kk = new G4int[Nrz]; // start indices along contour
910 G4int kfree = 1; // current free position in array of vertices pV
911
912 // set start indices, set vertices for nodes with r == 0
913 for(G4int i = 0; i < Nrz; ++i)
914 {
915 kk[i] = kfree;
916 if (r[i] == 0.) pV[kfree++] = G4Point3D(0, 0, z[i]);
917 if (r[i] != 0.) kfree += nVphi;
918 }
919
920 // set vertices by rotating r
921 for(G4int j = 0; j < nVphi; ++j)
922 {
923 G4double cosPhi = std::cos(phi + j*delPhi/nSphi);
924 G4double sinPhi = std::sin(phi + j*delPhi/nSphi);
925 for(G4int i = 0; i < Nrz; ++i)
926 {
927 if (r[i] != 0.)
928 pV[kk[i] + j] = G4Point3D(r[i]*cosPhi, r[i]*sinPhi, z[i]);
929 }
930 }
931
932 // S E T F A C E S
933
934 kfree = 1; // current free position in array of faces pF
935 for(G4int i = 0; i < Nrz; ++i)
936 {
937 G4int i1 = (i < Nrz - 1) ? i + 1 : 0; // inverse order if area > 0
938 G4int i2 = i;
939 if (area < 0.) std::swap(i1, i2);
940 RotateEdge(kk[i1], kk[i2], r[i1], r[i2], nodeVis, nodeVis,
941 edgeVis, ifWholeCircle, nSphi, kfree);
942 }
943
944 // S E T P H I _ W E D G E F A C E S
945
946 if (!ifWholeCircle)
947 {
948 std::vector<G4int> triangles;
949 TriangulatePolygon(rz, triangles);
950
951 G4int ii[4], vv[4];
952 G4int ntria = G4int(triangles.size()/3);
953 for (G4int i = 0; i < ntria; ++i)
954 {
955 G4int i1 = triangles[0 + i*3];
956 G4int i2 = triangles[1 + i*3];
957 G4int i3 = triangles[2 + i*3];
958 if (area < 0.) std::swap(i1, i3);
959 G4int v1 = (std::abs(i2-i1) == 1 || std::abs(i2-i1) == Nrz-1) ? 1 : -1;
960 G4int v2 = (std::abs(i3-i2) == 1 || std::abs(i3-i2) == Nrz-1) ? 1 : -1;
961 G4int v3 = (std::abs(i1-i3) == 1 || std::abs(i1-i3) == Nrz-1) ? 1 : -1;
962 ii[0] = i1; ii[1] = i2; ii[2] = i2; ii[3] = i3;
963 vv[0] = v1; vv[1] = -1; vv[2] = v2; vv[3] = v3;
964 SetSideFacets(ii, vv, kk, r, delPhi, nSphi, kfree);
965 }
966 }
967
968 // free memory
969 delete [] r;
970 delete [] z;
971 delete [] kk;
972
973 // final check
974 if (kfree - 1 != nface)
975 {
976 std::cerr
977 << "HepPolyhedron::RotateContourAroundZ: number of generated faces ("
978 << kfree-1 << ") is not equal to the number of allocated faces ("
979 << nface << ")"
980 << std::endl;
981 }
982}
983
984G4bool
985HepPolyhedron::TriangulatePolygon(const std::vector<G4TwoVector> &polygon,
986 std::vector<G4int> &result)
987/***********************************************************************
988 * *
989 * Name: HepPolyhedron::TriangulatePolygon Date: 12.05.21 *
990 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
991 * *
992 * Function: Simple implementation of "ear clipping" algorithm for *
993 * triangulation of a simple contour/polygon, it places *
994 * the result in a std::vector as triplets of vertex indices *
995 * *
996 * If triangulation is sucsessfull then the function *
997 * returns true, otherwise false *
998 * *
999 * Remark: It's a copy of G4GeomTools::TriangulatePolygon() *
1000 * *
1001 ***********************************************************************/
1002{
1003 result.resize(0);
1004 G4int n = (G4int)polygon.size();
1005 if (n < 3) return false;
1006
1007 // calculate area
1008 //
1009 G4double area = 0.;
1010 for(G4int i = 0; i < n; ++i)
1011 {
1012 G4int k = (i == 0) ? n - 1 : i - 1;
1013 area += polygon[k].x()*polygon[i].y() - polygon[i].x()*polygon[k].y();
1014 }
1015
1016 // allocate and initialize list of Vertices
1017 // we want a counter-clockwise polygon in V
1018 //
1019 auto V = new G4int[n];
1020 if (area > 0.)
1021 for (G4int i = 0; i < n; ++i) V[i] = i;
1022 else
1023 for (G4int i = 0; i < n; ++i) V[i] = (n - 1) - i;
1024
1025 // Triangulation: remove nv-2 Vertices, creating 1 triangle every time
1026 //
1027 G4int nv = n;
1028 G4int count = 2*nv; // error detection counter
1029 for(G4int b = nv - 1; nv > 2; )
1030 {
1031 // ERROR: if we loop, it is probably a non-simple polygon
1032 if ((count--) <= 0)
1033 {
1034 delete [] V;
1035 if (area < 0.) std::reverse(result.begin(),result.end());
1036 return false;
1037 }
1038
1039 // three consecutive vertices in current polygon, <a,b,c>
1040 G4int a = (b < nv) ? b : 0; // previous
1041 b = (a+1 < nv) ? a+1 : 0; // current
1042 G4int c = (b+1 < nv) ? b+1 : 0; // next
1043
1044 if (CheckSnip(polygon, a,b,c, nv,V))
1045 {
1046 // output Triangle
1047 result.push_back(V[a]);
1048 result.push_back(V[b]);
1049 result.push_back(V[c]);
1050
1051 // remove vertex b from remaining polygon
1052 nv--;
1053 for(G4int i = b; i < nv; ++i) V[i] = V[i+1];
1054
1055 count = 2*nv; // resest error detection counter
1056 }
1057 }
1058 delete [] V;
1059 if (area < 0.) std::reverse(result.begin(),result.end());
1060 return true;
1061}
1062
1063G4bool HepPolyhedron::CheckSnip(const std::vector<G4TwoVector> &contour,
1064 G4int a, G4int b, G4int c,
1065 G4int n, const G4int* V)
1066/***********************************************************************
1067 * *
1068 * Name: HepPolyhedron::CheckSnip Date: 12.05.21 *
1069 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
1070 * *
1071 * Function: Check for a valid snip, *
1072 * it is a helper functionfor TriangulatePolygon() *
1073 * *
1074 ***********************************************************************/
1075{
1076 static const G4double kCarTolerance = 1.e-9;
1077
1078 // check orientation of Triangle
1079 G4double Ax = contour[V[a]].x(), Ay = contour[V[a]].y();
1080 G4double Bx = contour[V[b]].x(), By = contour[V[b]].y();
1081 G4double Cx = contour[V[c]].x(), Cy = contour[V[c]].y();
1082 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) < kCarTolerance) return false;
1083
1084 // check that there is no point inside Triangle
1085 G4double xmin = std::min(std::min(Ax,Bx),Cx);
1086 G4double xmax = std::max(std::max(Ax,Bx),Cx);
1087 G4double ymin = std::min(std::min(Ay,By),Cy);
1088 G4double ymax = std::max(std::max(Ay,By),Cy);
1089
1090 for (G4int i=0; i<n; ++i)
1091 {
1092 if((i == a) || (i == b) || (i == c)) continue;
1093 G4double Px = contour[V[i]].x();
1094 if (Px < xmin || Px > xmax) continue;
1095 G4double Py = contour[V[i]].y();
1096 if (Py < ymin || Py > ymax) continue;
1097 // if (PointInTriangle(Ax,Ay,Bx,By,Cx,Cy,Px,Py)) return false;
1098 if ((Bx-Ax)*(Cy-Ay) - (By-Ay)*(Cx-Ax) > 0.)
1099 {
1100 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) < 0.) continue;
1101 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) < 0.) continue;
1102 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) < 0.) continue;
1103 }
1104 else
1105 {
1106 if ((Ax-Cx)*(Py-Cy) - (Ay-Cy)*(Px-Cx) > 0.) continue;
1107 if ((Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax) > 0.) continue;
1108 if ((Cx-Bx)*(Py-By) - (Cy-By)*(Px-Bx) > 0.) continue;
1109 }
1110 return false;
1111 }
1112 return true;
1113}
1114
1115void HepPolyhedron::SetReferences()
1116/***********************************************************************
1117 * *
1118 * Name: HepPolyhedron::SetReferences Date: 04.12.96 *
1119 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1120 * *
1121 * Function: For each edge set reference to neighbouring facet *
1122 * *
1123 ***********************************************************************/
1124{
1125 if (nface <= 0) return;
1126
1127 struct edgeListMember {
1128 edgeListMember *next;
1129 G4int v2;
1130 G4int iface;
1131 G4int iedge;
1132 } *edgeList, *freeList, **headList;
1133
1134
1135 // A L L O C A T E A N D I N I T I A T E L I S T S
1136
1137 edgeList = new edgeListMember[2*nface];
1138 headList = new edgeListMember*[nvert];
1139
1140 G4int i;
1141 for (i=0; i<nvert; i++) {
1142 headList[i] = nullptr;
1143 }
1144 freeList = edgeList;
1145 for (i=0; i<2*nface-1; i++) {
1146 edgeList[i].next = &edgeList[i+1];
1147 }
1148 edgeList[2*nface-1].next = nullptr;
1149
1150 // L O O P A L O N G E D G E S
1151
1152 G4int iface, iedge, nedge, i1, i2, k1, k2;
1153 edgeListMember *prev, *cur;
1154
1155 for(iface=1; iface<=nface; iface++) {
1156 nedge = (pF[iface].edge[3].v == 0) ? 3 : 4;
1157 for (iedge=0; iedge<nedge; iedge++) {
1158 i1 = iedge;
1159 i2 = (iedge < nedge-1) ? iedge+1 : 0;
1160 i1 = std::abs(pF[iface].edge[i1].v);
1161 i2 = std::abs(pF[iface].edge[i2].v);
1162 k1 = (i1 < i2) ? i1 : i2; // k1 = ::min(i1,i2);
1163 k2 = (i1 > i2) ? i1 : i2; // k2 = ::max(i1,i2);
1164
1165 // check head of the List corresponding to k1
1166 cur = headList[k1];
1167 if (cur == nullptr) {
1168 headList[k1] = freeList;
1169 if (freeList == nullptr) {
1170 std::cerr
1171 << "Polyhedron::SetReferences: bad link "
1172 << std::endl;
1173 break;
1174 }
1175 freeList = freeList->next;
1176 cur = headList[k1];
1177 cur->next = nullptr;
1178 cur->v2 = k2;
1179 cur->iface = iface;
1180 cur->iedge = iedge;
1181 continue;
1182 }
1183
1184 if (cur->v2 == k2) {
1185 headList[k1] = cur->next;
1186 cur->next = freeList;
1187 freeList = cur;
1188 pF[iface].edge[iedge].f = cur->iface;
1189 pF[cur->iface].edge[cur->iedge].f = iface;
1190 i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
1191 i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
1192 if (i1 != i2) {
1193 std::cerr
1194 << "Polyhedron::SetReferences: different edge visibility "
1195 << iface << "/" << iedge << "/"
1196 << pF[iface].edge[iedge].v << " and "
1197 << cur->iface << "/" << cur->iedge << "/"
1198 << pF[cur->iface].edge[cur->iedge].v
1199 << std::endl;
1200 }
1201 continue;
1202 }
1203
1204 // check List itself
1205 for (;;) {
1206 prev = cur;
1207 cur = prev->next;
1208 if (cur == nullptr) {
1209 prev->next = freeList;
1210 if (freeList == nullptr) {
1211 std::cerr
1212 << "Polyhedron::SetReferences: bad link "
1213 << std::endl;
1214 break;
1215 }
1216 freeList = freeList->next;
1217 cur = prev->next;
1218 cur->next = nullptr;
1219 cur->v2 = k2;
1220 cur->iface = iface;
1221 cur->iedge = iedge;
1222 break;
1223 }
1224
1225 if (cur->v2 == k2) {
1226 prev->next = cur->next;
1227 cur->next = freeList;
1228 freeList = cur;
1229 pF[iface].edge[iedge].f = cur->iface;
1230 pF[cur->iface].edge[cur->iedge].f = iface;
1231 i1 = (pF[iface].edge[iedge].v < 0) ? -1 : 1;
1232 i2 = (pF[cur->iface].edge[cur->iedge].v < 0) ? -1 : 1;
1233 if (i1 != i2) {
1234 std::cerr
1235 << "Polyhedron::SetReferences: different edge visibility "
1236 << iface << "/" << iedge << "/"
1237 << pF[iface].edge[iedge].v << " and "
1238 << cur->iface << "/" << cur->iedge << "/"
1239 << pF[cur->iface].edge[cur->iedge].v
1240 << std::endl;
1241 }
1242 break;
1243 }
1244 }
1245 }
1246 }
1247
1248 // C H E C K T H A T A L L L I S T S A R E E M P T Y
1249
1250 for (i=0; i<nvert; i++) {
1251 if (headList[i] != nullptr) {
1252 std::cerr
1253 << "Polyhedron::SetReferences: List " << i << " is not empty"
1254 << std::endl;
1255 }
1256 }
1257
1258 // F R E E M E M O R Y
1259
1260 delete [] edgeList;
1261 delete [] headList;
1262}
1263
1264void HepPolyhedron::JoinCoplanarFacets(G4double tolerance)
1265/***********************************************************************
1266 * *
1267 * Name: HepPolyhedron::JoinCoplanarFacets Date: 04.04.22 *
1268 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
1269 * *
1270 * Function: Join couples of triangular facets to quadrangular facets *
1271 * where it is possible *
1272 * *
1273 ***********************************************************************/
1274{
1275 G4int njoin = 0;
1276 for (G4int icur = 1; icur <= nface; ++icur)
1277 {
1278 // skip if already joined or quadrangle
1279 if (pF[icur].edge[0].v == 0) continue;
1280 if (pF[icur].edge[3].v != 0) continue;
1281 // skip if all references point to already checked facets
1282 if (pF[icur].edge[0].f < icur &&
1283 pF[icur].edge[1].f < icur &&
1284 pF[icur].edge[2].f < icur) continue;
1285 // compute plane equation
1286 G4Normal3D norm = GetUnitNormal(icur);
1287 G4double dd = norm.dot(pV[pF[icur].edge[0].v]);
1288 G4int vcur0 = std::abs(pF[icur].edge[0].v);
1289 G4int vcur1 = std::abs(pF[icur].edge[1].v);
1290 G4int vcur2 = std::abs(pF[icur].edge[2].v);
1291 // select neighbouring facet
1292 G4int kcheck = 0, icheck = 0, vcheck = 0;
1293 G4double dist = DBL_MAX;
1294 for (G4int k = 0; k < 3; ++k)
1295 {
1296 G4int itmp = pF[icur].edge[k].f;
1297 // skip if already checked, joined or quadrangle
1298 if (itmp < icur) continue;
1299 if (pF[itmp].edge[0].v == 0 ||
1300 pF[itmp].edge[3].v != 0) continue;
1301 // get candidate vertex
1302 G4int vtmp = 0;
1303 for (G4int j = 0; j < 3; ++j)
1304 {
1305 vtmp = std::abs(pF[itmp].edge[j].v);
1306 if (vtmp != vcur0 && vtmp != vcur1 && vtmp != vcur2) break;
1307 }
1308 // check distance to the plane
1309 G4double dtmp = std::abs(norm.dot(pV[vtmp]) - dd);
1310 if (dtmp > tolerance || dtmp >= dist) continue;
1311 dist = dtmp;
1312 kcheck = k;
1313 icheck = itmp;
1314 vcheck = vtmp;
1315 }
1316 if (icheck == 0) continue; // no facet selected
1317 // join facets
1318 njoin++;
1319 pF[icheck].edge[0].v = 0; // mark facet as joined
1320 if (kcheck == 0)
1321 {
1322 pF[icur].edge[3].v = pF[icur].edge[2].v;
1323 pF[icur].edge[2].v = pF[icur].edge[1].v;
1324 pF[icur].edge[1].v = vcheck;
1325 }
1326 else if (kcheck == 1)
1327 {
1328 pF[icur].edge[3].v = pF[icur].edge[2].v;
1329 pF[icur].edge[2].v = vcheck;
1330 }
1331 else
1332 {
1333 pF[icur].edge[3].v = vcheck;
1334 }
1335 }
1336 if (njoin == 0) return; // no joined facets
1337
1338 // restructure facets
1339 G4int nnew = 0;
1340 for (G4int icur = 1; icur <= nface; ++icur)
1341 {
1342 if (pF[icur].edge[0].v == 0) continue;
1343 nnew++;
1344 pF[nnew].edge[0].v = pF[icur].edge[0].v;
1345 pF[nnew].edge[1].v = pF[icur].edge[1].v;
1346 pF[nnew].edge[2].v = pF[icur].edge[2].v;
1347 pF[nnew].edge[3].v = pF[icur].edge[3].v;
1348 }
1349 nface = nnew;
1350 SetReferences();
1351}
1352
1353void HepPolyhedron::InvertFacets()
1354/***********************************************************************
1355 * *
1356 * Name: HepPolyhedron::InvertFacets Date: 01.12.99 *
1357 * Author: E.Chernyaev Revised: *
1358 * *
1359 * Function: Invert the order of the nodes in the facets *
1360 * *
1361 ***********************************************************************/
1362{
1363 if (nface <= 0) return;
1364 G4int i, k, nnode, v[4],f[4];
1365 for (i=1; i<=nface; i++) {
1366 nnode = (pF[i].edge[3].v == 0) ? 3 : 4;
1367 for (k=0; k<nnode; k++) {
1368 v[k] = (k+1 == nnode) ? pF[i].edge[0].v : pF[i].edge[k+1].v;
1369 if (v[k] * pF[i].edge[k].v < 0) v[k] = -v[k];
1370 f[k] = pF[i].edge[k].f;
1371 }
1372 for (k=0; k<nnode; k++) {
1373 pF[i].edge[nnode-1-k].v = v[k];
1374 pF[i].edge[nnode-1-k].f = f[k];
1375 }
1376 }
1377}
1378
1379HepPolyhedron & HepPolyhedron::Transform(const G4Transform3D &t)
1380/***********************************************************************
1381 * *
1382 * Name: HepPolyhedron::Transform Date: 01.12.99 *
1383 * Author: E.Chernyaev Revised: *
1384 * *
1385 * Function: Make transformation of the polyhedron *
1386 * *
1387 ***********************************************************************/
1388{
1389 if (nvert > 0) {
1390 for (G4int i=1; i<=nvert; i++) { pV[i] = t * pV[i]; }
1391
1392 // C H E C K D E T E R M I N A N T A N D
1393 // I N V E R T F A C E T S I F I T I S N E G A T I V E
1394
1395 G4Vector3D d = t * G4Vector3D(0,0,0);
1396 G4Vector3D x = t * G4Vector3D(1,0,0) - d;
1397 G4Vector3D y = t * G4Vector3D(0,1,0) - d;
1398 G4Vector3D z = t * G4Vector3D(0,0,1) - d;
1399 if ((x.cross(y))*z < 0) InvertFacets();
1400 }
1401 return *this;
1402}
1403
1404G4bool HepPolyhedron::GetNextVertexIndex(G4int &index, G4int &edgeFlag) const
1405/***********************************************************************
1406 * *
1407 * Name: HepPolyhedron::GetNextVertexIndex Date: 03.09.96 *
1408 * Author: Yasuhide Sawada Revised: *
1409 * *
1410 * Function: *
1411 * *
1412 ***********************************************************************/
1413{
1414 static G4ThreadLocal G4int iFace = 1;
1415 static G4ThreadLocal G4int iQVertex = 0;
1416 G4int vIndex = pF[iFace].edge[iQVertex].v;
1417
1418 edgeFlag = (vIndex > 0) ? 1 : 0;
1419 index = std::abs(vIndex);
1420
1421 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1422 iQVertex = 0;
1423 if (++iFace > nface) iFace = 1;
1424 return false; // Last Edge
1425 }
1426
1427 ++iQVertex;
1428 return true; // not Last Edge
1429}
1430
1431G4Point3D HepPolyhedron::GetVertex(G4int index) const
1432/***********************************************************************
1433 * *
1434 * Name: HepPolyhedron::GetVertex Date: 03.09.96 *
1435 * Author: Yasuhide Sawada Revised: 17.11.99 *
1436 * *
1437 * Function: Get vertex of the index. *
1438 * *
1439 ***********************************************************************/
1440{
1441 if (index <= 0 || index > nvert) {
1442 std::cerr
1443 << "HepPolyhedron::GetVertex: irrelevant index " << index
1444 << std::endl;
1445 return G4Point3D();
1446 }
1447 return pV[index];
1448}
1449
1450G4bool
1451HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag) const
1452/***********************************************************************
1453 * *
1454 * Name: HepPolyhedron::GetNextVertex Date: 22.07.96 *
1455 * Author: John Allison Revised: *
1456 * *
1457 * Function: Get vertices of the quadrilaterals in order for each *
1458 * face in face order. Returns false when finished each *
1459 * face. *
1460 * *
1461 ***********************************************************************/
1462{
1463 G4int index;
1464 G4bool rep = GetNextVertexIndex(index, edgeFlag);
1465 vertex = pV[index];
1466 return rep;
1467}
1468
1469G4bool HepPolyhedron::GetNextVertex(G4Point3D &vertex, G4int &edgeFlag,
1470 G4Normal3D &normal) const
1471/***********************************************************************
1472 * *
1473 * Name: HepPolyhedron::GetNextVertex Date: 26.11.99 *
1474 * Author: E.Chernyaev Revised: *
1475 * *
1476 * Function: Get vertices with normals of the quadrilaterals in order *
1477 * for each face in face order. *
1478 * Returns false when finished each face. *
1479 * *
1480 ***********************************************************************/
1481{
1482 static G4ThreadLocal G4int iFace = 1;
1483 static G4ThreadLocal G4int iNode = 0;
1484
1485 if (nface == 0) return false; // empty polyhedron
1486
1487 G4int k = pF[iFace].edge[iNode].v;
1488 if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
1489 vertex = pV[k];
1490 normal = FindNodeNormal(iFace,k);
1491 if (iNode >= 3 || pF[iFace].edge[iNode+1].v == 0) {
1492 iNode = 0;
1493 if (++iFace > nface) iFace = 1;
1494 return false; // last node
1495 }
1496 ++iNode;
1497 return true; // not last node
1498}
1499
1500G4bool HepPolyhedron::GetNextEdgeIndices(G4int &i1, G4int &i2, G4int &edgeFlag,
1501 G4int &iface1, G4int &iface2) const
1502/***********************************************************************
1503 * *
1504 * Name: HepPolyhedron::GetNextEdgeIndices Date: 30.09.96 *
1505 * Author: E.Chernyaev Revised: 17.11.99 *
1506 * *
1507 * Function: Get indices of the next edge together with indices of *
1508 * of the faces which share the edge. *
1509 * Returns false when the last edge. *
1510 * *
1511 ***********************************************************************/
1512{
1513 static G4ThreadLocal G4int iFace = 1;
1514 static G4ThreadLocal G4int iQVertex = 0;
1515 static G4ThreadLocal G4int iOrder = 1;
1516 G4int k1, k2, kflag, kface1, kface2;
1517
1518 if (iFace == 1 && iQVertex == 0) {
1519 k2 = pF[nface].edge[0].v;
1520 k1 = pF[nface].edge[3].v;
1521 if (k1 == 0) k1 = pF[nface].edge[2].v;
1522 if (std::abs(k1) > std::abs(k2)) iOrder = -1;
1523 }
1524
1525 do {
1526 k1 = pF[iFace].edge[iQVertex].v;
1527 kflag = k1;
1528 k1 = std::abs(k1);
1529 kface1 = iFace;
1530 kface2 = pF[iFace].edge[iQVertex].f;
1531 if (iQVertex >= 3 || pF[iFace].edge[iQVertex+1].v == 0) {
1532 iQVertex = 0;
1533 k2 = std::abs(pF[iFace].edge[iQVertex].v);
1534 iFace++;
1535 }else{
1536 iQVertex++;
1537 k2 = std::abs(pF[iFace].edge[iQVertex].v);
1538 }
1539 } while (iOrder*k1 > iOrder*k2);
1540
1541 i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1542 iface1 = kface1; iface2 = kface2;
1543
1544 if (iFace > nface) {
1545 iFace = 1; iOrder = 1;
1546 return false;
1547 }
1548
1549 return true;
1550}
1551
1552G4bool
1553HepPolyhedron::GetNextEdgeIndices(G4int &i1, G4int &i2, G4int &edgeFlag) const
1554/***********************************************************************
1555 * *
1556 * Name: HepPolyhedron::GetNextEdgeIndices Date: 17.11.99 *
1557 * Author: E.Chernyaev Revised: *
1558 * *
1559 * Function: Get indices of the next edge. *
1560 * Returns false when the last edge. *
1561 * *
1562 ***********************************************************************/
1563{
1564 G4int kface1, kface2;
1565 return GetNextEdgeIndices(i1, i2, edgeFlag, kface1, kface2);
1566}
1567
1568G4bool
1569HepPolyhedron::GetNextEdge(G4Point3D &p1,
1570 G4Point3D &p2,
1571 G4int &edgeFlag) const
1572/***********************************************************************
1573 * *
1574 * Name: HepPolyhedron::GetNextEdge Date: 30.09.96 *
1575 * Author: E.Chernyaev Revised: *
1576 * *
1577 * Function: Get next edge. *
1578 * Returns false when the last edge. *
1579 * *
1580 ***********************************************************************/
1581{
1582 G4int i1,i2;
1583 G4bool rep = GetNextEdgeIndices(i1,i2,edgeFlag);
1584 p1 = pV[i1];
1585 p2 = pV[i2];
1586 return rep;
1587}
1588
1589G4bool
1590HepPolyhedron::GetNextEdge(G4Point3D &p1, G4Point3D &p2,
1591 G4int &edgeFlag, G4int &iface1, G4int &iface2) const
1592/***********************************************************************
1593 * *
1594 * Name: HepPolyhedron::GetNextEdge Date: 17.11.99 *
1595 * Author: E.Chernyaev Revised: *
1596 * *
1597 * Function: Get next edge with indices of the faces which share *
1598 * the edge. *
1599 * Returns false when the last edge. *
1600 * *
1601 ***********************************************************************/
1602{
1603 G4int i1,i2;
1604 G4bool rep = GetNextEdgeIndices(i1,i2,edgeFlag,iface1,iface2);
1605 p1 = pV[i1];
1606 p2 = pV[i2];
1607 return rep;
1608}
1609
1610void HepPolyhedron::GetFacet(G4int iFace, G4int &n, G4int *iNodes,
1611 G4int *edgeFlags, G4int *iFaces) const
1612/***********************************************************************
1613 * *
1614 * Name: HepPolyhedron::GetFacet Date: 15.12.99 *
1615 * Author: E.Chernyaev Revised: *
1616 * *
1617 * Function: Get face by index *
1618 * *
1619 ***********************************************************************/
1620{
1621 if (iFace < 1 || iFace > nface) {
1622 std::cerr
1623 << "HepPolyhedron::GetFacet: irrelevant index " << iFace
1624 << std::endl;
1625 n = 0;
1626 }else{
1627 G4int i, k;
1628 for (i=0; i<4; i++) {
1629 k = pF[iFace].edge[i].v;
1630 if (k == 0) break;
1631 if (iFaces != nullptr) iFaces[i] = pF[iFace].edge[i].f;
1632 if (k > 0) {
1633 iNodes[i] = k;
1634 if (edgeFlags != nullptr) edgeFlags[i] = 1;
1635 }else{
1636 iNodes[i] = -k;
1637 if (edgeFlags != nullptr) edgeFlags[i] = -1;
1638 }
1639 }
1640 n = i;
1641 }
1642}
1643
1644void HepPolyhedron::GetFacet(G4int index, G4int &n, G4Point3D *nodes,
1645 G4int *edgeFlags, G4Normal3D *normals) const
1646/***********************************************************************
1647 * *
1648 * Name: HepPolyhedron::GetFacet Date: 17.11.99 *
1649 * Author: E.Chernyaev Revised: *
1650 * *
1651 * Function: Get face by index *
1652 * *
1653 ***********************************************************************/
1654{
1655 G4int iNodes[4];
1656 GetFacet(index, n, iNodes, edgeFlags);
1657 if (n != 0) {
1658 for (G4int i=0; i<n; i++) {
1659 nodes[i] = pV[iNodes[i]];
1660 if (normals != nullptr) normals[i] = FindNodeNormal(index,iNodes[i]);
1661 }
1662 }
1663}
1664
1665G4bool
1666HepPolyhedron::GetNextFacet(G4int &n, G4Point3D *nodes,
1667 G4int *edgeFlags, G4Normal3D *normals) const
1668/***********************************************************************
1669 * *
1670 * Name: HepPolyhedron::GetNextFacet Date: 19.11.99 *
1671 * Author: E.Chernyaev Revised: *
1672 * *
1673 * Function: Get next face with normals of unit length at the nodes. *
1674 * Returns false when finished all faces. *
1675 * *
1676 ***********************************************************************/
1677{
1678 static G4ThreadLocal G4int iFace = 1;
1679
1680 if (edgeFlags == nullptr) {
1681 GetFacet(iFace, n, nodes);
1682 }else if (normals == nullptr) {
1683 GetFacet(iFace, n, nodes, edgeFlags);
1684 }else{
1685 GetFacet(iFace, n, nodes, edgeFlags, normals);
1686 }
1687
1688 if (++iFace > nface) {
1689 iFace = 1;
1690 return false;
1691 }
1692
1693 return true;
1694}
1695
1696G4Normal3D HepPolyhedron::GetNormal(G4int iFace) const
1697/***********************************************************************
1698 * *
1699 * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1700 * Author: E.Chernyaev Revised: *
1701 * *
1702 * Function: Get normal of the face given by index *
1703 * *
1704 ***********************************************************************/
1705{
1706 if (iFace < 1 || iFace > nface) {
1707 std::cerr
1708 << "HepPolyhedron::GetNormal: irrelevant index " << iFace
1709 << std::endl;
1710 return G4Normal3D();
1711 }
1712
1713 G4int i0 = std::abs(pF[iFace].edge[0].v);
1714 G4int i1 = std::abs(pF[iFace].edge[1].v);
1715 G4int i2 = std::abs(pF[iFace].edge[2].v);
1716 G4int i3 = std::abs(pF[iFace].edge[3].v);
1717 if (i3 == 0) i3 = i0;
1718 return (pV[i2] - pV[i0]).cross(pV[i3] - pV[i1]);
1719}
1720
1721G4Normal3D HepPolyhedron::GetUnitNormal(G4int iFace) const
1722/***********************************************************************
1723 * *
1724 * Name: HepPolyhedron::GetNormal Date: 19.11.99 *
1725 * Author: E.Chernyaev Revised: *
1726 * *
1727 * Function: Get unit normal of the face given by index *
1728 * *
1729 ***********************************************************************/
1730{
1731 if (iFace < 1 || iFace > nface) {
1732 std::cerr
1733 << "HepPolyhedron::GetUnitNormal: irrelevant index " << iFace
1734 << std::endl;
1735 return G4Normal3D();
1736 }
1737
1738 G4int i0 = std::abs(pF[iFace].edge[0].v);
1739 G4int i1 = std::abs(pF[iFace].edge[1].v);
1740 G4int i2 = std::abs(pF[iFace].edge[2].v);
1741 G4int i3 = std::abs(pF[iFace].edge[3].v);
1742 if (i3 == 0) i3 = i0;
1743 return ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).unit();
1744}
1745
1746G4bool HepPolyhedron::GetNextNormal(G4Normal3D &normal) const
1747/***********************************************************************
1748 * *
1749 * Name: HepPolyhedron::GetNextNormal Date: 22.07.96 *
1750 * Author: John Allison Revised: 19.11.99 *
1751 * *
1752 * Function: Get normals of each face in face order. Returns false *
1753 * when finished all faces. *
1754 * *
1755 ***********************************************************************/
1756{
1757 static G4ThreadLocal G4int iFace = 1;
1758 normal = GetNormal(iFace);
1759 if (++iFace > nface) {
1760 iFace = 1;
1761 return false;
1762 }
1763 return true;
1764}
1765
1766G4bool HepPolyhedron::GetNextUnitNormal(G4Normal3D &normal) const
1767/***********************************************************************
1768 * *
1769 * Name: HepPolyhedron::GetNextUnitNormal Date: 16.09.96 *
1770 * Author: E.Chernyaev Revised: *
1771 * *
1772 * Function: Get normals of unit length of each face in face order. *
1773 * Returns false when finished all faces. *
1774 * *
1775 ***********************************************************************/
1776{
1777 G4bool rep = GetNextNormal(normal);
1778 normal = normal.unit();
1779 return rep;
1780}
1781
1782G4double HepPolyhedron::GetSurfaceArea() const
1783/***********************************************************************
1784 * *
1785 * Name: HepPolyhedron::GetSurfaceArea Date: 25.05.01 *
1786 * Author: E.Chernyaev Revised: *
1787 * *
1788 * Function: Returns area of the surface of the polyhedron. *
1789 * *
1790 ***********************************************************************/
1791{
1792 G4double srf = 0.;
1793 for (G4int iFace=1; iFace<=nface; iFace++) {
1794 G4int i0 = std::abs(pF[iFace].edge[0].v);
1795 G4int i1 = std::abs(pF[iFace].edge[1].v);
1796 G4int i2 = std::abs(pF[iFace].edge[2].v);
1797 G4int i3 = std::abs(pF[iFace].edge[3].v);
1798 if (i3 == 0) i3 = i0;
1799 srf += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).mag();
1800 }
1801 return srf/2.;
1802}
1803
1804G4double HepPolyhedron::GetVolume() const
1805/***********************************************************************
1806 * *
1807 * Name: HepPolyhedron::GetVolume Date: 25.05.01 *
1808 * Author: E.Chernyaev Revised: *
1809 * *
1810 * Function: Returns volume of the polyhedron. *
1811 * *
1812 ***********************************************************************/
1813{
1814 G4double v = 0.;
1815 for (G4int iFace=1; iFace<=nface; iFace++) {
1816 G4int i0 = std::abs(pF[iFace].edge[0].v);
1817 G4int i1 = std::abs(pF[iFace].edge[1].v);
1818 G4int i2 = std::abs(pF[iFace].edge[2].v);
1819 G4int i3 = std::abs(pF[iFace].edge[3].v);
1820 G4Point3D pt;
1821 if (i3 == 0) {
1822 i3 = i0;
1823 pt = (pV[i0]+pV[i1]+pV[i2]) * (1./3.);
1824 }else{
1825 pt = (pV[i0]+pV[i1]+pV[i2]+pV[i3]) * 0.25;
1826 }
1827 v += ((pV[i2] - pV[i0]).cross(pV[i3] - pV[i1])).dot(pt);
1828 }
1829 return v/6.;
1830}
1831
1832G4int
1833HepPolyhedron::createTwistedTrap(G4double Dz,
1834 const G4double xy1[][2],
1835 const G4double xy2[][2])
1836/***********************************************************************
1837 * *
1838 * Name: createTwistedTrap Date: 05.11.02 *
1839 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1840 * *
1841 * Function: Creates polyhedron for twisted trapezoid *
1842 * *
1843 * Input: Dz - half-length along Z 8----7 *
1844 * xy1[2,4] - quadrilateral at Z=-Dz 5----6 ! *
1845 * xy2[2,4] - quadrilateral at Z=+Dz ! 4-!--3 *
1846 * 1----2 *
1847 * *
1848 ***********************************************************************/
1849{
1850 AllocateMemory(12,18);
1851
1852 pV[ 1] = G4Point3D(xy1[0][0],xy1[0][1],-Dz);
1853 pV[ 2] = G4Point3D(xy1[1][0],xy1[1][1],-Dz);
1854 pV[ 3] = G4Point3D(xy1[2][0],xy1[2][1],-Dz);
1855 pV[ 4] = G4Point3D(xy1[3][0],xy1[3][1],-Dz);
1856
1857 pV[ 5] = G4Point3D(xy2[0][0],xy2[0][1], Dz);
1858 pV[ 6] = G4Point3D(xy2[1][0],xy2[1][1], Dz);
1859 pV[ 7] = G4Point3D(xy2[2][0],xy2[2][1], Dz);
1860 pV[ 8] = G4Point3D(xy2[3][0],xy2[3][1], Dz);
1861
1862 pV[ 9] = (pV[1]+pV[2]+pV[5]+pV[6])/4.;
1863 pV[10] = (pV[2]+pV[3]+pV[6]+pV[7])/4.;
1864 pV[11] = (pV[3]+pV[4]+pV[7]+pV[8])/4.;
1865 pV[12] = (pV[4]+pV[1]+pV[8]+pV[5])/4.;
1866
1867 enum {DUMMY, BOTTOM,
1868 LEFT_BOTTOM, LEFT_FRONT, LEFT_TOP, LEFT_BACK,
1869 BACK_BOTTOM, BACK_LEFT, BACK_TOP, BACK_RIGHT,
1870 RIGHT_BOTTOM, RIGHT_BACK, RIGHT_TOP, RIGHT_FRONT,
1871 FRONT_BOTTOM, FRONT_RIGHT, FRONT_TOP, FRONT_LEFT,
1872 TOP};
1873
1874 pF[ 1]=G4Facet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
1875
1876 pF[ 2]=G4Facet(4,BOTTOM, -1,LEFT_FRONT, -12,LEFT_BACK, 0,0);
1877 pF[ 3]=G4Facet(1,FRONT_LEFT, -5,LEFT_TOP, -12,LEFT_BOTTOM, 0,0);
1878 pF[ 4]=G4Facet(5,TOP, -8,LEFT_BACK, -12,LEFT_FRONT, 0,0);
1879 pF[ 5]=G4Facet(8,BACK_LEFT, -4,LEFT_BOTTOM, -12,LEFT_TOP, 0,0);
1880
1881 pF[ 6]=G4Facet(3,BOTTOM, -4,BACK_LEFT, -11,BACK_RIGHT, 0,0);
1882 pF[ 7]=G4Facet(4,LEFT_BACK, -8,BACK_TOP, -11,BACK_BOTTOM, 0,0);
1883 pF[ 8]=G4Facet(8,TOP, -7,BACK_RIGHT, -11,BACK_LEFT, 0,0);
1884 pF[ 9]=G4Facet(7,RIGHT_BACK, -3,BACK_BOTTOM, -11,BACK_TOP, 0,0);
1885
1886 pF[10]=G4Facet(2,BOTTOM, -3,RIGHT_BACK, -10,RIGHT_FRONT, 0,0);
1887 pF[11]=G4Facet(3,BACK_RIGHT, -7,RIGHT_TOP, -10,RIGHT_BOTTOM, 0,0);
1888 pF[12]=G4Facet(7,TOP, -6,RIGHT_FRONT, -10,RIGHT_BACK, 0,0);
1889 pF[13]=G4Facet(6,FRONT_RIGHT,-2,RIGHT_BOTTOM,-10,RIGHT_TOP, 0,0);
1890
1891 pF[14]=G4Facet(1,BOTTOM, -2,FRONT_RIGHT, -9,FRONT_LEFT, 0,0);
1892 pF[15]=G4Facet(2,RIGHT_FRONT,-6,FRONT_TOP, -9,FRONT_BOTTOM, 0,0);
1893 pF[16]=G4Facet(6,TOP, -5,FRONT_LEFT, -9,FRONT_RIGHT, 0,0);
1894 pF[17]=G4Facet(5,LEFT_FRONT, -1,FRONT_BOTTOM, -9,FRONT_TOP, 0,0);
1895
1896 pF[18]=G4Facet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
1897
1898 return 0;
1899}
1900
1901G4int
1902HepPolyhedron::createPolyhedron(G4int Nnodes, G4int Nfaces,
1903 const G4double xyz[][3],
1904 const G4int faces[][4])
1905/***********************************************************************
1906 * *
1907 * Name: createPolyhedron Date: 05.11.02 *
1908 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1909 * *
1910 * Function: Creates user defined polyhedron *
1911 * *
1912 * Input: Nnodes - number of nodes *
1913 * Nfaces - number of faces *
1914 * nodes[][3] - node coordinates *
1915 * faces[][4] - faces *
1916 * *
1917 ***********************************************************************/
1918{
1919 AllocateMemory(Nnodes, Nfaces);
1920 if (nvert == 0) return 1;
1921
1922 for (G4int i=0; i<Nnodes; i++) {
1923 pV[i+1] = G4Point3D(xyz[i][0], xyz[i][1], xyz[i][2]);
1924 }
1925 for (G4int k=0; k<Nfaces; k++) {
1926 pF[k+1] = G4Facet(faces[k][0],0,faces[k][1],0,faces[k][2],0,faces[k][3],0);
1927 }
1928 SetReferences();
1929 return 0;
1930}
1931
1932HepPolyhedronTrd2::HepPolyhedronTrd2(G4double Dx1, G4double Dx2,
1933 G4double Dy1, G4double Dy2,
1934 G4double Dz)
1935/***********************************************************************
1936 * *
1937 * Name: HepPolyhedronTrd2 Date: 22.07.96 *
1938 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1939 * *
1940 * Function: Create GEANT4 TRD2-trapezoid *
1941 * *
1942 * Input: Dx1 - half-length along X at -Dz 8----7 *
1943 * Dx2 - half-length along X ay +Dz 5----6 ! *
1944 * Dy1 - half-length along Y ay -Dz ! 4-!--3 *
1945 * Dy2 - half-length along Y ay +Dz 1----2 *
1946 * Dz - half-length along Z *
1947 * *
1948 ***********************************************************************/
1949{
1950 AllocateMemory(8,6);
1951
1952 pV[1] = G4Point3D(-Dx1,-Dy1,-Dz);
1953 pV[2] = G4Point3D( Dx1,-Dy1,-Dz);
1954 pV[3] = G4Point3D( Dx1, Dy1,-Dz);
1955 pV[4] = G4Point3D(-Dx1, Dy1,-Dz);
1956 pV[5] = G4Point3D(-Dx2,-Dy2, Dz);
1957 pV[6] = G4Point3D( Dx2,-Dy2, Dz);
1958 pV[7] = G4Point3D( Dx2, Dy2, Dz);
1959 pV[8] = G4Point3D(-Dx2, Dy2, Dz);
1960
1961 CreatePrism();
1962}
1963
1964HepPolyhedronTrd2::~HepPolyhedronTrd2() = default;
1965
1966HepPolyhedronTrd1::HepPolyhedronTrd1(G4double Dx1, G4double Dx2,
1967 G4double Dy, G4double Dz)
1968 : HepPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
1969
1970HepPolyhedronTrd1::~HepPolyhedronTrd1() = default;
1971
1972HepPolyhedronBox::HepPolyhedronBox(G4double Dx, G4double Dy, G4double Dz)
1973 : HepPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1974
1975HepPolyhedronBox::~HepPolyhedronBox() = default;
1976
1977HepPolyhedronTrap::HepPolyhedronTrap(G4double Dz,
1978 G4double Theta,
1979 G4double Phi,
1980 G4double Dy1,
1981 G4double Dx1,
1982 G4double Dx2,
1983 G4double Alp1,
1984 G4double Dy2,
1985 G4double Dx3,
1986 G4double Dx4,
1987 G4double Alp2)
1988/***********************************************************************
1989 * *
1990 * Name: HepPolyhedronTrap Date: 20.11.96 *
1991 * Author: E.Chernyaev Revised: *
1992 * *
1993 * Function: Create GEANT4 TRAP-trapezoid *
1994 * *
1995 * Input: DZ - half-length in Z *
1996 * Theta,Phi - polar angles of the line joining centres of the *
1997 * faces at Z=-Dz and Z=+Dz *
1998 * Dy1 - half-length in Y of the face at Z=-Dz *
1999 * Dx1 - half-length in X of low edge of the face at Z=-Dz *
2000 * Dx2 - half-length in X of top edge of the face at Z=-Dz *
2001 * Alp1 - angle between Y-axis and the median joining top and *
2002 * low edges of the face at Z=-Dz *
2003 * Dy2 - half-length in Y of the face at Z=+Dz *
2004 * Dx3 - half-length in X of low edge of the face at Z=+Dz *
2005 * Dx4 - half-length in X of top edge of the face at Z=+Dz *
2006 * Alp2 - angle between Y-axis and the median joining top and *
2007 * low edges of the face at Z=+Dz *
2008 * *
2009 ***********************************************************************/
2010{
2011 G4double DzTthetaCphi = Dz*std::tan(Theta)*std::cos(Phi);
2012 G4double DzTthetaSphi = Dz*std::tan(Theta)*std::sin(Phi);
2013 G4double Dy1Talp1 = Dy1*std::tan(Alp1);
2014 G4double Dy2Talp2 = Dy2*std::tan(Alp2);
2015
2016 AllocateMemory(8,6);
2017
2018 pV[1] = G4Point3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
2019 pV[2] = G4Point3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
2020 pV[3] = G4Point3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
2021 pV[4] = G4Point3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
2022 pV[5] = G4Point3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
2023 pV[6] = G4Point3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
2024 pV[7] = G4Point3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
2025 pV[8] = G4Point3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
2026
2027 CreatePrism();
2028}
2029
2030HepPolyhedronTrap::~HepPolyhedronTrap() = default;
2031
2032HepPolyhedronPara::HepPolyhedronPara(G4double Dx, G4double Dy, G4double Dz,
2033 G4double Alpha, G4double Theta,
2034 G4double Phi)
2035 : HepPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
2036
2037HepPolyhedronPara::~HepPolyhedronPara() = default;
2038
2039HepPolyhedronParaboloid::HepPolyhedronParaboloid(G4double r1,
2040 G4double r2,
2041 G4double dz,
2042 G4double sPhi,
2043 G4double dPhi)
2044/***********************************************************************
2045 * *
2046 * Name: HepPolyhedronParaboloid Date: 28.06.07 *
2047 * Author: L.Lindroos, T.Nikitina (CERN), July 2007 Revised: 28.06.07 *
2048 * *
2049 * Function: Constructor for paraboloid *
2050 * *
2051 * Input: r1 - inside and outside radiuses at -Dz *
2052 * r2 - inside and outside radiuses at +Dz *
2053 * dz - half length in Z *
2054 * sPhi - starting angle of the segment *
2055 * dPhi - segment range *
2056 * *
2057 ***********************************************************************/
2058{
2059 static const G4double wholeCircle=twopi;
2060
2061 // C H E C K I N P U T P A R A M E T E R S
2062
2063 G4int k = 0;
2064 if (r1 < 0. || r2 <= 0.) k = 1;
2065
2066 if (dz <= 0.) k += 2;
2067
2068 G4double phi1, phi2, dphi;
2069
2070 if(dPhi < 0.)
2071 {
2072 phi2 = sPhi; phi1 = phi2 + dPhi;
2073 }
2074 else if(dPhi == 0.)
2075 {
2076 phi1 = sPhi; phi2 = phi1 + wholeCircle;
2077 }
2078 else
2079 {
2080 phi1 = sPhi; phi2 = phi1 + dPhi;
2081 }
2082 dphi = phi2 - phi1;
2083
2084 if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
2085 if (dphi > wholeCircle) k += 4;
2086
2087 if (k != 0) {
2088 std::cerr << "HepPolyhedronParaboloid: error in input parameters";
2089 if ((k & 1) != 0) std::cerr << " (radiuses)";
2090 if ((k & 2) != 0) std::cerr << " (half-length)";
2091 if ((k & 4) != 0) std::cerr << " (angles)";
2092 std::cerr << std::endl;
2093 std::cerr << " r1=" << r1;
2094 std::cerr << " r2=" << r2;
2095 std::cerr << " dz=" << dz << " sPhi=" << sPhi << " dPhi=" << dPhi
2096 << std::endl;
2097 return;
2098 }
2099
2100 // P R E P A R E T W O P O L Y L I N E S
2101
2102 G4int n = GetNumberOfRotationSteps();
2103 G4double dl = (r2 - r1) / n;
2104 G4double k1 = (r2*r2 - r1*r1) / 2 / dz;
2105 G4double k2 = (r2*r2 + r1*r1) / 2;
2106
2107 auto zz = new G4double[n + 2], rr = new G4double[n + 2];
2108
2109 zz[0] = dz;
2110 rr[0] = r2;
2111
2112 for(G4int i = 1; i < n - 1; i++)
2113 {
2114 rr[i] = rr[i-1] - dl;
2115 zz[i] = (rr[i]*rr[i] - k2) / k1;
2116 if(rr[i] < 0)
2117 {
2118 rr[i] = 0;
2119 zz[i] = 0;
2120 }
2121 }
2122
2123 zz[n-1] = -dz;
2124 rr[n-1] = r1;
2125
2126 zz[n] = dz;
2127 rr[n] = 0;
2128
2129 zz[n+1] = -dz;
2130 rr[n+1] = 0;
2131
2132 // R O T A T E P O L Y L I N E S
2133
2134 RotateAroundZ(0, phi1, dphi, n, 2, zz, rr, -1, -1);
2135 SetReferences();
2136
2137 delete [] zz;
2138 delete [] rr;
2139}
2140
2141HepPolyhedronParaboloid::~HepPolyhedronParaboloid() = default;
2142
2143HepPolyhedronHype::HepPolyhedronHype(G4double r1,
2144 G4double r2,
2145 G4double sqrtan1,
2146 G4double sqrtan2,
2147 G4double halfZ)
2148/***********************************************************************
2149 * *
2150 * Name: HepPolyhedronHype Date: 14.04.08 *
2151 * Author: Tatiana Nikitina (CERN) Revised: 14.04.08 *
2152 * Evgueni Tcherniaev 01.12.20 *
2153 * *
2154 * Function: Constructor for Hype *
2155 * *
2156 * Input: r1 - inside radius at z=0 *
2157 * r2 - outside radiuses at z=0 *
2158 * sqrtan1 - sqr of tan of Inner Stereo Angle *
2159 * sqrtan2 - sqr of tan of Outer Stereo Angle *
2160 * halfZ - half length in Z *
2161 * *
2162 ***********************************************************************/
2163{
2164 static const G4double wholeCircle = twopi;
2165
2166 // C H E C K I N P U T P A R A M E T E R S
2167
2168 G4int k = 0;
2169 if (r1 < 0. || r2 < 0. || r1 >= r2) k = 1;
2170 if (halfZ <= 0.) k += 2;
2171 if (sqrtan1 < 0.|| sqrtan2 < 0.) k += 4;
2172
2173 if (k != 0)
2174 {
2175 std::cerr << "HepPolyhedronHype: error in input parameters";
2176 if ((k & 1) != 0) std::cerr << " (radiuses)";
2177 if ((k & 2) != 0) std::cerr << " (half-length)";
2178 if ((k & 4) != 0) std::cerr << " (angles)";
2179 std::cerr << std::endl;
2180 std::cerr << " r1=" << r1 << " r2=" << r2;
2181 std::cerr << " halfZ=" << halfZ << " sqrTan1=" << sqrtan1
2182 << " sqrTan2=" << sqrtan2
2183 << std::endl;
2184 return;
2185 }
2186
2187 // P R E P A R E T W O P O L Y L I N E S
2188
2189 G4int ns = std::max(3, GetNumberOfRotationSteps()/4);
2190 G4int nz1 = (sqrtan1 == 0.) ? 2 : ns + 1;
2191 G4int nz2 = (sqrtan2 == 0.) ? 2 : ns + 1;
2192 auto zz = new G4double[nz1 + nz2];
2193 auto rr = new G4double[nz1 + nz2];
2194
2195 // external polyline
2196 G4double dz2 = 2.*halfZ/(nz2 - 1);
2197 for(G4int i = 0; i < nz2; ++i)
2198 {
2199 zz[i] = halfZ - dz2*i;
2200 rr[i] = std::sqrt(sqrtan2*zz[i]*zz[i] + r2*r2);
2201 }
2202
2203 // internal polyline
2204 G4double dz1 = 2.*halfZ/(nz1 - 1);
2205 for(G4int i = 0; i < nz1; ++i)
2206 {
2207 G4int j = nz2 + i;
2208 zz[j] = halfZ - dz1*i;
2209 rr[j] = std::sqrt(sqrtan1*zz[j]*zz[j] + r1*r1);
2210 }
2211
2212 // R O T A T E P O L Y L I N E S
2213
2214 RotateAroundZ(0, 0., wholeCircle, nz2, nz1, zz, rr, -1, -1);
2215 SetReferences();
2216
2217 delete [] zz;
2218 delete [] rr;
2219}
2220
2221HepPolyhedronHype::~HepPolyhedronHype() = default;
2222
2223HepPolyhedronCons::HepPolyhedronCons(G4double Rmn1,
2224 G4double Rmx1,
2225 G4double Rmn2,
2226 G4double Rmx2,
2227 G4double Dz,
2228 G4double Phi1,
2229 G4double Dphi)
2230/***********************************************************************
2231 * *
2232 * Name: HepPolyhedronCons::HepPolyhedronCons Date: 15.12.96 *
2233 * Author: E.Chernyaev (IHEP/Protvino) Revised: 15.12.96 *
2234 * *
2235 * Function: Constructor for CONS, TUBS, CONE, TUBE *
2236 * *
2237 * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz *
2238 * Rmn2, Rmx2 - inside and outside radiuses at +Dz *
2239 * Dz - half length in Z *
2240 * Phi1 - starting angle of the segment *
2241 * Dphi - segment range *
2242 * *
2243 ***********************************************************************/
2244{
2245 static const G4double wholeCircle=twopi;
2246
2247 // C H E C K I N P U T P A R A M E T E R S
2248
2249 G4int k = 0;
2250 if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
2251 if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
2252 if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
2253
2254 if (Dz <= 0.) k += 2;
2255
2256 G4double phi1, phi2, dphi;
2257 if (Dphi < 0.) {
2258 phi2 = Phi1; phi1 = phi2 - Dphi;
2259 }else if (Dphi == 0.) {
2260 phi1 = Phi1; phi2 = phi1 + wholeCircle;
2261 }else{
2262 phi1 = Phi1; phi2 = phi1 + Dphi;
2263 }
2264 dphi = phi2 - phi1;
2265 if (std::abs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
2266 if (dphi > wholeCircle) k += 4;
2267
2268 if (k != 0) {
2269 std::cerr << "HepPolyhedronCone(s)/Tube(s): error in input parameters";
2270 if ((k & 1) != 0) std::cerr << " (radiuses)";
2271 if ((k & 2) != 0) std::cerr << " (half-length)";
2272 if ((k & 4) != 0) std::cerr << " (angles)";
2273 std::cerr << std::endl;
2274 std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
2275 std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
2276 std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
2277 << std::endl;
2278 return;
2279 }
2280
2281 // P R E P A R E T W O P O L Y L I N E S
2282
2283 G4double zz[4], rr[4];
2284 zz[0] = Dz;
2285 zz[1] = -Dz;
2286 zz[2] = Dz;
2287 zz[3] = -Dz;
2288 rr[0] = Rmx2;
2289 rr[1] = Rmx1;
2290 rr[2] = Rmn2;
2291 rr[3] = Rmn1;
2292
2293 // R O T A T E P O L Y L I N E S
2294
2295 RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1);
2296 SetReferences();
2297}
2298
2299HepPolyhedronCons::~HepPolyhedronCons() = default;
2300
2301HepPolyhedronCone::HepPolyhedronCone(G4double Rmn1, G4double Rmx1,
2302 G4double Rmn2, G4double Rmx2,
2303 G4double Dz) :
2304 HepPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
2305
2306HepPolyhedronCone::~HepPolyhedronCone() = default;
2307
2308HepPolyhedronTubs::HepPolyhedronTubs(G4double Rmin, G4double Rmax,
2309 G4double Dz,
2310 G4double Phi1, G4double Dphi)
2311 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
2312
2313HepPolyhedronTubs::~HepPolyhedronTubs() = default;
2314
2315HepPolyhedronTube::HepPolyhedronTube (G4double Rmin, G4double Rmax,
2316 G4double Dz)
2317 : HepPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
2318
2319HepPolyhedronTube::~HepPolyhedronTube () = default;
2320
2321HepPolyhedronPgon::HepPolyhedronPgon(G4double phi,
2322 G4double dphi,
2323 G4int npdv,
2324 G4int nz,
2325 const G4double *z,
2326 const G4double *rmin,
2327 const G4double *rmax)
2328/***********************************************************************
2329 * *
2330 * Name: HepPolyhedronPgon Date: 09.12.96 *
2331 * Author: E.Chernyaev Revised: *
2332 * *
2333 * Function: Constructor of polyhedron for PGON, PCON *
2334 * *
2335 * Input: phi - initial phi *
2336 * dphi - delta phi *
2337 * npdv - number of steps along phi *
2338 * nz - number of z-planes (at least two) *
2339 * z[] - z coordinates of the slices *
2340 * rmin[] - smaller r at the slices *
2341 * rmax[] - bigger r at the slices *
2342 * *
2343 ***********************************************************************/
2344{
2345 // C H E C K I N P U T P A R A M E T E R S
2346
2347 if (dphi <= 0. || dphi > twopi) {
2348 std::cerr
2349 << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
2350 << std::endl;
2351 return;
2352 }
2353
2354 if (nz < 2) {
2355 std::cerr
2356 << "HepPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
2357 << std::endl;
2358 return;
2359 }
2360
2361 if (npdv < 0) {
2362 std::cerr
2363 << "HepPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
2364 << std::endl;
2365 return;
2366 }
2367
2368 G4int i;
2369 for (i=0; i<nz; i++) {
2370 if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
2371 std::cerr
2372 << "HepPolyhedronPgon: error in radiuses rmin[" << i << "]="
2373 << rmin[i] << " rmax[" << i << "]=" << rmax[i]
2374 << std::endl;
2375 return;
2376 }
2377 }
2378
2379 // P R E P A R E T W O P O L Y L I N E S
2380
2381 G4double *zz, *rr;
2382 zz = new G4double[2*nz];
2383 rr = new G4double[2*nz];
2384
2385 if (z[0] > z[nz-1]) {
2386 for (i=0; i<nz; i++) {
2387 zz[i] = z[i];
2388 rr[i] = rmax[i];
2389 zz[i+nz] = z[i];
2390 rr[i+nz] = rmin[i];
2391 }
2392 }else{
2393 for (i=0; i<nz; i++) {
2394 zz[i] = z[nz-i-1];
2395 rr[i] = rmax[nz-i-1];
2396 zz[i+nz] = z[nz-i-1];
2397 rr[i+nz] = rmin[nz-i-1];
2398 }
2399 }
2400
2401 // R O T A T E P O L Y L I N E S
2402
2403 G4int nodeVis = 1;
2404 G4int edgeVis = (npdv == 0) ? -1 : 1;
2405 RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, nodeVis, edgeVis);
2406 SetReferences();
2407
2408 delete [] zz;
2409 delete [] rr;
2410}
2411
2412HepPolyhedronPgon::HepPolyhedronPgon(G4double phi,
2413 G4double dphi,
2414 G4int npdv,
2415 const std::vector<G4TwoVector> &rz)
2416/***********************************************************************
2417 * *
2418 * Name: HepPolyhedronPgon Date: 12.05.21 *
2419 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
2420 * *
2421 * Function: Constructor of polyhedron for PGON, PCON *
2422 * *
2423 * Input: phi - initial phi *
2424 * dphi - delta phi *
2425 * npdv - number of steps along phi *
2426 * rz - rz-contour *
2427 * *
2428 ***********************************************************************/
2429{
2430 // C H E C K I N P U T P A R A M E T E R S
2431
2432 if (dphi <= 0. || dphi > twopi) {
2433 std::cerr
2434 << "HepPolyhedronPgon/Pcon: wrong delta phi = " << dphi
2435 << std::endl;
2436 return;
2437 }
2438
2439 if (npdv < 0) {
2440 std::cerr
2441 << "HepPolyhedronPgon/Pcon: error in number of phi-steps = " << npdv
2442 << std::endl;
2443 return;
2444 }
2445
2446 G4int nrz = (G4int)rz.size();
2447 if (nrz < 3) {
2448 std::cerr
2449 << "HepPolyhedronPgon/Pcon: invalid number of nodes in rz-contour = " << nrz
2450 << std::endl;
2451 return;
2452 }
2453
2454 // R O T A T E P O L Y L I N E
2455
2456 G4int nodeVis = 1;
2457 G4int edgeVis = (npdv == 0) ? -1 : 1;
2458 RotateContourAroundZ(npdv, phi, dphi, rz, nodeVis, edgeVis);
2459 SetReferences();
2460}
2461
2462HepPolyhedronPgon::~HepPolyhedronPgon() = default;
2463
2464HepPolyhedronPcon::HepPolyhedronPcon(G4double phi, G4double dphi, G4int nz,
2465 const G4double *z,
2466 const G4double *rmin,
2467 const G4double *rmax)
2468 : HepPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
2469
2470HepPolyhedronPcon::HepPolyhedronPcon(G4double phi, G4double dphi,
2471 const std::vector<G4TwoVector> &rz)
2472 : HepPolyhedronPgon(phi, dphi, 0, rz) {}
2473
2474HepPolyhedronPcon::~HepPolyhedronPcon() = default;
2475
2476HepPolyhedronSphere::HepPolyhedronSphere(G4double rmin, G4double rmax,
2477 G4double phi, G4double dphi,
2478 G4double the, G4double dthe)
2479/***********************************************************************
2480 * *
2481 * Name: HepPolyhedronSphere Date: 11.12.96 *
2482 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
2483 * *
2484 * Function: Constructor of polyhedron for SPHERE *
2485 * *
2486 * Input: rmin - internal radius *
2487 * rmax - external radius *
2488 * phi - initial phi *
2489 * dphi - delta phi *
2490 * the - initial theta *
2491 * dthe - delta theta *
2492 * *
2493 ***********************************************************************/
2494{
2495 // C H E C K I N P U T P A R A M E T E R S
2496
2497 if (dphi <= 0. || dphi > twopi) {
2498 std::cerr
2499 << "HepPolyhedronSphere: wrong delta phi = " << dphi
2500 << std::endl;
2501 return;
2502 }
2503
2504 if (the < 0. || the > pi) {
2505 std::cerr
2506 << "HepPolyhedronSphere: wrong theta = " << the
2507 << std::endl;
2508 return;
2509 }
2510
2511 if (dthe <= 0. || dthe > pi) {
2512 std::cerr
2513 << "HepPolyhedronSphere: wrong delta theta = " << dthe
2514 << std::endl;
2515 return;
2516 }
2517
2518 if (the+dthe > pi) {
2519 std::cerr
2520 << "HepPolyhedronSphere: wrong theta + delta theta = "
2521 << the << " " << dthe
2522 << std::endl;
2523 return;
2524 }
2525
2526 if (rmin < 0. || rmin >= rmax) {
2527 std::cerr
2528 << "HepPolyhedronSphere: error in radiuses"
2529 << " rmin=" << rmin << " rmax=" << rmax
2530 << std::endl;
2531 return;
2532 }
2533
2534 // P R E P A R E T W O P O L Y L I N E S
2535
2536 G4int nds = (GetNumberOfRotationSteps() + 1) / 2;
2537 G4int np1 = G4int(dthe*nds/pi+.5) + 1;
2538 if (np1 <= 1) np1 = 2;
2539 G4int np2 = rmin < spatialTolerance ? 1 : np1;
2540
2541 G4double *zz, *rr;
2542 zz = new G4double[np1+np2];
2543 rr = new G4double[np1+np2];
2544
2545 G4double a = dthe/(np1-1);
2546 G4double cosa, sina;
2547 for (G4int i=0; i<np1; i++) {
2548 cosa = std::cos(the+i*a);
2549 sina = std::sin(the+i*a);
2550 zz[i] = rmax*cosa;
2551 rr[i] = rmax*sina;
2552 if (np2 > 1) {
2553 zz[i+np1] = rmin*cosa;
2554 rr[i+np1] = rmin*sina;
2555 }
2556 }
2557 if (np2 == 1) {
2558 zz[np1] = 0.;
2559 rr[np1] = 0.;
2560 }
2561
2562 // R O T A T E P O L Y L I N E S
2563
2564 RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
2565 SetReferences();
2566
2567 delete [] zz;
2568 delete [] rr;
2569}
2570
2571HepPolyhedronSphere::~HepPolyhedronSphere() = default;
2572
2573HepPolyhedronTorus::HepPolyhedronTorus(G4double rmin,
2574 G4double rmax,
2575 G4double rtor,
2576 G4double phi,
2577 G4double dphi)
2578/***********************************************************************
2579 * *
2580 * Name: HepPolyhedronTorus Date: 11.12.96 *
2581 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
2582 * *
2583 * Function: Constructor of polyhedron for TORUS *
2584 * *
2585 * Input: rmin - internal radius *
2586 * rmax - external radius *
2587 * rtor - radius of torus *
2588 * phi - initial phi *
2589 * dphi - delta phi *
2590 * *
2591 ***********************************************************************/
2592{
2593 // C H E C K I N P U T P A R A M E T E R S
2594
2595 if (dphi <= 0. || dphi > twopi) {
2596 std::cerr
2597 << "HepPolyhedronTorus: wrong delta phi = " << dphi
2598 << std::endl;
2599 return;
2600 }
2601
2602 if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
2603 std::cerr
2604 << "HepPolyhedronTorus: error in radiuses"
2605 << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
2606 << std::endl;
2607 return;
2608 }
2609
2610 // P R E P A R E T W O P O L Y L I N E S
2611
2612 G4int np1 = GetNumberOfRotationSteps();
2613 G4int np2 = rmin < spatialTolerance ? 1 : np1;
2614
2615 G4double *zz, *rr;
2616 zz = new G4double[np1+np2];
2617 rr = new G4double[np1+np2];
2618
2619 G4double a = twopi/np1;
2620 G4double cosa, sina;
2621 for (G4int i=0; i<np1; i++) {
2622 cosa = std::cos(i*a);
2623 sina = std::sin(i*a);
2624 zz[i] = rmax*cosa;
2625 rr[i] = rtor+rmax*sina;
2626 if (np2 > 1) {
2627 zz[i+np1] = rmin*cosa;
2628 rr[i+np1] = rtor+rmin*sina;
2629 }
2630 }
2631 if (np2 == 1) {
2632 zz[np1] = 0.;
2633 rr[np1] = rtor;
2634 np2 = -1;
2635 }
2636
2637 // R O T A T E P O L Y L I N E S
2638
2639 RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1);
2640 SetReferences();
2641
2642 delete [] zz;
2643 delete [] rr;
2644}
2645
2646HepPolyhedronTorus::~HepPolyhedronTorus() = default;
2647
2648HepPolyhedronTet::HepPolyhedronTet(const G4double p0[3],
2649 const G4double p1[3],
2650 const G4double p2[3],
2651 const G4double p3[3])
2652/***********************************************************************
2653 * *
2654 * Name: HepPolyhedronTet Date: 21.02.2020 *
2655 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
2656 * *
2657 * Function: Constructor of polyhedron for TETrahedron *
2658 * *
2659 * Input: p0,p1,p2,p3 - vertices *
2660 * *
2661 ***********************************************************************/
2662{
2663 AllocateMemory(4,4);
2664
2665 pV[1].set(p0[0], p0[1], p0[2]);
2666 pV[2].set(p1[0], p1[1], p1[2]);
2667 pV[3].set(p2[0], p2[1], p2[2]);
2668 pV[4].set(p3[0], p3[1], p3[2]);
2669
2670 G4Vector3D v1(pV[2] - pV[1]);
2671 G4Vector3D v2(pV[3] - pV[1]);
2672 G4Vector3D v3(pV[4] - pV[1]);
2673
2674 if (v1.cross(v2).dot(v3) < 0.)
2675 {
2676 pV[3].set(p3[0], p3[1], p3[2]);
2677 pV[4].set(p2[0], p2[1], p2[2]);
2678 }
2679
2680 pF[1] = G4Facet(1,2, 3,4, 2,3);
2681 pF[2] = G4Facet(1,3, 4,4, 3,1);
2682 pF[3] = G4Facet(1,1, 2,4, 4,2);
2683 pF[4] = G4Facet(2,1, 3,2, 4,3);
2684}
2685
2686HepPolyhedronTet::~HepPolyhedronTet() = default;
2687
2688HepPolyhedronEllipsoid::HepPolyhedronEllipsoid(G4double ax, G4double by,
2689 G4double cz, G4double zCut1,
2690 G4double zCut2)
2691/***********************************************************************
2692 * *
2693 * Name: HepPolyhedronEllipsoid Date: 25.02.05 *
2694 * Author: G.Guerrieri Revised: *
2695 * Evgueni Tcherniaev 20.01.22 *
2696 * *
2697 * Function: Constructor of polyhedron for ELLIPSOID *
2698 * *
2699 * Input: ax - semiaxis x *
2700 * by - semiaxis y *
2701 * cz - semiaxis z *
2702 * zCut1 - lower cut plane level (solid lies above this plane) *
2703 * zCut2 - upper cut plane level (solid lies below this plane) *
2704 * *
2705 ***********************************************************************/
2706{
2707 // C H E C K I N P U T P A R A M E T E R S
2708
2709 if (zCut1 >= cz || zCut2 <= -cz || zCut1 > zCut2) {
2710 std::cerr << "HepPolyhedronEllipsoid: wrong zCut1 = " << zCut1
2711 << " zCut2 = " << zCut2
2712 << " for given cz = " << cz << std::endl;
2713 return;
2714 }
2715 if (cz <= 0.0) {
2716 std::cerr << "HepPolyhedronEllipsoid: bad z semi-axis: cz = " << cz
2717 << std::endl;
2718 return;
2719 }
2720
2721 // P R E P A R E T W O P O L Y L I N E S
2722 // generate sphere of radius cz first, then rescale x and y later
2723
2724 G4double sthe = std::acos(zCut2/cz);
2725 G4double dthe = std::acos(zCut1/cz) - sthe;
2726 G4int nds = (GetNumberOfRotationSteps() + 1)/2;
2727 G4int np1 = G4int(dthe*nds/pi + 0.5) + 1;
2728 if (np1 <= 1) np1 = 2;
2729 G4int np2 = 2;
2730
2731 G4double *zz, *rr;
2732 zz = new G4double[np1 + np2];
2733 rr = new G4double[np1 + np2];
2734 if ((zz == nullptr) || (rr == nullptr))
2735 {
2736 G4Exception("HepPolyhedronEllipsoid::HepPolyhedronEllipsoid",
2737 "greps1002", FatalException, "Out of memory");
2738 }
2739
2740 G4double a = dthe/(np1 - 1);
2741 G4double cosa, sina;
2742 for (G4int i = 0; i < np1; ++i)
2743 {
2744 cosa = std::cos(sthe + i*a);
2745 sina = std::sin(sthe + i*a);
2746 zz[i] = cz*cosa;
2747 rr[i] = cz*sina;
2748 }
2749 zz[np1 + 0] = zCut2;
2750 rr[np1 + 0] = 0.;
2751 zz[np1 + 1] = zCut1;
2752 rr[np1 + 1] = 0.;
2753
2754 // R O T A T E P O L Y L I N E S
2755
2756 RotateAroundZ(0, 0., twopi, np1, np2, zz, rr, -1, -1);
2757 SetReferences();
2758
2759 delete [] zz;
2760 delete [] rr;
2761
2762 // rescale x and y vertex coordinates
2763 G4double kx = ax/cz;
2764 G4double ky = by/cz;
2765 G4Point3D* p = pV;
2766 for (G4int i = 0; i < nvert; ++i, ++p)
2767 {
2768 p->setX(p->x()*kx);
2769 p->setY(p->y()*ky);
2770 }
2771}
2772
2773HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid() = default;
2774
2775HepPolyhedronEllipticalCone::HepPolyhedronEllipticalCone(G4double ax,
2776 G4double ay,
2777 G4double h,
2778 G4double zTopCut)
2779/***********************************************************************
2780 * *
2781 * Name: HepPolyhedronEllipticalCone Date: 8.9.2005 *
2782 * Author: D.Anninos Revised: 9.9.2005 *
2783 * *
2784 * Function: Constructor for EllipticalCone *
2785 * *
2786 * Input: ax, ay - X & Y semi axes at z = 0 *
2787 * h - height of full cone *
2788 * zTopCut - Top Cut in Z Axis *
2789 * *
2790 ***********************************************************************/
2791{
2792 // C H E C K I N P U T P A R A M E T E R S
2793
2794 G4int k = 0;
2795 if ( (ax <= 0.) || (ay <= 0.) || (h <= 0.) || (zTopCut <= 0.) ) { k = 1; }
2796
2797 if (k != 0) {
2798 std::cerr << "HepPolyhedronCone: error in input parameters";
2799 std::cerr << std::endl;
2800 return;
2801 }
2802
2803 // P R E P A R E T W O P O L Y L I N E S
2804
2805 zTopCut = (h >= zTopCut ? zTopCut : h);
2806
2807 G4double *zz, *rr;
2808 zz = new G4double[4];
2809 rr = new G4double[4];
2810 zz[0] = zTopCut;
2811 zz[1] = -zTopCut;
2812 zz[2] = zTopCut;
2813 zz[3] = -zTopCut;
2814 rr[0] = (h-zTopCut);
2815 rr[1] = (h+zTopCut);
2816 rr[2] = 0.;
2817 rr[3] = 0.;
2818
2819 // R O T A T E P O L Y L I N E S
2820
2821 RotateAroundZ(0, 0., twopi, 2, 2, zz, rr, -1, -1);
2822 SetReferences();
2823
2824 delete [] zz;
2825 delete [] rr;
2826
2827 // rescale x and y vertex coordinates
2828 {
2829 G4Point3D * p= pV;
2830 for (G4int i=0; i<nvert; i++, p++) {
2831 p->setX( p->x() * ax );
2832 p->setY( p->y() * ay );
2833 }
2834 }
2835}
2836
2837HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone() = default;
2838
2839HepPolyhedronHyperbolicMirror::HepPolyhedronHyperbolicMirror(G4double a,
2840 G4double h,
2841 G4double r)
2842/***********************************************************************
2843 * *
2844 * Name: HepPolyhedronHyperbolicMirror Date: 22.02.2020 *
2845 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
2846 * *
2847 * Function: Create polyhedron for Hyperbolic mirror *
2848 * *
2849 * Input: a - half-separation *
2850 * h - height *
2851 * r - radius *
2852 * *
2853 ***********************************************************************/
2854{
2855 G4double H = std::abs(h);
2856 G4double R = std::abs(r);
2857 G4double A = std::abs(a);
2858 G4double B = A*R/std::sqrt(2*A*H + H*H);
2859
2860 // P R E P A R E T W O P O L Y L I N E S
2861
2862 G4int np1 = (A == 0.) ? 2 : std::max(3, GetNumberOfRotationSteps()/4) + 1;
2863 G4int np2 = 2;
2864 G4double maxAng = (A == 0.) ? 0. : std::acosh(1. + H/A);
2865 G4double delAng = maxAng/(np1 - 1);
2866
2867 auto zz = new G4double[np1 + np2];
2868 auto rr = new G4double[np1 + np2];
2869
2870 // 1st polyline
2871 zz[0] = H;
2872 rr[0] = R;
2873 for (G4int iz = 1; iz < np1 - 1; ++iz)
2874 {
2875 G4double ang = maxAng - iz*delAng;
2876 zz[iz] = A*std::cosh(ang) - A;
2877 rr[iz] = B*std::sinh(ang);
2878 }
2879 zz[np1 - 1] = 0.;
2880 rr[np1 - 1] = 0.;
2881
2882 // 2nd polyline
2883 zz[np1] = H;
2884 rr[np1] = 0.;
2885 zz[np1 + 1] = 0.;
2886 rr[np1 + 1] = 0.;
2887
2888 // R O T A T E P O L Y L I N E S
2889
2890 G4double phi = 0.;
2891 G4double dphi = CLHEP::twopi;
2892 RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
2893 SetReferences();
2894
2895 delete [] zz;
2896 delete [] rr;
2897}
2898
2899HepPolyhedronHyperbolicMirror::~HepPolyhedronHyperbolicMirror() = default;
2900
2901HepPolyhedronTetMesh::
2902HepPolyhedronTetMesh(const std::vector<G4ThreeVector>& tetrahedra)
2903/***********************************************************************
2904 * *
2905 * Name: HepPolyhedronTetMesh Date: 26.03.2022 *
2906 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
2907 * *
2908 * Function: Create polyhedron for tetrahedron mesh *
2909 * *
2910 * Input: tetrahedra - array of tetrahedron vertices, four vertices *
2911 * per tetrahedron *
2912 * *
2913 ***********************************************************************/
2914{
2915 // Check size of input vector
2916 G4int nnodes = (G4int)tetrahedra.size();
2917 if (nnodes == 0)
2918 {
2919 std::cerr
2920 << "HepPolyhedronTetMesh: Empty tetrahedron mesh" << std::endl;
2921 return;
2922 }
2923 G4int ntet = nnodes/4;
2924 if (nnodes != ntet*4)
2925 {
2926 std::cerr << "HepPolyhedronTetMesh: Number of nodes = " << nnodes
2927 << " in tetrahedron mesh is NOT multiple of 4"
2928 << std::endl;
2929 return;
2930 }
2931
2932 // Find coincident vertices using hash table techniques.
2933 // This could be done using std::unordered_map, but the code
2934 // below runs faster.
2935 std::vector<G4int> iheads(nnodes, -1);
2936 std::vector<std::pair<G4int,G4int>> ipairs(nnodes,std::pair(-1,-1));
2937 for (G4int i = 0; i < nnodes; ++i)
2938 {
2939 // Generate hash key
2940 G4ThreeVector point = tetrahedra[i];
2941 auto key = std::hash<G4double>()(point.x());
2942 key ^= std::hash<G4double>()(point.y());
2943 key ^= std::hash<G4double>()(point.z());
2944 key %= nnodes;
2945 // Check head of the list
2946 if (iheads[key] < 0)
2947 {
2948 iheads[key] = i;
2949 ipairs[i].first = i;
2950 continue;
2951 }
2952 // Loop along the list
2953 for (G4int icur = iheads[key], iprev = 0;;)
2954 {
2955 G4int icheck = ipairs[icur].first;
2956 if (tetrahedra[icheck] == point)
2957 {
2958 ipairs[i].first = icheck; // coincident vertex
2959 break;
2960 }
2961 iprev = icur;
2962 icur = ipairs[icur].second;
2963 // Append vertex to the list
2964 if (icur < 0)
2965 {
2966 ipairs[i].first = i;
2967 ipairs[iprev].second = i;
2968 break;
2969 }
2970 }
2971 }
2972
2973 // Create vector of original facets
2974 struct facet
2975 {
2976 G4int i1, i2, i3;
2977 facet() : i1(0), i2(0), i3(0) {};
2978 facet(G4int k1, G4int k2, G4int k3) : i1(k1), i2(k2), i3(k3) {};
2979 };
2980 G4int nfacets = nnodes;
2981 std::vector<facet> ifacets(nfacets);
2982 for (G4int i = 0; i < nfacets; i += 4)
2983 {
2984 G4int i0 = ipairs[i + 0].first;
2985 G4int i1 = ipairs[i + 1].first;
2986 G4int i2 = ipairs[i + 2].first;
2987 G4int i3 = ipairs[i + 3].first;
2988 if (i0 > i1) std::swap(i0, i1);
2989 if (i0 > i2) std::swap(i0, i2);
2990 if (i0 > i3) std::swap(i0, i3);
2991 if (i1 > i2) std::swap(i1, i2);
2992 if (i1 > i3) std::swap(i1, i3);
2993 G4ThreeVector e1 = tetrahedra[i1] - tetrahedra[i0];
2994 G4ThreeVector e2 = tetrahedra[i2] - tetrahedra[i0];
2995 G4ThreeVector e3 = tetrahedra[i3] - tetrahedra[i0];
2996 G4double volume = (e1.cross(e2)).dot(e3);
2997 if (volume > 0.) std::swap(i2, i3);
2998 ifacets[i + 0] = facet(i0, i1, i2);
2999 ifacets[i + 1] = facet(i0, i2, i3);
3000 ifacets[i + 2] = facet(i0, i3, i1);
3001 ifacets[i + 3] = facet(i1, i3, i2);
3002 }
3003
3004 // Find shared facets
3005 std::fill(iheads.begin(), iheads.end(), -1);
3006 std::fill(ipairs.begin(), ipairs.end(), std::pair(-1,-1));
3007 for (G4int i = 0; i < nfacets; ++i)
3008 {
3009 // Check head of the list
3010 G4int key = ifacets[i].i1;
3011 if (iheads[key] < 0)
3012 {
3013 iheads[key] = i;
3014 ipairs[i].first = i;
3015 continue;
3016 }
3017 // Loop along the list
3018 G4int i2 = ifacets[i].i2, i3 = ifacets[i].i3;
3019 for (G4int icur = iheads[key], iprev = -1;;)
3020 {
3021 G4int icheck = ipairs[icur].first;
3022 if (ifacets[icheck].i2 == i3 && ifacets[icheck].i3 == i2)
3023 {
3024 if (iprev < 0)
3025 {
3026 iheads[key] = ipairs[icur].second;
3027 }
3028 else
3029 {
3030 ipairs[iprev].second = ipairs[icur].second;
3031 }
3032 ipairs[icur].first = -1; // shared facet
3033 ipairs[icur].second = -1;
3034 break;
3035 }
3036 iprev = icur;
3037 icur = ipairs[icur].second;
3038 // Append facet to the list
3039 if (icur < 0)
3040 {
3041 ipairs[i].first = i;
3042 ipairs[iprev].second = i;
3043 break;
3044 }
3045 }
3046 }
3047
3048 // Count vertices and facets skipping shared facets
3049 std::fill(iheads.begin(), iheads.end(), -1);
3050 G4int nver = 0, nfac = 0;
3051 for (G4int i = 0; i < nfacets; ++i)
3052 {
3053 if (ipairs[i].first < 0) continue;
3054 G4int i1 = ifacets[i].i1;
3055 G4int i2 = ifacets[i].i2;
3056 G4int i3 = ifacets[i].i3;
3057 if (iheads[i1] < 0) iheads[i1] = nver++;
3058 if (iheads[i2] < 0) iheads[i2] = nver++;
3059 if (iheads[i3] < 0) iheads[i3] = nver++;
3060 nfac++;
3061 }
3062
3063 // Construct polyhedron
3064 AllocateMemory(nver, nfac);
3065 for (G4int i = 0; i < nnodes; ++i)
3066 {
3067 G4int k = iheads[i];
3068 if (k >= 0) SetVertex(k + 1, tetrahedra[i]);
3069 }
3070 for (G4int i = 0, k = 0; i < nfacets; ++i)
3071 {
3072 if (ipairs[i].first < 0) continue;
3073 G4int i1 = iheads[ifacets[i].i1] + 1;
3074 G4int i2 = iheads[ifacets[i].i2] + 1;
3075 G4int i3 = iheads[ifacets[i].i3] + 1;
3076 SetFacet(++k, i1, i2, i3);
3077 }
3078 SetReferences();
3079}
3080
3081HepPolyhedronTetMesh::~HepPolyhedronTetMesh() = default;
3082
3083HepPolyhedronBoxMesh::
3084HepPolyhedronBoxMesh(G4double sizeX, G4double sizeY, G4double sizeZ,
3085 const std::vector<G4ThreeVector>& positions)
3086/***********************************************************************
3087 * *
3088 * Name: HepPolyhedronBoxMesh Date: 07.04.2022 *
3089 * Author: E.Tcherniaev (E.Chernyaev) Revised: *
3090 * *
3091 * Function: Create polyhedron for box mesh *
3092 * *
3093 * Input: sizeX, sizeY, sizeZ - dimensions of the mesh cell *
3094 * positions - vector of cell centres *
3095 * *
3096 ***********************************************************************/
3097{
3098 G4int nbox = (G4int)positions.size();
3099 if (nbox == 0)
3100 {
3101 std::cerr << "HepPolyhedronBoxMesh: Empty box mesh" << std::endl;
3102 return;
3103 }
3104 // compute inverse dimensions
3105 G4double invx = 1./sizeX, invy = 1./sizeY, invz = 1./sizeZ;
3106 // find mesh bounding box
3107 G4ThreeVector pmin = positions[0], pmax = positions[0];
3108 for (const auto& p: positions)
3109 {
3110 if (pmin.x() > p.x()) pmin.setX(p.x());
3111 if (pmin.y() > p.y()) pmin.setY(p.y());
3112 if (pmin.z() > p.z()) pmin.setZ(p.z());
3113 if (pmax.x() < p.x()) pmax.setX(p.x());
3114 if (pmax.y() < p.y()) pmax.setY(p.y());
3115 if (pmax.z() < p.z()) pmax.setZ(p.z());
3116 }
3117 // find number of voxels
3118 G4int nx = (pmax.x() - pmin.x())*invx + 1.5;
3119 G4int ny = (pmax.y() - pmin.y())*invy + 1.5;
3120 G4int nz = (pmax.z() - pmin.z())*invz + 1.5;
3121 // create structures for voxels and node indices
3122 std::vector<char> voxels(nx*ny*nz, 0);
3123 std::vector<G4int> indices((nx+1)*(ny+1)*(nz+1), 0);
3124 // mark voxels listed in positions
3125 G4int kx = ny*nz, ky = nz;
3126 for (const auto& p: positions)
3127 {
3128 G4int ix = (p.x() - pmin.x())*invx + 0.5;
3129 G4int iy = (p.y() - pmin.y())*invy + 0.5;
3130 G4int iz = (p.z() - pmin.z())*invz + 0.5;
3131 G4int i = ix*kx + iy*ky + iz;
3132 voxels[i] = 1;
3133 }
3134 // count number of vertices and facets
3135 // set indices
3136 G4int kvx = (ny + 1)*(nz + 1), kvy = nz + 1;
3137 G4int nver = 0, nfac = 0;
3138 for (const auto& p: positions)
3139 {
3140 G4int ix = (p.x() - pmin.x())*invx + 0.5;
3141 G4int iy = (p.y() - pmin.y())*invy + 0.5;
3142 G4int iz = (p.z() - pmin.z())*invz + 0.5;
3143 //
3144 // 011 111
3145 // +---–---+
3146 // | 001 | 101
3147 // | +---–---+
3148 // | | | |
3149 // +---|---+ |
3150 // 010 | 110 |
3151 // +-------+
3152 // 000 100
3153 //
3154 G4int vcheck = 0;
3155 // check (ix - 1) side
3156 vcheck = (ix == 0) ? 0 : voxels[(ix-1)*kx + iy*ky + iz];
3157 if (vcheck == 0)
3158 {
3159 nfac++;
3160 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3161 G4int i2 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3162 G4int i3 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3163 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3164 if (indices[i1] == 0) indices[i1] = ++nver;
3165 if (indices[i2] == 0) indices[i2] = ++nver;
3166 if (indices[i3] == 0) indices[i3] = ++nver;
3167 if (indices[i4] == 0) indices[i4] = ++nver;
3168 }
3169 // check (ix + 1) side
3170 vcheck = (ix == nx - 1) ? 0 : voxels[(ix+1)*kx + iy*ky + iz];
3171 if (vcheck == 0)
3172 {
3173 nfac++;
3174 G4int i1 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3175 G4int i2 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3176 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3177 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3178 if (indices[i1] == 0) indices[i1] = ++nver;
3179 if (indices[i2] == 0) indices[i2] = ++nver;
3180 if (indices[i3] == 0) indices[i3] = ++nver;
3181 if (indices[i4] == 0) indices[i4] = ++nver;
3182 }
3183 // check (iy - 1) side
3184 vcheck = (iy == 0) ? 0 : voxels[ix*kx + (iy-1)*ky + iz];
3185 if (vcheck == 0)
3186 {
3187 nfac++;
3188 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3189 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3190 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3191 G4int i4 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3192 if (indices[i1] == 0) indices[i1] = ++nver;
3193 if (indices[i2] == 0) indices[i2] = ++nver;
3194 if (indices[i3] == 0) indices[i3] = ++nver;
3195 if (indices[i4] == 0) indices[i4] = ++nver;
3196 }
3197 // check (iy + 1) side
3198 vcheck = (iy == ny - 1) ? 0 : voxels[ix*kx + (iy+1)*ky + iz];
3199 if (vcheck == 0)
3200 {
3201 nfac++;
3202 G4int i1 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3203 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3204 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3205 G4int i4 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3206 if (indices[i1] == 0) indices[i1] = ++nver;
3207 if (indices[i2] == 0) indices[i2] = ++nver;
3208 if (indices[i3] == 0) indices[i3] = ++nver;
3209 if (indices[i4] == 0) indices[i4] = ++nver;
3210 }
3211 // check (iz - 1) side
3212 vcheck = (iz == 0) ? 0 : voxels[ix*kx + iy*ky + (iz-1)];
3213 if (vcheck == 0)
3214 {
3215 nfac++;
3216 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3217 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3218 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3219 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3220 if (indices[i1] == 0) indices[i1] = ++nver;
3221 if (indices[i2] == 0) indices[i2] = ++nver;
3222 if (indices[i3] == 0) indices[i3] = ++nver;
3223 if (indices[i4] == 0) indices[i4] = ++nver;
3224 }
3225 // check (iz + 1) side
3226 vcheck = (iz == nz - 1) ? 0 : voxels[ix*kx + iy*ky + (iz+1)];
3227 if (vcheck == 0)
3228 {
3229 nfac++;
3230 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3231 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3232 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3233 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3234 if (indices[i1] == 0) indices[i1] = ++nver;
3235 if (indices[i2] == 0) indices[i2] = ++nver;
3236 if (indices[i3] == 0) indices[i3] = ++nver;
3237 if (indices[i4] == 0) indices[i4] = ++nver;
3238 }
3239 }
3240 // Construct polyhedron
3241 AllocateMemory(nver, nfac);
3242 G4ThreeVector p0(pmin.x() - 0.5*sizeX, pmin.y() - 0.5*sizeY, pmin.z() - 0.5*sizeZ);
3243 for (G4int ix = 0; ix <= nx; ++ix)
3244 {
3245 for (G4int iy = 0; iy <= ny; ++iy)
3246 {
3247 for (G4int iz = 0; iz <= nz; ++iz)
3248 {
3249 G4int i = ix*kvx + iy*kvy + iz;
3250 if (indices[i] == 0) continue;
3251 SetVertex(indices[i], p0 + G4ThreeVector(ix*sizeX, iy*sizeY, iz*sizeZ));
3252 }
3253 }
3254 }
3255 nfac = 0;
3256 for (const auto& p: positions)
3257 {
3258 G4int ix = (p.x() - pmin.x())*invx + 0.5;
3259 G4int iy = (p.y() - pmin.y())*invy + 0.5;
3260 G4int iz = (p.z() - pmin.z())*invz + 0.5;
3261 G4int vcheck = 0;
3262 // check (ix - 1) side
3263 vcheck = (ix == 0) ? 0 : voxels[(ix-1)*kx + iy*ky + iz];
3264 if (vcheck == 0)
3265 {
3266 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3267 G4int i2 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3268 G4int i3 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3269 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3270 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3271 }
3272 // check (ix + 1) side
3273 vcheck = (ix == nx - 1) ? 0 : voxels[(ix+1)*kx + iy*ky + iz];
3274 if (vcheck == 0)
3275 {
3276 G4int i1 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3277 G4int i2 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3278 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3279 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3280 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3281
3282 }
3283 // check (iy - 1) side
3284 vcheck = (iy == 0) ? 0 : voxels[ix*kx + (iy-1)*ky + iz];
3285 if (vcheck == 0)
3286 {
3287 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3288 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3289 G4int i3 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3290 G4int i4 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3291 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3292 }
3293 // check (iy + 1) side
3294 vcheck = (iy == ny - 1) ? 0 : voxels[ix*kx + (iy+1)*ky + iz];
3295 if (vcheck == 0)
3296 {
3297 G4int i1 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3298 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3299 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3300 G4int i4 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3301 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3302 }
3303 // check (iz - 1) side
3304 vcheck = (iz == 0) ? 0 : voxels[ix*kx + iy*ky + (iz-1)];
3305 if (vcheck == 0)
3306 {
3307 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+0); // 000
3308 G4int i2 = (ix+0)*kvx + (iy+1)*kvy + (iz+0); // 010
3309 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+0); // 110
3310 G4int i4 = (ix+1)*kvx + (iy+0)*kvy + (iz+0); // 100
3311 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3312 }
3313 // check (iz + 1) side
3314 vcheck = (iz == nz - 1) ? 0 : voxels[ix*kx + iy*ky + (iz+1)];
3315 if (vcheck == 0)
3316 {
3317 G4int i1 = (ix+0)*kvx + (iy+0)*kvy + (iz+1); // 001
3318 G4int i2 = (ix+1)*kvx + (iy+0)*kvy + (iz+1); // 101
3319 G4int i3 = (ix+1)*kvx + (iy+1)*kvy + (iz+1); // 111
3320 G4int i4 = (ix+0)*kvx + (iy+1)*kvy + (iz+1); // 011
3321 SetFacet(++nfac, indices[i1], indices[i2], indices[i3], indices[i4]);
3322 }
3323 }
3324 SetReferences();
3325}
3326
3327HepPolyhedronBoxMesh::~HepPolyhedronBoxMesh() = default;
3328
3329G4ThreadLocal
3330G4int HepPolyhedron::fNumberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
3331/***********************************************************************
3332 * *
3333 * Name: HepPolyhedron::fNumberOfRotationSteps Date: 24.06.97 *
3334 * Author: J.Allison (Manchester University) Revised: *
3335 * *
3336 * Function: Number of steps for whole circle *
3337 * *
3338 ***********************************************************************/
3339
3340#include "BooleanProcessor.src"
3341
3342HepPolyhedron HepPolyhedron::add(const HepPolyhedron & p) const
3343/***********************************************************************
3344 * *
3345 * Name: HepPolyhedron::add Date: 19.03.00 *
3346 * Author: E.Chernyaev Revised: *
3347 * *
3348 * Function: Boolean "union" of two polyhedra *
3349 * *
3350 ***********************************************************************/
3351{
3352 G4int ierr;
3353 BooleanProcessor processor;
3354 return processor.execute(OP_UNION, *this, p,ierr);
3355}
3356
3357HepPolyhedron HepPolyhedron::intersect(const HepPolyhedron & p) const
3358/***********************************************************************
3359 * *
3360 * Name: HepPolyhedron::intersect Date: 19.03.00 *
3361 * Author: E.Chernyaev Revised: *
3362 * *
3363 * Function: Boolean "intersection" of two polyhedra *
3364 * *
3365 ***********************************************************************/
3366{
3367 G4int ierr;
3368 BooleanProcessor processor;
3369 return processor.execute(OP_INTERSECTION, *this, p,ierr);
3370}
3371
3372HepPolyhedron HepPolyhedron::subtract(const HepPolyhedron & p) const
3373/***********************************************************************
3374 * *
3375 * Name: HepPolyhedron::add Date: 19.03.00 *
3376 * Author: E.Chernyaev Revised: *
3377 * *
3378 * Function: Boolean "subtraction" of "p" from "this" *
3379 * *
3380 ***********************************************************************/
3381{
3382 G4int ierr;
3383 BooleanProcessor processor;
3384 return processor.execute(OP_SUBTRACTION, *this, p,ierr);
3385}
3386
3387//NOTE : include the code of HepPolyhedronProcessor here
3388// since there is no BooleanProcessor.h
3389
3390#undef INTERSECTION
3391
3392#include "HepPolyhedronProcessor.src"
kCarTolerance
const G4double kCarTolerance
Definition: G4BoundingEnvelope.cc:38
B
G4double B(G4double temperature)
Definition: G4DNAElectronHoleRecombination.cc:68
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
G4Normal3D
HepGeom::Normal3D< G4double > G4Normal3D
Definition: G4Normal3D.hh:34
G4Point3D
HepGeom::Point3D< G4double > G4Point3D
Definition: G4Point3D.hh:34
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4Vector3D
HepGeom::Vector3D< G4double > G4Vector3D
Definition: G4Vector3D.hh:34
A
const G4double A[17]
Definition: G4WaterStopping.cc:53
operator<<
std::ostream & operator<<(std::ostream &ostr, const G4Facet &facet)
Definition: HepPolyhedron.cc:98
spatialTolerance
const G4double spatialTolerance
Definition: HepPolyhedron.cc:88
DEFAULT_NUMBER_OF_STEPS
#define DEFAULT_NUMBER_OF_STEPS
Definition: HepPolyhedron.h:213
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::setY
void setY(double)
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::Hep3Vector::setX
void setX(double)
HepGeom::BasicVector3D::cross
BasicVector3D< T > cross(const BasicVector3D< T > &v) const
Definition: BasicVector3D.h:281
HepGeom::BasicVector3D::unit
BasicVector3D< T > unit() const
Definition: BasicVector3D.h:310
HepGeom::BasicVector3D::set
void set(T x1, T y1, T z1)
Definition: BasicVector3D.h:167
HepGeom::BasicVector3D::dot
T dot(const BasicVector3D< T > &v) const
Definition: BasicVector3D.h:275
HepPolyhedronBoxMesh::~HepPolyhedronBoxMesh
~HepPolyhedronBoxMesh() override
HepPolyhedronBoxMesh::HepPolyhedronBoxMesh
HepPolyhedronBoxMesh(G4double sizeX, G4double sizeY, G4double sizeZ, const std::vector< G4ThreeVector > &positions)
Definition: HepPolyhedron.cc:3084
HepPolyhedronBox::HepPolyhedronBox
HepPolyhedronBox(G4double Dx, G4double Dy, G4double Dz)
Definition: HepPolyhedron.cc:1972
HepPolyhedronBox::~HepPolyhedronBox
~HepPolyhedronBox() override
HepPolyhedronCone::~HepPolyhedronCone
~HepPolyhedronCone() override
HepPolyhedronCone::HepPolyhedronCone
HepPolyhedronCone(G4double Rmn1, G4double Rmx1, G4double Rmn2, G4double Rmx2, G4double Dz)
Definition: HepPolyhedron.cc:2301
HepPolyhedronCons::~HepPolyhedronCons
~HepPolyhedronCons() override
HepPolyhedronCons::HepPolyhedronCons
HepPolyhedronCons(G4double Rmn1, G4double Rmx1, G4double Rmn2, G4double Rmx2, G4double Dz, G4double Phi1, G4double Dphi)
Definition: HepPolyhedron.cc:2223
HepPolyhedronEllipsoid::~HepPolyhedronEllipsoid
~HepPolyhedronEllipsoid() override
HepPolyhedronEllipsoid::HepPolyhedronEllipsoid
HepPolyhedronEllipsoid(G4double dx, G4double dy, G4double dz, G4double zcut1, G4double zcut2)
Definition: HepPolyhedron.cc:2688
HepPolyhedronEllipticalCone::HepPolyhedronEllipticalCone
HepPolyhedronEllipticalCone(G4double dx, G4double dy, G4double z, G4double zcut1)
Definition: HepPolyhedron.cc:2775
HepPolyhedronEllipticalCone::~HepPolyhedronEllipticalCone
~HepPolyhedronEllipticalCone() override
HepPolyhedronHype::HepPolyhedronHype
HepPolyhedronHype(G4double r1, G4double r2, G4double tan1, G4double tan2, G4double halfZ)
Definition: HepPolyhedron.cc:2143
HepPolyhedronHype::~HepPolyhedronHype
~HepPolyhedronHype() override
HepPolyhedronHyperbolicMirror::~HepPolyhedronHyperbolicMirror
~HepPolyhedronHyperbolicMirror() override
HepPolyhedronHyperbolicMirror::HepPolyhedronHyperbolicMirror
HepPolyhedronHyperbolicMirror(G4double a, G4double h, G4double r)
Definition: HepPolyhedron.cc:2839
HepPolyhedronPara::~HepPolyhedronPara
~HepPolyhedronPara() override
HepPolyhedronPara::HepPolyhedronPara
HepPolyhedronPara(G4double Dx, G4double Dy, G4double Dz, G4double Alpha, G4double Theta, G4double Phi)
Definition: HepPolyhedron.cc:2032
HepPolyhedronParaboloid::~HepPolyhedronParaboloid
~HepPolyhedronParaboloid() override
HepPolyhedronParaboloid::HepPolyhedronParaboloid
HepPolyhedronParaboloid(G4double r1, G4double r2, G4double dz, G4double Phi1, G4double Dphi)
Definition: HepPolyhedron.cc:2039
HepPolyhedronPcon::~HepPolyhedronPcon
~HepPolyhedronPcon() override
HepPolyhedronPcon::HepPolyhedronPcon
HepPolyhedronPcon(G4double phi, G4double dphi, G4int nz, const G4double *z, const G4double *rmin, const G4double *rmax)
Definition: HepPolyhedron.cc:2464
HepPolyhedronPgon::~HepPolyhedronPgon
~HepPolyhedronPgon() override
HepPolyhedronPgon::HepPolyhedronPgon
HepPolyhedronPgon(G4double phi, G4double dphi, G4int npdv, G4int nz, const G4double *z, const G4double *rmin, const G4double *rmax)
Definition: HepPolyhedron.cc:2321
HepPolyhedronSphere::~HepPolyhedronSphere
~HepPolyhedronSphere() override
HepPolyhedronSphere::HepPolyhedronSphere
HepPolyhedronSphere(G4double rmin, G4double rmax, G4double phi, G4double dphi, G4double the, G4double dthe)
Definition: HepPolyhedron.cc:2476
HepPolyhedronTetMesh::~HepPolyhedronTetMesh
~HepPolyhedronTetMesh() override
HepPolyhedronTetMesh::HepPolyhedronTetMesh
HepPolyhedronTetMesh(const std::vector< G4ThreeVector > &tetrahedra)
Definition: HepPolyhedron.cc:2902
HepPolyhedronTet::HepPolyhedronTet
HepPolyhedronTet(const G4double p0[3], const G4double p1[3], const G4double p2[3], const G4double p3[3])
Definition: HepPolyhedron.cc:2648
HepPolyhedronTet::~HepPolyhedronTet
~HepPolyhedronTet() override
HepPolyhedronTorus::HepPolyhedronTorus
HepPolyhedronTorus(G4double rmin, G4double rmax, G4double rtor, G4double phi, G4double dphi)
Definition: HepPolyhedron.cc:2573
HepPolyhedronTorus::~HepPolyhedronTorus
~HepPolyhedronTorus() override
HepPolyhedronTrap::~HepPolyhedronTrap
~HepPolyhedronTrap() override
HepPolyhedronTrap::HepPolyhedronTrap
HepPolyhedronTrap(G4double Dz, G4double Theta, G4double Phi, G4double Dy1, G4double Dx1, G4double Dx2, G4double Alp1, G4double Dy2, G4double Dx3, G4double Dx4, G4double Alp2)
Definition: HepPolyhedron.cc:1977
HepPolyhedronTrd1::HepPolyhedronTrd1
HepPolyhedronTrd1(G4double Dx1, G4double Dx2, G4double Dy, G4double Dz)
Definition: HepPolyhedron.cc:1966
HepPolyhedronTrd1::~HepPolyhedronTrd1
~HepPolyhedronTrd1() override
HepPolyhedronTrd2::HepPolyhedronTrd2
HepPolyhedronTrd2(G4double Dx1, G4double Dx2, G4double Dy1, G4double Dy2, G4double Dz)
Definition: HepPolyhedron.cc:1932
HepPolyhedronTrd2::~HepPolyhedronTrd2
~HepPolyhedronTrd2() override
HepPolyhedronTube::~HepPolyhedronTube
~HepPolyhedronTube() override
HepPolyhedronTube::HepPolyhedronTube
HepPolyhedronTube(G4double Rmin, G4double Rmax, G4double Dz)
Definition: HepPolyhedron.cc:2315
HepPolyhedronTubs::HepPolyhedronTubs
HepPolyhedronTubs(G4double Rmin, G4double Rmax, G4double Dz, G4double Phi1, G4double Dphi)
Definition: HepPolyhedron.cc:2308
HepPolyhedronTubs::~HepPolyhedronTubs
~HepPolyhedronTubs() override
HepPolyhedron::SetNumberOfRotationSteps
static void SetNumberOfRotationSteps(G4int n)
Definition: HepPolyhedron.cc:339
HepPolyhedron::RotateAroundZ
void RotateAroundZ(G4int nstep, G4double phi, G4double dphi, G4int np1, G4int np2, const G4double *z, G4double *r, G4int nodeVis, G4int edgeVis)
Definition: HepPolyhedron.cc:553
HepPolyhedron::GetNextEdgeIndices
G4bool GetNextEdgeIndices(G4int &i1, G4int &i2, G4int &edgeFlag, G4int &iface1, G4int &iface2) const
Definition: HepPolyhedron.cc:1500
HepPolyhedron::GetUnitNormal
G4Normal3D GetUnitNormal(G4int iFace) const
Definition: HepPolyhedron.cc:1721
HepPolyhedron::createTwistedTrap
G4int createTwistedTrap(G4double Dz, const G4double xy1[][2], const G4double xy2[][2])
Definition: HepPolyhedron.cc:1833
HepPolyhedron::pV
G4Point3D * pV
Definition: HepPolyhedron.h:237
HepPolyhedron::GetNextNormal
G4bool GetNextNormal(G4Normal3D &normal) const
Definition: HepPolyhedron.cc:1746
HepPolyhedron::operator=
HepPolyhedron & operator=(const HepPolyhedron &from)
Definition: HepPolyhedron.cc:168
HepPolyhedron::FindNodeNormal
G4Normal3D FindNodeNormal(G4int iFace, G4int iNode) const
Definition: HepPolyhedron.cc:242
HepPolyhedron::GetNumberOfRotationSteps
static G4int GetNumberOfRotationSteps()
Definition: HepPolyhedron.cc:269
HepPolyhedron::GetNormal
G4Normal3D GetNormal(G4int iFace) const
Definition: HepPolyhedron.cc:1696
HepPolyhedron::FindNeighbour
G4int FindNeighbour(G4int iFace, G4int iNode, G4int iOrder) const
Definition: HepPolyhedron.cc:214
HepPolyhedron::TriangulatePolygon
G4bool TriangulatePolygon(const std::vector< G4TwoVector > &polygon, std::vector< G4int > &result)
Definition: HepPolyhedron.cc:985
HepPolyhedron::SetVertex
void SetVertex(G4int index, const G4Point3D &v)
Definition: HepPolyhedron.cc:282
HepPolyhedron::RotateContourAroundZ
void RotateContourAroundZ(G4int nstep, G4double phi, G4double dphi, const std::vector< G4TwoVector > &rz, G4int nodeVis, G4int edgeVis)
Definition: HepPolyhedron.cc:828
HepPolyhedron::SetFacet
void SetFacet(G4int index, G4int iv1, G4int iv2, G4int iv3, G4int iv4=0)
Definition: HepPolyhedron.cc:305
HepPolyhedron::GetVertex
G4Point3D GetVertex(G4int index) const
Definition: HepPolyhedron.cc:1431
HepPolyhedron::GetSurfaceArea
G4double GetSurfaceArea() const
Definition: HepPolyhedron.cc:1782
HepPolyhedron::subtract
HepPolyhedron subtract(const HepPolyhedron &p) const
Definition: HepPolyhedron.cc:3372
HepPolyhedron::GetNextVertex
G4bool GetNextVertex(G4Point3D &vertex, G4int &edgeFlag) const
Definition: HepPolyhedron.cc:1451
HepPolyhedron::GetFacet
void GetFacet(G4int iFace, G4int &n, G4int *iNodes, G4int *edgeFlags=nullptr, G4int *iFaces=nullptr) const
Definition: HepPolyhedron.cc:1610
HepPolyhedron::GetVolume
G4double GetVolume() const
Definition: HepPolyhedron.cc:1804
HepPolyhedron::intersect
HepPolyhedron intersect(const HepPolyhedron &p) const
Definition: HepPolyhedron.cc:3357
HepPolyhedron::Transform
HepPolyhedron & Transform(const G4Transform3D &t)
Definition: HepPolyhedron.cc:1379
HepPolyhedron::createPolyhedron
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
Definition: HepPolyhedron.cc:1902
HepPolyhedron::AllocateMemory
void AllocateMemory(G4int Nvert, G4int Nface)
Definition: HepPolyhedron.cc:374
HepPolyhedron::GetNextUnitNormal
G4bool GetNextUnitNormal(G4Normal3D &normal) const
Definition: HepPolyhedron.cc:1766
HepPolyhedron::ResetNumberOfRotationSteps
static void ResetNumberOfRotationSteps()
Definition: HepPolyhedron.cc:361
HepPolyhedron::CheckSnip
G4bool CheckSnip(const std::vector< G4TwoVector > &contour, G4int a, G4int b, G4int c, G4int n, const G4int *V)
Definition: HepPolyhedron.cc:1063
HepPolyhedron::RotateEdge
void RotateEdge(G4int k1, G4int k2, G4double r1, G4double r2, G4int v1, G4int v2, G4int vEdge, G4bool ifWholeCircle, G4int ns, G4int &kface)
Definition: HepPolyhedron.cc:420
HepPolyhedron::fNumberOfRotationSteps
static G4ThreadLocal G4int fNumberOfRotationSteps
Definition: HepPolyhedron.h:235
HepPolyhedron::GetNextVertexIndex
G4bool GetNextVertexIndex(G4int &index, G4int &edgeFlag) const
Definition: HepPolyhedron.cc:1404
HepPolyhedron::pF
G4Facet * pF
Definition: HepPolyhedron.h:238
HepPolyhedron::GetNextFacet
G4bool GetNextFacet(G4int &n, G4Point3D *nodes, G4int *edgeFlags=nullptr, G4Normal3D *normals=nullptr) const
Definition: HepPolyhedron.cc:1666
HepPolyhedron::JoinCoplanarFacets
void JoinCoplanarFacets(G4double tolerance)
Definition: HepPolyhedron.cc:1264
HepPolyhedron::add
HepPolyhedron add(const HepPolyhedron &p) const
Definition: HepPolyhedron.cc:3342
HepPolyhedron::GetNextEdge
G4bool GetNextEdge(G4Point3D &p1, G4Point3D &p2, G4int &edgeFlag) const
Definition: HepPolyhedron.cc:1569
HepPolyhedron::SetSideFacets
void SetSideFacets(G4int ii[4], G4int vv[4], G4int *kk, G4double *r, G4double dphi, G4int ns, G4int &kface)
Definition: HepPolyhedron.cc:482
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:77
ns
#define ns(x)
Definition: xmltok.c:1649