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
16#define perMillion 0.000001
17#define deg (M_PI/180.0)
18
19#include <float.h> //G.Barrand : to have DBL_EPSILON on Windows.
20
21// G.Barrand : introduce iabs to avoid a mess with cmath and some compiler.
22inline int iabs(int a) {
23 return a < 0 ? -a : a;
24}
25inline float ffabs(float a) {
26 return a < 0.0f ? -a : a;
27}
28
29
30// rbianchi
31HVPoint3D::HVPoint3D(): SbVec3d(0,0,0){}
32HVPoint3D::HVPoint3D(double x,double y,double z): SbVec3d(x,y,z){}
33HVPoint3D::HVPoint3D(const HVPoint3D& v): SbVec3d(v){}
34HVPoint3D::HVPoint3D(const SbVec3d& v): SbVec3d(v){}
35HVPoint3D& HVPoint3D::operator=(const HVPoint3D& v) {
36 setValue(v[0],v[1],v[2]);
37 return *this;
38}
39HVPoint3D& HVPoint3D::operator=(const SbVec3d& v) {
40 setValue(v[0],v[1],v[2]);
41 return *this;
42}
43HVPoint3D operator +(const HVPoint3D& v1,const HVPoint3D& v2) {
44 return HVPoint3D(v1[0] + v2[0],v1[1] + v2[1],v1[2] + v2[2]);
45}
46//---
47
48
49
50//--------------------------------------------------------------------//
51// JFB: //
52// SbPolyhedron was SbPolyhedron, retrofitted to Open Inventor //
53// infrastructure: //
54//--------------------------------------------------------------------//
55
56
57//
58// ********************************************************************
59// * DISCLAIMER *
60// * *
61// * The following disclaimer summarizes all the specific disclaimers *
62// * of contributors to this software. The specific disclaimers,which *
63// * govern, are listed with their locations in: *
64// * http://cern.ch/geant4/license *
65// * *
66// * Neither the authors of this software system, nor their employing *
67// * institutes,nor the agencies providing financial support for this *
68// * work make any representation or warranty, express or implied, *
69// * regarding this software system or assume any liability for its *
70// * use. *
71// * *
72// * This code implementation is the intellectual property of the *
73// * GEANT4 collaboration. *
74// * By copying, distributing or modifying the Program (or any work *
75// * based on the Program) you indicate your acceptance of this *
76// * statement, and all its terms. *
77// ********************************************************************
78//
79//
80// $Id: SbPolyhedron.cxx,v 1.3 2008-09-14 19:04:40 tkittel Exp $
81// GEANT4 tag $Name: not supported by cvs2svn $
82//
83//
84//
85// G4 Polyhedron library
86//
87// History:
88// 23.07.96 E.Chernyaev <Evgueni.Tcherniaev@cern.ch> - initial version
89//
90// 30.09.96 E.Chernyaev
91// - added GetNextVertexIndex, GetVertex by Yasuhide Sawada
92// - added GetNextUnitNormal, GetNextEdgeIndeces, GetNextEdge
93//
94// 15.12.96 E.Chernyaev
95// - added GetNumberOfRotationSteps, RotateEdge, RotateAroundZ, SetReferences
96// - rewritten G4PolyhedronCons;
97// - added G4PolyhedronPara, ...Trap, ...Pgon, ...Pcon, ...Sphere, ...Torus
98//
99// 01.06.97 E.Chernyaev
100// - modified RotateAroundZ, added SetSideFacets
101//
102// 19.03.00 E.Chernyaev
103// - implemented boolean operations (add, subtract, intersect) on polyhedra;
104//
105// 25.05.01 E.Chernyaev
106// - added GetSurfaceArea() and GetVolume();
107//
108
109
110/***********************************************************************
111 * *
112 * Name: SbPolyhedron operator << Date: 09.05.96 *
113 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
114 * *
115 * Function: Print contents of G4 polyhedron *
116 * *
117 ***********************************************************************/
118std::ostream & operator<<(std::ostream & ostr, const SbFacet & facet) {
119 for (int k=0; k<4; k++) {
120 ostr << " " << facet.m_edge[k].v << "/" << facet.m_edge[k].f;
121 }
122 return ostr;
123}
124
125std::ostream & operator<<(std::ostream & ostr, const SbPolyhedron & ph) {
126 ostr << std::endl;
127 ostr << "Nverteces=" << ph.m_nvert << ", Nfacets=" << ph.m_nface << std::endl;
128 int i;
129 for (i=1; i<=ph.m_nvert; i++) {
130 ostr << "xyz(" << i << ")="
131 << ph.m_pV[i][0] << ' ' << ph.m_pV[i][1] << ' ' << ph.m_pV[i][2]
132 << std::endl;
133 }
134 for (i=1; i<=ph.m_nface; i++) {
135 ostr << "face(" << i << ")=" << ph.m_pF[i] << std::endl;
136 }
137 return ostr;
138}
139
140SbPolyhedron::SbPolyhedron(const SbPolyhedron &from)
141/***********************************************************************
142 * *
143 * Name: SbPolyhedron copy constructor Date: 23.07.96 *
144 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
145 * *
146 ***********************************************************************/
147{
148 if (from.m_nvert > 0 && from.m_nface > 0) {
149 m_nvert = from.m_nvert;
150 m_nface = from.m_nface;
151 m_pV = new HVPoint3D[m_nvert + 1];
152 m_pF = new SbFacet[m_nface + 1];
153 int i;
154 for (i=1; i<=m_nvert; i++) m_pV[i] = from.m_pV[i];
155 for (i=1; i<=m_nface; i++) m_pF[i] = from.m_pF[i];
156 }else{
157 m_nvert = 0; m_nface = 0; m_pV = 0; m_pF = 0;
158 }
159}
160
161SbPolyhedron & SbPolyhedron::operator=(const SbPolyhedron &from)
162/***********************************************************************
163 * *
164 * Name: SbPolyhedron operator = Date: 23.07.96 *
165 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
166 * *
167 * Function: Copy contents of one GEANT4 polyhedron to another *
168 * *
169 ***********************************************************************/
170{
171 if (this == &from) return *this;
172 delete [] m_pV;
173 delete [] m_pF;
174 if (from.m_nvert > 0 && from.m_nface > 0) {
175 m_nvert = from.m_nvert;
176 m_nface = from.m_nface;
177 m_pV = new HVPoint3D[m_nvert + 1];
178 m_pF = new SbFacet[m_nface + 1];
179 int i;
180 for (i=1; i<=m_nvert; i++) m_pV[i] = from.m_pV[i];
181 for (i=1; i<=m_nface; i++) m_pF[i] = from.m_pF[i];
182 }else{
183 m_nvert = 0; m_nface = 0; m_pV = 0; m_pF = 0;
184 }
185 return *this;
186}
187
188int
189SbPolyhedron::FindNeighbour(int iFace, int iNode, int iOrder) const
190/***********************************************************************
191 * *
192 * Name: SbPolyhedron::FindNeighbour Date: 22.11.99 *
193 * Author: E.Chernyaev Revised: *
194 * *
195 * Function: Find neighbouring face *
196 * *
197 ***********************************************************************/
198{
199 int i;
200 for (i=0; i<4; i++) {
201 if (iNode == iabs(m_pF[iFace].m_edge[i].v)) break;
202 }
203 if (i == 4) {
204 std::cerr
205 << "SbPolyhedron::FindNeighbour: face " << iFace
206 << " has no node " << iNode
207 << std::endl;
208 return 0;
209 }
210 if (iOrder < 0) {
211 if ( --i < 0) i = 3;
212 if (m_pF[iFace].m_edge[i].v == 0) i = 2;
213 }
214 return (m_pF[iFace].m_edge[i].v > 0) ? 0 : m_pF[iFace].m_edge[i].f;
215}
216
217HVNormal3D SbPolyhedron::FindNodeNormal(int iFace, int iNode) const
218/***********************************************************************
219 * *
220 * Name: SbPolyhedron::FindNodeNormal Date: 22.11.99 *
221 * Author: E.Chernyaev Revised: *
222 * *
223 * Function: Find normal at given node *
224 * *
225 ***********************************************************************/
226{
227 HVNormal3D normal = GetUnitNormal(iFace);
228 int k = iFace, iOrder = 1;
229
230 for(;;) {
231 k = FindNeighbour(k, iNode, iOrder);
232 if (k == iFace) break;
233 if (k > 0) {
234 normal += GetUnitNormal(k);
235 }else{
236 if (iOrder < 0) break;
237 k = iFace;
238 iOrder = -iOrder;
239 }
240 }
241 normal.normalize();
242 return normal;
243}
244
245void SbPolyhedron::SetNumberOfRotationSteps(int n)
246/***********************************************************************
247 * *
248 * Name: SbPolyhedron::SetNumberOfRotationSteps Date: 24.06.97 *
249 * Author: J.Allison (Manchester University) Revised: *
250 * *
251 * Function: Set number of steps for whole circle *
252 * *
253 ***********************************************************************/
254{
255 const int nMin = 3;
256 if (n < nMin) {
257 std::cerr
258 << "SbPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
259 << "number of steps per circle < " << nMin << "; forced to " << nMin
260 << std::endl;
261 s_numberOfRotationSteps = nMin;
262 }else{
263 s_numberOfRotationSteps = n;
264 }
265}
266
267void SbPolyhedron::AllocateMemory(int Nvert, int Nface)
268/***********************************************************************
269 * *
270 * Name: SbPolyhedron::AllocateMemory Date: 19.06.96 *
271 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
272 * *
273 * Function: Allocate memory for GEANT4 polyhedron *
274 * *
275 * Input: Nvert - number of nodes *
276 * Nface - number of faces *
277 * *
278 ***********************************************************************/
279{
280 m_nvert = Nvert;
281 m_nface = Nface;
282 m_pV = new HVPoint3D[m_nvert+1];
283 m_pF = new SbFacet[m_nface+1];
284}
285
286void SbPolyhedron::CreatePrism()
287/***********************************************************************
288 * *
289 * Name: SbPolyhedron::CreatePrism Date: 15.07.96 *
290 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
291 * *
292 * Function: Set facets for a prism *
293 * *
294 ***********************************************************************/
295{
296 enum {DUMMY, BOTTOM, LEFT, BACK, RIGHT, FRONT, TOP};
297
298 m_pF[1] = SbFacet(1,LEFT, 4,BACK, 3,RIGHT, 2,FRONT);
299 m_pF[2] = SbFacet(5,TOP, 8,BACK, 4,BOTTOM, 1,FRONT);
300 m_pF[3] = SbFacet(8,TOP, 7,RIGHT, 3,BOTTOM, 4,LEFT);
301 m_pF[4] = SbFacet(7,TOP, 6,FRONT, 2,BOTTOM, 3,BACK);
302 m_pF[5] = SbFacet(6,TOP, 5,LEFT, 1,BOTTOM, 2,RIGHT);
303 m_pF[6] = SbFacet(5,FRONT, 6,RIGHT, 7,BACK, 8,LEFT);
304}
305
306void SbPolyhedron::RotateEdge(int k1, int k2, double r1, double r2,
307 int v1, int v2, int vEdge,
308 bool ifWholeCircle, int ns, int &kface)
309/***********************************************************************
310 * *
311 * Name: SbPolyhedron::RotateEdge Date: 05.12.96 *
312 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
313 * *
314 * Function: Create set of facets by rotation of an edge around Z-axis *
315 * *
316 * Input: k1, k2 - end vertices of the edge *
317 * r1, r2 - radiuses of the end vertices *
318 * v1, v2 - visibility of edges produced by rotation of the end *
319 * vertices *
320 * vEdge - visibility of the edge *
321 * ifWholeCircle - is true in case of whole circle rotation *
322 * ns - number of discrete steps *
323 * r[] - r-coordinates *
324 * kface - current free cell in the m_pF array *
325 * *
326 ***********************************************************************/
327{
328 if (r1 == 0. && r2 == 0) return;
329
330 int i;
331 int i1 = k1;
332 int i2 = k2;
333 int ii1 = ifWholeCircle ? i1 : i1+ns;
334 int ii2 = ifWholeCircle ? i2 : i2+ns;
335 int vv = ifWholeCircle ? vEdge : 1;
336
337 if (ns == 1) {
338 if (r1 == 0.) {
339 m_pF[kface++] = SbFacet(i1,0, v2*i2,0, (i2+1),0);
340 }else if (r2 == 0.) {
341 m_pF[kface++] = SbFacet(i1,0, i2,0, v1*(i1+1),0);
342 }else{
343 m_pF[kface++] = SbFacet(i1,0, v2*i2,0, (i2+1),0, v1*(i1+1),0);
344 }
345 }else{
346 if (r1 == 0.) {
347 m_pF[kface++] = SbFacet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0);
348 for (i2++,i=1; i<ns-1; i2++,i++) {
349 m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0);
350 }
351 m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vv*ii2,0);
352 }else if (r2 == 0.) {
353 m_pF[kface++] = SbFacet(vv*i1,0, vEdge*i2,0, v1*(i1+1),0);
354 for (i1++,i=1; i<ns-1; i1++,i++) {
355 m_pF[kface++] = SbFacet(vEdge*i1,0, vEdge*i2,0, v1*(i1+1),0);
356 }
357 m_pF[kface++] = SbFacet(vEdge*i1,0, vv*i2,0, v1*ii1,0);
358 }else{
359 m_pF[kface++] = SbFacet(vv*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
360 for (i1++,i2++,i=1; i<ns-1; i1++,i2++,i++) {
361 m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vEdge*(i2+1),0,v1*(i1+1),0);
362 }
363 m_pF[kface++] = SbFacet(vEdge*i1,0, v2*i2,0, vv*ii2,0, v1*ii1,0);
364 }
365 }
366}
367
368void SbPolyhedron::SetSideFacets(int ii[4], int vv[4],
369 int *kk, double *r,
370 double dphi, int ns, int &kface)
371/***********************************************************************
372 * *
373 * Name: SbPolyhedron::SetSideFacets Date: 20.05.97 *
374 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
375 * *
376 * Function: Set side facets for the case of incomplete rotation *
377 * *
378 * Input: ii[4] - indeces of original verteces *
379 * vv[4] - visibility of edges *
380 * kk[] - indeces of nodes *
381 * r[] - radiuses *
382 * dphi - delta phi *
383 * ns - number of discrete steps *
384 * kface - current free cell in the m_pF array *
385 * *
386 ***********************************************************************/
387{
388 int k1, k2, k3, k4;
389 if (fabs(dphi-M_PI) < perMillion) { // half a circle
390 for (int i=0; i<4; i++) {
391 k1 = ii[i];
392 k2 = (i == 3) ? ii[0] : ii[i+1];
393 if (r[k1] == 0. && r[k2] == 0.) vv[i] = -1;
394 }
395 }
396
397 if (ii[1] == ii[2]) {
398 k1 = kk[ii[0]];
399 k2 = kk[ii[2]];
400 k3 = kk[ii[3]];
401 m_pF[kface++] = SbFacet(vv[0]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
402 if (r[ii[0]] != 0.) k1 += ns;
403 if (r[ii[2]] != 0.) k2 += ns;
404 if (r[ii[3]] != 0.) k3 += ns;
405 m_pF[kface++] = SbFacet(vv[2]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
406 }else if (kk[ii[0]] == kk[ii[1]]) {
407 k1 = kk[ii[0]];
408 k2 = kk[ii[2]];
409 k3 = kk[ii[3]];
410 m_pF[kface++] = SbFacet(vv[1]*k1,0, vv[2]*k2,0, vv[3]*k3,0);
411 if (r[ii[0]] != 0.) k1 += ns;
412 if (r[ii[2]] != 0.) k2 += ns;
413 if (r[ii[3]] != 0.) k3 += ns;
414 m_pF[kface++] = SbFacet(vv[2]*k3,0, vv[1]*k2,0, vv[3]*k1,0);
415 }else if (kk[ii[2]] == kk[ii[3]]) {
416 k1 = kk[ii[0]];
417 k2 = kk[ii[1]];
418 k3 = kk[ii[2]];
419 m_pF[kface++] = SbFacet(vv[0]*k1,0, vv[1]*k2,0, vv[3]*k3,0);
420 if (r[ii[0]] != 0.) k1 += ns;
421 if (r[ii[1]] != 0.) k2 += ns;
422 if (r[ii[2]] != 0.) k3 += ns;
423 m_pF[kface++] = SbFacet(vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
424 }else{
425 k1 = kk[ii[0]];
426 k2 = kk[ii[1]];
427 k3 = kk[ii[2]];
428 k4 = kk[ii[3]];
429 m_pF[kface++] = SbFacet(vv[0]*k1,0, vv[1]*k2,0, vv[2]*k3,0, vv[3]*k4,0);
430 if (r[ii[0]] != 0.) k1 += ns;
431 if (r[ii[1]] != 0.) k2 += ns;
432 if (r[ii[2]] != 0.) k3 += ns;
433 if (r[ii[3]] != 0.) k4 += ns;
434 m_pF[kface++] = SbFacet(vv[2]*k4,0, vv[1]*k3,0, vv[0]*k2,0, vv[3]*k1,0);
435 }
436}
437
438void SbPolyhedron::RotateAroundZ(int nstep, double phi, double dphi,
439 int np1, int np2,
440 const double *z, double *r,
441 int nodeVis, int edgeVis)
442/***********************************************************************
443 * *
444 * Name: SbPolyhedron::RotateAroundZ Date: 27.11.96 *
445 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
446 * *
447 * Function: Create SbPolyhedron for a solid produced by rotation of *
448 * two polylines around Z-axis *
449 * *
450 * Input: nstep - number of discrete steps, if 0 then default *
451 * phi - starting phi angle *
452 * dphi - delta phi *
453 * np1 - number of points in external polyline *
454 * (must be negative in case of closed polyline) *
455 * np2 - number of points in internal polyline (may be 1) *
456 * z[] - z-coordinates (+z >>> -z for both polylines) *
457 * r[] - r-coordinates *
458 * nodeVis - how to Draw edges joing consecutive positions of *
459 * node during rotation *
460 * edgeVis - how to Draw edges *
461 * *
462 ***********************************************************************/
463{
464 static double wholeCircle = 2*M_PI;
465
466 // S E T R O T A T I O N P A R A M E T E R S
467
468 bool ifWholeCircle = (fabs(dphi-wholeCircle) < perMillion) ?
469 true : false;
470 double delPhi = ifWholeCircle ? wholeCircle : dphi;
471 int nSphi = (nstep > 0) ?
472 nstep : int(delPhi*GetNumberOfRotationSteps()/wholeCircle+.5);
473 if (nSphi == 0) nSphi = 1;
474 int nVphi = ifWholeCircle ? nSphi : nSphi+1;
475 bool ifClosed = np1 > 0 ? false : true;
476
477 // C O U N T V E R T E C E S
478
479 int absNp1 = iabs(np1);
480 int absNp2 = iabs(np2);
481 int i1beg = 0;
482 int i1end = absNp1-1;
483 int i2beg = absNp1;
484 int i2end = absNp1+absNp2-1;
485 int i, j, k;
486
487 for(i=i1beg; i<=i2end; i++) {
488 if (fabs(r[i]) < perMillion) r[i] = 0.;
489 }
490
491 j = 0; // external nodes
492 for (i=i1beg; i<=i1end; i++) {
493 j += (r[i] == 0.) ? 1 : nVphi;
494 }
495
496 bool ifSide1 = false; // internal nodes
497 bool ifSide2 = false;
498
499 if (r[i2beg] != r[i1beg] || z[i2beg] != z[i1beg]) {
500 j += (r[i2beg] == 0.) ? 1 : nVphi;
501 ifSide1 = true;
502 }
503
504 for(i=i2beg+1; i<i2end; i++) {
505 j += (r[i] == 0.) ? 1 : nVphi;
506 }
507
508 if (r[i2end] != r[i1end] || z[i2end] != z[i1end]) {
509 if (absNp2 > 1) j += (r[i2end] == 0.) ? 1 : nVphi;
510 ifSide2 = true;
511 }
512
513 // C O U N T F A C E S
514
515 k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi; // external faces
516
517 if (absNp2 > 1) { // internal faces
518 for(i=i2beg; i<i2end; i++) {
519 if (r[i] > 0. || r[i+1] > 0.) k += nSphi;
520 }
521
522 if (ifClosed) {
523 if (r[i2end] > 0. || r[i2beg] > 0.) k += nSphi;
524 }
525 }
526
527 if (!ifClosed) { // side faces
528 if (ifSide1 && (r[i1beg] > 0. || r[i2beg] > 0.)) k += nSphi;
529 if (ifSide2 && (r[i1end] > 0. || r[i2end] > 0.)) k += nSphi;
530 }
531
532 if (!ifWholeCircle) { // phi_side faces
533 k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
534 }
535
536 // A L L O C A T E M E M O R Y
537
538 AllocateMemory(j, k);
539
540 // G E N E R A T E V E R T E C E S
541
542 int *kk;
543 kk = new int[absNp1+absNp2];
544
545 k = 1;
546 for(i=i1beg; i<=i1end; i++) {
547 kk[i] = k;
548 if (r[i] == 0.) { m_pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
549 }
550
551 i = i2beg;
552 if (ifSide1) {
553 kk[i] = k;
554 if (r[i] == 0.) { m_pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
555 }else{
556 kk[i] = kk[i1beg];
557 }
558
559 for(i=i2beg+1; i<i2end; i++) {
560 kk[i] = k;
561 if (r[i] == 0.) { m_pV[k++] = HVPoint3D(0, 0, z[i]); } else { k += nVphi; }
562 }
563
564 if (absNp2 > 1) {
565 i = i2end;
566 if (ifSide2) {
567 kk[i] = k;
568 if (r[i] == 0.) m_pV[k] = HVPoint3D(0, 0, z[i]);
569 }else{
570 kk[i] = kk[i1end];
571 }
572 }
573
574 double cosPhi, sinPhi;
575
576 for(j=0; j<nVphi; j++) {
577 cosPhi = cos(phi+j*delPhi/nSphi);
578 sinPhi = sin(phi+j*delPhi/nSphi);
579 for(i=i1beg; i<=i2end; i++) {
580 if (r[i] != 0.) m_pV[kk[i]+j] = HVPoint3D(r[i]*cosPhi,r[i]*sinPhi,z[i]);
581 }
582 }
583
584 // G E N E R A T E E X T E R N A L F A C E S
585
586 int v1,v2;
587
588 k = 1;
589 v2 = ifClosed ? nodeVis : 1;
590 for(i=i1beg; i<i1end; i++) {
591 v1 = v2;
592 if (!ifClosed && i == i1end-1) {
593 v2 = 1;
594 }else{
595 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
596 }
597 RotateEdge(kk[i], kk[i+1], r[i], r[i+1], v1, v2,
598 edgeVis, ifWholeCircle, nSphi, k);
599 }
600 if (ifClosed) {
601 RotateEdge(kk[i1end], kk[i1beg], r[i1end],r[i1beg], nodeVis, nodeVis,
602 edgeVis, ifWholeCircle, nSphi, k);
603 }
604
605 // G E N E R A T E I N T E R N A L F A C E S
606
607 if (absNp2 > 1) {
608 v2 = ifClosed ? nodeVis : 1;
609 for(i=i2beg; i<i2end; i++) {
610 v1 = v2;
611 if (!ifClosed && i==i2end-1) {
612 v2 = 1;
613 }else{
614 v2 = (r[i] == r[i+1] && r[i+1] == r[i+2]) ? -1 : nodeVis;
615 }
616 RotateEdge(kk[i+1], kk[i], r[i+1], r[i], v2, v1,
617 edgeVis, ifWholeCircle, nSphi, k);
618 }
619 if (ifClosed) {
620 RotateEdge(kk[i2beg], kk[i2end], r[i2beg], r[i2end], nodeVis, nodeVis,
621 edgeVis, ifWholeCircle, nSphi, k);
622 }
623 }
624
625 // G E N E R A T E S I D E F A C E S
626
627 if (!ifClosed) {
628 if (ifSide1) {
629 RotateEdge(kk[i2beg], kk[i1beg], r[i2beg], r[i1beg], 1, 1,
630 -1, ifWholeCircle, nSphi, k);
631 }
632 if (ifSide2) {
633 RotateEdge(kk[i1end], kk[i2end], r[i1end], r[i2end], 1, 1,
634 -1, ifWholeCircle, nSphi, k);
635 }
636 }
637
638 // G E N E R A T E S I D E F A C E S for the case of incomplete circle
639
640 if (!ifWholeCircle) {
641
642 int ii[4], vv[4];
643
644 if (ifClosed) {
645 for (i=i1beg; i<=i1end; i++) {
646 ii[0] = i;
647 ii[3] = (i == i1end) ? i1beg : i+1;
648 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
649 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
650 vv[0] = -1;
651 vv[1] = 1;
652 vv[2] = -1;
653 vv[3] = 1;
654 SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
655 }
656 }else{
657 for (i=i1beg; i<i1end; i++) {
658 ii[0] = i;
659 ii[3] = i+1;
660 ii[1] = (absNp2 == 1) ? i2beg : ii[0]+absNp1;
661 ii[2] = (absNp2 == 1) ? i2beg : ii[3]+absNp1;
662 vv[0] = (i == i1beg) ? 1 : -1;
663 vv[1] = 1;
664 vv[2] = (i == i1end-1) ? 1 : -1;
665 vv[3] = 1;
666 SetSideFacets(ii, vv, kk, r, dphi, nSphi, k);
667 }
668 }
669 }
670
671 delete [] kk;
672
673 if (k-1 != m_nface) {
674 std::cerr
675 << "Polyhedron::RotateAroundZ: number of generated faces ("
676 << k-1 << ") is not equal to the number of allocated faces ("
677 << m_nface << ")"
678 << std::endl;
679 }
680}
681
682void SbPolyhedron::SetReferences()
683/***********************************************************************
684 * *
685 * Name: SbPolyhedron::SetReferences Date: 04.12.96 *
686 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
687 * *
688 * Function: For each edge set reference to neighbouring facet *
689 * *
690 ***********************************************************************/
691{
692 if (m_nface <= 0) return;
693
694 struct edgeListMember {
695 edgeListMember *next;
696 int v2;
697 int iface;
698 int iedge;
699 } *edgeList, *freeList, **headList;
700
701
702 // A L L O C A T E A N D I N I T I A T E L I S T S
703
704 edgeList = new edgeListMember[2*m_nface];
705 headList = new edgeListMember*[m_nvert];
706
707 int i;
708 for (i=0; i<m_nvert; i++) {
709 headList[i] = 0;
710 }
711 freeList = edgeList;
712 for (i=0; i<2*m_nface-1; i++) {
713 edgeList[i].next = &edgeList[i+1];
714 }
715 edgeList[2*m_nface-1].next = 0;
716
717 // L O O P A L O N G E D G E S
718
719 int iface, iedge, nedge, i1, i2, k1, k2;
720 edgeListMember *prev, *cur;
721
722 for(iface=1; iface<=m_nface; iface++) {
723 nedge = (m_pF[iface].m_edge[3].v == 0) ? 3 : 4;
724 for (iedge=0; iedge<nedge; iedge++) {
725 i1 = iedge;
726 i2 = (iedge < nedge-1) ? iedge+1 : 0;
727 i1 = iabs(m_pF[iface].m_edge[i1].v);
728 i2 = iabs(m_pF[iface].m_edge[i2].v);
729 k1 = (i1 < i2) ? i1 : i2; // k1 = ::min(i1,i2);
730 k2 = (i1 > i2) ? i1 : i2; // k2 = ::max(i1,i2);
731
732 // check head of the List corresponding to k1
733 cur = headList[k1];
734 if (cur == 0) {
735 headList[k1] = freeList;
736 freeList = freeList->next;
737 cur = headList[k1];
738 cur->next = 0;
739 cur->v2 = k2;
740 cur->iface = iface;
741 cur->iedge = iedge;
742 continue;
743 }
744
745 if (cur->v2 == k2) {
746 headList[k1] = cur->next;
747 cur->next = freeList;
748 freeList = cur;
749 m_pF[iface].m_edge[iedge].f = cur->iface;
750 m_pF[cur->iface].m_edge[cur->iedge].f = iface;
751 i1 = (m_pF[iface].m_edge[iedge].v < 0) ? -1 : 1;
752 i2 = (m_pF[cur->iface].m_edge[cur->iedge].v < 0) ? -1 : 1;
753 if (i1 != i2) {
754 std::cerr
755 << "Polyhedron::SetReferences: different edge visibility "
756 << iface << "/" << iedge << "/"
757 << m_pF[iface].m_edge[iedge].v << " and "
758 << cur->iface << "/" << cur->iedge << "/"
759 << m_pF[cur->iface].m_edge[cur->iedge].v
760 << std::endl;
761 }
762 continue;
763 }
764
765 // check List itself
766 for (;;) {
767 prev = cur;
768 cur = prev->next;
769 if (cur == 0) {
770 prev->next = freeList;
771 freeList = freeList->next;
772 cur = prev->next;
773 cur->next = 0;
774 cur->v2 = k2;
775 cur->iface = iface;
776 cur->iedge = iedge;
777 break;
778 }
779
780 if (cur->v2 == k2) {
781 prev->next = cur->next;
782 cur->next = freeList;
783 freeList = cur;
784 m_pF[iface].m_edge[iedge].f = cur->iface;
785 m_pF[cur->iface].m_edge[cur->iedge].f = iface;
786 i1 = (m_pF[iface].m_edge[iedge].v < 0) ? -1 : 1;
787 i2 = (m_pF[cur->iface].m_edge[cur->iedge].v < 0) ? -1 : 1;
788 if (i1 != i2) {
789 std::cerr
790 << "Polyhedron::SetReferences: different edge visibility "
791 << iface << "/" << iedge << "/"
792 << m_pF[iface].m_edge[iedge].v << " and "
793 << cur->iface << "/" << cur->iedge << "/"
794 << m_pF[cur->iface].m_edge[cur->iedge].v
795 << std::endl;
796 }
797 break;
798 }
799 }
800 }
801 }
802
803 // 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
804
805 for (i=0; i<m_nvert; i++) {
806 if (headList[i] != 0) {
807 std::cerr
808 << "Polyhedron::SetReferences: List " << i << " is not empty"
809 << std::endl;
810 }
811 }
812
813 // F R E E M E M O R Y
814
815 delete [] edgeList;
816 delete [] headList;
817}
818
819void SbPolyhedron::InvertFacets()
820/***********************************************************************
821 * *
822 * Name: SbPolyhedron::InvertFacets Date: 01.12.99 *
823 * Author: E.Chernyaev Revised: *
824 * *
825 * Function: Invert the order of the nodes in the facets *
826 * *
827 ***********************************************************************/
828{
829 if (m_nface <= 0) return;
830 int i, k, nnode, v[4],f[4];
831 for (i=1; i<=m_nface; i++) {
832 nnode = (m_pF[i].m_edge[3].v == 0) ? 3 : 4;
833 for (k=0; k<nnode; k++) {
834 v[k] = (k+1 == nnode) ? m_pF[i].m_edge[0].v : m_pF[i].m_edge[k+1].v;
835 if (v[k] * m_pF[i].m_edge[k].v < 0) v[k] = -v[k];
836 f[k] = m_pF[i].m_edge[k].f;
837 }
838 for (k=0; k<nnode; k++) {
839 m_pF[i].m_edge[nnode-1-k].v = v[k];
840 m_pF[i].m_edge[nnode-1-k].f = f[k];
841 }
842 }
843}
844
845
846// rbianchi change
847//- original transform method
848/*
849SbPolyhedron & SbPolyhedron::Transform(const HVRotation & rotation, const HVVector3D & translation)
850// ***********************************************************************
851// * *
852// * Name: SbPolyhedron::Transform Date: 01.12.99 *
853// * Author: E.Chernyaev Revised: *
854// * *
855// * Function: Make transformation of the polyhedron *
856// * *
857// ***********************************************************************
858{
859 if (m_nvert > 0) {
860 for (int i=1; i<=m_nvert; i++) {
861 HVVector3D tmp;
862 rotation.multVec(m_pV[i],tmp);
863 m_pV[i] = tmp+translation;
864 }
865
866 // C H E C K D E T E R M I N A N T A N D
867 // 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
868
869 HVVector3D x; rotation.multVec(HVVector3D(1,0,0),x);
870 HVVector3D y; rotation.multVec(HVVector3D(0,1,0),y);
871 HVVector3D z; rotation.multVec(HVVector3D(0,0,1),z);
872 if ((x.cross(y)).dot(z) < 0) InvertFacets();
873 }
874 return *this;
875}
876*/
877//- new transform method from the newest OpenScientist/HEPVis
878SbPolyhedron & SbPolyhedron::Transform(
879 const HEPVis::SbRotation& rotation
880,const SbVec3d& translation
881)
882{
883 if (m_nvert > 0) {
884 SbVec3d tmp;
885 for (int i=1; i<=m_nvert; i++) {
886 rotation.multVec(m_pV[i],tmp);
887 m_pV[i] = tmp+translation;
888 }
889 SbVec3d x; rotation.multVec(SbVec3d(1,0,0),x);
890 SbVec3d y; rotation.multVec(SbVec3d(0,1,0),y);
891 SbVec3d z; rotation.multVec(SbVec3d(0,0,1),z);
892 if ((x.cross(y)).dot(z) < 0) InvertFacets();
893 }
894 return *this;
895}
896//---
897
898
899
900
901
902bool SbPolyhedron::GetNextVertexIndex(int &index, int &edgeFlag) const
903/***********************************************************************
904 * *
905 * Name: SbPolyhedron::GetNextVertexIndex Date: 03.09.96 *
906 * Author: Yasuhide Sawada Revised: *
907 * *
908 * Function: *
909 * *
910 ***********************************************************************/
911{
912 static int iFace = 1;
913 static int iQVertex = 0;
914 int vIndex = m_pF[iFace].m_edge[iQVertex].v;
915
916 edgeFlag = (vIndex > 0) ? 1 : 0;
917 index = iabs(vIndex);
918
919 if (iQVertex >= 3 || m_pF[iFace].m_edge[iQVertex+1].v == 0) {
920 iQVertex = 0;
921 if (++iFace > m_nface) iFace = 1;
922 return false; // Last Edge
923 }else{
924 ++iQVertex;
925 return true; // not Last Edge
926 }
927}
928
929HVPoint3D SbPolyhedron::GetVertex(int index) const
930/***********************************************************************
931 * *
932 * Name: SbPolyhedron::GetVertex Date: 03.09.96 *
933 * Author: Yasuhide Sawada Revised: 17.11.99 *
934 * *
935 * Function: Get vertex of the index. *
936 * *
937 ***********************************************************************/
938{
939 if (index <= 0 || index > m_nvert) {
940 std::cerr
941 << "SbPolyhedron::GetVertex: irrelevant index " << index
942 << std::endl;
943 return HVPoint3D();
944 }
945 return m_pV[index];
946}
947
948
949// rbianchi - 14.12.2012
950const HVPoint3D& SbPolyhedron::GetVertexFast(int index) const { //G.Barrand
951 return m_pV[index];
952}
953//---
954
955
956bool
957SbPolyhedron::GetNextVertex(HVPoint3D &vertex, int &edgeFlag) const
958/***********************************************************************
959 * *
960 * Name: SbPolyhedron::GetNextVertex Date: 22.07.96 *
961 * Author: John Allison Revised: *
962 * *
963 * Function: Get vertices of the quadrilaterals in order for each *
964 * face in face order. Returns false when finished each *
965 * face. *
966 * *
967 ***********************************************************************/
968{
969 int index;
970 bool rep = GetNextVertexIndex(index, edgeFlag);
971 vertex = m_pV[index];
972 return rep;
973}
974
975bool SbPolyhedron::GetNextVertex(HVPoint3D &vertex, int &edgeFlag,
976 HVNormal3D &normal) const
977/***********************************************************************
978 * *
979 * Name: SbPolyhedron::GetNextVertex Date: 26.11.99 *
980 * Author: E.Chernyaev Revised: *
981 * *
982 * Function: Get vertices with normals of the quadrilaterals in order *
983 * for each face in face order. *
984 * Returns false when finished each face. *
985 * *
986 ***********************************************************************/
987{
988 static int iFace = 1;
989 static int iNode = 0;
990
991 if (m_nface == 0) return false; // empty polyhedron
992
993 int k = m_pF[iFace].m_edge[iNode].v;
994 if (k > 0) { edgeFlag = 1; } else { edgeFlag = -1; k = -k; }
995 vertex = m_pV[k];
996 normal = FindNodeNormal(iFace,k);
997 if (iNode >= 3 || m_pF[iFace].m_edge[iNode+1].v == 0) {
998 iNode = 0;
999 if (++iFace > m_nface) iFace = 1;
1000 return false; // last node
1001 }else{
1002 ++iNode;
1003 return true; // not last node
1004 }
1005}
1006
1007bool SbPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag,
1008 int &iface1, int &iface2) const
1009/***********************************************************************
1010 * *
1011 * Name: SbPolyhedron::GetNextEdgeIndeces Date: 30.09.96 *
1012 * Author: E.Chernyaev Revised: 17.11.99 *
1013 * *
1014 * Function: Get indeces of the next edge together with indeces of *
1015 * of the faces which share the edge. *
1016 * Returns false when the last edge. *
1017 * *
1018 ***********************************************************************/
1019{
1020 static int iFace = 1;
1021 static int iQVertex = 0;
1022 static int iOrder = 1;
1023 int k1, k2, kflag, kface1, kface2;
1024
1025 if (iFace == 1 && iQVertex == 0) {
1026 k2 = m_pF[m_nface].m_edge[0].v;
1027 k1 = m_pF[m_nface].m_edge[3].v;
1028 if (k1 == 0) k1 = m_pF[m_nface].m_edge[2].v;
1029 if (iabs(k1) > iabs(k2)) iOrder = -1;
1030 }
1031
1032 do {
1033 k1 = m_pF[iFace].m_edge[iQVertex].v;
1034 kflag = k1;
1035 k1 = iabs(k1);
1036 kface1 = iFace;
1037 kface2 = m_pF[iFace].m_edge[iQVertex].f;
1038 if (iQVertex >= 3 || m_pF[iFace].m_edge[iQVertex+1].v == 0) {
1039 iQVertex = 0;
1040 k2 = iabs(m_pF[iFace].m_edge[iQVertex].v);
1041 iFace++;
1042 }else{
1043 iQVertex++;
1044 k2 = iabs(m_pF[iFace].m_edge[iQVertex].v);
1045 }
1046 } while (iOrder*k1 > iOrder*k2);
1047
1048 i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1049 iface1 = kface1; iface2 = kface2;
1050
1051 if (iFace > m_nface) {
1052 iFace = 1; iOrder = 1;
1053 return false;
1054 }else{
1055 return true;
1056 }
1057}
1058
1059bool
1060SbPolyhedron::GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag) const
1061/***********************************************************************
1062 * *
1063 * Name: SbPolyhedron::GetNextEdgeIndeces Date: 17.11.99 *
1064 * Author: E.Chernyaev Revised: *
1065 * *
1066 * Function: Get indeces of the next edge. *
1067 * Returns false when the last edge. *
1068 * *
1069 ***********************************************************************/
1070{
1071 int kface1, kface2;
1072 return GetNextEdgeIndeces(i1, i2, edgeFlag, kface1, kface2);
1073}
1074
1075bool
1076SbPolyhedron::GetNextEdge(HVPoint3D &p1,
1077 HVPoint3D &p2,
1078 int &edgeFlag) const
1079/***********************************************************************
1080 * *
1081 * Name: SbPolyhedron::GetNextEdge Date: 30.09.96 *
1082 * Author: E.Chernyaev Revised: *
1083 * *
1084 * Function: Get next edge. *
1085 * Returns false when the last edge. *
1086 * *
1087 ***********************************************************************/
1088{
1089 int i1,i2;
1090 bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag);
1091 p1 = m_pV[i1];
1092 p2 = m_pV[i2];
1093 return rep;
1094}
1095
1096bool
1097SbPolyhedron::GetNextEdge(HVPoint3D &p1, HVPoint3D &p2,
1098 int &edgeFlag, int &iface1, int &iface2) const
1099/***********************************************************************
1100 * *
1101 * Name: SbPolyhedron::GetNextEdge Date: 17.11.99 *
1102 * Author: E.Chernyaev Revised: *
1103 * *
1104 * Function: Get next edge with indeces of the faces which share *
1105 * the edge. *
1106 * Returns false when the last edge. *
1107 * *
1108 ***********************************************************************/
1109{
1110 int i1,i2;
1111 bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag,iface1,iface2);
1112 p1 = m_pV[i1];
1113 p2 = m_pV[i2];
1114 return rep;
1115}
1116
1117void SbPolyhedron::GetFacet(int iFace, int &n, int *iNodes,
1118 int *edgeFlags, int *iFaces) const
1119/***********************************************************************
1120 * *
1121 * Name: SbPolyhedron::GetFacet Date: 15.12.99 *
1122 * Author: E.Chernyaev Revised: *
1123 * *
1124 * Function: Get face by index *
1125 * *
1126 ***********************************************************************/
1127{
1128 if (iFace < 1 || iFace > m_nface) {
1129 std::cerr
1130 << "SbPolyhedron::GetFacet: irrelevant index " << iFace
1131 << std::endl;
1132 n = 0;
1133 }else{
1134 int i, k;
1135 for (i=0; i<4; i++) {
1136 k = m_pF[iFace].m_edge[i].v;
1137 if (k == 0) break;
1138 if (iFaces != 0) iFaces[i] = m_pF[iFace].m_edge[i].f;
1139 if (k > 0) {
1140 iNodes[i] = k;
1141 if (edgeFlags != 0) edgeFlags[i] = 1;
1142 }else{
1143 iNodes[i] = -k;
1144 if (edgeFlags != 0) edgeFlags[i] = -1;
1145 }
1146 }
1147 n = i;
1148 }
1149}
1150
1151void SbPolyhedron::GetFacet(int index, int &n, HVPoint3D *nodes,
1152 int *edgeFlags, HVNormal3D *normals) const
1153/***********************************************************************
1154 * *
1155 * Name: SbPolyhedron::GetFacet Date: 17.11.99 *
1156 * Author: E.Chernyaev Revised: *
1157 * *
1158 * Function: Get face by index *
1159 * *
1160 ***********************************************************************/
1161{
1162 int iNodes[4];
1163 GetFacet(index, n, iNodes, edgeFlags);
1164 if (n != 0) {
1165 for (int i=0; i<4; i++) {
1166 nodes[i] = m_pV[iNodes[i]];
1167 if (normals != 0) normals[i] = FindNodeNormal(index,iNodes[i]);
1168 }
1169 }
1170}
1171
1172bool
1173SbPolyhedron::GetNextFacet(int &n, HVPoint3D *nodes,
1174 int *edgeFlags, HVNormal3D *normals) const
1175/***********************************************************************
1176 * *
1177 * Name: SbPolyhedron::GetNextFacet Date: 19.11.99 *
1178 * Author: E.Chernyaev Revised: *
1179 * *
1180 * Function: Get next face with normals of unit length at the nodes. *
1181 * Returns false when finished all faces. *
1182 * *
1183 ***********************************************************************/
1184{
1185 static int iFace = 1;
1186
1187 if (edgeFlags == 0) {
1188 GetFacet(iFace, n, nodes);
1189 }else if (normals == 0) {
1190 GetFacet(iFace, n, nodes, edgeFlags);
1191 }else{
1192 GetFacet(iFace, n, nodes, edgeFlags, normals);
1193 }
1194
1195 if (++iFace > m_nface) {
1196 iFace = 1;
1197 return false;
1198 }else{
1199 return true;
1200 }
1201}
1202
1203HVNormal3D SbPolyhedron::GetNormal(int iFace) const
1204/***********************************************************************
1205 * *
1206 * Name: SbPolyhedron::GetNormal Date: 19.11.99 *
1207 * Author: E.Chernyaev Revised: *
1208 * *
1209 * Function: Get normal of the face given by index *
1210 * *
1211 ***********************************************************************/
1212{
1213 if (iFace < 1 || iFace > m_nface) {
1214 std::cerr
1215 << "SbPolyhedron::GetNormal: irrelevant index " << iFace
1216 << std::endl;
1217 return HVNormal3D();
1218 }
1219
1220 int i0 = iabs(m_pF[iFace].m_edge[0].v);
1221 int i1 = iabs(m_pF[iFace].m_edge[1].v);
1222 int i2 = iabs(m_pF[iFace].m_edge[2].v);
1223 int i3 = iabs(m_pF[iFace].m_edge[3].v);
1224 if (i3 == 0) i3 = i0;
1225 return (m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1]);
1226}
1227
1228HVNormal3D SbPolyhedron::GetUnitNormal(int iFace) const
1229/***********************************************************************
1230 * *
1231 * Name: SbPolyhedron::GetNormal Date: 19.11.99 *
1232 * Author: E.Chernyaev Revised: *
1233 * *
1234 * Function: Get unit normal of the face given by index *
1235 * *
1236 ***********************************************************************/
1237{
1238 if (iFace < 1 || iFace > m_nface) {
1239 std::cerr
1240 << "SbPolyhedron::GetUnitNormal: irrelevant index " << iFace
1241 << std::endl;
1242 return HVNormal3D();
1243 }
1244
1245 int i0 = iabs(m_pF[iFace].m_edge[0].v);
1246 int i1 = iabs(m_pF[iFace].m_edge[1].v);
1247 int i2 = iabs(m_pF[iFace].m_edge[2].v);
1248 int i3 = iabs(m_pF[iFace].m_edge[3].v);
1249 if (i3 == 0) i3 = i0;
1250 HVNormal3D nm = (m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1]);
1251 nm.normalize();
1252 return nm;
1253}
1254
1255bool SbPolyhedron::GetNextNormal(HVNormal3D &normal) const
1256/***********************************************************************
1257 * *
1258 * Name: SbPolyhedron::GetNextNormal Date: 22.07.96 *
1259 * Author: John Allison Revised: 19.11.99 *
1260 * *
1261 * Function: Get normals of each face in face order. Returns false *
1262 * when finished all faces. *
1263 * *
1264 ***********************************************************************/
1265{
1266 static int iFace = 1;
1267 normal = GetNormal(iFace);
1268 if (++iFace > m_nface) {
1269 iFace = 1;
1270 return false;
1271 }else{
1272 return true;
1273 }
1274}
1275
1276bool SbPolyhedron::GetNextUnitNormal(HVNormal3D &normal) const
1277/***********************************************************************
1278 * *
1279 * Name: SbPolyhedron::GetNextUnitNormal Date: 16.09.96 *
1280 * Author: E.Chernyaev Revised: *
1281 * *
1282 * Function: Get normals of unit length of each face in face order. *
1283 * Returns false when finished all faces. *
1284 * *
1285 ***********************************************************************/
1286{
1287 bool rep = GetNextNormal(normal);
1288 normal.normalize();
1289 return rep;
1290}
1291
1292double SbPolyhedron::GetSurfaceArea() const
1293/***********************************************************************
1294 * *
1295 * Name: SbPolyhedron::GetSurfaceArea Date: 25.05.01 *
1296 * Author: E.Chernyaev Revised: *
1297 * *
1298 * Function: Returns area of the surface of the polyhedron. *
1299 * *
1300 ***********************************************************************/
1301{
1302 double s = 0.;
1303 for (int iFace=1; iFace<=m_nface; iFace++) {
1304 int i0 = iabs(m_pF[iFace].m_edge[0].v);
1305 int i1 = iabs(m_pF[iFace].m_edge[1].v);
1306 int i2 = iabs(m_pF[iFace].m_edge[2].v);
1307 int i3 = iabs(m_pF[iFace].m_edge[3].v);
1308 if (i3 == 0) i3 = i0;
1309 s += ((m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1])).length();
1310 }
1311 return s/2.;
1312}
1313
1314double SbPolyhedron::GetVolume() const
1315/***********************************************************************
1316 * *
1317 * Name: SbPolyhedron::GetVolume Date: 25.05.01 *
1318 * Author: E.Chernyaev Revised: *
1319 * *
1320 * Function: Returns volume of the polyhedron. *
1321 * *
1322 ***********************************************************************/
1323{
1324 double v = 0.;
1325 for (int iFace=1; iFace<=m_nface; iFace++) {
1326 int i0 = iabs(m_pF[iFace].m_edge[0].v);
1327 int i1 = iabs(m_pF[iFace].m_edge[1].v);
1328 int i2 = iabs(m_pF[iFace].m_edge[2].v);
1329 int i3 = iabs(m_pF[iFace].m_edge[3].v);
1330 HVPoint3D g;
1331 if (i3 == 0) {
1332 i3 = i0;
1333 g = (m_pV[i0]+m_pV[i1]+m_pV[i2]) * (1.0f/3.0f);
1334 }else{
1335 g = (m_pV[i0]+m_pV[i1]+m_pV[i2]+m_pV[i3]) * 0.25f;
1336 }
1337 v += ((m_pV[i2] - m_pV[i0]).cross(m_pV[i3] - m_pV[i1])).dot(g);
1338 }
1339 return v/6.;
1340}
1341
1342SbPolyhedronTrd2::SbPolyhedronTrd2(double Dx1, double Dx2,
1343 double Dy1, double Dy2,
1344 double Dz)
1345/***********************************************************************
1346 * *
1347 * Name: SbPolyhedronTrd2 Date: 22.07.96 *
1348 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1349 * *
1350 * Function: Create GEANT4 TRD2-trapezoid *
1351 * *
1352 * Input: Dx1 - half-length along X at -Dz 8----7 *
1353 * Dx2 - half-length along X ay +Dz 5----6 ! *
1354 * Dy1 - half-length along Y ay -Dz ! 4-!--3 *
1355 * Dy2 - half-length along Y ay +Dz 1----2 *
1356 * Dz - half-length along Z *
1357 * *
1358 ***********************************************************************/
1359{
1360 AllocateMemory(8,6);
1361
1362 m_pV[1] = HVPoint3D(-Dx1,-Dy1,-Dz);
1363 m_pV[2] = HVPoint3D( Dx1,-Dy1,-Dz);
1364 m_pV[3] = HVPoint3D( Dx1, Dy1,-Dz);
1365 m_pV[4] = HVPoint3D(-Dx1, Dy1,-Dz);
1366 m_pV[5] = HVPoint3D(-Dx2,-Dy2, Dz);
1367 m_pV[6] = HVPoint3D( Dx2,-Dy2, Dz);
1368 m_pV[7] = HVPoint3D( Dx2, Dy2, Dz);
1369 m_pV[8] = HVPoint3D(-Dx2, Dy2, Dz);
1370
1371 CreatePrism();
1372}
1373
1374SbPolyhedronTrd2::~SbPolyhedronTrd2() {}
1375
1376SbPolyhedronTrd1::SbPolyhedronTrd1(double Dx1, double Dx2,
1377 double Dy, double Dz)
1378 : SbPolyhedronTrd2(Dx1, Dx2, Dy, Dy, Dz) {}
1379
1380SbPolyhedronTrd1::~SbPolyhedronTrd1() {}
1381
1382SbPolyhedronBox::SbPolyhedronBox(double Dx, double Dy, double Dz)
1383 : SbPolyhedronTrd2(Dx, Dx, Dy, Dy, Dz) {}
1384
1385SbPolyhedronBox::~SbPolyhedronBox() {}
1386
1387SbPolyhedronTrap::SbPolyhedronTrap(double Dz,
1388 double Theta,
1389 double Phi,
1390 double Dy1,
1391 double Dx1,
1392 double Dx2,
1393 double Alp1,
1394 double Dy2,
1395 double Dx3,
1396 double Dx4,
1397 double Alp2)
1398/***********************************************************************
1399 * *
1400 * Name: SbPolyhedronTrap Date: 20.11.96 *
1401 * Author: E.Chernyaev Revised: *
1402 * *
1403 * Function: Create GEANT4 TRAP-trapezoid *
1404 * *
1405 * Input: DZ - half-length in Z *
1406 * Theta,Phi - polar angles of the line joining centres of the *
1407 * faces at Z=-Dz and Z=+Dz *
1408 * Dy1 - half-length in Y of the face at Z=-Dz *
1409 * Dx1 - half-length in X of low edge of the face at Z=-Dz *
1410 * Dx2 - half-length in X of top edge of the face at Z=-Dz *
1411 * Alp1 - angle between Y-axis and the median joining top and *
1412 * low edges of the face at Z=-Dz *
1413 * Dy2 - half-length in Y of the face at Z=+Dz *
1414 * Dx3 - half-length in X of low edge of the face at Z=+Dz *
1415 * Dx4 - half-length in X of top edge of the face at Z=+Dz *
1416 * Alp2 - angle between Y-axis and the median joining top and *
1417 * low edges of the face at Z=+Dz *
1418 * *
1419 ***********************************************************************/
1420{
1421 double DzTthetaCphi = Dz*tan(Theta)*cos(Phi);
1422 double DzTthetaSphi = Dz*tan(Theta)*sin(Phi);
1423 double Dy1Talp1 = Dy1*tan(Alp1);
1424 double Dy2Talp2 = Dy2*tan(Alp2);
1425
1426 AllocateMemory(8,6);
1427
1428 m_pV[1] = HVPoint3D(-DzTthetaCphi-Dy1Talp1-Dx1,-DzTthetaSphi-Dy1,-Dz);
1429 m_pV[2] = HVPoint3D(-DzTthetaCphi-Dy1Talp1+Dx1,-DzTthetaSphi-Dy1,-Dz);
1430 m_pV[3] = HVPoint3D(-DzTthetaCphi+Dy1Talp1+Dx2,-DzTthetaSphi+Dy1,-Dz);
1431 m_pV[4] = HVPoint3D(-DzTthetaCphi+Dy1Talp1-Dx2,-DzTthetaSphi+Dy1,-Dz);
1432 m_pV[5] = HVPoint3D( DzTthetaCphi-Dy2Talp2-Dx3, DzTthetaSphi-Dy2, Dz);
1433 m_pV[6] = HVPoint3D( DzTthetaCphi-Dy2Talp2+Dx3, DzTthetaSphi-Dy2, Dz);
1434 m_pV[7] = HVPoint3D( DzTthetaCphi+Dy2Talp2+Dx4, DzTthetaSphi+Dy2, Dz);
1435 m_pV[8] = HVPoint3D( DzTthetaCphi+Dy2Talp2-Dx4, DzTthetaSphi+Dy2, Dz);
1436
1437 CreatePrism();
1438}
1439
1440SbPolyhedronTrap::~SbPolyhedronTrap() {}
1441
1442SbPolyhedronPara::SbPolyhedronPara(double Dx, double Dy, double Dz,
1443 double Alpha, double Theta,
1444 double Phi)
1445 : SbPolyhedronTrap(Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
1446
1447SbPolyhedronPara::~SbPolyhedronPara() {}
1448
1449SbPolyhedronCons::SbPolyhedronCons(double Rmn1,
1450 double Rmx1,
1451 double Rmn2,
1452 double Rmx2,
1453 double Dz,
1454 double Phi1,
1455 double Dphi)
1456/***********************************************************************
1457 * *
1458 * Name: SbPolyhedronCons::SbPolyhedronCons Date: 15.12.96 *
1459 * Author: E.Chernyaev (IHEP/Protvino) Revised: 15.12.96 *
1460 * *
1461 * Function: Constructor for CONS, TUBS, CONE, TUBE *
1462 * *
1463 * Input: Rmn1, Rmx1 - inside and outside radiuses at -Dz *
1464 * Rmn2, Rmx2 - inside and outside radiuses at +Dz *
1465 * Dz - half length in Z *
1466 * Phi1 - starting angle of the segment *
1467 * Dphi - segment range *
1468 * *
1469 ***********************************************************************/
1470{
1471 static double wholeCircle=2*M_PI;
1472
1473 // C H E C K I N P U T P A R A M E T E R S
1474
1475 int k = 0;
1476 if (Rmn1 < 0. || Rmx1 < 0. || Rmn2 < 0. || Rmx2 < 0.) k = 1;
1477 if (Rmn1 > Rmx1 || Rmn2 > Rmx2) k = 1;
1478 if (Rmn1 == Rmx1 && Rmn2 == Rmx2) k = 1;
1479
1480 // We can get this from the tracking geometry. Don't complain.
1481 if (Dz == 0) return;
1482
1483 if (Dz < 0.) k += 2;
1484
1485 double phi1, phi2, dphi;
1486 if (Dphi < 0.) {
1487 phi2 = Phi1; phi1 = phi2 - Dphi;
1488 }else if (Dphi == 0.) {
1489 phi1 = Phi1; phi2 = phi1 + wholeCircle;
1490 }else{
1491 phi1 = Phi1; phi2 = phi1 + Dphi;
1492 }
1493 dphi = phi2 - phi1;
1494 if (fabs(dphi-wholeCircle) < perMillion) dphi = wholeCircle;
1495 if (dphi > wholeCircle) k += 4;
1496
1497 if (k != 0) {
1498 std::cerr << "SbPolyhedronCone(s)/Tube(s): error in input parameters";
1499 if ((k & 1) != 0) std::cerr << " (radiuses)";
1500 if ((k & 2) != 0) std::cerr << " (half-length)";
1501 if ((k & 4) != 0) std::cerr << " (angles)";
1502 std::cerr << std::endl;
1503 std::cerr << " Rmn1=" << Rmn1 << " Rmx1=" << Rmx1;
1504 std::cerr << " Rmn2=" << Rmn2 << " Rmx2=" << Rmx2;
1505 std::cerr << " Dz=" << Dz << " Phi1=" << Phi1 << " Dphi=" << Dphi
1506 << std::endl;
1507 return;
1508 }
1509
1510 // P R E P A R E T W O P O L Y L I N E S
1511
1512 double zz[4], rr[4];
1513 zz[0] = Dz;
1514 zz[1] = -Dz;
1515 zz[2] = Dz;
1516 zz[3] = -Dz;
1517 rr[0] = Rmx2;
1518 rr[1] = Rmx1;
1519 rr[2] = Rmn2;
1520 rr[3] = Rmn1;
1521
1522 // R O T A T E P O L Y L I N E S
1523
1524 RotateAroundZ(0, phi1, dphi, 2, 2, zz, rr, -1, -1);
1525 SetReferences();
1526}
1527
1528SbPolyhedronCons::~SbPolyhedronCons() {}
1529
1530SbPolyhedronCone::SbPolyhedronCone(double Rmn1, double Rmx1,
1531 double Rmn2, double Rmx2,
1532 double Dz) :
1533 SbPolyhedronCons(Rmn1, Rmx1, Rmn2, Rmx2, Dz, 0*deg, 360*deg) {}
1534
1535SbPolyhedronCone::~SbPolyhedronCone() {}
1536
1537SbPolyhedronTubs::SbPolyhedronTubs(double Rmin, double Rmax,
1538 double Dz,
1539 double Phi1, double Dphi)
1540 : SbPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, Phi1, Dphi) {}
1541
1542SbPolyhedronTubs::~SbPolyhedronTubs() {}
1543
1544SbPolyhedronTube::SbPolyhedronTube (double Rmin, double Rmax,
1545 double Dz)
1546 : SbPolyhedronCons(Rmin, Rmax, Rmin, Rmax, Dz, 0*deg, 360*deg) {}
1547
1548SbPolyhedronTube::~SbPolyhedronTube () {}
1549
1550SbPolyhedronPgon::SbPolyhedronPgon(double phi,
1551 double dphi,
1552 int npdv,
1553 int nz,
1554 const double *z,
1555 const double *rmin,
1556 const double *rmax)
1557/***********************************************************************
1558 * *
1559 * Name: SbPolyhedronPgon Date: 09.12.96 *
1560 * Author: E.Chernyaev Revised: *
1561 * *
1562 * Function: Constructor of polyhedron for PGON, PCON *
1563 * *
1564 * Input: phi - initial phi *
1565 * dphi - delta phi *
1566 * npdv - number of steps along phi *
1567 * nz - number of z-planes (at least two) *
1568 * z[] - z coordinates of the slices *
1569 * rmin[] - smaller r at the slices *
1570 * rmax[] - bigger r at the slices *
1571 * *
1572 ***********************************************************************/
1573{
1574 // C H E C K I N P U T P A R A M E T E R S
1575
1576 if (dphi <= 0. || dphi > 2*M_PI) {
1577 std::cerr
1578 << "SbPolyhedronPgon/Pcon: wrong delta phi = " << dphi
1579 << std::endl;
1580 return;
1581 }
1582
1583 if (nz < 2) {
1584 std::cerr
1585 << "SbPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
1586 << std::endl;
1587 return;
1588 }
1589
1590 if (npdv < 0) {
1591 std::cerr
1592 << "SbPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
1593 << std::endl;
1594 return;
1595 }
1596
1597 int i;
1598 for (i=0; i<nz; i++) {
1599 if (rmin[i] < 0. || rmax[i] < 0. || rmin[i] > rmax[i]) {
1600 std::cerr
1601 << "SbPolyhedronPgon: error in radiuses rmin[" << i << "]="
1602 << rmin[i] << " rmax[" << i << "]=" << rmax[i]
1603 << std::endl;
1604 return;
1605 }
1606 }
1607
1608 // P R E P A R E T W O P O L Y L I N E S
1609
1610 double *zz, *rr;
1611 zz = new double[2*nz];
1612 rr = new double[2*nz];
1613
1614 if (z[0] > z[nz-1]) {
1615 for (i=0; i<nz; i++) {
1616 zz[i] = z[i];
1617 rr[i] = rmax[i];
1618 zz[i+nz] = z[i];
1619 rr[i+nz] = rmin[i];
1620 }
1621 }else{
1622 for (i=0; i<nz; i++) {
1623 zz[i] = z[nz-i-1];
1624 rr[i] = rmax[nz-i-1];
1625 zz[i+nz] = z[nz-i-1];
1626 rr[i+nz] = rmin[nz-i-1];
1627 }
1628 }
1629
1630 // R O T A T E P O L Y L I N E S
1631
1632 RotateAroundZ(npdv, phi, dphi, nz, nz, zz, rr, -1, (npdv == 0) ? -1 : 1);
1633 SetReferences();
1634
1635 delete [] zz;
1636 delete [] rr;
1637}
1638
1639SbPolyhedronPgon::~SbPolyhedronPgon() {}
1640
1641SbPolyhedronPcon::SbPolyhedronPcon(double phi, double dphi, int nz,
1642 const double *z,
1643 const double *rmin,
1644 const double *rmax)
1645 : SbPolyhedronPgon(phi, dphi, 0, nz, z, rmin, rmax) {}
1646
1647SbPolyhedronPcon::~SbPolyhedronPcon() {}
1648
1649SbPolyhedronSphere::SbPolyhedronSphere(double rmin, double rmax,
1650 double phi, double dphi,
1651 double the, double dthe)
1652/***********************************************************************
1653 * *
1654 * Name: SbPolyhedronSphere Date: 11.12.96 *
1655 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1656 * *
1657 * Function: Constructor of polyhedron for SPHERE *
1658 * *
1659 * Input: rmin - internal radius *
1660 * rmax - external radius *
1661 * phi - initial phi *
1662 * dphi - delta phi *
1663 * the - initial theta *
1664 * dthe - delta theta *
1665 * *
1666 ***********************************************************************/
1667{
1668 // C H E C K I N P U T P A R A M E T E R S
1669
1670 if (dphi <= 0. || dphi > 2*M_PI) {
1671 std::cerr
1672 << "SbPolyhedronSphere: wrong delta phi = " << dphi
1673 << std::endl;
1674 return;
1675 }
1676
1677 if (the < 0. || the > M_PI) {
1678 std::cerr
1679 << "SbPolyhedronSphere: wrong theta = " << the
1680 << std::endl;
1681 return;
1682 }
1683
1684 if (dthe <= 0. || dthe > M_PI) {
1685 std::cerr
1686 << "SbPolyhedronSphere: wrong delta theta = " << dthe
1687 << std::endl;
1688 return;
1689 }
1690
1691 if ( (the+dthe >= M_PI) && (the+dthe < M_PI + 2*DBL_EPSILON) )
1692 dthe = M_PI - the; //G.Barrand : coming from LHCb/S.Ponce.
1693
1694 if (the+dthe > M_PI) {
1695 std::cerr
1696 << "SbPolyhedronSphere: wrong theta + delta theta = "
1697 << the << " " << dthe
1698 << std::endl;
1699 return;
1700 }
1701
1702 if (rmin < 0. || rmin >= rmax) {
1703 std::cerr
1704 << "SbPolyhedronSphere: error in radiuses"
1705 << " rmin=" << rmin << " rmax=" << rmax
1706 << std::endl;
1707 return;
1708 }
1709
1710 // P R E P A R E T W O P O L Y L I N E S
1711
1712 int ns = (GetNumberOfRotationSteps() + 1) / 2;
1713 int np1 = int(dthe*ns/M_PI+.5) + 1;
1714 if (np1 <= 1) np1 = 2;
1715 int np2 = rmin < perMillion ? 1 : np1;
1716
1717 double *zz, *rr;
1718 zz = new double[np1+np2];
1719 rr = new double[np1+np2];
1720
1721 double a = dthe/(np1-1);
1722 double cosa, sina;
1723 for (int i=0; i<np1; i++) {
1724 cosa = cos(the+i*a);
1725 sina = sin(the+i*a);
1726 zz[i] = rmax*cosa;
1727 rr[i] = rmax*sina;
1728 if (np2 > 1) {
1729 zz[i+np1] = rmin*cosa;
1730 rr[i+np1] = rmin*sina;
1731 }
1732 }
1733 if (np2 == 1) {
1734 zz[np1] = 0.;
1735 rr[np1] = 0.;
1736 }
1737
1738 // R O T A T E P O L Y L I N E S
1739
1740 RotateAroundZ(0, phi, dphi, np1, np2, zz, rr, -1, -1);
1741 SetReferences();
1742
1743 delete [] zz;
1744 delete [] rr;
1745}
1746
1747SbPolyhedronSphere::~SbPolyhedronSphere() {}
1748
1749SbPolyhedronTorus::SbPolyhedronTorus(double rmin,
1750 double rmax,
1751 double rtor,
1752 double phi,
1753 double dphi)
1754/***********************************************************************
1755 * *
1756 * Name: SbPolyhedronTorus Date: 11.12.96 *
1757 * Author: E.Chernyaev (IHEP/Protvino) Revised: *
1758 * *
1759 * Function: Constructor of polyhedron for TORUS *
1760 * *
1761 * Input: rmin - internal radius *
1762 * rmax - external radius *
1763 * rtor - radius of torus *
1764 * phi - initial phi *
1765 * dphi - delta phi *
1766 * *
1767 ***********************************************************************/
1768{
1769 // C H E C K I N P U T P A R A M E T E R S
1770
1771 if (dphi <= 0. || dphi > 2*M_PI) {
1772 std::cerr
1773 << "SbPolyhedronTorus: wrong delta phi = " << dphi
1774 << std::endl;
1775 return;
1776 }
1777
1778 if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
1779 std::cerr
1780 << "SbPolyhedronTorus: error in radiuses"
1781 << " rmin=" << rmin << " rmax=" << rmax << " rtorus=" << rtor
1782 << std::endl;
1783 return;
1784 }
1785
1786 // P R E P A R E T W O P O L Y L I N E S
1787
1788 int np1 = GetNumberOfRotationSteps();
1789 int np2 = rmin < perMillion ? 1 : np1;
1790
1791 //double *zz(nullptr), *rr;
1792 const auto totNumPts = np1+np2;
1793 auto zz = new double[totNumPts]{};
1794 auto rr = new double[totNumPts]{};
1795 for (unsigned int i(0);i!=(unsigned int) totNumPts;++i){
1796 rr[i]=0.0;
1797 zz[i]=0.0;
1798 }
1799
1800 double a = 2*M_PI/np1;
1801 double cosa, sina;
1802 for (int i=0; i<np1; i++) {
1803 cosa = cos(i*a);
1804 sina = sin(i*a);
1805 zz[i] = rmax*cosa;
1806 rr[i] = rtor+rmax*sina;
1807 if (np2 > 1) {
1808 zz[i+np1] = rmin*cosa;
1809 rr[i+np1] = rtor+rmin*sina;
1810 }
1811 }
1812 if (np2 == 1) {
1813 zz[np1] = 0.;
1814 rr[np1] = rtor;
1815 np2 = -1;
1816 }
1817
1818 // R O T A T E P O L Y L I N E S
1819
1820 RotateAroundZ(0, phi, dphi, -np1, -np2, zz, rr, -1,-1);
1821 SetReferences();
1822
1823 delete [] zz;
1824 delete [] rr;
1825}
1826
1827SbPolyhedronTorus::~SbPolyhedronTorus() {}
1828
1829int SbPolyhedron::s_numberOfRotationSteps = DEFAULT_NUMBER_OF_STEPS;
1830/***********************************************************************
1831 * *
1832 * Name: SbPolyhedron::s_numberOfRotationSteps Date: 24.06.97 *
1833 * Author: J.Allison (Manchester University) Revised: *
1834 * *
1835 * Function: Number of steps for whole circle *
1836 * *
1837 ***********************************************************************/
1838
1839#include "BooleanProcessor.h"
1840static HEPVis_BooleanProcessor processor;
1841
1842SbPolyhedron SbPolyhedron::add(const SbPolyhedron & p) const
1843/***********************************************************************
1844 * *
1845 * Name: SbPolyhedron::add Date: 19.03.00 *
1846 * Author: E.Chernyaev Revised: *
1847 * *
1848 * Function: Boolean "union" of two polyhedra *
1849 * *
1850 ***********************************************************************/
1851{
1852 // rbianchi change
1853 //return processor.execute(OP_UNION, *this, p);
1854 HEPVis_BooleanProcessor processor; //G.Barrand
1855 int err;
1856 return processor.execute(OP_UNION, *this, p, err);
1857 //---
1858}
1859
1860SbPolyhedron SbPolyhedron::intersect(const SbPolyhedron & p) const
1861/***********************************************************************
1862 * *
1863 * Name: SbPolyhedron::intersect Date: 19.03.00 *
1864 * Author: E.Chernyaev Revised: *
1865 * *
1866 * Function: Boolean "intersection" of two polyhedra *
1867 * *
1868 ***********************************************************************/
1869{
1870 // rbianchi change
1871 //return processor.execute(OP_INTERSECTION, *this, p);
1872 HEPVis_BooleanProcessor processor; //G.Barrand
1873 int err;
1874 return processor.execute(OP_INTERSECTION, *this, p, err);
1875 //---
1876}
1877
1878SbPolyhedron SbPolyhedron::subtract(const SbPolyhedron & p) const
1879/***********************************************************************
1880 * *
1881 * Name: SbPolyhedron::add Date: 19.03.00 *
1882 * Author: E.Chernyaev Revised: *
1883 * *
1884 * Function: Boolean "subtraction" of "p" from "this" *
1885 * *
1886 ***********************************************************************/
1887{
1888 // rbianchi change
1889 //return processor.execute(OP_SUBTRACTION, *this, p);
1890 HEPVis_BooleanProcessor processor; //G.Barrand
1891 int err;
1892 return processor.execute(OP_SUBTRACTION, *this, p, err);
1893 //---
1894}
1895
1896
1901
1902#include "PolygonTriangulator.h"
1903#include <map>
1904#include <set>
1905
1906class SbPolyhedronPolygonXSect::Internals {
1907public:
1908
1909 typedef std::pair<int,int> Edge;
1910 typedef PolygonTriangulator::Triangle Triangle;
1911 typedef PolygonTriangulator::Triangles Triangles;
1912
1913 // Methods:
1917 ~Internals() { delete poly; }
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 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//_______________________________________________________________________________________________________________________
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:43
perMillion
#define perMillion
Definition SbPolyhedron.cxx:16
iabs
int iabs(int a)
Definition SbPolyhedron.cxx:22
ffabs
float ffabs(float a)
Definition SbPolyhedron.cxx:25
deg
#define deg
Definition SbPolyhedron.cxx:17
processor
static HEPVis_BooleanProcessor processor
Definition SbPolyhedron.cxx:1840
operator<<
std::ostream & operator<<(std::ostream &ostr, const SbFacet &facet)
Definition SbPolyhedron.cxx:118
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:35
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:1382
SbPolyhedronBox::~SbPolyhedronBox
virtual ~SbPolyhedronBox()
Definition SbPolyhedron.cxx:1385
SbPolyhedronCone::~SbPolyhedronCone
virtual ~SbPolyhedronCone()
Definition SbPolyhedron.cxx:1535
SbPolyhedronCone::SbPolyhedronCone
SbPolyhedronCone(double Rmn1, double Rmx1, double Rmn2, double Rmx2, double Dz)
Definition SbPolyhedron.cxx:1530
SbPolyhedronCons::SbPolyhedronCons
SbPolyhedronCons(double Rmn1, double Rmx1, double Rmn2, double Rmx2, double Dz, double Phi1, double Dphi)
Definition SbPolyhedron.cxx:1449
SbPolyhedronCons::~SbPolyhedronCons
virtual ~SbPolyhedronCons()
Definition SbPolyhedron.cxx:1528
SbPolyhedronGenericTrap::SbPolyhedronGenericTrap
SbPolyhedronGenericTrap(double Dz, const std::vector< std::pair< double, double > > &Vertices)
Definition SbPolyhedron.cxx:2411
SbPolyhedronPara::~SbPolyhedronPara
virtual ~SbPolyhedronPara()
Definition SbPolyhedron.cxx:1447
SbPolyhedronPara::SbPolyhedronPara
SbPolyhedronPara(double Dx, double Dy, double Dz, double Alpha, double Theta, double Phi)
Definition SbPolyhedron.cxx:1442
SbPolyhedronPcon::SbPolyhedronPcon
SbPolyhedronPcon(double phi, double dphi, int nz, const double *z, const double *rmin, const double *rmax)
Definition SbPolyhedron.cxx:1641
SbPolyhedronPcon::~SbPolyhedronPcon
virtual ~SbPolyhedronPcon()
Definition SbPolyhedron.cxx:1647
SbPolyhedronPgon::~SbPolyhedronPgon
virtual ~SbPolyhedronPgon()
Definition SbPolyhedron.cxx:1639
SbPolyhedronPgon::SbPolyhedronPgon
SbPolyhedronPgon(double phi, double dphi, int npdv, int nz, const double *z, const double *rmin, const double *rmax)
Definition SbPolyhedron.cxx:1550
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:1910
SbPolyhedronPolygonXSect::Internals::Triangles
PolygonTriangulator::Triangles Triangles
Definition SbPolyhedron.cxx:1911
SbPolyhedronPolygonXSect::Internals::y
const std::vector< double > * y
Definition SbPolyhedron.cxx:1928
SbPolyhedronPolygonXSect::Internals::Internals
Internals(SbPolyhedronPolygonXSect *sbp)
Definition SbPolyhedron.cxx:1914
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::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:1747
SbPolyhedronSphere::SbPolyhedronSphere
SbPolyhedronSphere(double rmin, double rmax, double phi, double dphi, double the, double dthe)
Definition SbPolyhedron.cxx:1649
SbPolyhedronTorus::SbPolyhedronTorus
SbPolyhedronTorus(double rmin, double rmax, double rtor, double phi, double dphi)
Definition SbPolyhedron.cxx:1749
SbPolyhedronTorus::~SbPolyhedronTorus
virtual ~SbPolyhedronTorus()
Definition SbPolyhedron.cxx:1827
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:1387
SbPolyhedronTrap::~SbPolyhedronTrap
virtual ~SbPolyhedronTrap()
Definition SbPolyhedron.cxx:1440
SbPolyhedronTrd1::~SbPolyhedronTrd1
virtual ~SbPolyhedronTrd1()
Definition SbPolyhedron.cxx:1380
SbPolyhedronTrd1::SbPolyhedronTrd1
SbPolyhedronTrd1(double Dx1, double Dx2, double Dy, double Dz)
Definition SbPolyhedron.cxx:1376
SbPolyhedronTrd2::SbPolyhedronTrd2
SbPolyhedronTrd2(double Dx1, double Dx2, double Dy1, double Dy2, double Dz)
Definition SbPolyhedron.cxx:1342
SbPolyhedronTrd2::~SbPolyhedronTrd2
virtual ~SbPolyhedronTrd2()
Definition SbPolyhedron.cxx:1374
SbPolyhedronTube::SbPolyhedronTube
SbPolyhedronTube(double Rmin, double Rmax, double Dz)
Definition SbPolyhedron.cxx:1544
SbPolyhedronTube::~SbPolyhedronTube
virtual ~SbPolyhedronTube()
Definition SbPolyhedron.cxx:1548
SbPolyhedronTubs::~SbPolyhedronTubs
virtual ~SbPolyhedronTubs()
Definition SbPolyhedron.cxx:1542
SbPolyhedronTubs::SbPolyhedronTubs
SbPolyhedronTubs(double Rmin, double Rmax, double Dz, double Phi1, double Dphi)
Definition SbPolyhedron.cxx:1537
SbPolyhedron::GetNextVertex
bool GetNextVertex(HVPoint3D &vertex, int &edgeFlag) const
Definition SbPolyhedron.cxx:957
SbPolyhedron::GetNextEdge
bool GetNextEdge(HVPoint3D &p1, HVPoint3D &p2, int &edgeFlag) const
Definition SbPolyhedron.cxx:1076
SbPolyhedron::intersect
SbPolyhedron intersect(const SbPolyhedron &p) const
Definition SbPolyhedron.cxx:1860
SbPolyhedron::add
SbPolyhedron add(const SbPolyhedron &p) const
Definition SbPolyhedron.cxx:1842
SbPolyhedron::FindNeighbour
int FindNeighbour(int iFace, int iNode, int iOrder) const
Definition SbPolyhedron.cxx:189
SbPolyhedron::operator=
virtual SbPolyhedron & operator=(const SbPolyhedron &from)
Definition SbPolyhedron.cxx:161
SbPolyhedron::GetNextNormal
bool GetNextNormal(HVNormal3D &normal) const
Definition SbPolyhedron.cxx:1255
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:438
SbPolyhedron::GetNextEdgeIndeces
bool GetNextEdgeIndeces(int &i1, int &i2, int &edgeFlag, int &iface1, int &iface2) const
Definition SbPolyhedron.cxx:1007
SbPolyhedron::GetSurfaceArea
double GetSurfaceArea() const
Definition SbPolyhedron.cxx:1292
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:267
SbPolyhedron::GetVertexFast
const HVPoint3D & GetVertexFast(int index) const
Definition SbPolyhedron.cxx:950
SbPolyhedron::GetUnitNormal
HVNormal3D GetUnitNormal(int iFace) const
Definition SbPolyhedron.cxx:1228
SbPolyhedron::GetNextUnitNormal
bool GetNextUnitNormal(HVNormal3D &normal) const
Definition SbPolyhedron.cxx:1276
SbPolyhedron::subtract
SbPolyhedron subtract(const SbPolyhedron &p) const
Definition SbPolyhedron.cxx:1878
SbPolyhedron::GetVolume
double GetVolume() const
Definition SbPolyhedron.cxx:1314
SbPolyhedron::SetSideFacets
void SetSideFacets(int ii[4], int vv[4], int *kk, double *r, double dphi, int ns, int &kface)
Definition SbPolyhedron.cxx:368
SbPolyhedron::GetNextVertexIndex
bool GetNextVertexIndex(int &index, int &edgeFlag) const
Definition SbPolyhedron.cxx:902
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:306
SbPolyhedron::GetVertex
HVPoint3D GetVertex(int index) const
Definition SbPolyhedron.cxx:929
SbPolyhedron::Transform
SbPolyhedron & Transform(const HEPVis::SbRotation &rot, const SbVec3d &trans)
Definition SbPolyhedron.cxx:878
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:217
SbPolyhedron::GetFacet
void GetFacet(int iFace, int &n, int *iNodes, int *edgeFlags=0, int *iFaces=0) const
Definition SbPolyhedron.cxx:1117
SbPolyhedron::GetNextFacet
bool GetNextFacet(int &n, HVPoint3D *nodes, int *edgeFlags=0, HVNormal3D *normals=0) const
Definition SbPolyhedron.cxx:1173
SbPolyhedron::GetNormal
HVNormal3D GetNormal(int iFace) const
Definition SbPolyhedron.cxx:1203
SbPolyhedron::SetNumberOfRotationSteps
static void SetNumberOfRotationSteps(int n)
Definition SbPolyhedron.cxx:245
r
int r
Definition globals.cxx:22
const
index
Definition index.py:1
Generated on for ATLAS Offline Software by doxygen 1.14.0