ATLAS Offline Software
Loading...
Searching...
No Matches
SbPolyhedron.cxx
Go to the documentation of this file.
1
2// //
3// Source file for class SbPolyhedron //
4// //
5// Update:
6// Riccardo-Maria BIANCHI (rbianchi@cern.ch) //
7// 13.12.2012 //
8// //
10
11
12// this :
13#include <VP1HEPVis/SbPolyhedron.h>
14#include <cassert>
15#include <memory>
16
17#define perMillion 0.000001
18#define deg (M_PI/180.0)
19
20#include <float.h> //G.Barrand : to have DBL_EPSILON on Windows.
21
22// G.Barrand : introduce iabs to avoid a mess with cmath and some compiler.
23inline int iabs(int a) {
24 return a < 0 ? -a : a;
25}
26inline float ffabs(float a) {
27 return a < 0.0f ? -a : a;
28}
29
30
31// rbianchi
32HVPoint3D::HVPoint3D(): SbVec3d(0,0,0){}
33HVPoint3D::HVPoint3D(double x,double y,double z): SbVec3d(x,y,z){}
34HVPoint3D::HVPoint3D(const HVPoint3D& v): SbVec3d(v){}
35HVPoint3D::HVPoint3D(const SbVec3d& v): SbVec3d(v){}
36HVPoint3D& HVPoint3D::operator=(const HVPoint3D& v) {
37 setValue(v[0],v[1],v[2]);
38 return *this;
39}
40HVPoint3D& HVPoint3D::operator=(const SbVec3d& v) {
41 setValue(v[0],v[1],v[2]);
42 return *this;
43}
44HVPoint3D operator +(const HVPoint3D& v1,const HVPoint3D& v2) {
45 return HVPoint3D(v1[0] + v2[0],v1[1] + v2[1],v1[2] + v2[2]);
46}
47//---
48
49
50
51//--------------------------------------------------------------------//
52// JFB: //
53// SbPolyhedron was SbPolyhedron, retrofitted to Open Inventor //
54// infrastructure: //
55//--------------------------------------------------------------------//
56
57
58//
59// ********************************************************************
60// * DISCLAIMER *
61// * *
62// * The following disclaimer summarizes all the specific disclaimers *
63// * of contributors to this software. The specific disclaimers,which *
64// * govern, are listed with their locations in: *
65// * http://cern.ch/geant4/license *
66// * *
67// * Neither the authors of this software system, nor their employing *
68// * institutes,nor the agencies providing financial support for this *
69// * work make any representation or warranty, express or implied, *
70// * regarding this software system or assume any liability for its *
71// * use. *
72// * *
73// * This code implementation is the intellectual property of the *
74// * GEANT4 collaboration. *
75// * By copying, distributing or modifying the Program (or any work *
76// * based on the Program) you indicate your acceptance of this *
77// * statement, and all its terms. *
78// ********************************************************************
79//
80//
81// $Id: SbPolyhedron.cxx,v 1.3 2008-09-14 19:04:40 tkittel Exp $
82// GEANT4 tag $Name: not supported by cvs2svn $
83//
84//
85//
86// G4 Polyhedron library
87//
88// History:
89// 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@cern.ch> - initial version
90//
91// 30.09.96 E.Chernyaev
92// - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
93// - added GetNextUnitNormal, GetNextEdgeIndeces, GetNextEdge
94//
95// 15.12.96 E.Chernyaev
96// - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
97// - rewritten G4PolyhedronCons;
98// - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
99//
100// 01.06.97 E.Chernyaev
101// - modified RotateAroundZ, added SetSideFacets
102//
103// 19.03.00 E.Chernyaev
104// - implemented boolean operations (add, subtract, intersect) on polyhedra;
105//
106// 25.05.01 E.Chernyaev
107// - added GetSurfaceArea() and GetVolume();
108//
109
110
111/***********************************************************************
112 * *
113 * Name: SbPolyhedron operator << Date: 09.05.96 *
114 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
115 * *
116 * Function: Print contents of G4 polyhedron *
117 * *
118 ***********************************************************************/
119std::ostream & operator<<(std::ostream & ostr, const SbFacet & facet) {
120 for (int k=0; k<4; k++) {
121 ostr << " " << facet.m_edge[k].v << "/" << facet.m_edge[k].f;
122 }
123 return ostr;
124}
125
126std::ostream & operator<<(std::ostream & ostr, const SbPolyhedron & ph) {
127 ostr << std::endl;
128 ostr << "Nverteces=" << ph.m_nvert << ", Nfacets=" << ph.m_nface << std::endl;
129 int i;
130 for (i=1; i<=ph.m_nvert; i++) {
131 ostr << "xyz(" << i << ")="
132 << ph.m_pV[i][0] << ' ' << ph.m_pV[i][1] << ' ' << ph.m_pV[i][2]
133 << std::endl;
134 }
135 for (i=1; i<=ph.m_nface; i++) {
136 ostr << "face(" << i << ")=" << ph.m_pF[i] << std::endl;
137 }
138 return ostr;
139}
140
141SbPolyhedron::SbPolyhedron(const SbPolyhedron &from)
142/***********************************************************************
143 * *
144 * Name: SbPolyhedron copy constructor Date: 23.07.96 *
145 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
146 * *
147 ***********************************************************************/
148{
149 if (from.m_nvert > 0 && from.m_nface > 0) {
150 m_nvert = from.m_nvert;
151 m_nface = from.m_nface;
152 m_pV = new HVPoint3D[m_nvert + 1];
153 m_pF = new SbFacet[m_nface + 1];
154 int i;
155 for (i=1; i<=m_nvert; i++) m_pV[i] = from.m_pV[i];
156 for (i=1; i<=m_nface; i++) m_pF[i] = from.m_pF[i];
157 }else{
158 m_nvert = 0; m_nface = 0; m_pV = 0; m_pF = 0;
159 }
160}
161
162SbPolyhedron & SbPolyhedron::operator=(const SbPolyhedron &from)
163/***********************************************************************
164 * *
165 * Name: SbPolyhedron operator = Date: 23.07.96 *
166 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
167 * *
168 * Function: Copy contents of one GEANT4 polyhedron to another *
169 * *
170 ***********************************************************************/
171{
172 if (this == &from) return *this;
173 delete [] m_pV;
174 delete [] m_pF;
175 if (from.m_nvert > 0 && from.m_nface > 0) {
176 m_nvert = from.m_nvert;
177 m_nface = from.m_nface;
178 m_pV = new HVPoint3D[m_nvert + 1];
179 m_pF = new SbFacet[m_nface + 1];
180 int i;
181 for (i=1; i<=m_nvert; i++) m_pV[i] = from.m_pV[i];
182 for (i=1; i<=m_nface; i++) m_pF[i] = from.m_pF[i];
183 }else{
184 m_nvert = 0; m_nface = 0; m_pV = 0; m_pF = 0;
185 }
186 return *this;
187}
188
189int
190SbPolyhedron::FindNeighbour(int iFace, int iNode, int iOrder) const
191/***********************************************************************
192 * *
193 * Name: SbPolyhedron::FindNeighbour Date: 22.11.99 *
194 * Author: E.Chernyaev Revised: *
195 * *
196 * Function: Find neighbouring face *
197 * *
198 ***********************************************************************/
199{
200 int i;
201 for (i=0; i<4; i++) {
202 if (iNode == iabs(m_pF[iFace].m_edge[i].v)) break;
203 }
204 if (i == 4) {
205 std::cerr
206 << "SbPolyhedron::FindNeighbour: face " << iFace
207 << " has no node " << iNode
208 << std::endl;
209 return 0;
210 }
211 if (iOrder < 0) {
212 if ( --i < 0) i = 3;
213 if (m_pF[iFace].m_edge[i].v == 0) i = 2;
214 }
215 return (m_pF[iFace].m_edge[i].v > 0) ? 0 : m_pF[iFace].m_edge[i].f;
216}
217
218HVNormal3D SbPolyhedron::FindNodeNormal(int iFace, int iNode) const
219/***********************************************************************
220 * *
221 * Name: SbPolyhedron::FindNodeNormal Date: 22.11.99 *
222 * Author: E.Chernyaev Revised: *
223 * *
224 * Function: Find normal at given node *
225 * *
226 ***********************************************************************/
227{
228 HVNormal3D normal = GetUnitNormal(iFace);
229 int k = iFace, iOrder = 1;
230
231 for(;;) {
232 k = FindNeighbour(k, iNode, iOrder);
233 if (k == iFace) break;
234 if (k > 0) {
235 normal += GetUnitNormal(k);
236 }else{
237 if (iOrder < 0) break;
238 k = iFace;
239 iOrder = -iOrder;
240 }
241 }
242 normal.normalize();
243 return normal;
244}
245
246void SbPolyhedron::SetNumberOfRotationSteps(int n)
247/***********************************************************************
248 * *
249 * Name: SbPolyhedron::SetNumberOfRotationSteps Date: 24.06.97 *
250 * Author: J.Allison (Manchester University) Revised: *
251 * *
252 * Function: Set number of steps for whole circle *
253 * *
254 ***********************************************************************/
255{
256 const int nMin = 3;
257 if (n < nMin) {
258 std::cerr
259 << "SbPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
260 << "number of steps per circle < " << nMin << "; forced to " << nMin
261 << std::endl;
262 s_numberOfRotationSteps = nMin;
263 }else{
264 s_numberOfRotationSteps = n;
265 }
266}
267
268void SbPolyhedron::AllocateMemory(int Nvert, int Nface)
269/***********************************************************************
270 * *
271 * Name: SbPolyhedron::AllocateMemory Date: 19.06.96 *
272 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
273 * *
274 * Function: Allocate memory for GEANT4 polyhedron *
275 * *
276 * Input: Nvert - number of nodes *
277 * Nface - number of faces *
278 * *
279 ***********************************************************************/
280{
281 m_nvert = Nvert;
282 m_nface = Nface;
283 m_pV = new HVPoint3D[m_nvert+1];
284 m_pF = new SbFacet[m_nface+1];
285}
286
287void SbPolyhedron::CreatePrism()
288/***********************************************************************
289 * *
290 * Name: SbPolyhedron::CreatePrism Date: 15.07.96 *
291 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
292 * *
293 * Function: Set facets for a prism *
294 * *
295 ***********************************************************************/
296{
297 enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
298
299 m_pF[1] = SbFacet(1,LEFT, 4,BACK, 3,RIGHT, 2,FRONT);
300 m_pF[2] = SbFacet(5,TOP, 8,BACK, 4,BOTTOM, 1,FRONT);
301 m_pF[3] = SbFacet(8,TOP, 7,RIGHT, 3,BOTTOM, 4,LEFT);
302 m_pF[4] = SbFacet(7,TOP, 6,FRONT, 2,BOTTOM, 3,BACK);
303 m_pF[5] = SbFacet(6,TOP, 5,LEFT, 1,BOTTOM, 2,RIGHT);
304 m_pF[6] = SbFacet(5,FRONT, 6,RIGHT, 7,BACK, 8,LEFT);
305}
306
307void SbPolyhedron::RotateEdge(int k1, int k2, double r1, double r2,
308 int v1, int v2, int vEdge,
309 bool ifWholeCircle, int ns, int &kface)
310/***********************************************************************
311 * *
312 * Name: SbPolyhedron::RotateEdge Date: 05.12.96 *
313 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
314 * *
315 * Function: Create set of facets by rotation of an edge around Z-axis *
316 * *
317 * Input: k1, k2 - end vertices of the edge *
318 * r1, r2 - radiuses of the end vertices *
319 * v1, v2 - visibility of edges produced by rotation of the end *
320 * vertices *
321 * vEdge - visibility of the edge *
322 * ifWholeCircle - is true in case of whole circle rotation *
323 * ns - number of discrete steps *
324 * r[] - r-coordinates *
325 * kface - current free cell in the m_pF array *
326 * *
327 ***********************************************************************/
328{
329 if (r1 == 0. && r2 == 0) return;
330
331 int i;
332 int i1 = k1;
333 int i2 = k2;
334 int ii1 = ifWholeCircle ? i1 : i1+ns;
335 int ii2 = ifWholeCircle ? i2 : i2+ns;
336 int vv = ifWholeCircle ? vEdge : 1;
337
338 if (ns == 1) {
339 if (r1 == 0.) {
340 m_pF[kface++] = SbFacet(i1,0, v2*i2,0, (i2+1),0);
341 }else if (r2 == 0.) {
342 m_pF[kface++] = SbFacet(i1,0, i2,0, v1*(i1+1),0);
343 }else{
344 m_pF[kface++] = SbFacet(i1,0, v2*i2,0, (i2+1),0, v1*(i1+1),0);
345 }
346 }else{
347 if (r1 == 0.) {
348 m_pF[kface++] = SbFacet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0);
349 for (i2++,i=1; i<ns-1; i2++,i++) {
350 m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
351 }
352 m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
353 }else if (r2 == 0.) {
354 m_pF[kface++] = SbFacet(vv*i1,0, vEdge*i2,0, v1*(i1+1),0);
355 for (i1++,i=1; i<ns-1; i1++,i++) {
356 m_pF[kface++] = SbFacet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
357 }
358 m_pF[kface++] = SbFacet(vEdge*i1,0, vv*i2,0, v1*ii1,0);
359 }else{
360 m_pF[kface++] = SbFacet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
361 for (i1++,i2++,i=1; i<ns-1; i1++,i2++,i++) {
362 m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
363 }
364 m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vv*ii2,0, v1*ii1,0);
365 }
366 }
367}
368
369void SbPolyhedron::SetSideFacets(int ii[4], int vv[4],
370 int *kk, double *r,
371 double dphi, int ns, int &kface)
372/***********************************************************************
373 * *
374 * Name: SbPolyhedron::SetSideFacets Date: 20.05.97 *
375 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
376 * *
377 * Function: Set side facets for the case of incomplete rotation *
378 * *
379 * Input: ii[4] - indeces of original verteces *
380 * vv[4] - visibility of edges *
381 * kk[] - indeces of nodes *
382 * r[] - radiuses *
383 * dphi - delta phi *
384 * ns - number of discrete steps *
385 * kface - current free cell in the m_pF array *
386 * *
387 ***********************************************************************/
388{
389 int k1, k2, k3, k4;
390 if (fabs(dphi-M_PI) < perMillion) { // half a circle
391 for (int i=0; i<4; i++) {
392 k1 = ii[i];
393 k2 = (i == 3) ? ii[0] : ii[i+1];
394 if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;
395 }
396 }
397
398 if (ii[1] == ii[2]) {
399 k1 = kk[ii[0]];
400 k2 = kk[ii[2]];
401 k3 = kk[ii[3]];
402 m_pF[kface++] = SbFacet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
403 if (r[ii[0]] != 0.) k1 += ns;
404 if (r[ii[2]] != 0.) k2 += ns;
405 if (r[ii[3]] != 0.) k3 += ns;
406 m_pF[kface++] = SbFacet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
407 }else if (kk[ii[0]] == kk[ii[1]]) {
408 k1 = kk[ii[0]];
409 k2 = kk[ii[2]];
410 k3 = kk[ii[3]];
411 m_pF[kface++] = SbFacet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
412 if (r[ii[0]] != 0.) k1 += ns;
413 if (r[ii[2]] != 0.) k2 += ns;
414 if (r[ii[3]] != 0.) k3 += ns;
415 m_pF[kface++] = SbFacet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
416 }else if (kk[ii[2]] == kk[ii[3]]) {
417 k1 = kk[ii[0]];
418 k2 = kk[ii[1]];
419 k3 = kk[ii[2]];
420 m_pF[kface++] = SbFacet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
421 if (r[ii[0]] != 0.) k1 += ns;
422 if (r[ii[1]] != 0.) k2 += ns;
423 if (r[ii[2]] != 0.) k3 += ns;
424 m_pF[kface++] = SbFacet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
425 }else{
426 k1 = kk[ii[0]];
427 k2 = kk[ii[1]];
428 k3 = kk[ii[2]];
429 k4 = kk[ii[3]];
430 m_pF[kface++] = SbFacet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
431 if (r[ii[0]] != 0.) k1 += ns;
432 if (r[ii[1]] != 0.) k2 += ns;
433 if (r[ii[2]] != 0.) k3 += ns;
434 if (r[ii[3]] != 0.) k4 += ns;
435 m_pF[kface++] = SbFacet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
436 }
437}
438
439void SbPolyhedron::RotateAroundZ(int nstep, double phi, double dphi,
440 int np1, int np2,
441 const double *z, double *r,
442 int nodeVis, int edgeVis)
443/***********************************************************************
444 * *
445 * Name: SbPolyhedron::RotateAroundZ Date: 27.11.96 *
446 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
447 * *
448 * Function: Create SbPolyhedron for a solid produced by rotation of *
449 * two polylines around Z-axis *
450 * *
451 * Input: nstep - number of discrete steps, if 0 then default *
452 * phi - starting phi angle *
453 * dphi - delta phi *
454 * np1 - number of points in external polyline *
455 * (must be negative in case of closed polyline) *
456 * np2 - number of points in internal polyline (may be 1) *
457 * z[] - z-coordinates (+z >>> -z for both polylines) *
458 * r[] - r-coordinates *
459 * nodeVis - how to Draw edges joing consecutive positions of *
460 * node during rotation *
461 * edgeVis - how to Draw edges *
462 * *
463 ***********************************************************************/
464{
465 static double wholeCircle = 2*M_PI;
466
467 // S E T R O T A T I O N P A R A M E T E R S
468
469 bool ifWholeCircle = (fabs(dphi-wholeCircle) < perMillion) ?
470 true : false;
471 double delPhi = ifWholeCircle ? wholeCircle : dphi;
472 int nSphi = (nstep > 0) ?
473 nstep : int(delPhi*GetNumberOfRotationSteps()/wholeCircle+.5);
474 if (nSphi == 0) nSphi = 1;
475 int nVphi = ifWholeCircle ? nSphi : nSphi+1;
476 bool ifClosed = np1 > 0 ? false : true;
477
478 // C O U N T V E R T E C E S
479
480 int absNp1 = iabs(np1);
481 int absNp2 = iabs(np2);
482 int i1beg = 0;
483 int i1end = absNp1-1;
484 int i2beg = absNp1;
485 int i2end = absNp1+absNp2-1;
486 int i, j, k;
487
488 for(i=i1beg; i<=i2end; i++) {
489 if (fabs(r[i]) < perMillion) r[i] = 0.;
490 }
491
492 j = 0; // external nodes
493 for (i=i1beg; i<=i1end; i++) {
494 j += (r[i] == 0.) ? 1 : nVphi;
495 }
496
497 bool ifSide1 = false; // internal nodes
498 bool ifSide2 = false;
499
500 if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) {
501 j += (r[i2beg] == 0.) ? 1 : nVphi;
502 ifSide1 = true;
503 }
504
505 for(i=i2beg+1; i<i2end; i++) {
506 j += (r[i] == 0.) ? 1 : nVphi;
507 }
508
509 if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) {
510 if (absNp2 > 1) j += (r[i2end] == 0.) ? 1 : nVphi;
511 ifSide2 = true;
512 }
513
514 // C O U N T F A C E S
515
516 k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi; // external faces
517
518 if (absNp2 > 1) { // internal faces
519 for(i=i2beg; i<i2end; i++) {
520 if (r[i] > 0. || r[i+1] > 0.) k += nSphi;
521 }
522
523 if (ifClosed) {
524 if (r[i2end] > 0. || r[i2beg] > 0.) k += nSphi;
525 }
526 }
527
528 if (!ifClosed) { // side faces
529 if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) k += nSphi;
530 if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) k += nSphi;
531 }
532
533 if (!ifWholeCircle) { // phi_side faces
534 k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
535 }
536
537 // A L L O C A T E M E M O R Y
538
539 AllocateMemory(j, k);
540
541 // G E N E R A T E V E R T E C E S
542
543 int *kk;
544 kk = new int[absNp1+absNp2];
545
546 k = 1;
547 for(i=i1beg; i<=i1end; i++) {
548 kk[i] = k;
549 if (r[i] == 0.) { m_pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
550 }
551
552 i = i2beg;
553 if (ifSide1) {
554 kk[i] = k;
555 if (r[i] == 0.) { m_pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
556 }else{
557 kk[i] = kk[i1beg];
558 }
559
560 for(i=i2beg+1; i<i2end; i++) {
561 kk[i] = k;
562 if (r[i] == 0.) { m_pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
563 }
564
565 if (absNp2 > 1) {
566 i = i2end;
567 if (ifSide2) {
568 kk[i] = k;
569 if (r[i] == 0.) m_pV[k] = HVPoint3D(0, 0, z[i]);
570 }else{
571 kk[i] = kk[i1end];
572 }
573 }
574
575 double cosPhi, sinPhi;
576
577 for(j=0; j<nVphi; j++) {
578 cosPhi = cos(phi+j*delPhi/nSphi);
579 sinPhi = sin(phi+j*delPhi/nSphi);
580 for(i=i1beg; i<=i2end; i++) {
581 if (r[i] != 0.) m_pV[kk[i]+j] = HVPoint3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
582 }
583 }
584
585 // G E N E R A T E E X T E R N A L F A C E S
586
587 int v1,v2;
588
589 k = 1;
590 v2 = ifClosed ? nodeVis : 1;
591 for(i=i1beg; i<i1end; i++) {
592 v1 = v2;
593 if (!ifClosed && i == i1end-1) {
594 v2 = 1;
595 }else{
596 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
597 }
598 RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
599 edgeVis, ifWholeCircle, nSphi, k);
600 }
601 if (ifClosed) {
602 RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
603 edgeVis, ifWholeCircle, nSphi, k);
604 }
605
606 // G E N E R A T E I N T E R N A L F A C E S
607
608 if (absNp2 > 1) {
609 v2 = ifClosed ? nodeVis : 1;
610 for(i=i2beg; i<i2end; i++) {
611 v1 = v2;
612 if (!ifClosed && i==i2end-1) {
613 v2 = 1;
614 }else{
615 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
616 }
617 RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
618 edgeVis, ifWholeCircle, nSphi, k);
619 }
620 if (ifClosed) {
621 RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
622 edgeVis, ifWholeCircle, nSphi, k);
623 }
624 }
625
626 // G E N E R A T E S I D E F A C E S
627
628 if (!ifClosed) {
629 if (ifSide1) {
630 RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
631 -1, ifWholeCircle, nSphi, k);
632 }
633 if (ifSide2) {
634 RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
635 -1, ifWholeCircle, nSphi, k);
636 }
637 }
638
639 // G E N E R A T E S I D E F A C E S for the case of incomplete circle
640
641 if (!ifWholeCircle) {
642
643 int ii[4], vv[4];
644
645 if (ifClosed) {
646 for (i=i1beg; i<=i1end; i++) {
647 ii[0] = i;
648 ii[3] = (i == i1end) ? i1beg : i+1;
649 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
650 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
651 vv[0] = -1;
652 vv[1] = 1;
653 vv[2] = -1;
654 vv[3] = 1;
655 SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
656 }
657 }else{
658 for (i=i1beg; i<i1end; i++) {
659 ii[0] = i;
660 ii[3] = i+1;
661 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
662 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
663 vv[0] = (i == i1beg) ? 1 : -1;
664 vv[1] = 1;
665 vv[2] = (i == i1end-1) ? 1 : -1;
666 vv[3] = 1;
667 SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
668 }
669 }
670 }
671
672 delete [] kk;
673
674 if (k-1 != m_nface) {
675 std::cerr
676 << "Polyhedron::RotateAroundZ: number of generated faces ("
677 << k-1 << ") is not equal to the number of allocated faces ("
678 << m_nface << ")"
679 << std::endl;
680 }
681}
682
683void SbPolyhedron::SetReferences()
684/***********************************************************************
685 * *
686 * Name: SbPolyhedron::SetReferences Date: 04.12.96 *
687 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
688 * *
689 * Function: For each edge set reference to neighbouring facet *
690 * *
691 ***********************************************************************/
692{
693 if (m_nface <= 0) return;
694
695 struct edgeListMember {
696 edgeListMember *next;
697 int v2;
698 int iface;
699 int iedge;
700 } *edgeList, *freeList, **headList;
701
702
703 // A L L O C A T E A N D I N I T I A T E L I S T S
704
705 edgeList = new edgeListMember[2*m_nface];
706 headList = new edgeListMember*[m_nvert];
707
708 int i;
709 for (i=0; i<m_nvert; i++) {
710 headList[i] = 0;
711 }
712 freeList = edgeList;
713 for (i=0; i<2*m_nface-1; i++) {
714 edgeList[i].next = &edgeList[i+1];
715 }
716 edgeList[2*m_nface-1].next = 0;
717
718 // L O O P A L O N G E D G E S
719
720 int iface, iedge, nedge, i1, i2, k1, k2;
721 edgeListMember *prev, *cur;
722
723 for(iface=1; iface<=m_nface; iface++) {
724 nedge = (m_pF[iface].m_edge[3].v == 0) ? 3 : 4;
725 for (iedge=0; iedge<nedge; iedge++) {
726 i1 = iedge;
727 i2 = (iedge < nedge-1) ? iedge+1 : 0;
728 i1 = iabs(m_pF[iface].m_edge[i1].v);
729 i2 = iabs(m_pF[iface].m_edge[i2].v);
730 k1 = (i1 < i2) ? i1 : i2; // k1 = ::min(i1,i2);
731 k2 = (i1 > i2) ? i1 : i2; // k2 = ::max(i1,i2);
732
733 // check head of the List corresponding to k1
734 cur = headList[k1];
735 if (cur == 0) {
736 headList[k1] = freeList;
737 freeList = freeList->next;
738 cur = headList[k1];
739 cur->next = 0;
740 cur->v2 = k2;
741 cur->iface = iface;
742 cur->iedge = iedge;
743 continue;
744 }
745
746 if (cur->v2 == k2) {
747 headList[k1] = cur->next;
748 cur->next = freeList;
749 freeList = cur;
750 m_pF[iface].m_edge[iedge].f = cur->iface;
751 m_pF[cur->iface].m_edge[cur->iedge].f = iface;
752 i1 = (m_pF[iface].m_edge[iedge].v < 0) ? -1 : 1;
753 i2 = (m_pF[cur->iface].m_edge[cur->iedge].v < 0) ? -1 : 1;
754 if (i1 != i2) {
755 std::cerr
756 << "Polyhedron::SetReferences: different edge visibility "
757 << iface << "/" << iedge << "/"
758 << m_pF[iface].m_edge[iedge].v << " and "
759 << cur->iface << "/" << cur->iedge << "/"
760 << m_pF[cur->iface].m_edge[cur->iedge].v
761 << std::endl;
762 }
763 continue;
764 }
765
766 // check List itself
767 for (;;) {
768 prev = cur;
769 cur = prev->next;
770 if (cur == 0) {
771 prev->next = freeList;
772 freeList = freeList->next;
773 cur = prev->next;
774 cur->next = 0;
775 cur->v2 = k2;
776 cur->iface = iface;
777 cur->iedge = iedge;
778 break;
779 }
780
781 if (cur->v2 == k2) {
782 prev->next = cur->next;
783 cur->next = freeList;
784 freeList = cur;
785 m_pF[iface].m_edge[iedge].f = cur->iface;
786 m_pF[cur->iface].m_edge[cur->iedge].f = iface;
787 i1 = (m_pF[iface].m_edge[iedge].v < 0) ? -1 : 1;
788 i2 = (m_pF[cur->iface].m_edge[cur->iedge].v < 0) ? -1 : 1;
789 if (i1 != i2) {
790 std::cerr
791 << "Polyhedron::SetReferences: different edge visibility "
792 << iface << "/" << iedge << "/"
793 << m_pF[iface].m_edge[iedge].v << " and "
794 << cur->iface << "/" << cur->iedge << "/"
795 << m_pF[cur->iface].m_edge[cur->iedge].v
796 << std::endl;
797 }
798 break;
799 }
800 }
801 }
802 }
803
804 // 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
805
806 for (i=0; i<m_nvert; i++) {
807 if (headList[i] != 0) {
808 std::cerr
809 << "Polyhedron::SetReferences: List " << i << " is not empty"
810 << std::endl;
811 }
812 }
813
814 // F R E E M E M O R Y
815
816 delete [] edgeList;
817 delete [] headList;
818}
819
820void SbPolyhedron::InvertFacets()
821/***********************************************************************
822 * *
823 * Name: SbPolyhedron::InvertFacets Date: 01.12.99 *
824 * Author: E.Chernyaev Revised: *
825 * *
826 * Function: Invert the order of the nodes in the facets *
827 * *
828 ***********************************************************************/
829{
830 if (m_nface <= 0) return;
831 int i, k, nnode, v[4],f[4];
832 for (i=1; i<=m_nface; i++) {
833 nnode = (m_pF[i].m_edge[3].v == 0) ? 3 : 4;
834 for (k=0; k<nnode; k++) {
835 v[k] = (k+1 == nnode) ? m_pF[i].m_edge[0].v : m_pF[i].m_edge[k+1].v;
836 if (v[k] * m_pF[i].m_edge[k].v < 0) v[k] = -v[k];
837 f[k] = m_pF[i].m_edge[k].f;
838 }
839 for (k=0; k<nnode; k++) {
840 m_pF[i].m_edge[nnode-1-k].v = v[k];
841 m_pF[i].m_edge[nnode-1-k].f = f[k];
842 }
843 }
844}
845
846
847// rbianchi change
848//- original transform method
849/*
850SbPolyhedron & SbPolyhedron::Transform(const HVRotation & rotation, const HVVector3D & translation)
851// ***********************************************************************
852// * *
853// * Name: SbPolyhedron::Transform Date: 01.12.99 *
854// * Author: E.Chernyaev Revised: *
855// * *
856// * Function: Make transformation of the polyhedron *
857// * *
858// ***********************************************************************
859{
860 if (m_nvert > 0) {
861 for (int i=1; i<=m_nvert; i++) {
862 HVVector3D tmp;
863 rotation.multVec(m_pV[i],tmp);
864 m_pV[i] = tmp+translation;
865 }
866
867 // C H E C K D E T E R M I N A N T A N D
868 // 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
869
870 HVVector3D x; rotation.multVec(HVVector3D(1,0,0),x);
871 HVVector3D y; rotation.multVec(HVVector3D(0,1,0),y);
872 HVVector3D z; rotation.multVec(HVVector3D(0,0,1),z);
873 if ((x.cross(y)).dot(z) < 0) InvertFacets();
874 }
875 return *this;
876}
877*/
878//- new transform method from the newest OpenScientist/HEPVis
879SbPolyhedron & SbPolyhedron::Transform(
880 const HEPVis::SbRotation& rotation
881,const SbVec3d& translation
882)
883{
884 if (m_nvert > 0) {
885 SbVec3d tmp;
886 for (int i=1; i<=m_nvert; i++) {
887 rotation.multVec(m_pV[i],tmp);
888 m_pV[i] = tmp+translation;
889 }
890 SbVec3d x; rotation.multVec(SbVec3d(1,0,0),x);
891 SbVec3d y; rotation.multVec(SbVec3d(0,1,0),y);
892 SbVec3d z; rotation.multVec(SbVec3d(0,0,1),z);
893 if ((x.cross(y)).dot(z) < 0) InvertFacets();
894 }
895 return *this;
896}
897//---
898
899
900
901
902
903bool SbPolyhedron::GetNextVertexIndex(int &index, int &edgeFlag) const
904/***********************************************************************
905 * *
906 * Name: SbPolyhedron::GetNextVertexIndex Date: 03.09.96 *
907 * Author: Yasuhide Sawada Revised: *
908 * *
909 * Function: *
910 * *
911 ***********************************************************************/
912{
913 static int iFace = 1;
914 static int iQVertex = 0;
915 int vIndex = m_pF[iFace].m_edge[iQVertex].v;
916
917 edgeFlag = (vIndex > 0) ? 1 : 0;
918 index = iabs(vIndex);
919
920 if (iQVertex >= 3 || m_pF[iFace].m_edge[iQVertex+1].v == 0) {
921 iQVertex = 0;
922 if (++iFace > m_nface) iFace = 1;
923 return false; // Last Edge
924 }else{
925 ++iQVertex;
926 return true; // not Last Edge
927 }
928}
929
930HVPoint3D SbPolyhedron::GetVertex(int index) const
931/***********************************************************************
932 * *
933 * Name: SbPolyhedron::GetVertex Date: 03.09.96 *
934 * Author: Yasuhide Sawada Revised: 17.11.99 *
935 * *
936 * Function: Get vertex of the index. *
937 * *
938 ***********************************************************************/
939{
940 if (index <= 0 || index > m_nvert) {
941 std::cerr
942 << "SbPolyhedron::GetVertex: irrelevant index " << index
943 << std::endl;
944 return HVPoint3D();
945 }
946 return m_pV[index];
947}
948
949
950// rbianchi - 14.12.2012
951const HVPoint3D& SbPolyhedron::GetVertexFast(int index) const { //G.Barrand
952 return m_pV[index];
953}
954//---
955
956
957bool
958SbPolyhedron::GetNextVertex(HVPoint3D &vertex, int &edgeFlag) const
959/***********************************************************************
960 * *
961 * Name: SbPolyhedron::GetNextVertex Date: 22.07.96 *
962 * Author: John Allison Revised: *
963 * *
964 * Function: Get vertices of the quadrilaterals in order for each *
965 * face in face order. Returns false when finished each *
966 * face. *
967 * *
968 ***********************************************************************/
969{
970 int index;
971 bool rep = GetNextVertexIndex(index, edgeFlag);
972 vertex = m_pV[index];
973 return rep;
974}
975
976bool SbPolyhedron::GetNextVertex(HVPoint3D &vertex, int &edgeFlag,
977 HVNormal3D &normal) const
978/***********************************************************************
979 * *
980 * Name: SbPolyhedron::GetNextVertex Date: 26.11.99 *
981 * Author: E.Chernyaev Revised: *
982 * *
983 * Function: Get vertices with normals of the quadrilaterals in order *
984 * for each face in face order. *
985 * Returns false when finished each face. *
986 * *
987 ***********************************************************************/
988{
989 static int iFace = 1;
990 static int iNode = 0;
991
992 if (m_nface == 0) return false; // empty polyhedron
993
994 int k = m_pF[iFace].m_edge[iNode].v;
995 if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
996 vertex = m_pV[k];
997 normal = FindNodeNormal(iFace,k);
998 if (iNode >= 3 || m_pF[iFace].m_edge[iNode+1].v == 0) {
999 iNode = 0;
1000 if (++iFace > m_nface) iFace = 1;
1001 return false; // last node
1002 }else{
1003 ++iNode;
1004 return true; // not last node
1005 }
1006}
1007
1008bool SbPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag,
1009 int &iface1, int &iface2) const
1010/***********************************************************************
1011 * *
1012 * Name: SbPolyhedron::GetNextEdgeIndeces Date: 30.09.96 *
1013 * Author: E.Chernyaev Revised: 17.11.99 *
1014 * *
1015 * Function: Get indeces of the next edge together with indeces of *
1016 * of the faces which share the edge. *
1017 * Returns false when the last edge. *
1018 * *
1019 ***********************************************************************/
1020{
1021 static int iFace = 1;
1022 static int iQVertex = 0;
1023 static int iOrder = 1;
1024 int k1, k2, kflag, kface1, kface2;
1025
1026 if (iFace == 1 && iQVertex == 0) {
1027 k2 = m_pF[m_nface].m_edge[0].v;
1028 k1 = m_pF[m_nface].m_edge[3].v;
1029 if (k1 == 0) k1 = m_pF[m_nface].m_edge[2].v;
1030 if (iabs(k1) > iabs(k2)) iOrder = -1;
1031 }
1032
1033 do {
1034 k1 = m_pF[iFace].m_edge[iQVertex].v;
1035 kflag = k1;
1036 k1 = iabs(k1);
1037 kface1 = iFace;
1038 kface2 = m_pF[iFace].m_edge[iQVertex].f;
1039 if (iQVertex >= 3 || m_pF[iFace].m_edge[iQVertex+1].v == 0) {
1040 iQVertex = 0;
1041 k2 = iabs(m_pF[iFace].m_edge[iQVertex].v);
1042 iFace++;
1043 }else{
1044 iQVertex++;
1045 k2 = iabs(m_pF[iFace].m_edge[iQVertex].v);
1046 }
1047 } while (iOrder*k1 > iOrder*k2);
1048
1049 i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1050 iface1 = kface1; iface2 = kface2;
1051
1052 if (iFace > m_nface) {
1053 iFace = 1; iOrder = 1;
1054 return false;
1055 }else{
1056 return true;
1057 }
1058}
1059
1060bool
1061SbPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag) const
1062/***********************************************************************
1063 * *
1064 * Name: SbPolyhedron::GetNextEdgeIndeces Date: 17.11.99 *
1065 * Author: E.Chernyaev Revised: *
1066 * *
1067 * Function: Get indeces of the next edge. *
1068 * Returns false when the last edge. *
1069 * *
1070 ***********************************************************************/
1071{
1072 int kface1, kface2;
1073 return GetNextEdgeIndeces(i1, i2, edgeFlag, kface1, kface2);
1074}
1075
1076bool
1077SbPolyhedron::GetNextEdge(HVPoint3D &p1,
1078 HVPoint3D &p2,
1079 int &edgeFlag) const
1080/***********************************************************************
1081 * *
1082 * Name: SbPolyhedron::GetNextEdge Date: 30.09.96 *
1083 * Author: E.Chernyaev Revised: *
1084 * *
1085 * Function: Get next edge. *
1086 * Returns false when the last edge. *
1087 * *
1088 ***********************************************************************/
1089{
1090 int i1,i2;
1091 bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag);
1092 p1 = m_pV[i1];
1093 p2 = m_pV[i2];
1094 return rep;
1095}
1096
1097bool
1098SbPolyhedron::GetNextEdge(HVPoint3D &p1, HVPoint3D &p2,
1099 int &edgeFlag, int &iface1, int &iface2) const
1100/***********************************************************************
1101 * *
1102 * Name: SbPolyhedron::GetNextEdge Date: 17.11.99 *
1103 * Author: E.Chernyaev Revised: *
1104 * *
1105 * Function: Get next edge with indeces of the faces which share *
1106 * the edge. *
1107 * Returns false when the last edge. *
1108 * *
1109 ***********************************************************************/
1110{
1111 int i1,i2;
1112 bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag,iface1,iface2);
1113 p1 = m_pV[i1];
1114 p2 = m_pV[i2];
1115 return rep;
1116}
1117
1118void SbPolyhedron::GetFacet(int iFace, int &n, int *iNodes,
1119 int *edgeFlags, int *iFaces) const
1120/***********************************************************************
1121 * *
1122 * Name: SbPolyhedron::GetFacet Date: 15.12.99 *
1123 * Author: E.Chernyaev Revised: *
1124 * *
1125 * Function: Get face by index *
1126 * *
1127 ***********************************************************************/
1128{
1129 if (iFace < 1 || iFace > m_nface) {
1130 std::cerr
1131 << "SbPolyhedron::GetFacet: irrelevant index " << iFace
1132 << std::endl;
1133 n = 0;
1134 }else{
1135 int i, k;
1136 for (i=0; i<4; i++) {
1137 k = m_pF[iFace].m_edge[i].v;
1138 if (k == 0) break;
1139 if (iFaces != 0) iFaces[i] = m_pF[iFace].m_edge[i].f;
1140 if (k > 0) {
1141 iNodes[i] = k;
1142 if (edgeFlags != 0) edgeFlags[i] = 1;
1143 }else{
1144 iNodes[i] = -k;
1145 if (edgeFlags != 0) edgeFlags[i] = -1;
1146 }
1147 }
1148 n = i;
1149 }
1150}
1151
1152void SbPolyhedron::GetFacet(int index, int &n, HVPoint3D *nodes,
1153 int *edgeFlags, HVNormal3D *normals) const
1154/***********************************************************************
1155 * *
1156 * Name: SbPolyhedron::GetFacet Date: 17.11.99 *
1157 * Author: E.Chernyaev Revised: *
1158 * *
1159 * Function: Get face by index *
1160 * *
1161 ***********************************************************************/
1162{
1163 int iNodes[4];
1164 GetFacet(index, n, iNodes, edgeFlags);
1165 if (n != 0) {
1166 for (int i=0; i<4; i++) {
1167 nodes[i] = m_pV[iNodes[i]];
1168 if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
1169 }
1170 }
1171}
1172
1173bool
1174SbPolyhedron::GetNextFacet(int &n, HVPoint3D *nodes,
1175 int *edgeFlags, HVNormal3D *normals) const
1176/***********************************************************************
1177 * *
1178 * Name: SbPolyhedron::GetNextFacet Date: 19.11.99 *
1179 * Author: E.Chernyaev Revised: *
1180 * *
1181 * Function: Get next face with normals of unit length at the nodes. *
1182 * Returns false when finished all faces. *
1183 * *
1184 ***********************************************************************/
1185{
1186 static int iFace = 1;
1187
1188 if (edgeFlags == 0) {
1189 GetFacet(iFace, n, nodes);
1190 }else if (normals == 0) {
1191 GetFacet(iFace, n, nodes, edgeFlags);
1192 }else{
1193 GetFacet(iFace, n, nodes, edgeFlags, normals);
1194 }
1195
1196 if (++iFace > m_nface) {
1197 iFace = 1;
1198 return false;
1199 }else{
1200 return true;
1201 }
1202}
1203
1204HVNormal3D SbPolyhedron::GetNormal(int iFace) const
1205/***********************************************************************
1206 * *
1207 * Name: SbPolyhedron::GetNormal Date: 19.11.99 *
1208 * Author: E.Chernyaev Revised: *
1209 * *
1210 * Function: Get normal of the face given by index *
1211 * *
1212 ***********************************************************************/
1213{
1214 if (iFace < 1 || iFace > m_nface) {
1215 std::cerr
1216 << "SbPolyhedron::GetNormal: irrelevant index " << iFace
1217 << std::endl;
1218 return HVNormal3D();
1219 }
1220
1221 int i0 = iabs(m_pF[iFace].m_edge[0].v);
1222 int i1 = iabs(m_pF[iFace].m_edge[1].v);
1223 int i2 = iabs(m_pF[iFace].m_edge[2].v);
1224 int i3 = iabs(m_pF[iFace].m_edge[3].v);
1225 if (i3 == 0) i3 = i0;
1226 return (m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1]);
1227}
1228
1229HVNormal3D SbPolyhedron::GetUnitNormal(int iFace) const
1230/***********************************************************************
1231 * *
1232 * Name: SbPolyhedron::GetNormal Date: 19.11.99 *
1233 * Author: E.Chernyaev Revised: *
1234 * *
1235 * Function: Get unit normal of the face given by index *
1236 * *
1237 ***********************************************************************/
1238{
1239 if (iFace < 1 || iFace > m_nface) {
1240 std::cerr
1241 << "SbPolyhedron::GetUnitNormal: irrelevant index " << iFace
1242 << std::endl;
1243 return HVNormal3D();
1244 }
1245
1246 int i0 = iabs(m_pF[iFace].m_edge[0].v);
1247 int i1 = iabs(m_pF[iFace].m_edge[1].v);
1248 int i2 = iabs(m_pF[iFace].m_edge[2].v);
1249 int i3 = iabs(m_pF[iFace].m_edge[3].v);
1250 if (i3 == 0) i3 = i0;
1251 HVNormal3D nm = (m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1]);
1252 nm.normalize();
1253 return nm;
1254}
1255
1256bool SbPolyhedron::GetNextNormal(HVNormal3D &normal) const
1257/***********************************************************************
1258 * *
1259 * Name: SbPolyhedron::GetNextNormal Date: 22.07.96 *
1260 * Author: John Allison Revised: 19.11.99 *
1261 * *
1262 * Function: Get normals of each face in face order. Returns false *
1263 * when finished all faces. *
1264 * *
1265 ***********************************************************************/
1266{
1267 static int iFace = 1;
1268 normal = GetNormal(iFace);
1269 if (++iFace > m_nface) {
1270 iFace = 1;
1271 return false;
1272 }else{
1273 return true;
1274 }
1275}
1276
1277bool SbPolyhedron::GetNextUnitNormal(HVNormal3D &normal) const
1278/***********************************************************************
1279 * *
1280 * Name: SbPolyhedron::GetNextUnitNormal Date: 16.09.96 *
1281 * Author: E.Chernyaev Revised: *
1282 * *
1283 * Function: Get normals of unit length of each face in face order. *
1284 * Returns false when finished all faces. *
1285 * *
1286 ***********************************************************************/
1287{
1288 bool rep = GetNextNormal(normal);
1289 normal.normalize();
1290 return rep;
1291}
1292
1293double SbPolyhedron::GetSurfaceArea() const
1294/***********************************************************************
1295 * *
1296 * Name: SbPolyhedron::GetSurfaceArea Date: 25.05.01 *
1297 * Author: E.Chernyaev Revised: *
1298 * *
1299 * Function: Returns area of the surface of the polyhedron. *
1300 * *
1301 ***********************************************************************/
1302{
1303 double s = 0.;
1304 for (int iFace=1; iFace<=m_nface; iFace++) {
1305 int i0 = iabs(m_pF[iFace].m_edge[0].v);
1306 int i1 = iabs(m_pF[iFace].m_edge[1].v);
1307 int i2 = iabs(m_pF[iFace].m_edge[2].v);
1308 int i3 = iabs(m_pF[iFace].m_edge[3].v);
1309 if (i3 == 0) i3 = i0;
1310 s += ((m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1])).length();
1311 }
1312 return s/2.;
1313}
1314
1315double SbPolyhedron::GetVolume() const
1316/***********************************************************************
1317 * *
1318 * Name: SbPolyhedron::GetVolume Date: 25.05.01 *
1319 * Author: E.Chernyaev Revised: *
1320 * *
1321 * Function: Returns volume of the polyhedron. *
1322 * *
1323 ***********************************************************************/
1324{
1325 double v = 0.;
1326 for (int iFace=1; iFace<=m_nface; iFace++) {
1327 int i0 = iabs(m_pF[iFace].m_edge[0].v);
1328 int i1 = iabs(m_pF[iFace].m_edge[1].v);
1329 int i2 = iabs(m_pF[iFace].m_edge[2].v);
1330 int i3 = iabs(m_pF[iFace].m_edge[3].v);
1331 HVPoint3D g;
1332 if (i3 == 0) {
1333 i3 = i0;
1334 g = (m_pV[i0]+m_pV[i1]+m_pV[i2]) * (1.0f/3.0f);
1335 }else{
1336 g = (m_pV[i0]+m_pV[i1]+m_pV[i2]+m_pV[i3]) * 0.25f;
1337 }
1338 v += ((m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1])).dot(g);
1339 }
1340 return v/6.;
1341}
1342
1343SbPolyhedronTrd2::SbPolyhedronTrd2(double Dx1, double Dx2,
1344 double Dy1, double Dy2,
1345 double Dz)
1346/***********************************************************************
1347 * *
1348 * Name: SbPolyhedronTrd2 Date: 22.07.96 *
1349 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1350 * *
1351 * Function: Create GEANT4 TRD2-trapezoid *
1352 * *
1353 * Input: Dx1 - half-length along X at -Dz 8----7 *
1354 * Dx2 - half-length along X ay +Dz 5----6 ! *
1355 * Dy1 - half-length along Y ay -Dz ! 4-!--3 *
1356 * Dy2 - half-length along Y ay +Dz 1----2 *
1357 * Dz - half-length along Z *
1358 * *
1359 ***********************************************************************/
1360{
1361 AllocateMemory(8,6);
1362
1363 m_pV[1] = HVPoint3D(-Dx1,-Dy1,-Dz);
1364 m_pV[2] = HVPoint3D( Dx1,-Dy1,-Dz);
1365 m_pV[3] = HVPoint3D( Dx1, Dy1,-Dz);
1366 m_pV[4] = HVPoint3D(-Dx1, Dy1,-Dz);
1367 m_pV[5] = HVPoint3D(-Dx2,-Dy2, Dz);
1368 m_pV[6] = HVPoint3D( Dx2,-Dy2, Dz);
1369 m_pV[7] = HVPoint3D( Dx2, Dy2, Dz);
1370 m_pV[8] = HVPoint3D(-Dx2, Dy2, Dz);
1371
1372 CreatePrism();
1373}
1374
1375SbPolyhedronTrd2::~SbPolyhedronTrd2() {}
1376
1377SbPolyhedronTrd1::SbPolyhedronTrd1(double Dx1, double Dx2,
1378 double Dy, double Dz)
1379 : SbPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
1380
1381SbPolyhedronTrd1::~SbPolyhedronTrd1() {}
1382
1383SbPolyhedronBox::SbPolyhedronBox(double Dx, double Dy, double Dz)
1384 : SbPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1385
1386SbPolyhedronBox::~SbPolyhedronBox() {}
1387
1388SbPolyhedronTrap::SbPolyhedronTrap(double Dz,
1389 double Theta,
1390 double Phi,
1391 double Dy1,
1392 double Dx1,
1393 double Dx2,
1394 double Alp1,
1395 double Dy2,
1396 double Dx3,
1397 double Dx4,
1398 double Alp2)
1399/***********************************************************************
1400 * *
1401 * Name: SbPolyhedronTrap Date: 20.11.96 *
1402 * Author: E.Chernyaev Revised: *
1403 * *
1404 * Function: Create GEANT4 TRAP-trapezoid *
1405 * *
1406 * Input: DZ - half-length in Z *
1407 * Theta,Phi - polar angles of the line joining centres of the *
1408 * faces at Z=-Dz and Z=+Dz *
1409 * Dy1 - half-length in Y of the face at Z=-Dz *
1410 * Dx1 - half-length in X of low edge of the face at Z=-Dz *
1411 * Dx2 - half-length in X of top edge of the face at Z=-Dz *
1412 * Alp1 - angle between Y-axis and the median joining top and *
1413 * low edges of the face at Z=-Dz *
1414 * Dy2 - half-length in Y of the face at Z=+Dz *
1415 * Dx3 - half-length in X of low edge of the face at Z=+Dz *
1416 * Dx4 - half-length in X of top edge of the face at Z=+Dz *
1417 * Alp2 - angle between Y-axis and the median joining top and *
1418 * low edges of the face at Z=+Dz *
1419 * *
1420 ***********************************************************************/
1421{
1422 double DzTthetaCphi = Dz*tan(Theta)*cos(Phi);
1423 double DzTthetaSphi = Dz*tan(Theta)*sin(Phi);
1424 double Dy1Talp1 = Dy1*tan(Alp1);
1425 double Dy2Talp2 = Dy2*tan(Alp2);
1426
1427 AllocateMemory(8,6);
1428
1429 m_pV[1] = HVPoint3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
1430 m_pV[2] = HVPoint3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
1431 m_pV[3] = HVPoint3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
1432 m_pV[4] = HVPoint3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
1433 m_pV[5] = HVPoint3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
1434 m_pV[6] = HVPoint3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
1435 m_pV[7] = HVPoint3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
1436 m_pV[8] = HVPoint3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
1437
1438 CreatePrism();
1439}
1440
1441SbPolyhedronTrap::~SbPolyhedronTrap() {}
1442
1443SbPolyhedronPara::SbPolyhedronPara(double Dx, double Dy, double Dz,
1444 double Alpha, double Theta,
1445 double Phi)
1446 : SbPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
1447
1448SbPolyhedronPara::~SbPolyhedronPara() {}
1449
1450SbPolyhedronCons::SbPolyhedronCons(double Rmn1,
1451 double Rmx1,
1452 double Rmn2,
1453 double Rmx2,
1454 double Dz,
1455 double Phi1,
1456 double Dphi)
1457/***********************************************************************
1458 * *
1459 * Name: SbPolyhedronCons::SbPolyhedronCons Date: 15.12.96 *
1460 * Author: E.Chernyaev (IHEP/Protvino) Revised: 15.12.96 *
1461 * *
1462 * Function: Constructor for CONS, TUBS, CONE, TUBE *
1463 * *
1464 * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz *
1465 * Rmn2, Rmx2 - inside and outside radiuses at +Dz *
1466 * Dz - half length in Z *
1467 * Phi1 - starting angle of the segment *
1468 * Dphi - segment range *
1469 * *
1470 ***********************************************************************/
1471{
1472 static double wholeCircle=2*M_PI;
1473
1474 // C H E C K I N P U T P A R A M E T E R S
1475
1476 int k = 0;
1477 if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
1478 if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
1479 if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
1480
1481 // We can get this from the tracking geometry. Don't complain.
1482 if (Dz == 0) return;
1483
1484 if (Dz < 0.) k += 2;
1485
1486 double phi1, phi2, dphi;
1487 if (Dphi < 0.) {
1488 phi2 = Phi1; phi1 = phi2 - Dphi;
1489 }else if (Dphi == 0.) {
1490 phi1 = Phi1; phi2 = phi1 + wholeCircle;
1491 }else{
1492 phi1 = Phi1; phi2 = phi1 + Dphi;
1493 }
1494 dphi = phi2 - phi1;
1495 if (fabs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
1496 if (dphi > wholeCircle) k += 4;
1497
1498 if (k != 0) {
1499 std::cerr << "SbPolyhedronCone(s)/Tube(s): error in input parameters";
1500 if ((k & 1) != 0) std::cerr << " (radiuses)";
1501 if ((k & 2) != 0) std::cerr << " (half-length)";
1502 if ((k & 4) != 0) std::cerr << " (angles)";
1503 std::cerr << std::endl;
1504 std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
1505 std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
1506 std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
1507 << std::endl;
1508 return;
1509 }
1510
1511 // P R E P A R E T W O P O L Y L I N E S
1512
1513 double zz[4], rr[4];
1514 zz[0] = Dz;
1515 zz[1] = -Dz;
1516 zz[2] = Dz;
1517 zz[3] = -Dz;
1518 rr[0] = Rmx2;
1519 rr[1] = Rmx1;
1520 rr[2] = Rmn2;
1521 rr[3] = Rmn1;
1522
1523 // R O T A T E P O L Y L I N E S
1524
1525 RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1);
1526 SetReferences();
1527}
1528
1529SbPolyhedronCons::~SbPolyhedronCons() {}
1530
1531SbPolyhedronCone::SbPolyhedronCone(double Rmn1, double Rmx1,
1532 double Rmn2, double Rmx2,
1533 double Dz) :
1534 SbPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
1535
1536SbPolyhedronCone::~SbPolyhedronCone() {}
1537
1538SbPolyhedronTubs::SbPolyhedronTubs(double Rmin, double Rmax,
1539 double Dz,
1540 double Phi1, double Dphi)
1541 : SbPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
1542
1543SbPolyhedronTubs::~SbPolyhedronTubs() {}
1544
1545SbPolyhedronTube::SbPolyhedronTube (double Rmin, double Rmax,
1546 double Dz)
1547 : SbPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
1548
1549SbPolyhedronTube::~SbPolyhedronTube () {}
1550
1551SbPolyhedronPgon::SbPolyhedronPgon(double phi,
1552 double dphi,
1553 int npdv,
1554 int nz,
1555 const double *z,
1556 const double *rmin,
1557 const double *rmax)
1558/***********************************************************************
1559 * *
1560 * Name: SbPolyhedronPgon Date: 09.12.96 *
1561 * Author: E.Chernyaev Revised: *
1562 * *
1563 * Function: Constructor of polyhedron for PGON, PCON *
1564 * *
1565 * Input: phi - initial phi *
1566 * dphi - delta phi *
1567 * npdv - number of steps along phi *
1568 * nz - number of z-planes (at least two) *
1569 * z[] - z coordinates of the slices *
1570 * rmin[] - smaller r at the slices *
1571 * rmax[] - bigger r at the slices *
1572 * *
1573 ***********************************************************************/
1574{
1575 // C H E C K I N P U T P A R A M E T E R S
1576
1577 if (dphi <= 0. || dphi > 2*M_PI) {
1578 std::cerr
1579 << "SbPolyhedronPgon/Pcon: wrong delta phi = " << dphi
1580 << std::endl;
1581 return;
1582 }
1583
1584 if (nz < 2) {
1585 std::cerr
1586 << "SbPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
1587 << std::endl;
1588 return;
1589 }
1590
1591 if (npdv < 0) {
1592 std::cerr
1593 << "SbPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
1594 << std::endl;
1595 return;
1596 }
1597
1598 int i;
1599 for (i=0; i<nz; i++) {
1600 if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
1601 std::cerr
1602 << "SbPolyhedronPgon: error in radiuses rmin[" << i << "]="
1603 << rmin[i] << " rmax[" << i << "]=" << rmax[i]
1604 << std::endl;
1605 return;
1606 }
1607 }
1608
1609 // P R E P A R E T W O P O L Y L I N E S
1610
1611 double *zz, *rr;
1612 zz = new double[2*nz];
1613 rr = new double[2*nz];
1614
1615 if (z[0] > z[nz-1]) {
1616 for (i=0; i<nz; i++) {
1617 zz[i] = z[i];
1618 rr[i] = rmax[i];
1619 zz[i+nz] = z[i];
1620 rr[i+nz] = rmin[i];
1621 }
1622 }else{
1623 for (i=0; i<nz; i++) {
1624 zz[i] = z[nz-i-1];
1625 rr[i] = rmax[nz-i-1];
1626 zz[i+nz] = z[nz-i-1];
1627 rr[i+nz] = rmin[nz-i-1];
1628 }
1629 }
1630
1631 // R O T A T E P O L Y L I N E S
1632
1633 RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1);
1634 SetReferences();
1635
1636 delete [] zz;
1637 delete [] rr;
1638}
1639
1640SbPolyhedronPgon::~SbPolyhedronPgon() {}
1641
1642SbPolyhedronPcon::SbPolyhedronPcon(double phi, double dphi, int nz,
1643 const double *z,
1644 const double *rmin,
1645 const double *rmax)
1646 : SbPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
1647
1648SbPolyhedronPcon::~SbPolyhedronPcon() {}
1649
1650SbPolyhedronSphere::SbPolyhedronSphere(double rmin, double rmax,
1651 double phi, double dphi,
1652 double the, double dthe)
1653/***********************************************************************
1654 * *
1655 * Name: SbPolyhedronSphere Date: 11.12.96 *
1656 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1657 * *
1658 * Function: Constructor of polyhedron for SPHERE *
1659 * *
1660 * Input: rmin - internal radius *
1661 * rmax - external radius *
1662 * phi - initial phi *
1663 * dphi - delta phi *
1664 * the - initial theta *
1665 * dthe - delta theta *
1666 * *
1667 ***********************************************************************/
1668{
1669 // C H E C K I N P U T P A R A M E T E R S
1670
1671 if (dphi <= 0. || dphi > 2*M_PI) {
1672 std::cerr
1673 << "SbPolyhedronSphere: wrong delta phi = " << dphi
1674 << std::endl;
1675 return;
1676 }
1677
1678 if (the < 0. || the > M_PI) {
1679 std::cerr
1680 << "SbPolyhedronSphere: wrong theta = " << the
1681 << std::endl;
1682 return;
1683 }
1684
1685 if (dthe <= 0. || dthe > M_PI) {
1686 std::cerr
1687 << "SbPolyhedronSphere: wrong delta theta = " << dthe
1688 << std::endl;
1689 return;
1690 }
1691
1692 if ( (the+dthe >= M_PI) && (the+dthe < M_PI + 2*DBL_EPSILON) )
1693 dthe = M_PI - the; //G.Barrand : coming from LHCb/S.Ponce.
1694
1695 if (the+dthe > M_PI) {
1696 std::cerr
1697 << "SbPolyhedronSphere: wrong theta + delta theta = "
1698 << the << " " << dthe
1699 << std::endl;
1700 return;
1701 }
1702
1703 if (rmin < 0. || rmin >= rmax) {
1704 std::cerr
1705 << "SbPolyhedronSphere: error in radiuses"
1706 << " rmin=" << rmin << " rmax=" << rmax
1707 << std::endl;
1708 return;
1709 }
1710
1711 // P R E P A R E T W O P O L Y L I N E S
1712
1713 int ns = (GetNumberOfRotationSteps() + 1) / 2;
1714 int np1 = int(dthe*ns/M_PI+.5) + 1;
1715 if (np1 <= 1) np1 = 2;
1716 int np2 = rmin < perMillion ? 1 : np1;
1717
1718 double *zz, *rr;
1719 zz = new double[np1+np2];
1720 rr = new double[np1+np2];
1721
1722 double a = dthe/(np1-1);
1723 double cosa, sina;
1724 for (int i=0; i<np1; i++) {
1725 cosa = cos(the+i*a);
1726 sina = sin(the+i*a);
1727 zz[i] = rmax*cosa;
1728 rr[i] = rmax*sina;
1729 if (np2 > 1) {
1730 zz[i+np1] = rmin*cosa;
1731 rr[i+np1] = rmin*sina;
1732 }
1733 }
1734 if (np2 == 1) {
1735 zz[np1] = 0.;
1736 rr[np1] = 0.;
1737 }
1738
1739 // R O T A T E P O L Y L I N E S
1740
1741 RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
1742 SetReferences();
1743
1744 delete [] zz;
1745 delete [] rr;
1746}
1747
1748SbPolyhedronSphere::~SbPolyhedronSphere() {}
1749
1750SbPolyhedronTorus::SbPolyhedronTorus(double rmin,
1751 double rmax,
1752 double rtor,
1753 double phi,
1754 double dphi)
1755/***********************************************************************
1756 * *
1757 * Name: SbPolyhedronTorus Date: 11.12.96 *
1758 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1759 * *
1760 * Function: Constructor of polyhedron for TORUS *
1761 * *
1762 * Input: rmin - internal radius *
1763 * rmax - external radius *
1764 * rtor - radius of torus *
1765 * phi - initial phi *
1766 * dphi - delta phi *
1767 * *
1768 ***********************************************************************/
1769{
1770 // C H E C K I N P U T P A R A M E T E R S
1771
1772 if (dphi <= 0. || dphi > 2*M_PI) {
1773 std::cerr
1774 << "SbPolyhedronTorus: wrong delta phi = " << dphi
1775 << std::endl;
1776 return;
1777 }
1778
1779 if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
1780 std::cerr
1781 << "SbPolyhedronTorus: error in radiuses"
1782 << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
1783 << std::endl;
1784 return;
1785 }
1786
1787 // P R E P A R E T W O P O L Y L I N E S
1788
1789 int np1 = GetNumberOfRotationSteps();
1790 int np2 = rmin < perMillion ? 1 : np1;
1791
1792 //double *zz(nullptr), *rr;
1793 const auto totNumPts = np1+np2;
1794 auto zz = new double[totNumPts]{};
1795 auto rr = new double[totNumPts]{};
1796 for (unsigned int i(0);i!=(unsigned int) totNumPts;++i){
1797 rr[i]=0.0;
1798 zz[i]=0.0;
1799 }
1800
1801 double a = 2*M_PI/np1;
1802 double cosa, sina;
1803 for (int i=0; i<np1; i++) {
1804 cosa = cos(i*a);
1805 sina = sin(i*a);
1806 zz[i] = rmax*cosa;
1807 rr[i] = rtor+rmax*sina;
1808 if (np2 > 1) {
1809 zz[i+np1] = rmin*cosa;
1810 rr[i+np1] = rtor+rmin*sina;
1811 }
1812 }
1813 if (np2 == 1) {
1814 zz[np1] = 0.;
1815 rr[np1] = rtor;
1816 np2 = -1;
1817 }
1818
1819 // R O T A T E P O L Y L I N E S
1820
1821 RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1);
1822 SetReferences();
1823
1824 delete [] zz;
1825 delete [] rr;
1826}
1827
1828SbPolyhedronTorus::~SbPolyhedronTorus() {}
1829
1830int SbPolyhedron::s_numberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
1831/***********************************************************************
1832 * *
1833 * Name: SbPolyhedron::s_numberOfRotationSteps Date: 24.06.97 *
1834 * Author: J.Allison (Manchester University) Revised: *
1835 * *
1836 * Function: Number of steps for whole circle *
1837 * *
1838 ***********************************************************************/
1839
1840#include "BooleanProcessor.h"
1841static HEPVis_BooleanProcessor processor;
1842
1843SbPolyhedron SbPolyhedron::add(const SbPolyhedron & p) const
1844/***********************************************************************
1845 * *
1846 * Name: SbPolyhedron::add Date: 19.03.00 *
1847 * Author: E.Chernyaev Revised: *
1848 * *
1849 * Function: Boolean "union" of two polyhedra *
1850 * *
1851 ***********************************************************************/
1852{
1853 // rbianchi change
1854 //return processor.execute(OP_UNION, *this, p);
1855 HEPVis_BooleanProcessor processor; //G.Barrand
1856 int err;
1857 return processor.execute(OP_UNION, *this, p, err);
1858 //---
1859}
1860
1861SbPolyhedron SbPolyhedron::intersect(const SbPolyhedron & p) const
1862/***********************************************************************
1863 * *
1864 * Name: SbPolyhedron::intersect Date: 19.03.00 *
1865 * Author: E.Chernyaev Revised: *
1866 * *
1867 * Function: Boolean "intersection" of two polyhedra *
1868 * *
1869 ***********************************************************************/
1870{
1871 // rbianchi change
1872 //return processor.execute(OP_INTERSECTION, *this, p);
1873 HEPVis_BooleanProcessor processor; //G.Barrand
1874 int err;
1875 return processor.execute(OP_INTERSECTION, *this, p, err);
1876 //---
1877}
1878
1879SbPolyhedron SbPolyhedron::subtract(const SbPolyhedron & p) const
1880/***********************************************************************
1881 * *
1882 * Name: SbPolyhedron::add Date: 19.03.00 *
1883 * Author: E.Chernyaev Revised: *
1884 * *
1885 * Function: Boolean "subtraction" of "p" from "this" *
1886 * *
1887 ***********************************************************************/
1888{
1889 // rbianchi change
1890 //return processor.execute(OP_SUBTRACTION, *this, p);
1891 HEPVis_BooleanProcessor processor; //G.Barrand
1892 int err;
1893 return processor.execute(OP_SUBTRACTION, *this, p, err);
1894 //---
1895}
1896
1897
1902
1903#include "PolygonTriangulator.h"
1904#include <map>
1905#include <set>
1906
1907class SbPolyhedronPolygonXSect::Internals {
1908public:
1909
1910 typedef std::pair<int,int> Edge;
1911 typedef PolygonTriangulator::Triangle Triangle;
1912 typedef PolygonTriangulator::Triangles Triangles;
1913
1914 // Methods:
1918
1919 void setData(const std::vector<double> * xx,const std::vector<double>* yy, const double& dz);
1920 void triangulate();
1921 void initEdgeClassificationsAndNeighbours();
1922 void addExtraVertices();
1923 void allocateMemoryAndDefineVertexCoordinates();
1924 void defineFacesTopology();
1925 //Data:
1926 double dz;
1927 const std::vector<double> * x;
1928 const std::vector<double> * y;
1929 unsigned n;
1930 unsigned ntriangles;
1931 SbPolyhedronPolygonXSect* sbpolyhedron;
1932
1933 //The triangulation:
1934 std::unique_ptr<PolygonTriangulator> poly;
1935 //Various helper maps, sets, ...:
1936 std::set<Edge> edges_internal;
1937 std::set<Edge> edges_external;
1938 std::map<Edge,const Triangle*> neighbourmap;//In this map the Edges are oriented!
1939 std::set<Edge> edges_with_extra_vertex;
1940 unsigned nextraexternalvertices;
1941 unsigned nextrainternalvertices;
1942 std::map<Edge,unsigned> edgewithextravertex_2_id;
1943 std::map<Edge,unsigned> externaledgewithextravertex_2_2ndedgeid;
1944
1945 //Helper methods:
1946
1947 bool isinternaledge(const Edge& edge, const std::map<Edge,std::vector<Triangles::const_iterator> >& edge2triangles_map) const;
1948 Edge GetEdge(const Triangle * tr, const unsigned& i, const bool& oriented=false) const;//If oriented, 2-4 and 4-2 are different. If not, 4-2 is returned as 2-4.
1949 void setupface(const unsigned& face_id, const std::vector<unsigned>&v) const;
1950 void setupExternalFace(const unsigned& id_extedge,const unsigned& prev_edgeface,const unsigned& next_edgeface,
1951 const unsigned& idfirstvertex,const unsigned& idsecondvertex,
1952 const unsigned& id_triangleface,const unsigned& n, std::vector<unsigned>& faceinfo) const;
1953
1954 void getSurroundingValues(const std::vector<unsigned>& numbercycle, const unsigned& centralValue,
1955 unsigned& prevValue, unsigned& nextValue) const;
1956
1957 unsigned top2bottomvertexid(const unsigned& topid) const;
1958 unsigned top2bottomfaceid(const unsigned& topid) const;
1959};
1960
1961
1962//_______________________________________________________________________________________________________________________
1963SbPolyhedronPolygonXSect::SbPolyhedronPolygonXSect(const std::vector<double>& x,const std::vector<double>& y, const double& dz)
1964{
1965
1966 Internals * internals = new Internals(this);
1967 internals->setData(&x,&y,dz);
1968 internals->triangulate();
1969
1970 internals->initEdgeClassificationsAndNeighbours();
1971 internals->addExtraVertices();
1972 internals->allocateMemoryAndDefineVertexCoordinates();
1973 internals->defineFacesTopology();
1974
1975 delete internals;
1976}
1977
1978//_______________________________________________________________________________________________________________________
1979SbPolyhedronPolygonXSect::~SbPolyhedronPolygonXSect() {}
1980
1981//_______________________________________________________________________________________________________________________
1982void SbPolyhedronPolygonXSect::Internals::setData(const std::vector<double> * xx,const std::vector<double>* yy, const double& the_dz)
1983{
1984 n = xx->size();
1985 ntriangles = n-2;
1986 assert (n==yy->size()&&n>2);//fixme n>2, and special code for n==3.
1987 dz = the_dz;
1988 x = xx;
1989 y = yy;
1990 nextraexternalvertices = nextrainternalvertices = 0;
1991}
1992
1993//_______________________________________________________________________________________________________________________
1994void SbPolyhedronPolygonXSect::Internals::triangulate()
1995{
1996 poly = std::make_unique<PolygonTriangulator>(*x,*y);
1997}
1998
1999//_______________________________________________________________________________________________________________________
2000void SbPolyhedronPolygonXSect::Internals::initEdgeClassificationsAndNeighbours() {
2001 //Goal: classify edges as external/internal, get ability to find neighbours.
2002
2003 //As a help, we create a map of edges->triangles. We define an edge
2004 //as a pair of int's, where the first should always be lower than
2005 //the second (otherwise the edges 2-4 and 4-2 would seem different):
2006 std::map<Edge,std::vector<PolygonTriangulator::Triangles::const_iterator> > edge2triangles_map;
2007 std::set<Edge> alledges;
2008 Triangles::const_iterator itt=poly->triangles()->begin();
2009 Triangles::const_iterator ittE=poly->triangles()->end();
2010 for(; itt!=ittE; ++itt)
2011 for (unsigned iedge=0;iedge<3;++iedge) {
2012 Edge eun = GetEdge(&(*itt),iedge,false/*unoriented!*/);
2013 if (edge2triangles_map.find(eun)==edge2triangles_map.end()) {
2014 alledges.insert(eun);
2015 }
2016 edge2triangles_map[eun].push_back(itt);
2017 }
2018
2019
2020 //set of external/internal edges:
2021 std::set<Edge>::const_iterator itedges = alledges.begin();
2022 std::set<Edge>::const_iterator itedgesE = alledges.end();
2023 for (;itedges!=itedgesE;++itedges) {
2024 if (isinternaledge((*itedges),edge2triangles_map))
2025 edges_internal.insert(*itedges);
2026 else
2027 edges_external.insert(*itedges);
2028 }
2029
2030 //Neighbourmap:
2031 itt=poly->triangles()->begin();
2032 for(; itt!=ittE; ++itt)
2033 for (unsigned iedge=0;iedge<3;++iedge) {
2034 Edge e = GetEdge(&(*itt),iedge,true/*oriented!*/);
2035 Edge eun = GetEdge(&(*itt),iedge,false/*unoriented!*/);
2036 assert(edge2triangles_map.find(eun)!=edge2triangles_map.end());//fixme
2037 if (edges_external.find(eun)!=edges_external.end()) {
2038 neighbourmap[e]=0;
2039 neighbourmap[eun]=0;
2040 } else {
2041 assert(edge2triangles_map[eun].size()==2);//fixme
2042 const Triangle* triangle = static_cast<const Triangle*>( &(*itt) );
2043 const Triangle* triangle1 = static_cast<const Triangle*>( &(*(edge2triangles_map[eun].front())) );
2044 const Triangle* triangle2 = static_cast<const Triangle*>( &(*(edge2triangles_map[eun].back())) );
2045 neighbourmap[e]=(triangle==triangle1?triangle2:triangle1);
2046 }
2047 }
2048}
2049
2050//_______________________________________________________________________________________________________________________
2051void SbPolyhedronPolygonXSect::Internals::addExtraVertices() {
2052
2062
2063 std::set<const Triangle*> triangles_with_extra_vertex;//Eventually, all triangles should end up here.
2064
2066 Triangles::const_iterator itt=poly->triangles()->begin();
2067 Triangles::const_iterator ittE=poly->triangles()->end();
2068 for(; itt!=ittE; ++itt) {
2069 unsigned ninternaledges(0);
2070 for (unsigned iedge=0;iedge<3;++iedge) {
2071 Edge eun = GetEdge(&(*itt),iedge,false/*unoriented!*/);
2072 if (edges_internal.find(eun)!=edges_internal.end())
2073 ++ninternaledges;
2074 }
2075 if (ninternaledges<3)
2076 continue;
2077
2078
2079 //Check the three edges to find a neighbour that has not already become a quadrilateral.
2080 unsigned iedge(0);
2081 for (;iedge<3;++iedge) {
2082 Edge e = GetEdge(&(*itt),iedge,true/*oriented!*/);
2083 const Triangle* trneighbour = neighbourmap[e];
2084 assert(trneighbour);//fixme
2085 if (triangles_with_extra_vertex.insert(trneighbour).second) {
2086 triangles_with_extra_vertex.insert(&(*itt));
2087 Edge eun = GetEdge(&(*itt),iedge,false/*unoriented*/);
2088 edges_with_extra_vertex.insert(eun);
2089 ++nextrainternalvertices;
2090 break;
2091 }
2092 }
2093 assert(iedge<3);//If we cannot add a vertex to a triangle with 3 internal edges, then we have a problem...
2094 }
2095
2097 itt=poly->triangles()->begin();
2098 for(; itt!=ittE; ++itt) {
2099 Triangle* triangle = const_cast<Triangle*>( &(*itt) );
2100 if (triangles_with_extra_vertex.find(triangle)!=triangles_with_extra_vertex.end())
2101 continue;
2102 for (unsigned iedge=0;iedge<3;++iedge) {
2103 Edge eun = GetEdge(&(*itt),iedge,false/*unoriented!*/);
2104 //only look at external edges
2105 if (edges_external.find(eun)==edges_external.end())
2106 continue;
2107 assert(edges_with_extra_vertex.find(eun)==edges_with_extra_vertex.end());//fixme
2108 triangles_with_extra_vertex.insert(triangle);
2109 edges_with_extra_vertex.insert(eun);
2110 ++nextraexternalvertices;
2111 break;
2112 }
2113 }
2114 assert(ntriangles==triangles_with_extra_vertex.size());
2115}
2116
2117//_______________________________________________________________________________________________________________________
2118void SbPolyhedronPolygonXSect::Internals::allocateMemoryAndDefineVertexCoordinates() {
2119
2120 const unsigned nfaces = 2*ntriangles+n+nextraexternalvertices;
2121 const unsigned nvertices = 2*(n+nextraexternalvertices+nextrainternalvertices);
2122 sbpolyhedron->AllocateMemory(nvertices,nfaces);
2123
2125
2126 //Make edgewithextravertex_2_id map.
2127 std::set<Edge>::const_iterator itl = edges_with_extra_vertex.begin();
2128 std::set<Edge>::const_iterator itlE = edges_with_extra_vertex.end();
2129 unsigned il = 2*n;
2130 unsigned ielev = 3*n-4;
2131 for (;itl!=itlE;++itl) {
2132 edgewithextravertex_2_id[*itl]=++il;
2133 if (edges_external.find(*itl)!=edges_external.end()) {
2134 externaledgewithextravertex_2_2ndedgeid[*itl]=++ielev;
2135 }
2136 }
2137 //Now place the vertices:
2138
2139 //First the original n vertices at +dz
2140 for (unsigned i = 0; i<n; ++i)
2141 sbpolyhedron->m_pV[i+1] = HVPoint3D(x->at(i),y->at(i),dz);
2142 //Then the original n vertices at -dz
2143 for (unsigned i = 0; i<n; ++i)
2144 sbpolyhedron->m_pV[i+1+n] = HVPoint3D(x->at(i),y->at(i),-dz);
2145 //Extra vertices at +dz.
2146 std::map<Edge,unsigned>::const_iterator itlid = edgewithextravertex_2_id.begin();
2147 std::map<Edge,unsigned>::const_iterator itlidE = edgewithextravertex_2_id.end();
2148 for (; itlid!=itlidE; ++itlid) {
2149 double xx = 0.5*(x->at(itlid->first.first-1)+x->at(itlid->first.second-1));
2150 double yy = 0.5*(y->at(itlid->first.first-1)+y->at(itlid->first.second-1));
2151 sbpolyhedron->m_pV[itlid->second] = HVPoint3D(xx,yy,dz);
2152 }
2153 //Extra vertices at -dz.
2154 itlid = edgewithextravertex_2_id.begin();
2155 for (; itlid!=itlidE; ++itlid) {
2156 double xx = 0.5*(x->at(itlid->first.first-1)+x->at(itlid->first.second-1));
2157 double yy = 0.5*(y->at(itlid->first.first-1)+y->at(itlid->first.second-1));
2158 sbpolyhedron->m_pV[itlid->second+edgewithextravertex_2_id.size()] = HVPoint3D(xx,yy,-dz);
2159 }
2160
2161}
2162
2163//_______________________________________________________________________________________________________________________
2164bool SbPolyhedronPolygonXSect::Internals::isinternaledge(const Edge& edge,
2165 const std::map<Edge,std::vector<Triangles::const_iterator> >& edge2triangles_map) const {
2166
2167 Edge reverse_edge(edge.second,edge.first);
2168 unsigned k = 0;
2169 std::map<Edge,std::vector<Triangles::const_iterator> >::const_iterator it = edge2triangles_map.find(edge);
2170 if (it!=edge2triangles_map.end()) k += it->second.size();
2171 it = edge2triangles_map.find(reverse_edge);
2172 if (it!=edge2triangles_map.end()) k += it->second.size();
2173 return k==2;
2174
2175}
2176
2177//_______________________________________________________________________________________________________________________
2178SbPolyhedronPolygonXSect::Internals::Edge SbPolyhedronPolygonXSect::Internals::GetEdge(const Triangle * tr, const unsigned& i, const bool& oriented ) const {
2179 assert(i<=2);
2180 if (i==2)
2181 return oriented ? Edge((*tr)[2],(*tr)[0]) : Edge((*tr)[0],(*tr)[2]);
2182 else
2183 return Edge((*tr)[i],(*tr)[i+1]);
2184}
2185
2186//____________________________________________________________________________________________________
2187void SbPolyhedronPolygonXSect::Internals::setupface(const unsigned& face_id, const std::vector<unsigned>&v) const {
2188 assert(v.size()==8);
2189 sbpolyhedron->m_pF[face_id] = SbFacet(v.at(0),v.at(1),v.at(2),v.at(3),v.at(4),v.at(5),v.at(6),v.at(7));
2190}
2191
2192
2193//____________________________________________________________________________________________________
2194void SbPolyhedronPolygonXSect::Internals::setupExternalFace(const unsigned& id_extedge,const unsigned& prev_edgeface,const unsigned& next_edgeface,
2195 const unsigned& idfirstvertex,const unsigned& idsecondvertex,
2196 const unsigned& id_triangleface,const unsigned& n, std::vector<unsigned>& faceinfo) const {
2197 faceinfo.clear();
2198 faceinfo.push_back(idfirstvertex);//point
2199 faceinfo.push_back(prev_edgeface);//neighbour face
2200 faceinfo.push_back(top2bottomvertexid(idfirstvertex));//point
2201 faceinfo.push_back(id_triangleface+n-2);//neighbour face
2202 faceinfo.push_back(top2bottomvertexid(idsecondvertex));//point
2203 faceinfo.push_back(next_edgeface);//neighbour face
2204 faceinfo.push_back(idsecondvertex);//point
2205 faceinfo.push_back(id_triangleface);//neighbour face
2206 setupface(id_extedge,faceinfo);
2207}
2208
2209
2210//____________________________________________________________________________________________________
2211void SbPolyhedronPolygonXSect::Internals::getSurroundingValues(const std::vector<unsigned>& numbercycle, const unsigned& centralValue,
2212 unsigned& prevValue, unsigned& nextValue) const
2213{
2214 //First, find the index of the centralValue
2215 unsigned int i;
2216 for ( i=0/*fixme, not O(n)*/;i<numbercycle.size();++i)
2217 if (numbercycle.at(i)==centralValue)
2218 break;
2219 //Sanity check that the central value was found
2220 assert(i<numbercycle.size());//fixme
2221 prevValue = ( i==0 ? numbercycle.back() : numbercycle.at(i-1) );
2222 nextValue = ( i==numbercycle.size()-1 ? numbercycle.front() : numbercycle.at(i+1) );
2223
2224}
2225
2226
2227//_______________________________________________________________________________________________________________________
2228unsigned SbPolyhedronPolygonXSect::Internals::top2bottomvertexid(const unsigned& topid) const {
2229 if (topid<=n)
2230 return topid+n;//original vertex
2231 return topid+nextraexternalvertices+nextrainternalvertices;//added extra vertex;
2232}
2233
2234//_______________________________________________________________________________________________________________________
2235unsigned SbPolyhedronPolygonXSect::Internals::top2bottomfaceid(const unsigned& topid) const {
2236 if (topid<=ntriangles)
2237 return topid+ntriangles;//triangle
2238 return topid;//side face
2239}
2240
2241
2242//_______________________________________________________________________________________________________________________
2243void SbPolyhedronPolygonXSect::Internals::defineFacesTopology() {
2244 //This is definitely the most tricky part to get right...
2245
2246 //Make triangle_2_id map.
2247 std::map<const Triangle*,unsigned> triangle_2_id;
2248 Triangles::const_iterator itt=poly->triangles()->begin();
2249 Triangles::const_iterator ittE=poly->triangles()->end();
2250 unsigned itri = 0;
2251 for(; itt!=ittE; ++itt) {
2252 Triangle* triangle = const_cast<Triangle*>( &(*itt) );
2253 triangle_2_id[triangle]=++itri;
2254 }
2255 assert(triangle_2_id.size()==ntriangles);
2256
2257 //Loop over the external edges, and get a list of edge face id's.
2258 std::vector<unsigned> edgeface_ids;
2259 std::map<Edge,unsigned> edge_2_firstid;
2260 std::map<Edge,unsigned> edge_2_secondid;
2261 for (unsigned i = 1; i<=n;++i) {
2262 Edge eun((i==n?1:i),(i==n?n:i+1));
2263 edgeface_ids.push_back(2*n-4+i);
2264 edge_2_firstid[eun]=edgeface_ids.back();
2265 bool hasextravertex=
2266 (edges_with_extra_vertex.find(eun)!=edges_with_extra_vertex.end());
2267 if (hasextravertex) {
2268 edgeface_ids.push_back(externaledgewithextravertex_2_2ndedgeid[eun]);
2269 edge_2_secondid[eun]=edgeface_ids.back();
2270 }
2271 }
2272
2273 std::vector<unsigned> faceinfo_top, faceinfo_bottom(8), faceinfo_sides;
2274 itt = poly->triangles()->begin();
2275 for(; itt!=ittE; ++itt) {
2276 unsigned triangleid = triangle_2_id[&(*itt)];
2277
2279 for (unsigned iedge=0;iedge<3;++iedge) {
2280 Edge e = GetEdge(&(*itt),iedge,true/*oriented!*/);
2281 Edge eun = GetEdge(&(*itt),iedge,false/*unoriented!*/);
2282
2283 bool hasextravertex(edges_with_extra_vertex.find(eun)!=edges_with_extra_vertex.end());
2284 bool isinternal(edges_external.find(eun)==edges_external.end());
2285 //add for the triangle.
2286 faceinfo_top.push_back(e.first);
2287 if (isinternal) {
2288 assert(neighbourmap[e]);//fixme
2289 assert(triangle_2_id.find(neighbourmap[e])!=triangle_2_id.end());
2290 faceinfo_top.push_back(triangle_2_id[neighbourmap[e]]);
2291 if (hasextravertex) {
2292 faceinfo_top.push_back(edgewithextravertex_2_id[eun]);
2293 faceinfo_top.push_back(triangle_2_id[neighbourmap[e]]);//NB: Same neighbourface again!
2294 // assert(false&&"need to implement extravertex for internal edges");
2295 //fixme!
2296 }
2297 } else {
2298 assert(edge_2_firstid.find(eun)!=edge_2_firstid.end());
2299 unsigned firstsideid = edge_2_firstid[eun];
2300 faceinfo_top.push_back(firstsideid);
2301 if (hasextravertex) {
2302 faceinfo_top.push_back(edgewithextravertex_2_id[eun]);
2303 assert(edge_2_secondid.find(eun)!=edge_2_secondid.end());
2304 faceinfo_top.push_back(edge_2_secondid[eun]);
2305 }
2307 unsigned prev_edgeface,next_edgeface;
2308 getSurroundingValues(edgeface_ids, firstsideid,prev_edgeface, next_edgeface);
2309 unsigned firstvertexid = e.first;
2310 unsigned secondvertexid = (hasextravertex? edgewithextravertex_2_id[eun] : e.second);
2311 setupExternalFace(firstsideid,prev_edgeface, next_edgeface,firstvertexid,secondvertexid,
2312 /*firstvertexid+n,secondvertexid + (hasextravertex?externaledgewithextravertex_2_2ndedgeid.size():n),*/
2313 triangleid,n,faceinfo_sides);
2314 if (hasextravertex) {
2315 firstsideid=next_edgeface;
2316 getSurroundingValues(edgeface_ids, firstsideid,prev_edgeface, next_edgeface);
2317 firstvertexid=secondvertexid;
2318 secondvertexid=e.second;
2319 setupExternalFace(firstsideid,prev_edgeface, next_edgeface,
2320 edgewithextravertex_2_id[eun],secondvertexid,
2321 /*edgewithextravertex_2_id[eun]+externaledgewithextravertex_2_2ndedgeid.size(),secondvertexid+n,*/
2322 triangleid,n,faceinfo_sides);
2323 }
2324
2325 }
2326 }//End edge loop
2327
2329 assert(faceinfo_bottom.size()==8);//fixme
2330 assert(faceinfo_top.size()==8);//fixme
2331 faceinfo_bottom[0]=top2bottomvertexid(faceinfo_top[0]);
2332 faceinfo_bottom[1]=top2bottomfaceid(faceinfo_top[7]);
2333 faceinfo_bottom[2]=top2bottomvertexid(faceinfo_top[6]);
2334 faceinfo_bottom[3]=top2bottomfaceid(faceinfo_top[5]);
2335 faceinfo_bottom[4]=top2bottomvertexid(faceinfo_top[4]);
2336 faceinfo_bottom[5]=top2bottomfaceid(faceinfo_top[3]);
2337 faceinfo_bottom[6]=top2bottomvertexid(faceinfo_top[2]);
2338 faceinfo_bottom[7]=top2bottomfaceid(faceinfo_top[1]);
2339
2340 //Actually setup triangle faces.
2341 setupface(triangleid,faceinfo_top);
2342 setupface(triangleid+ntriangles,faceinfo_bottom);
2343 faceinfo_top.clear();
2344 }//End triangle loop
2345
2346}
2347
2348// ************************************************************************************************ //
2349// * * //
2350// * Class SbPolyhedronArbitrary for Tessellated Solids * //
2351// * * //
2352// ************************************************************************************************ //
2353
2354SbPolyhedronArbitrary::SbPolyhedronArbitrary(const int nVertices, const int nFacets)
2355{
2356 AllocateMemory(nVertices, nFacets);
2357 m_nVertexCount = 0;
2358 m_nFacetCount = 0;
2359}
2360
2364
2365void SbPolyhedronArbitrary::AddVertex(const double v1, const double v2, const double v3)
2366{
2367 if(m_nVertexCount==m_nvert+1) {
2368 std::cerr <<"ERROR in SbPolyhedronArbitrary::AddVertex. Attempt to exceed maximum number of vertices: "
2369 << m_nVertexCount << std::endl;
2370 }
2371 else {
2372 m_nVertexCount++;
2373 m_pV[m_nVertexCount] = HVPoint3D(v1,v2,v3);
2374 }
2375}
2376
2377void SbPolyhedronArbitrary::AddFacet(const int iv1, const int iv2, const int iv3, const int iv4)
2378{
2379 if(m_nFacetCount==m_nface) {
2380 std::cerr <<"ERROR in SbPolyhedronArbitrary::AddFacet. Attempt to exceed maximum number of facets: "
2381 << m_nFacetCount << std::endl;
2382 }
2383 else if(iv1 < 1
2384 || iv1 > m_nvert
2385 || iv2 < 1
2386 || iv2 > m_nvert
2387 || iv3 < 1
2388 || iv3 > m_nvert
2389 || iv4 > m_nvert) {
2390 std::cerr <<"ERROR in SbPolyhedronArbitrary::AddFacet. Attempt to index vertex number which is out-of-range: ("
2391 << iv1 << "," << iv2 << "," << iv3 << "," << iv4 << ")" << std::endl;
2392 }
2393 else if(iv1 > m_nVertexCount
2394 || iv2 > m_nVertexCount
2395 || iv3 > m_nVertexCount
2396 || iv4 > m_nVertexCount) {
2397 std::cerr <<"ERROR in SbPolyhedronArbitrary::AddFacet. Vertex needs to be defined first : ("
2398 << iv1 << "," << iv2 << "," << iv3 << "," << iv4 << ")" << std::endl;
2399 }
2400 else {
2401 m_nFacetCount++;
2402 m_pF[m_nFacetCount] = SbFacet(iv1, 0, iv2, 0, iv3, 0, iv4, 0);
2403 }
2404}
2405
2410
2411SbPolyhedronGenericTrap::SbPolyhedronGenericTrap(double Dz, const std::vector<std::pair<double,double> >& Vertices)
2412{
2413 AllocateMemory(8,6);
2414
2415 m_pV[1] = HVPoint3D(Vertices[0].first,Vertices[0].second,-Dz);
2416 m_pV[2] = HVPoint3D(Vertices[1].first,Vertices[1].second,-Dz);
2417 m_pV[3] = HVPoint3D(Vertices[2].first,Vertices[2].second,-Dz);
2418 m_pV[4] = HVPoint3D(Vertices[3].first,Vertices[3].second,-Dz);
2419 m_pV[5] = HVPoint3D(Vertices[4].first,Vertices[4].second,Dz);
2420 m_pV[6] = HVPoint3D(Vertices[5].first,Vertices[5].second,Dz);
2421 m_pV[7] = HVPoint3D(Vertices[6].first,Vertices[6].second,Dz);
2422 m_pV[8] = HVPoint3D(Vertices[7].first,Vertices[7].second,Dz);
2423
2424 CreatePrism();
2425}
2426
rr
const boost::regex rr(r_r)
M_PI
#define M_PI
Definition ActiveFraction.h:11
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
OP_INTERSECTION
#define OP_INTERSECTION
Definition BooleanProcessor.h:90
OP_SUBTRACTION
#define OP_SUBTRACTION
Definition BooleanProcessor.h:91
OP_UNION
#define OP_UNION
Definition BooleanProcessor.h:89
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
operator+
HVPoint3D operator+(const HVPoint3D &v1, const HVPoint3D &v2)
Definition SbPolyhedron.cxx:44
perMillion
#define perMillion
Definition SbPolyhedron.cxx:17
iabs
int iabs(int a)
Definition SbPolyhedron.cxx:23
ffabs
float ffabs(float a)
Definition SbPolyhedron.cxx:26
deg
#define deg
Definition SbPolyhedron.cxx:18
processor
static HEPVis_BooleanProcessor processor
Definition SbPolyhedron.cxx:1841
operator<<
std::ostream & operator<<(std::ostream &ostr, const SbFacet &facet)
Definition SbPolyhedron.cxx:119
HVNormal3D
HVPoint3D HVNormal3D
Definition SbPolyhedron.h:199
DEFAULT_NUMBER_OF_STEPS
#define DEFAULT_NUMBER_OF_STEPS
Definition SbPolyhedron.h:214
y
#define y
x
#define x
z
#define z
HVPoint3D::operator=
HVPoint3D & operator=(const HVPoint3D &v)
Definition SbPolyhedron.cxx:36
PolygonTriangulator::Triangles
std::list< Triangle > Triangles
Definition PolygonTriangulator.h:36
PolygonTriangulator::Triangle
std::vector< unsigned > Triangle
Definition PolygonTriangulator.h:35
SbFacet::SbFacet
SbFacet(int v1=0, int f1=0, int v2=0, int f2=0, int v3=0, int f3=0, int v4=0, int f4=0)
Definition SbPolyhedron.h:225
SbFacet::m_edge
struct SbFacet::@146305307241173307306147201343177135163056232201 m_edge[4]
SbPolyhedronArbitrary::SbPolyhedronArbitrary
SbPolyhedronArbitrary(const int nVertices, const int nFacets)
Definition SbPolyhedron.cxx:2354
SbPolyhedronArbitrary::AddFacet
void AddFacet(const int iv1, const int iv2, const int iv3, const int iv4=0)
Definition SbPolyhedron.cxx:2377
SbPolyhedronArbitrary::AddVertex
void AddVertex(const double v1, const double v2, const double v3)
Definition SbPolyhedron.cxx:2365
SbPolyhedronBox::SbPolyhedronBox
SbPolyhedronBox(double Dx, double Dy, double Dz)
Definition SbPolyhedron.cxx:1383
SbPolyhedronBox::~SbPolyhedronBox
virtual ~SbPolyhedronBox()
Definition SbPolyhedron.cxx:1386
SbPolyhedronCone::~SbPolyhedronCone
virtual ~SbPolyhedronCone()
Definition SbPolyhedron.cxx:1536
SbPolyhedronCone::SbPolyhedronCone
SbPolyhedronCone(double Rmn1, double Rmx1, double Rmn2, double Rmx2, double Dz)
Definition SbPolyhedron.cxx:1531
SbPolyhedronCons::SbPolyhedronCons
SbPolyhedronCons(double Rmn1, double Rmx1, double Rmn2, double Rmx2, double Dz, double Phi1, double Dphi)
Definition SbPolyhedron.cxx:1450
SbPolyhedronCons::~SbPolyhedronCons
virtual ~SbPolyhedronCons()
Definition SbPolyhedron.cxx:1529
SbPolyhedronGenericTrap::SbPolyhedronGenericTrap
SbPolyhedronGenericTrap(double Dz, const std::vector< std::pair< double, double > > &Vertices)
Definition SbPolyhedron.cxx:2411
SbPolyhedronPara::~SbPolyhedronPara
virtual ~SbPolyhedronPara()
Definition SbPolyhedron.cxx:1448
SbPolyhedronPara::SbPolyhedronPara
SbPolyhedronPara(double Dx, double Dy, double Dz, double Alpha, double Theta, double Phi)
Definition SbPolyhedron.cxx:1443
SbPolyhedronPcon::SbPolyhedronPcon
SbPolyhedronPcon(double phi, double dphi, int nz, const double *z, const double *rmin, const double *rmax)
Definition SbPolyhedron.cxx:1642
SbPolyhedronPcon::~SbPolyhedronPcon
virtual ~SbPolyhedronPcon()
Definition SbPolyhedron.cxx:1648
SbPolyhedronPgon::~SbPolyhedronPgon
virtual ~SbPolyhedronPgon()
Definition SbPolyhedron.cxx:1640
SbPolyhedronPgon::SbPolyhedronPgon
SbPolyhedronPgon(double phi, double dphi, int npdv, int nz, const double *z, const double *rmin, const double *rmax)
Definition SbPolyhedron.cxx:1551
SbPolyhedronPolygonXSect::Internals::setupExternalFace
void setupExternalFace(const unsigned &id_extedge, const unsigned &prev_edgeface, const unsigned &next_edgeface, const unsigned &idfirstvertex, const unsigned &idsecondvertex, const unsigned &id_triangleface, const unsigned &n, std::vector< unsigned > &faceinfo) const
Definition SbPolyhedron.cxx:2194
SbPolyhedronPolygonXSect::Internals::Triangle
PolygonTriangulator::Triangle Triangle
Definition SbPolyhedron.cxx:1911
SbPolyhedronPolygonXSect::Internals::Triangles
PolygonTriangulator::Triangles Triangles
Definition SbPolyhedron.cxx:1912
SbPolyhedronPolygonXSect::Internals::y
const std::vector< double > * y
Definition SbPolyhedron.cxx:1928
SbPolyhedronPolygonXSect::Internals::Internals
Internals(SbPolyhedronPolygonXSect *sbp)
Definition SbPolyhedron.cxx:1915
SbPolyhedronPolygonXSect::Internals::GetEdge
Edge GetEdge(const Triangle *tr, const unsigned &i, const bool &oriented=false) const
Definition SbPolyhedron.cxx:2178
SbPolyhedronPolygonXSect::Internals::isinternaledge
bool isinternaledge(const Edge &edge, const std::map< Edge, std::vector< Triangles::const_iterator > > &edge2triangles_map) const
Definition SbPolyhedron.cxx:2164
SbPolyhedronPolygonXSect::Internals::externaledgewithextravertex_2_2ndedgeid
std::map< Edge, unsigned > externaledgewithextravertex_2_2ndedgeid
Definition SbPolyhedron.cxx:1943
SbPolyhedronPolygonXSect::Internals::edgewithextravertex_2_id
std::map< Edge, unsigned > edgewithextravertex_2_id
Definition SbPolyhedron.cxx:1942
SbPolyhedronPolygonXSect::Internals::getSurroundingValues
void getSurroundingValues(const std::vector< unsigned > &numbercycle, const unsigned &centralValue, unsigned &prevValue, unsigned &nextValue) const
Definition SbPolyhedron.cxx:2211
SbPolyhedronPolygonXSect::Internals::setupface
void setupface(const unsigned &face_id, const std::vector< unsigned > &v) const
Definition SbPolyhedron.cxx:2187
SbPolyhedronPolygonXSect::Internals::top2bottomfaceid
unsigned top2bottomfaceid(const unsigned &topid) const
Definition SbPolyhedron.cxx:2235
SbPolyhedronPolygonXSect::Internals::top2bottomvertexid
unsigned top2bottomvertexid(const unsigned &topid) const
Definition SbPolyhedron.cxx:2228
SbPolyhedronPolygonXSect::Internals::x
const std::vector< double > * x
Definition SbPolyhedron.cxx:1927
SbPolyhedronPolygonXSect::Internals::setData
void setData(const std::vector< double > *xx, const std::vector< double > *yy, const double &dz)
Definition SbPolyhedron.cxx:1982
SbPolyhedronPolygonXSect::Internals::neighbourmap
std::map< Edge, const Triangle * > neighbourmap
Definition SbPolyhedron.cxx:1938
SbPolyhedronPolygonXSect::Internals::sbpolyhedron
SbPolyhedronPolygonXSect * sbpolyhedron
Definition SbPolyhedron.cxx:1931
SbPolyhedronPolygonXSect::Internals::poly
std::unique_ptr< PolygonTriangulator > poly
Definition SbPolyhedron.cxx:1934
SbPolyhedronPolygonXSect::SbPolyhedronPolygonXSect
SbPolyhedronPolygonXSect(const std::vector< double > &x, const std::vector< double > &y, const double &dz)
Definition SbPolyhedron.cxx:1963
SbPolyhedronSphere::~SbPolyhedronSphere
virtual ~SbPolyhedronSphere()
Definition SbPolyhedron.cxx:1748
SbPolyhedronSphere::SbPolyhedronSphere
SbPolyhedronSphere(double rmin, double rmax, double phi, double dphi, double the, double dthe)
Definition SbPolyhedron.cxx:1650
SbPolyhedronTorus::SbPolyhedronTorus
SbPolyhedronTorus(double rmin, double rmax, double rtor, double phi, double dphi)
Definition SbPolyhedron.cxx:1750
SbPolyhedronTorus::~SbPolyhedronTorus
virtual ~SbPolyhedronTorus()
Definition SbPolyhedron.cxx:1828
SbPolyhedronTrap::SbPolyhedronTrap
SbPolyhedronTrap(double Dz, double Theta, double Phi, double Dy1, double Dx1, double Dx2, double Alp1, double Dy2, double Dx3, double Dx4, double Alp2)
Definition SbPolyhedron.cxx:1388
SbPolyhedronTrap::~SbPolyhedronTrap
virtual ~SbPolyhedronTrap()
Definition SbPolyhedron.cxx:1441
SbPolyhedronTrd1::~SbPolyhedronTrd1
virtual ~SbPolyhedronTrd1()
Definition SbPolyhedron.cxx:1381
SbPolyhedronTrd1::SbPolyhedronTrd1
SbPolyhedronTrd1(double Dx1, double Dx2, double Dy, double Dz)
Definition SbPolyhedron.cxx:1377
SbPolyhedronTrd2::SbPolyhedronTrd2
SbPolyhedronTrd2(double Dx1, double Dx2, double Dy1, double Dy2, double Dz)
Definition SbPolyhedron.cxx:1343
SbPolyhedronTrd2::~SbPolyhedronTrd2
virtual ~SbPolyhedronTrd2()
Definition SbPolyhedron.cxx:1375
SbPolyhedronTube::SbPolyhedronTube
SbPolyhedronTube(double Rmin, double Rmax, double Dz)
Definition SbPolyhedron.cxx:1545
SbPolyhedronTube::~SbPolyhedronTube
virtual ~SbPolyhedronTube()
Definition SbPolyhedron.cxx:1549
SbPolyhedronTubs::~SbPolyhedronTubs
virtual ~SbPolyhedronTubs()
Definition SbPolyhedron.cxx:1543
SbPolyhedronTubs::SbPolyhedronTubs
SbPolyhedronTubs(double Rmin, double Rmax, double Dz, double Phi1, double Dphi)
Definition SbPolyhedron.cxx:1538
SbPolyhedron::GetNextVertex
bool GetNextVertex(HVPoint3D &vertex, int &edgeFlag) const
Definition SbPolyhedron.cxx:958
SbPolyhedron::GetNextEdge
bool GetNextEdge(HVPoint3D &p1, HVPoint3D &p2, int &edgeFlag) const
Definition SbPolyhedron.cxx:1077
SbPolyhedron::intersect
SbPolyhedron intersect(const SbPolyhedron &p) const
Definition SbPolyhedron.cxx:1861
SbPolyhedron::add
SbPolyhedron add(const SbPolyhedron &p) const
Definition SbPolyhedron.cxx:1843
SbPolyhedron::FindNeighbour
int FindNeighbour(int iFace, int iNode, int iOrder) const
Definition SbPolyhedron.cxx:190
SbPolyhedron::operator=
virtual SbPolyhedron & operator=(const SbPolyhedron &from)
Definition SbPolyhedron.cxx:162
SbPolyhedron::GetNextNormal
bool GetNextNormal(HVNormal3D &normal) const
Definition SbPolyhedron.cxx:1256
SbPolyhedron::RotateAroundZ
void RotateAroundZ(int nstep, double phi, double dphi, int np1, int np2, const double *z, double *r, int nodeVis, int edgeVis)
Definition SbPolyhedron.cxx:439
SbPolyhedron::GetNextEdgeIndeces
bool GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag, int &iface1, int &iface2) const
Definition SbPolyhedron.cxx:1008
SbPolyhedron::GetSurfaceArea
double GetSurfaceArea() const
Definition SbPolyhedron.cxx:1293
SbPolyhedron::s_numberOfRotationSteps
static int s_numberOfRotationSteps
Definition SbPolyhedron.h:235
SbPolyhedron::m_pV
HVPoint3D * m_pV
Definition SbPolyhedron.h:239
SbPolyhedron::AllocateMemory
void AllocateMemory(int Nvert, int Nface)
Definition SbPolyhedron.cxx:268
SbPolyhedron::GetVertexFast
const HVPoint3D & GetVertexFast(int index) const
Definition SbPolyhedron.cxx:951
SbPolyhedron::GetUnitNormal
HVNormal3D GetUnitNormal(int iFace) const
Definition SbPolyhedron.cxx:1229
SbPolyhedron::GetNextUnitNormal
bool GetNextUnitNormal(HVNormal3D &normal) const
Definition SbPolyhedron.cxx:1277
SbPolyhedron::subtract
SbPolyhedron subtract(const SbPolyhedron &p) const
Definition SbPolyhedron.cxx:1879
SbPolyhedron::GetVolume
double GetVolume() const
Definition SbPolyhedron.cxx:1315
SbPolyhedron::SetSideFacets
void SetSideFacets(int ii[4], int vv[4], int *kk, double *r, double dphi, int ns, int &kface)
Definition SbPolyhedron.cxx:369
SbPolyhedron::GetNextVertexIndex
bool GetNextVertexIndex(int &index, int &edgeFlag) const
Definition SbPolyhedron.cxx:903
SbPolyhedron::RotateEdge
void RotateEdge(int k1, int k2, double r1, double r2, int v1, int v2, int vEdge, bool ifWholeCircle, int ns, int &kface)
Definition SbPolyhedron.cxx:307
SbPolyhedron::GetVertex
HVPoint3D GetVertex(int index) const
Definition SbPolyhedron.cxx:930
SbPolyhedron::Transform
SbPolyhedron & Transform(const HEPVis::SbRotation &rot, const SbVec3d &trans)
Definition SbPolyhedron.cxx:879
SbPolyhedron::SbPolyhedron
SbPolyhedron(int Nvert=0, int Nface=0)
Definition SbPolyhedron.h:278
SbPolyhedron::GetNumberOfRotationSteps
static int GetNumberOfRotationSteps()
Definition SbPolyhedron.h:375
SbPolyhedron::m_pF
SbFacet * m_pF
Definition SbPolyhedron.h:240
SbPolyhedron::FindNodeNormal
HVNormal3D FindNodeNormal(int iFace, int iNode) const
Definition SbPolyhedron.cxx:218
SbPolyhedron::GetFacet
void GetFacet(int iFace, int &n, int *iNodes, int *edgeFlags=0, int *iFaces=0) const
Definition SbPolyhedron.cxx:1118
SbPolyhedron::GetNextFacet
bool GetNextFacet(int &n, HVPoint3D *nodes, int *edgeFlags=0, HVNormal3D *normals=0) const
Definition SbPolyhedron.cxx:1174
SbPolyhedron::GetNormal
HVNormal3D GetNormal(int iFace) const
Definition SbPolyhedron.cxx:1204
SbPolyhedron::SetNumberOfRotationSteps
static void SetNumberOfRotationSteps(int n)
Definition SbPolyhedron.cxx:246
r
int r
Definition globals.cxx:22
const
index
Definition index.py:1
Generated on for ATLAS Offline Software by doxygen 1.14.0