16 #define perMillion 0.000001
17 #define deg (M_PI/180.0)
23 return a < 0 ? -
a :
a;
26 return a < 0.0f ? -
a :
a;
36 setValue(
v[0],
v[1],
v[2]);
40 setValue(
v[0],
v[1],
v[2]);
119 for (
int k=0;
k<4;
k++) {
127 ostr <<
"Nverteces=" << ph.
m_nvert <<
", Nfacets=" << ph.
m_nface << std::endl;
130 ostr <<
"xyz(" <<
i <<
")="
135 ostr <<
"face(" <<
i <<
")=" << ph.
m_pF[
i] << std::endl;
171 if (
this == &from)
return *
this;
200 for (
i=0;
i<4;
i++) {
201 if (iNode ==
iabs(m_pF[iFace].m_edge[
i].
v))
break;
205 <<
"SbPolyhedron::FindNeighbour: face " << iFace
206 <<
" has no node " << iNode
212 if (m_pF[iFace].m_edge[
i].
v == 0)
i = 2;
214 return (m_pF[iFace].m_edge[
i].
v > 0) ? 0 : m_pF[iFace].m_edge[
i].f;
228 int k = iFace, iOrder = 1;
231 k = FindNeighbour(
k, iNode, iOrder);
232 if (
k == iFace)
break;
234 normal += GetUnitNormal(
k);
236 if (iOrder < 0)
break;
258 <<
"SbPolyhedron::SetNumberOfRotationSteps: attempt to set the\n"
259 <<
"number of steps per circle < " << nMin <<
"; forced to " << nMin
299 m_pF[2] =
SbFacet(5,TOP, 8,BACK, 4,BOTTOM, 1,FRONT);
301 m_pF[4] =
SbFacet(7,TOP, 6,FRONT, 2,BOTTOM, 3,BACK);
307 int v1,
int v2,
int vEdge,
308 bool ifWholeCircle,
int ns,
int &kface)
328 if (r1 == 0. && r2 == 0)
return;
333 int ii1 = ifWholeCircle ? i1 : i1+
ns;
334 int ii2 = ifWholeCircle ? i2 : i2+
ns;
335 int vv = ifWholeCircle ? vEdge : 1;
340 }
else if (r2 == 0.) {
343 m_pF[kface++] =
SbFacet(i1,0,
v2*i2,0, (i2+1),0, v1*(i1+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);
352 }
else if (r2 == 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);
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);
370 double dphi,
int ns,
int &kface)
390 for (
int i=0;
i<4;
i++) {
392 k2 = (
i == 3) ? ii[0] : ii[
i+1];
393 if (
r[k1] == 0. &&
r[k2] == 0.)
vv[
i] = -1;
397 if (ii[1] == ii[2]) {
402 if (
r[ii[0]] != 0.) k1 +=
ns;
403 if (
r[ii[2]] != 0.) k2 +=
ns;
404 if (
r[ii[3]] != 0.) k3 +=
ns;
406 }
else if (
kk[ii[0]] ==
kk[ii[1]]) {
411 if (
r[ii[0]] != 0.) k1 +=
ns;
412 if (
r[ii[2]] != 0.) k2 +=
ns;
413 if (
r[ii[3]] != 0.) k3 +=
ns;
415 }
else if (
kk[ii[2]] ==
kk[ii[3]]) {
420 if (
r[ii[0]] != 0.) k1 +=
ns;
421 if (
r[ii[1]] != 0.) k2 +=
ns;
422 if (
r[ii[2]] != 0.) k3 +=
ns;
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;
440 const double *
z,
double *
r,
441 int nodeVis,
int edgeVis)
464 static double wholeCircle = 2*
M_PI;
468 bool ifWholeCircle = (fabs(dphi-wholeCircle) <
perMillion) ?
470 double delPhi = ifWholeCircle ? wholeCircle : dphi;
471 int nSphi = (nstep > 0) ?
473 if (nSphi == 0) nSphi = 1;
474 int nVphi = ifWholeCircle ? nSphi : nSphi+1;
475 bool ifClosed = np1 > 0 ? false :
true;
479 int absNp1 =
iabs(np1);
480 int absNp2 =
iabs(np2);
482 int i1end = absNp1-1;
484 int i2end = absNp1+absNp2-1;
487 for(
i=i1beg;
i<=i2end;
i++) {
492 for (
i=i1beg;
i<=i1end;
i++) {
493 j += (
r[
i] == 0.) ? 1 : nVphi;
496 bool ifSide1 =
false;
497 bool ifSide2 =
false;
499 if (
r[i2beg] !=
r[i1beg] ||
z[i2beg] !=
z[i1beg]) {
500 j += (
r[i2beg] == 0.) ? 1 : nVphi;
504 for(
i=i2beg+1;
i<i2end;
i++) {
505 j += (
r[
i] == 0.) ? 1 : nVphi;
508 if (
r[i2end] !=
r[i1end] ||
z[i2end] !=
z[i1end]) {
509 if (absNp2 > 1) j += (
r[i2end] == 0.) ? 1 : nVphi;
515 k = ifClosed ? absNp1*nSphi : (absNp1-1)*nSphi;
518 for(
i=i2beg;
i<i2end;
i++) {
519 if (
r[
i] > 0. ||
r[
i+1] > 0.)
k += nSphi;
523 if (
r[i2end] > 0. ||
r[i2beg] > 0.)
k += nSphi;
528 if (ifSide1 && (
r[i1beg] > 0. ||
r[i2beg] > 0.))
k += nSphi;
529 if (ifSide2 && (
r[i1end] > 0. ||
r[i2end] > 0.))
k += nSphi;
532 if (!ifWholeCircle) {
533 k += ifClosed ? 2*absNp1 : 2*(absNp1-1);
543 kk =
new int[absNp1+absNp2];
546 for(
i=i1beg;
i<=i1end;
i++) {
559 for(
i=i2beg+1;
i<i2end;
i++) {
574 double cosPhi, sinPhi;
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++) {
589 v2 = ifClosed ? nodeVis : 1;
590 for(
i=i1beg;
i<i1end;
i++) {
592 if (!ifClosed &&
i == i1end-1) {
595 v2 = (
r[
i] ==
r[
i+1] &&
r[
i+1] ==
r[
i+2]) ? -1 : nodeVis;
598 edgeVis, ifWholeCircle, nSphi,
k);
602 edgeVis, ifWholeCircle, nSphi,
k);
608 v2 = ifClosed ? nodeVis : 1;
609 for(
i=i2beg;
i<i2end;
i++) {
611 if (!ifClosed &&
i==i2end-1) {
614 v2 = (
r[
i] ==
r[
i+1] &&
r[
i+1] ==
r[
i+2]) ? -1 : nodeVis;
617 edgeVis, ifWholeCircle, nSphi,
k);
621 edgeVis, ifWholeCircle, nSphi,
k);
630 -1, ifWholeCircle, nSphi,
k);
634 -1, ifWholeCircle, nSphi,
k);
640 if (!ifWholeCircle) {
645 for (
i=i1beg;
i<=i1end;
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;
657 for (
i=i1beg;
i<i1end;
i++) {
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;
664 vv[2] = (
i == i1end-1) ? 1 : -1;
675 <<
"Polyhedron::RotateAroundZ: number of generated faces ("
676 <<
k-1 <<
") is not equal to the number of allocated faces ("
694 struct edgeListMember {
695 edgeListMember *
next;
699 } *edgeList, *freeList, **headList;
704 edgeList =
new edgeListMember[2*
m_nface];
705 headList =
new edgeListMember*[
m_nvert];
713 edgeList[
i].next = &edgeList[
i+1];
715 edgeList[2*
m_nface-1].next = 0;
719 int iface, iedge, nedge, i1, i2, k1, k2;
720 edgeListMember *prev, *
cur;
722 for(iface=1; iface<=
m_nface; iface++) {
724 for (iedge=0; iedge<nedge; iedge++) {
726 i2 = (iedge < nedge-1) ? iedge+1 : 0;
729 k1 = (i1 < i2) ? i1 : i2;
730 k2 = (i1 > i2) ? i1 : i2;
735 headList[k1] = freeList;
736 freeList = freeList->next;
746 headList[k1] =
cur->next;
747 cur->next = freeList;
755 <<
"Polyhedron::SetReferences: different edge visibility "
756 << iface <<
"/" << iedge <<
"/"
758 <<
cur->iface <<
"/" <<
cur->iedge <<
"/"
770 prev->next = freeList;
771 freeList = freeList->next;
781 prev->next =
cur->next;
782 cur->next = freeList;
790 <<
"Polyhedron::SetReferences: different edge visibility "
791 << iface <<
"/" << iedge <<
"/"
793 <<
cur->iface <<
"/" <<
cur->iedge <<
"/"
806 if (headList[
i] != 0) {
808 <<
"Polyhedron::SetReferences: List " <<
i <<
" is not empty"
830 int i,
k, nnode,
v[4],
f[4];
833 for (
k=0;
k<nnode;
k++) {
838 for (
k=0;
k<nnode;
k++) {
880 ,
const SbVec3d& translation
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);
912 static int iFace = 1;
913 static int iQVertex = 0;
914 int vIndex = m_pF[iFace].m_edge[iQVertex].v;
916 edgeFlag = (vIndex > 0) ? 1 : 0;
919 if (iQVertex >= 3 || m_pF[iFace].m_edge[iQVertex+1].
v == 0) {
921 if (++iFace > m_nface) iFace = 1;
939 if (index <= 0 || index > m_nvert) {
941 <<
"SbPolyhedron::GetVertex: irrelevant index " <<
index
970 bool rep = GetNextVertexIndex(
index, edgeFlag);
988 static int iFace = 1;
989 static int iNode = 0;
991 if (m_nface == 0)
return false;
993 int k = m_pF[iFace].m_edge[iNode].v;
994 if (
k > 0) { edgeFlag = 1; }
else { edgeFlag = -1;
k = -
k; }
996 normal = FindNodeNormal(iFace,
k);
997 if (iNode >= 3 || m_pF[iFace].m_edge[iNode+1].
v == 0) {
999 if (++iFace > m_nface) iFace = 1;
1008 int &iface1,
int &iface2)
const
1020 static int iFace = 1;
1021 static int iQVertex = 0;
1022 static int iOrder = 1;
1023 int k1, k2, kflag, kface1, kface2;
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;
1033 k1 = m_pF[iFace].m_edge[iQVertex].v;
1037 kface2 = m_pF[iFace].m_edge[iQVertex].f;
1038 if (iQVertex >= 3 || m_pF[iFace].m_edge[iQVertex+1].
v == 0) {
1040 k2 =
iabs(m_pF[iFace].m_edge[iQVertex].
v);
1044 k2 =
iabs(m_pF[iFace].m_edge[iQVertex].
v);
1046 }
while (iOrder*k1 > iOrder*k2);
1048 i1 = k1; i2 = k2; edgeFlag = (kflag > 0) ? 1 : 0;
1049 iface1 = kface1; iface2 = kface2;
1051 if (iFace > m_nface) {
1052 iFace = 1; iOrder = 1;
1072 return GetNextEdgeIndeces(i1, i2, edgeFlag, kface1, kface2);
1078 int &edgeFlag)
const
1090 bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag);
1098 int &edgeFlag,
int &iface1,
int &iface2)
const
1111 bool rep = GetNextEdgeIndeces(i1,i2,edgeFlag,iface1,iface2);
1118 int *edgeFlags,
int *iFaces)
const
1128 if (iFace < 1 || iFace > m_nface) {
1130 <<
"SbPolyhedron::GetFacet: irrelevant index " << iFace
1135 for (
i=0;
i<4;
i++) {
1136 k = m_pF[iFace].m_edge[
i].v;
1138 if (iFaces != 0) iFaces[
i] = m_pF[iFace].m_edge[
i].f;
1141 if (edgeFlags != 0) edgeFlags[
i] = 1;
1144 if (edgeFlags != 0) edgeFlags[
i] = -1;
1163 GetFacet(
index,
n, iNodes, edgeFlags);
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]);
1185 static int iFace = 1;
1187 if (edgeFlags == 0) {
1188 GetFacet(iFace,
n, nodes);
1189 }
else if (normals == 0) {
1190 GetFacet(iFace,
n, nodes, edgeFlags);
1192 GetFacet(iFace,
n, nodes, edgeFlags, normals);
1195 if (++iFace > m_nface) {
1213 if (iFace < 1 || iFace > m_nface) {
1215 <<
"SbPolyhedron::GetNormal: irrelevant index " << iFace
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]);
1238 if (iFace < 1 || iFace > m_nface) {
1240 <<
"SbPolyhedron::GetUnitNormal: irrelevant index " << iFace
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]);
1266 static int iFace = 1;
1267 normal = GetNormal(iFace);
1268 if (++iFace > m_nface) {
1287 bool rep = GetNextNormal(normal);
1303 for (
int iFace=1; iFace<=
m_nface; iFace++) {
1308 if (i3 == 0) i3 = i0;
1325 for (
int iFace=1; iFace<=
m_nface; iFace++) {
1343 double Dy1,
double Dy2,
1377 double Dy,
double Dz)
1423 double Dy1Talp1 = Dy1*
tan(Alp1);
1424 double Dy2Talp2 = Dy2*
tan(Alp2);
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);
1443 double Alpha,
double Theta,
1445 :
SbPolyhedronTrap(
Dz, Theta,
Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha) {}
1471 static double wholeCircle=2*
M_PI;
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;
1481 if (
Dz == 0)
return;
1483 if (
Dz < 0.)
k += 2;
1485 double phi1, phi2, dphi;
1487 phi2 = Phi1; phi1 = phi2 - Dphi;
1488 }
else if (Dphi == 0.) {
1489 phi1 = Phi1; phi2 = phi1 + wholeCircle;
1491 phi1 = Phi1; phi2 = phi1 + Dphi;
1494 if (fabs(dphi-wholeCircle) <
perMillion) dphi = wholeCircle;
1495 if (dphi > wholeCircle)
k += 4;
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
1512 double zz[4],
rr[4];
1531 double Rmn2,
double Rmx2,
1539 double Phi1,
double Dphi)
1576 if (dphi <= 0. || dphi > 2*
M_PI) {
1578 <<
"SbPolyhedronPgon/Pcon: wrong delta phi = " << dphi
1585 <<
"SbPolyhedronPgon/Pcon: number of z-planes less than two = " << nz
1592 <<
"SbPolyhedronPgon/Pcon: error in number of phi-steps =" << npdv
1598 for (
i=0;
i<nz;
i++) {
1599 if (rmin[
i] < 0. || rmax[
i] < 0. || rmin[
i] > rmax[
i]) {
1601 <<
"SbPolyhedronPgon: error in radiuses rmin[" <<
i <<
"]="
1602 << rmin[
i] <<
" rmax[" <<
i <<
"]=" << rmax[
i]
1611 zz =
new double[2*nz];
1612 rr =
new double[2*nz];
1614 if (
z[0] >
z[nz-1]) {
1615 for (
i=0;
i<nz;
i++) {
1622 for (
i=0;
i<nz;
i++) {
1624 rr[
i] = rmax[nz-
i-1];
1625 zz[
i+nz] =
z[nz-
i-1];
1626 rr[
i+nz] = rmin[nz-
i-1];
1650 double phi,
double dphi,
1651 double the,
double dthe)
1670 if (dphi <= 0. || dphi > 2*
M_PI) {
1672 <<
"SbPolyhedronSphere: wrong delta phi = " << dphi
1677 if (the < 0. || the >
M_PI) {
1679 <<
"SbPolyhedronSphere: wrong theta = " << the
1684 if (dthe <= 0. || dthe >
M_PI) {
1686 <<
"SbPolyhedronSphere: wrong delta theta = " << dthe
1691 if ( (the+dthe >=
M_PI) && (the+dthe <
M_PI + 2*DBL_EPSILON) )
1694 if (the+dthe >
M_PI) {
1696 <<
"SbPolyhedronSphere: wrong theta + delta theta = "
1697 << the <<
" " << dthe
1702 if (rmin < 0. || rmin >= rmax) {
1704 <<
"SbPolyhedronSphere: error in radiuses"
1705 <<
" rmin=" << rmin <<
" rmax=" << rmax
1714 if (np1 <= 1) np1 = 2;
1718 zz =
new double[np1+np2];
1719 rr =
new double[np1+np2];
1721 double a = dthe/(np1-1);
1723 for (
int i=0;
i<np1;
i++) {
1724 cosa =
cos(the+
i*
a);
1725 sina =
sin(the+
i*
a);
1729 zz[
i+np1] = rmin*cosa;
1730 rr[
i+np1] = rmin*sina;
1771 if (dphi <= 0. || dphi > 2*
M_PI) {
1773 <<
"SbPolyhedronTorus: wrong delta phi = " << dphi
1778 if (rmin < 0. || rmin >= rmax || rmax >= rtor) {
1780 <<
"SbPolyhedronTorus: error in radiuses"
1781 <<
" rmin=" << rmin <<
" rmax=" << rmax <<
" rtorus=" << rtor
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){
1800 double a = 2*
M_PI/np1;
1802 for (
int i=0;
i<np1;
i++) {
1806 rr[
i] = rtor+rmax*sina;
1808 zz[
i+np1] = rmin*cosa;
1809 rr[
i+np1] = rtor+rmin*sina;
1840 static HEPVis_BooleanProcessor processor;
1854 HEPVis_BooleanProcessor processor;
1872 HEPVis_BooleanProcessor processor;
1890 HEPVis_BooleanProcessor processor;
1919 void setData(
const std::vector<double> * xx,
const std::vector<double>*
yy,
const double&
dz);
1927 const std::vector<double> *
x;
1928 const std::vector<double> *
y;
1947 bool isinternaledge(
const Edge& edge,
const std::map<
Edge,std::vector<Triangles::const_iterator> >& edge2triangles_map)
const;
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;
1954 void getSurroundingValues(
const std::vector<unsigned>& numbercycle,
const unsigned& centralValue,
1955 unsigned& prevValue,
unsigned& nextValue)
const;
1986 assert (
n==
yy->size()&&
n>2);
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);
2013 if (edge2triangles_map.find(eun)==edge2triangles_map.end()) {
2014 alledges.insert(eun);
2016 edge2triangles_map[eun].push_back(itt);
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);
2027 edges_external.insert(*itedges);
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);
2035 Edge eun = GetEdge(&(*itt),iedge,
false);
2036 assert(edge2triangles_map.find(eun)!=edge2triangles_map.end());
2037 if (edges_external.find(eun)!=edges_external.end()) {
2039 neighbourmap[eun]=0;
2041 assert(edge2triangles_map[eun].
size()==2);
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);
2063 std::set<const Triangle*> triangles_with_extra_vertex;
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);
2072 if (edges_internal.find(eun)!=edges_internal.end())
2075 if (ninternaledges<3)
2081 for (;iedge<3;++iedge) {
2082 Edge e = GetEdge(&(*itt),iedge,
true);
2083 const Triangle* trneighbour = neighbourmap[
e];
2084 assert(trneighbour);
2085 if (triangles_with_extra_vertex.insert(trneighbour).second) {
2086 triangles_with_extra_vertex.insert(&(*itt));
2087 Edge eun = GetEdge(&(*itt),iedge,
false);
2088 edges_with_extra_vertex.insert(eun);
2089 ++nextrainternalvertices;
2097 itt=
poly->triangles()->begin();
2098 for(; itt!=ittE; ++itt) {
2100 if (triangles_with_extra_vertex.find(triangle)!=triangles_with_extra_vertex.end())
2102 for (
unsigned iedge=0;iedge<3;++iedge) {
2103 Edge eun = GetEdge(&(*itt),iedge,
false);
2105 if (edges_external.find(eun)==edges_external.end())
2107 assert(edges_with_extra_vertex.find(eun)==edges_with_extra_vertex.end());
2108 triangles_with_extra_vertex.insert(triangle);
2109 edges_with_extra_vertex.insert(eun);
2110 ++nextraexternalvertices;
2114 assert(ntriangles==triangles_with_extra_vertex.size());
2120 const unsigned nfaces = 2*ntriangles+
n+nextraexternalvertices;
2121 const unsigned nvertices = 2*(
n+nextraexternalvertices+nextrainternalvertices);
2122 sbpolyhedron->AllocateMemory(nvertices,nfaces);
2127 std::set<Edge>::const_iterator itl = edges_with_extra_vertex.begin();
2128 std::set<Edge>::const_iterator itlE = edges_with_extra_vertex.end();
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;
2140 for (
unsigned i = 0;
i<
n; ++
i)
2143 for (
unsigned i = 0;
i<
n; ++
i)
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);
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);
2165 const std::map<
Edge,std::vector<Triangles::const_iterator> >& edge2triangles_map)
const {
2167 Edge reverse_edge(edge.second,edge.first);
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();
2181 return oriented ?
Edge((*tr)[2],(*tr)[0]) :
Edge((*tr)[0],(*tr)[2]);
2183 return Edge((*tr)[
i],(*tr)[
i+1]);
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));
2195 const unsigned& idfirstvertex,
const unsigned& idsecondvertex,
2196 const unsigned& id_triangleface,
const unsigned&
n, std::vector<unsigned>& faceinfo)
const {
2198 faceinfo.push_back(idfirstvertex);
2199 faceinfo.push_back(prev_edgeface);
2200 faceinfo.push_back(top2bottomvertexid(idfirstvertex));
2201 faceinfo.push_back(id_triangleface+
n-2);
2202 faceinfo.push_back(top2bottomvertexid(idsecondvertex));
2203 faceinfo.push_back(next_edgeface);
2204 faceinfo.push_back(idsecondvertex);
2205 faceinfo.push_back(id_triangleface);
2206 setupface(id_extedge,faceinfo);
2212 unsigned& prevValue,
unsigned& nextValue)
const
2216 for (
i=0;
i<numbercycle.size();++
i)
2217 if (numbercycle.at(
i)==centralValue)
2220 assert(
i<numbercycle.size());
2221 prevValue = (
i==0 ? numbercycle.back() : numbercycle.at(
i-1) );
2222 nextValue = (
i==numbercycle.size()-1 ? numbercycle.front() : numbercycle.at(
i+1) );
2231 return topid+nextraexternalvertices+nextrainternalvertices;
2236 if (topid<=ntriangles)
2237 return topid+ntriangles;
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();
2251 for(; itt!=ittE; ++itt) {
2253 triangle_2_id[triangle]=++itri;
2255 assert(triangle_2_id.size()==ntriangles);
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) {
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();
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)];
2279 for (
unsigned iedge=0;iedge<3;++iedge) {
2280 Edge e = GetEdge(&(*itt),iedge,
true);
2281 Edge eun = GetEdge(&(*itt),iedge,
false);
2283 bool hasextravertex(edges_with_extra_vertex.find(eun)!=edges_with_extra_vertex.end());
2284 bool isinternal(edges_external.find(eun)==edges_external.end());
2286 faceinfo_top.push_back(
e.first);
2288 assert(neighbourmap[
e]);
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]]);
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]);
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,
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,
2322 triangleid,
n,faceinfo_sides);
2329 assert(faceinfo_bottom.size()==8);
2330 assert(faceinfo_top.size()==8);
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]);
2341 setupface(triangleid,faceinfo_top);
2342 setupface(triangleid+ntriangles,faceinfo_bottom);
2343 faceinfo_top.clear();
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;
2380 std::cerr <<
"ERROR in SbPolyhedronArbitrary::AddFacet. Attempt to exceed maximum number of facets: "
2381 << m_nFacetCount << std::endl;
2390 std::cerr <<
"ERROR in SbPolyhedronArbitrary::AddFacet. Attempt to index vertex number which is out-of-range: ("
2391 << iv1 <<
"," << iv2 <<
"," << iv3 <<
"," << iv4 <<
")" << std::endl;
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;
2402 m_pF[m_nFacetCount] =
SbFacet(iv1, 0, iv2, 0, iv3, 0, iv4, 0);