44 #ifdef BP_GEANT4 //G.Barrand
45 #include "G4Plane3D.hh"
46 #include "G4Point3D.hh"
47 #include "G4Normal3D.hh"
54 #define ExtNode HEPVis_ExtNode
55 #define ExtEdge HEPVis_ExtEdge
56 #define ExtFace HEPVis_ExtFace
57 #define FaceList HEPVis_FaceList
58 #define ExtPolyhedron HEPVis_ExtPolyhedron
59 #define BooleanProcessor HEPVis_BooleanProcessor
61 #define HepPolyhedron SbPolyhedron
62 #define G4Facet SbFacet
82 #define INITIAL_SIZE 200
83 #define CRAZY_POINT HVPoint3D(-10.e+6, -10.e+6, -10.e+6)
85 #define GRANULARITY 10.e+5 //G.Barrand : rm the trailing ;
87 #define SWAP(A,B) w = A; A = B; B = w
89 #define OP_UNION 0 // Operations
90 #define OP_INTERSECTION 1
91 #define OP_SUBTRACTION 2
93 #define OUT_OF_PLANE 0 // Face vs face statuses
95 #define INTERSECTION 2
97 #define NON_PLANAR_FACE 4
99 #define UNKNOWN_FACE 0 // Face statuses
100 #define ORIGINAL_FACE -1
102 #define UNSUITABLE_FACE -3
103 #define DEFECTIVE_FACE -4
160 if (&
node ==
this)
return *
this;
177 ExtEdge(
int k1=0,
int k2=0,
int kface1=0,
int kface2=0,
int kvis=0) :
187 if (&edge ==
this)
return *
this;
218 ExtFace(std::vector<ExtEdge>& a_edges,
int iedge)
237 if (&face ==
this)
return *
this;
307 return HepPolyhedron::operator = (from);
356 double &
a,
double &
b,
double &
c)
const;
360 int checkTriangle(
int iedge1,
int iedge2,
int ix,
int iy)
const;
383 for (
int i = 0;
i < 3;
i++) {
420 int iEprev, iEcur, iEnext;
422 iEprev = 0; iEcur =
iold;
433 iEprev = 0; iEcur =
inew;
444 #ifdef BP_GEANT4 //G.Barrand
452 double dx,
double dy,
double dz)
462 int i,
k, nnode, iNodes[5], iVis[4], iFaces[4];
463 int dnode =
m_nodes.size() - 1;
464 int dface =
m_faces.size() - 1;
473 for (
i=1;
i <=
p.GetNoVertices();
i++) {
474 #ifdef BP_GEANT4 //G.Barrand
475 ppp =
p.GetVertex(
i);
476 ppp.setX(ppp.x()+
dx);
477 ppp.setY(ppp.y()+
dy);
478 ppp.setZ(ppp.z()+dz);
480 ppp =
p.GetVertexFast(
i);
488 for (
int iface=1; iface <=
p.GetNoFacets(); iface++) {
493 p.GetFacet(iface, nnode, iNodes, iVis, iFaces);
494 for (
i=0;
i<nnode;
i++) {
498 if (iNodes[
i] < 1 || iNodes[
i] >
p.GetNoVertices()) {
502 <<
"BooleanProcessor::takePolyhedron : problem 1."
506 if (iFaces[
i] < 1 || iFaces[
i] >
p.GetNoFacets()) {
510 <<
"BooleanProcessor::takePolyhedron : problem 2."
520 iNodes[nnode] = iNodes[0];
522 for (
i=0;
i<nnode;
i++) {
525 iface+dface, iFaces[
i], iVis[
i]));
533 for (
i=0;
i<3;
i++) {
537 for (
i=1;
i<nnode;
i++) {
539 for (
k=0;
k<3;
k++) {
553 for (
i=0;
i<nnode;
i++) { point +=
m_nodes[iNodes[
i]].v; }
554 if (nnode > 1) point *= (1./nnode);
580 double rmin1[3], rmax1[3];
581 double rmin2[3], rmax2[3];
585 for (
i=0;
i<3;
i++) {
595 for (
i=0;
i<3;
i++) {
605 for (
i=0;
i<3;
i++) {
614 for (
i=0;
i<3;
i++) {
615 m_rmin[
i] = (rmin1[
i] > rmin2[
i]) ? rmin1[
i] : rmin2[
i];
616 m_rmax[
i] = (rmax1[
i] < rmax2[
i]) ? rmax1[
i] : rmax2[
i];
623 for (
i=0;
i<3;
i++) {
624 if ((rmax1[
i]-rmin1[
i]) > del1) del1 = rmax1[
i]-rmin1[
i];
625 if ((rmax2[
i]-rmin2[
i]) > del2) del2 = rmax2[
i]-rmin2[
i];
627 return ((del1 < del2) ? del1 : del2) /
GRANULARITY;
640 int i, outflag, iface = ifaces, *prev;
656 for (
i=0;
i<3;
i++) {
664 int npos = 0, nneg = 0;
666 for (
i=0;
i<8;
i++) {
671 if (npos == 8 || nneg == 8) outflag = 1;
700 int i, nnode, npos = 0, nneg = 0, nzer = 0;
705 nnode = (
m_faces[iface].iedges[3] == 0) ? 3 : 4;
706 for (
i=0;
i<nnode;
i++) {
719 if (npos == nnode || nneg == nnode)
return OUT_OF_PLANE;
727 int ie1 = 0, ie2 = 0,
s1 = 0,
s2 = 0,
status, nint = 0;
728 enum { PLUS_MINUS, MINUS_PLUS, ZERO_ZERO, ZERO_PLUS, ZERO_MINUS };
731 for (
i=0;
i<nnode;
i++) {
733 if (
dd[
i+1] >= 0)
continue;
735 }
else if (
dd[
i] < 0) {
736 if (
dd[
i+1] <= 0)
continue;
756 if (
s1 != ZERO_ZERO &&
s2 != ZERO_ZERO) {
758 int iedge, i1 = 0, i2 = 0, ii[2];
759 double d1 = 0.,
d2 = 0., d3 = 0.;
760 ii[0] = ie1; ii[1] = ie2;
761 for (
i=0;
i<2;
i++) {
762 iedge =
m_faces[iface].iedges[ii[
i]];
784 if (
s1 == MINUS_PLUS ||
s1 == ZERO_PLUS) {
794 edge.
inext = (
s1 == ZERO_ZERO) ? ie1+1 : ie2+1;
795 if (
s1 == ZERO_PLUS ||
s2 == ZERO_MINUS) {
820 if (i1 == i2)
return;
824 if (i1 == ilast) { i1 = i2;
m_nodes.pop_back();
return; }
825 if (i2 == ilast) { i2 = i1; }
826 if (i3 == ilast) { i3 = i1; }
827 if (i4 == ilast) { i4 = i1; }
844 for (
i=0;
i<3;
i++) {
847 if (
d >
dd) {
dd =
d; ii =
i; }
855 if (t3 <= t1+m_del || t4 >=
t2-
m_del)
return 0;
885 if (i1 < i2) { i2 = i1; }
886 else if (i1 > i2) { i1 = i2; }
887 else { iedges[1] = 0; }
888 if (iedges[0] == iedges[1])
return;
890 int ie1, ie2,
inode = i1;
892 for (
int i=0;
i<2;
i++) {
896 if ((ie1 = iedges[
i]) == 0)
continue;
965 <<
"BooleanProcessor::caseIE : unimplemented case"
983 <<
"BooleanProcessor::caseEE : unimplemented case"
1005 for (
int i=0;
i<3;
i++) {
1061 int iedge =
m_faces[iface].inew;
1100 #define INSERT_EDGE_TO_THE_LIST(A) \
1101 *ilink = A; ilink = &m_edges[A].inext; *ilink = 0
1106 if (face.
inew == 0)
break;
1122 if (edge_i.
i1 == edge_cur.
i2)
break;
1129 if (edge_i.
i1 == edge_cur.
i2) {
1140 if (
m_edges[icur].i2 == ifirst) {
break; }
else {
continue; }
1145 <<
"BooleanProcessor::assembleFace(" << iface <<
") : "
1146 <<
"could not find next edge of the contour"
1159 if (
what == 0 && ioldflag == 0 && iedge > 0) {
1161 if (
m_edges[iedge].inext > 0) {
1183 iface2 =
m_edges[iedge].iface2;
1202 if (
m_faces[iface].inew > 0) {
1327 int i, iedge, iface;
1336 iface =
m_edges[iedge].iface2;
1357 iface =
m_edges[iedge].iface2;
1400 double &
a,
double &
b,
double &
c)
const
1414 w = ::fabs(
a)+::fabs(
b);
1430 double a1, b1,
c1, a2, b2,
c2,
d1,
d2;
1434 findABC(
x[0],
y[0],
x[1],
y[1], a1, b1,
c1);
1435 findABC(
x[1],
y[1],
x[2],
y[2], a2, b2,
c2);
1436 d1 = a1*
x[4] + b1*
y[4] +
c1;
1437 d2 = a2*
x[4] + b2*
y[4] +
c2;
1438 if (
d1 <= m_del &&
d2 <= m_del)
return 1;
1439 if (! (
d1 > m_del &&
d2 > m_del)) {
1440 if ( a1*
x[2] + b1*
y[2] +
c1 >= -m_del)
return 1;
1445 findABC(
x[3],
y[3],
x[4],
y[4], a1, b1,
c1);
1446 findABC(
x[4],
y[4],
x[5],
y[5], a2, b2,
c2);
1447 d1 = a1*
x[1] + b1*
y[1] +
c1;
1448 d2 = a2*
x[1] + b2*
y[1] +
c2;
1449 if (
d1 <= m_del &&
d2 <= m_del)
return 1;
1450 if (!(
d1 > m_del &&
d2 > m_del)) {
1451 if ( a1*
x[5] + b1*
y[5] +
c1 >= -m_del)
return 1;
1469 x1 = m_nodes[i1].v[ix];
1470 y1 = m_nodes[i1].v[iy];
1471 x2 = m_nodes[i2].v[ix];
1472 y2 = m_nodes[i2].v[iy];
1477 int icontour, iedge, k1, k2;
1478 double x3, y3, x4, y4, a2, b2,
c2,
d1,
d2;
1479 for(icontour=0; icontour<(
int)m_external_contours.size(); icontour++) {
1480 iedge = m_external_contours[icontour];
1482 k1 = m_edges[iedge].i1;
1483 k2 = m_edges[iedge].i2;
1484 iedge = m_edges[iedge].inext;
1485 if (k1 == i1 || k2 == i1)
continue;
1486 if (k1 == i2 || k2 == i2)
continue;
1487 x3 = m_nodes[k1].v[ix];
1488 y3 = m_nodes[k1].v[iy];
1489 x4 = m_nodes[k2].v[ix];
1490 y4 = m_nodes[k2].v[iy];
1491 d1 = a1*x3 + b1*y3 +
c1;
1492 d2 = a1*x4 + b1*y4 +
c1;
1493 if (
d1 > m_del &&
d2 > m_del)
continue;
1494 if (
d1 < -m_del &&
d2 < -m_del)
continue;
1496 findABC(x3, y3, x4, y4, a2, b2,
c2);
1499 if (
d1 > m_del &&
d2 > m_del)
continue;
1500 if (
d1 < -m_del &&
d2 < -m_del)
continue;
1507 for(icontour=0; icontour<(
int)m_internal_contours.size(); icontour++) {
1508 iedge = m_internal_contours[icontour];
1510 k1 = m_edges[iedge].i1;
1511 k2 = m_edges[iedge].i2;
1512 iedge = m_edges[iedge].inext;
1513 if (k1 == i1 || k2 == i1)
continue;
1514 if (k1 == i2 || k2 == i2)
continue;
1515 x3 = m_nodes[k1].v[ix];
1516 y3 = m_nodes[k1].v[iy];
1517 x4 = m_nodes[k2].v[ix];
1518 y4 = m_nodes[k2].v[iy];
1519 d1 = a1*x3 + b1*y3 +
c1;
1520 d2 = a1*x4 + b1*y4 +
c1;
1521 if (
d1 > m_del &&
d2 > m_del)
continue;
1522 if (
d1 < -m_del &&
d2 < -m_del)
continue;
1524 findABC(x3, y3, x4, y4, a2, b2,
c2);
1527 if (
d1 > m_del &&
d2 > m_del)
continue;
1528 if (
d1 < -m_del &&
d2 < -m_del)
continue;
1545 int i1ext, i2ext, i1int, i2int,
i,
k[6];
1557 for (
i=0;
i<3;
i++) {
1571 for (
i=3;
i<6;
i++) {
1621 k[0] = m_edges[iedge1].i1;
1622 k[1] = m_edges[iedge1].i2;
1623 k[2] = m_edges[iedge2].i2;
1624 for (
int i=0;
i<3;
i++) {
1625 x[
i] = m_nodes[
k[
i]].v[ix];
1626 y[
i] = m_nodes[
k[
i]].v[iy];
1631 findABC(
x[2],
y[2],
x[0],
y[0], a1, b1,
c1);
1632 if (a1*
x[1]+b1*
y[1]+
c1 <= 0.1*m_del)
return 1;
1637 double a2, b2,
c2, a3, b3,
c3;
1639 findABC(
x[0],
y[0],
x[1],
y[1], a2, b2,
c2);
1640 findABC(
x[1],
y[1],
x[2],
y[2], a3, b3,
c3);
1643 iedge = m_edges[iedge].inext;
1644 if (m_edges[iedge].inext == iedge1)
return 0;
1645 inode = m_edges[iedge].i2;
1646 if (
inode ==
k[0])
continue;
1647 if (
inode ==
k[1])
continue;
1648 if (
inode ==
k[2])
continue;
1649 x[1] = m_nodes[
inode].v[ix];
1650 y[1] = m_nodes[
inode].v[iy];
1651 if (a1*
x[1]+b1*
y[1]+
c1 < -0.1*m_del)
continue;
1652 if (a2*
x[1]+b2*
y[1]+
c2 < -0.1*m_del)
continue;
1653 if (a3*
x[1]+b3*
y[1]+
c3 < -0.1*m_del)
continue;
1676 int ipnext = ihead, nnode = 1;
1678 if (
m_edges[ipnext].inext > 0) {
1679 ipnext =
m_edges[ipnext].inext;
1682 m_edges[ipnext].inext = ihead;
1692 int iedge1, iedge2, iedge3, istart = 0;
1694 iedge1 =
m_edges[ipnext].inext;
1695 iedge2 =
m_edges[iedge1].inext;
1707 iedge3 =
m_edges[iedge2].inext;
1720 }
else if (istart == iedge1) {
1724 <<
"BooleanProcessor::triangulateContour : "
1725 <<
"could not generate a triangle (infinite loop)"
1734 ipnext =
m_edges[ipnext].inext;
1742 int iface1 =
m_edges[iedge1].iface1;
1756 m_edges[iedge1].iface1 = iface2;
1757 m_edges[iedge2].iface1 = iface2;
1758 ipnext =
m_edges[ipnext].inext;
1777 int iedge =
m_faces[iface].iold;
1788 <<
"BooleanProcessor::modifyReference : could not find the edge, "
1789 <<
"iface=" << iface <<
", i1,i2=" << i1 <<
"," << i2 <<
", iref=" << iref
1808 #ifdef BP_GEANT4 //G.Barrand
1816 if (::fabs(normal[1]) > ::fabs(normal[iz])) iz = 1;
1817 if (::fabs(normal[2]) > ::fabs(normal[iz])) iz = 2;
1818 if (normal[iz] > 0) {
1819 ix = (iz+1)%3; iy = (ix+1)%3;
1821 iy = (iz+1)%3; ix = (iy+1)%3;
1829 int i1, i2, ifirst, iedge, icontour =
m_faces[iface].iold;
1830 while (icontour > 0) {
1840 z += node_1.
v[ix]*node_2.
v[iy]-node_2.
v[ix]*node_1.
v[iy];
1853 <<
"BooleanProcessor::triangulateFace : too small contour"
1857 icontour =
m_edges[iedge].inext;
1865 <<
"BooleanProcessor::triangulateFace : broken contour"
1886 <<
"BooleanProcessor::triangulateFace : "
1887 <<
"could not merge internal contour " << kint
1901 <<
"BooleanProcessor::triangulateFace : "
1902 <<
"triangulateContour failed."
1912 for (
int ifa=nface; ifa<(
int)
m_faces.size(); ifa++) {
1915 if (
m_edges[iedge].iface1 != ifa) {
1919 <<
"BooleanProcessor::triangulateFace : wrong reference to itself, "
1920 <<
"iface=" << ifa <<
", iface1=" <<
m_edges[iedge].iface1
1923 }
else if (
m_edges[iedge].iface2 > 0) {
1926 }
else if (
m_edges[iedge].iface2 < 0) {
1944 int i, iedge, nnode = 0, nface = 0;
1970 if (nface == 0)
return polyhedron;
1971 polyhedron.AllocateMemory(nnode, nface);
1981 int k,
v[4]={0},
f[4]={0};
1983 if (
m_faces[
i].inew == 0)
continue;
1984 v[3] =
f[3] =
k = 0;
1988 std::cerr <<
"BooleanProcessor::createPolyhedron : too many edges" << std::endl;
1998 std::cerr <<
"BooleanProcessor::createPolyhedron : "
1999 <<
"face has only " <<
k <<
" edges"
2062 <<
"BooleanProcessor: corrupted input polyhedron"
2074 <<
"BooleanProcessor: intersection with empty polyhedron"
2079 <<
"BooleanProcessor: subtraction from empty polyhedron"
2091 <<
"BooleanProcessor: intersection with empty polyhedron"
2111 #define PROCESSOR_ERROR(a_what)
2113 unsigned int try_count = 1;
2189 <<
"BooleanProcessor::execute : unknown faces !!!"
2212 <<
"BooleanProcessor::execute : had converged."
2235 <<
"BooleanProcessor::execute : try another tilt..."
2242 #undef PROCESSOR_ERROR //G.Barrand
2471 std::cerr <<
"nodes : " <<
number << std::endl;
2474 std::cerr <<
" " <<
index
2475 <<
" x = " <<
node.v[0]
2476 <<
" y = " <<
node.v[1]
2477 <<
" z = " <<
node.v[2]