ATLAS Offline Software
Loading...
Searching...
No Matches
BooleanProcessor.h
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3*/
4
5/*
6 * PLEASE NOTE: This file has been copied and pasted from GEANT4, rel. 4.9.6.cand00
7 *
8 * original source file:
9 * /afs/cern.ch/sw/geant4/releases/geant4.9.6.cand00/share/source/graphics_reps/src/BooleanProcessor.src
10 *
11 *
12 * Update: Riccardo-Maria BIANCHI <rbianchi@cern.ch> - 2012-11-16
13 *
14 *
15 * NOTES:
16 *
17 * ===================
18 * VP1 Customization:
19 *
20 * - look into the code for comments "// VP1 change"
21 *
22 *
23*/
24
25
26/***********************************************************************
27 * *
28 * Name: BooleanProcessor Date: 10.12.99 *
29 * Author: E.Chernyaev Revised: *
30 * *
31 * Function: Internal class for executing boolean operations *
32 * on Polyhedra *
33 * *
34 ***********************************************************************/
35
36
37
38
39/* // VP1 change
40//G.Barrand : begin
41#define BP_GEANT4
42*/
43
44#ifdef BP_GEANT4 //G.Barrand
45#include "G4Plane3D.hh"
46#include "G4Point3D.hh"
47#include "G4Normal3D.hh"
48typedef G4Plane3D HVPlane3D;
49typedef G4Point3D HVPoint3D;
50typedef G4Normal3D HVNormal3D;
51
52#else //BP_HEPVIS
53
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
60
61#define HepPolyhedron SbPolyhedron
62#define G4Facet SbFacet
63
64// VP1 change
65//#include <HEPVis/SbPlane.h>
66#include <VP1HEPVis/SbPlane.h>
68//---
69
70
71#endif
72
73
74
75
76//using namespace HepGeom;
77
78//#define BP_DEBUG
79
80//G.Barrand : end
81
82#define INITIAL_SIZE 200
83#define CRAZY_POINT HVPoint3D(-10.e+6, -10.e+6, -10.e+6)
84//#define GRANULARITY 10.e+5;
85#define GRANULARITY 10.e+5 //G.Barrand : rm the trailing ;
86
87#define SWAP(A,B) w = A; A = B; B = w
88
89#define OP_UNION 0 // Operations
90#define OP_INTERSECTION 1
91#define OP_SUBTRACTION 2
92
93#define OUT_OF_PLANE 0 // Face vs face statuses
94#define ON_PLANE 1
95#define INTERSECTION 2
96#define EDGE 3
97#define NON_PLANAR_FACE 4
98
99#define UNKNOWN_FACE 0 // Face statuses
100#define ORIGINAL_FACE -1
101#define NEW_FACE -2
102#define UNSUITABLE_FACE -3
103#define DEFECTIVE_FACE -4
104
105
106
107// -------------------------------------------- Simplified STL vector ---
108//G.Barrand : begin
109#include <vector>
110/*
111template<class T>
112class vector {
113 private:
114 int cur_size, max_size;
115 T * v;
116
117 public:
118 vector(): cur_size(0), max_size(INITIAL_SIZE), v(new T[INITIAL_SIZE]) {}
119 ~vector() { delete [] v; }
120
121 void clear() { cur_size = 0; }
122 int size() const { return cur_size; }
123 T & operator[](int i) { return v[i]; }
124 const T & operator[](int i) const { return v[i]; }
125 T & front() { return v[0]; }
126 const T & front() const { return v[0]; }
127 T & back() { return v[cur_size-1]; }
128 const T & back() const { return v[cur_size-1]; }
129 void pop_back() { cur_size--; }
130 void push_back(const T & a) {
131 if (cur_size == max_size) {
132 T * w = v;
133 max_size *= 2;
134 v = new T[max_size];
135 for (int i=0; i<cur_size; i++) v[i] = w[i];
136 v[cur_size++] = a;
137 delete [] w;
138 }else{
139 v[cur_size++] = a;
140 }
141 }
142};
143*/
144//G.Barrand : end
145
146// ---------------------------------------------------- Extended node ---
147class ExtNode {
148 public:
150 int s;
151
152 public:
153 ExtNode(HVPoint3D vertex=HVPoint3D(), int status=0)
154 : v(vertex), s(status) {}
156
157 ExtNode(const ExtNode & node) : v(node.v), s(node.s) {}
158
160 if (&node == this) return *this;
161 v = node.v;
162 s = node.s;
163 return *this;
164 }
165};
166
167// ---------------------------------------------------- Extended edge ---
168class ExtEdge {
169 public:
170 int i1, i2; // end points
171 int iface1; // native face
172 int iface2; // neighbouring face
173 int ivis; // visibility: +1 (visible), -1 (invisible)
174 int inext; // index of next edge
175
176 public:
177 ExtEdge(int k1=0, int k2=0, int kface1=0, int kface2=0, int kvis=0) :
178 i1(k1), i2(k2), iface1(kface1), iface2(kface2), ivis(kvis), inext(0) {}
179
181
182 ExtEdge(const ExtEdge & edge) :
183 i1(edge.i1), i2(edge.i2), iface1(edge.iface1), iface2(edge.iface2),
184 ivis(edge.ivis), inext(edge.inext) {}
185
186 ExtEdge & operator=(const ExtEdge & edge) {
187 if (&edge == this) return *this;
188 i1 = edge.i1;
189 i2 = edge.i2;
190 iface1 = edge.iface1;
191 iface2 = edge.iface2;
192 ivis = edge.ivis;
193 inext = edge.inext;
194 return *this;
195 }
196
197 void invert() {
198 int w;
199 SWAP(i1, i2);
200 }
201};
202
203// ---------------------------------------------------- Extended face ---
204class ExtFace {
205 private:
206 std::vector<ExtEdge>& m_edges; //G.Barrand // VP1 change
207 public:
208 int iedges[4]; // indices of original edges
209 HVPlane3D plane; // face plane
210 double rmin[3], rmax[3]; // bounding box
211 int iold; // head of the list of the original edges
212 int inew; // head of the list of the new edges
213 int iprev; // index of previous face
214 int inext; // index of next face
215
216 public:
217 //G.Barrand : ExtFace(int iedge=0) : iold(iedge), inew(0), iprev(iprev), inext(0) {}
218 ExtFace(std::vector<ExtEdge>& a_edges,int iedge)
219 : m_edges(a_edges), iold(iedge), inew(0), iprev(0), inext(0) {
220 //G.Barrand : initialize arrays to quiet valgrind.
221 {for (int i=0; i<4; i++) { iedges[i] = 0; }}
222 {for (int i=0; i<3; i++) { rmin[i] = 0; rmax[i] = 0; }}
223 }
225
226 ExtFace(const ExtFace & face) :
227 m_edges(face.m_edges), //G.Barrand
228 plane(face.plane), iold(face.iold), inew(face.inew),
229 iprev(face.iprev), inext(face.inext)
230 {
231 int i;
232 for (i=0; i<4; i++) { iedges[i] = face.iedges[i]; }
233 for (i=0; i<3; i++) { rmin[i] = face.rmin[i]; rmax[i] = face.rmax[i]; }
234 }
235
236 ExtFace & operator=(const ExtFace & face) {
237 if (&face == this) return *this;
238 //FIXME : edges(face.edges) ???? //G.Barrand
239 int i;
240 for (i=0; i<4; i++) { iedges[i] = face.iedges[i]; }
241 plane = face.plane;
242 for (i=0; i<3; i++) { rmin[i] = face.rmin[i]; rmax[i] = face.rmax[i]; }
243 iold = face.iold;
244 inew = face.inew;
245 iprev = face.iprev;
246 inext = face.inext;
247 return *this;
248 }
249
250 void invert();
251};
252
253// ---------------------------------------------------- Global arrays ---
254//G.Barrand : MacIntel : crash with g++-4.0.1 with -O on some subtract.
255// Anyway static of objects is proved to be not safe.
256// We put the below vector as members of BooleanProcessor.
257//GB static std::vector<ExtNode> nodes; // vector of nodes
258//GB static std::vector<ExtEdge> edges; // vector of edges
259//GB static std::vector<ExtFace> faces; // vector of faces
260
261// ---------------------------------------------------- List of faces ---
262class FaceList {
263 private:
264 std::vector<ExtFace>& m_faces; //G.Barrad : end // VP1 change
265 private:
268
269 public:
270 //G.Barrand : FaceList() : ihead(0), ilast(0) {}
271 FaceList(std::vector<ExtFace>& a_faces) : m_faces(a_faces),m_ihead(0),m_ilast(0) {}
273
274 void clean() { m_ihead = 0; m_ilast = 0; }
275 int front() { return m_ihead; }
276
277 void push_back(int i) {
278 if (m_ilast == 0) { m_ihead = i; } else { m_faces[m_ilast].inext = i; }
279 ExtFace& face = m_faces[i]; //G.Barrand : optimize.
280 face.iprev = m_ilast;
281 face.inext = 0;
282 m_ilast = i;
283 }
284
285 void remove(int i) {
286 ExtFace& face = m_faces[i]; //G.Barrand : optimize.
287 if (m_ihead == i) {
288 m_ihead = face.inext;
289 }else{
290 m_faces[face.iprev].inext = face.inext;
291 }
292 if (m_ilast == i) {
293 m_ilast = face.iprev;
294 }else{
295 m_faces[face.inext].iprev = face.iprev;
296 }
297 face.iprev = 0;
298 face.inext = 0;
299 }
300};
301
302// --------------------- Polyhedron with extended access to
303// its members from the BooleanProcessor class ---
305 friend class BooleanProcessor;
306 virtual HepPolyhedron& operator = (const HepPolyhedron& from) {
307 return HepPolyhedron::operator = (from);
308 }
309};
310
311// ----------------------------------------- Boolean processor class ---
313 private:
314 static int s_ishift; //G.Barrand // VP1 change
315 std::vector<ExtNode> m_nodes; // vector of nodes //G.Barrand // VP1 change
316 std::vector<ExtEdge> m_edges; // vector of edges //G.Barrand // VP1 change
317 std::vector<ExtFace> m_faces; // vector of faces //G.Barrand // VP1 change
318 private:
319 int m_processor_error; // is set in case of error // VP1 change
320 int m_operation; // 0 (union), 1 (intersection), 2 (subtraction) // VP1 change
321 int m_ifaces1, m_ifaces2; // lists of faces // VP1 change
322 int m_iout1, m_iout2; // lists of faces with status "out" // VP1 change
323 int m_iunk1, m_iunk2; // lists of faces with status "unknown" // VP1 change
324 double m_rmin[3], m_rmax[3]; // intersection of bounding boxes // VP1 change
325 double m_del; // precision (tolerance) // VP1 change
326
327 FaceList m_result_faces; // list of accepted faces // VP1 change
328 FaceList m_suitable_faces; // list of suitable faces // VP1 change
329 FaceList m_unsuitable_faces; // list of unsuitable faces // VP1 change
330 FaceList m_unknown_faces; // list of unknown faces // VP1 change
331
332 std::vector<int> m_external_contours; // heads of external contours // VP1 change
333 std::vector<int> m_internal_contours; // heads of internal contours // VP1 change
334
335 private:
336 void takePolyhedron(const HepPolyhedron & p, double, double, double);
337 double findMinMax();
338 void selectOutsideFaces(int & ifaces, int & iout);
339 int testFaceVsPlane(ExtEdge & edge);
340 void renumberNodes(int & i1, int & i2, int & i3, int & i4);
341 int testEdgeVsEdge(ExtEdge & edge1, ExtEdge & edge2);
342 void removeJunkNodes() { while(m_nodes.back().s != 0) m_nodes.pop_back(); }
343 void divideEdge(int & i1, int & i2);
344 void insertEdge(const ExtEdge & edge);
345 void caseII(ExtEdge & edge1, ExtEdge & edge2);
346 void caseIE(ExtEdge & edge1, ExtEdge & edge2);
347 void caseEE(ExtEdge & edge1, ExtEdge & edge2);
348 void testFaceVsFace(int iface1, int iface2);
349 void invertNewEdges(int iface);
350 void checkDoubleEdges(int iface);
351 void assembleFace(int what, int iface);
352 void assembleNewFaces(int what, int ihead);
353 void initiateLists();
354 void assemblePolyhedra();
355 void findABC(double x1, double y1, double x2, double y2,
356 double &a, double &b, double &c) const;
357 int checkDirection(double *x, double *y) const;
358 int checkIntersection(int ix, int iy, int i1, int i2) const;
359 void mergeContours(int ix, int iy, int kext, int kint);
360 int checkTriangle(int iedge1, int iedge2, int ix, int iy) const;
361 void triangulateContour(int ix, int iy, int ihead);
362 void modifyReference(int iface, int i1, int i2, int iref);
363 void triangulateFace(int iface);
364 HepPolyhedron createPolyhedron();
365
366 public:
367 //G.Barrand : BooleanProcessor() {}
368 BooleanProcessor() //G.Barrand
369 :m_processor_error(0) // The next few fields are work space, initialised
370 ,m_operation(0) // here to prevent Coverity warnings.
371 ,m_ifaces1(0) // "
372 ,m_ifaces2(0) // "
373 ,m_iout1(0) // "
374 ,m_iout2(0) // "
375 ,m_iunk1(0) // "
376 ,m_iunk2(0) // "
377 ,m_del(0.) // "
382 { // rmin, rmax also initialised here to prevent Coverity warnings.
383 for (int i = 0; i < 3; i++) {
384 m_rmin[i] = 0.;
385 m_rmax[i] = 0.;
386 }
387 }
388
390
391 HepPolyhedron execute(int op,
392 const HepPolyhedron &a,
393 const HepPolyhedron &b,
394 int& err);
395
396 void draw();
397 void draw_edge(int, int);
398 void draw_contour(int, int, int);
399 void draw_faces(int, int, int);
400 void print_face(int);
401 void print_edge(int);
403
404 void dump(); //G.Barrand
405 static int get_shift(); //G.Barrand
406 static void set_shift(int); //G.Barrand
407 static int get_num_shift(); //G.Barrand
408};
409
410inline void ExtFace::invert()
411/***********************************************************************
412 * *
413 * Name: ExtFace::invert() Date: 28.02.00 *
414 * Author: E.Chernyaev Revised: *
415 * *
416 * Function: Invert face *
417 * *
418 ***********************************************************************/
419{
420 int iEprev, iEcur, iEnext;
421
422 iEprev = 0; iEcur = iold;
423 while (iEcur > 0) {
424 ExtEdge& edge = m_edges[iEcur]; //G.Barrand : optimize.
425 edge.invert();
426 iEnext = edge.inext;
427 edge.inext = iEprev;
428 iEprev = iEcur;
429 iEcur = iEnext;
430 }
431 if (iold > 0) iold = iEprev;
432
433 iEprev = 0; iEcur = inew;
434 while (iEcur > 0) {
435 ExtEdge& edge = m_edges[iEcur]; //G.Barrand : optimize.
436 edge.invert();
437 iEnext = edge.inext;
438 edge.inext = iEprev;
439 iEprev = iEcur;
440 iEcur = iEnext;
441 }
442 if (inew > 0) inew = iEprev;
443
444#ifdef BP_GEANT4 //G.Barrand
445 plane = HVPlane3D(-plane.a(), -plane.b(), -plane.c(), -plane.d());
446#else
447 plane = HVPlane3D(-plane.getNormal(), -plane.getDistanceFromOrigin());
448#endif
449}
450
452 double dx, double dy, double dz)
453/***********************************************************************
454 * *
455 * Name: BooleanProcessor::takePolyhedron Date: 16.12.99 *
456 * Author: E.Chernyaev Revised: *
457 * *
458 * Function: Transfer Polyhedron to internal representation *
459 * *
460 ***********************************************************************/
461{
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;
465
466 // S E T N O D E S
467
468 // for (i=1; i <= p.GetNoVertices(); i++) {
469 // nodes.push_back(ExtNode(p.GetVertex(i)));
470 // }
471
472 HVPoint3D ppp;
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);
479#else
480 ppp = p.GetVertexFast(i);
481 ppp += HVPoint3D(dx,dy,dz);
482#endif
483 m_nodes.push_back(ExtNode(ppp));
484 }
485
486 // S E T F A C E S
487
488 for (int iface=1; iface <= p.GetNoFacets(); iface++) {
489 m_faces.push_back(ExtFace(m_edges,m_edges.size()));
490
491 // S E T F A C E N O D E S
492
493 p.GetFacet(iface, nnode, iNodes, iVis, iFaces);
494 for (i=0; i<nnode; i++) {
495 //if (iNodes[i] < 1 || iNodes[i] > p.GetNoVertices()) m_processor_error = 1;
496 //if (iFaces[i] < 1 || iFaces[i] > p.GetNoFacets()) m_processor_error = 1;
497
498 if (iNodes[i] < 1 || iNodes[i] > p.GetNoVertices()) { //G.Barrand
500#ifdef BP_DEBUG
501 G4cerr
502 << "BooleanProcessor::takePolyhedron : problem 1."
503 << G4endl;
504#endif
505 }
506 if (iFaces[i] < 1 || iFaces[i] > p.GetNoFacets()) { //G.Barrand
508#ifdef BP_DEBUG
509 G4cerr
510 << "BooleanProcessor::takePolyhedron : problem 2."
511 << G4endl;
512#endif
513 }
514 iNodes[i] += dnode;
515 iFaces[i] += dface;
516 }
517
518 // S E T E D G E S
519
520 iNodes[nnode] = iNodes[0];
521 m_faces.back().iedges[3] = 0;
522 for (i=0; i<nnode; i++) {
523 m_faces.back().iedges[i] = m_edges.size();
524 m_edges.push_back(ExtEdge(iNodes[i], iNodes[i+1],
525 iface+dface, iFaces[i], iVis[i]));
526 m_edges.back().inext = m_edges.size();
527 }
528 m_edges.back().inext = 0;
529
530 // S E T F A C E M I N - M A X
531
532 ExtFace& face = m_faces.back(); //G.Barrand : optimize.
533 for (i=0; i<3; i++) {
534 face.rmin[i] = m_nodes[iNodes[0]].v[i];
535 face.rmax[i] = m_nodes[iNodes[0]].v[i];
536 }
537 for (i=1; i<nnode; i++) {
538 ExtNode& node = m_nodes[iNodes[i]]; //G.Barrand : optimize.
539 for (k=0; k<3; k++) {
540 if (face.rmin[k] > node.v[k])
541 face.rmin[k] = node.v[k];
542 if (face.rmax[k] < node.v[k])
543 face.rmax[k] = node.v[k];
544 }
545 }
546
547 // S E T F A C E P L A N E
548
549 HVNormal3D n = (m_nodes[iNodes[2]].v-m_nodes[iNodes[0]].v).cross
550 (m_nodes[iNodes[3]].v-m_nodes[iNodes[1]].v);
551 HVPoint3D point(0,0,0);
552
553 for (i=0; i<nnode; i++) { point += m_nodes[iNodes[i]].v; }
554 if (nnode > 1) point *= (1./nnode);
555
556
557 //G.Barrand : faces.back().plane = HVPlane3D(n.unit(), point);
558 m_faces.back().plane = HVPlane3D(n, point); //G.Barrand
559
560 // S E T R E F E R E N C E T O T H E N E X T F A C E
561
562 m_faces.back().inext = m_faces.size();
563 }
564 m_faces.back().inext = 0;
565}
566
568/***********************************************************************
569 * *
570 * Name: BooleanProcessor::findMinMax Date: 16.12.99 *
571 * Author: E.Chernyaev Revised: *
572 * *
573 * Function: Find min-max (bounding) boxes for polyhedra *
574 * *
575 ***********************************************************************/
576{
577 if (m_ifaces1 == 0 || m_ifaces2 == 0) return 0;
578
579 int i, iface;
580 double rmin1[3], rmax1[3];
581 double rmin2[3], rmax2[3];
582
583 // F I N D B O U N D I N G B O X E S
584
585 for (i=0; i<3; i++) {
586 rmin1[i] = m_faces[m_ifaces1].rmin[i];
587 rmax1[i] = m_faces[m_ifaces1].rmax[i];
588 rmin2[i] = m_faces[m_ifaces2].rmin[i];
589 rmax2[i] = m_faces[m_ifaces2].rmax[i];
590 }
591
592 iface = m_faces[m_ifaces1].inext;
593 while(iface > 0) {
594 ExtFace& face = m_faces[iface]; //G.Barrand
595 for (i=0; i<3; i++) {
596 if (rmin1[i] > face.rmin[i]) rmin1[i] = face.rmin[i];
597 if (rmax1[i] < face.rmax[i]) rmax1[i] = face.rmax[i];
598 }
599 iface = face.inext;
600 }
601
602 iface = m_faces[m_ifaces2].inext;
603 while(iface > 0) {
604 ExtFace& face = m_faces[iface]; //G.Barrand
605 for (i=0; i<3; i++) {
606 if (rmin2[i] > face.rmin[i]) rmin2[i] = face.rmin[i];
607 if (rmax2[i] < face.rmax[i]) rmax2[i] = face.rmax[i];
608 }
609 iface = face.inext;
610 }
611
612 // F I N D I N T E R S E C T I O N O F B O U N D I N G B O X E S
613
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];
617 }
618
619 // F I N D T O L E R A N C E
620
621 double del1 = 0;
622 double del2 = 0;
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];
626 }
627 return ((del1 < del2) ? del1 : del2) / GRANULARITY;
628}
629
630void BooleanProcessor::selectOutsideFaces(int & ifaces, int & iout)
631/***********************************************************************
632 * *
633 * Name: BooleanProcessor::selectOutsideFaces Date: 10.01.00 *
634 * Author: E.Chernyaev Revised: *
635 * *
636 * Function: Preselection of outside faces *
637 * *
638 ***********************************************************************/
639{
640 int i, outflag, iface = ifaces, *prev;
641 HVPoint3D mmbox[8] = { HVPoint3D(m_rmin[0],m_rmin[1],m_rmin[2]),
642 HVPoint3D(m_rmax[0],m_rmin[1],m_rmin[2]),
643 HVPoint3D(m_rmin[0],m_rmax[1],m_rmin[2]),
644 HVPoint3D(m_rmax[0],m_rmax[1],m_rmin[2]),
645 HVPoint3D(m_rmin[0],m_rmin[1],m_rmax[2]),
646 HVPoint3D(m_rmax[0],m_rmin[1],m_rmax[2]),
647 HVPoint3D(m_rmin[0],m_rmax[1],m_rmax[2]),
648 HVPoint3D(m_rmax[0],m_rmax[1],m_rmax[2]) };
649 prev = &ifaces;
650 while (iface > 0) {
651
652 // B O U N D I N G B O X vs B O U N D I N G B O X
653
654 outflag = 0;
655 ExtFace& face = m_faces[iface]; //G.Barrand : optimize.
656 for (i=0; i<3; i++) {
657 if (face.rmin[i] > m_rmax[i] + m_del) { outflag = 1; break; }
658 if (face.rmax[i] < m_rmin[i] - m_del) { outflag = 1; break; }
659 }
660
661 // B O U N D I N G B O X vs P L A N E
662
663 if (outflag == 0) {
664 int npos = 0, nneg = 0;
665 double d;
666 for (i=0; i<8; i++) {
667 d = face.plane.distance(mmbox[i]); //G.Barrand : optimize
668 if (d > +m_del) npos++;
669 if (d < -m_del) nneg++;
670 }
671 if (npos == 8 || nneg == 8) outflag = 1;
672 }
673
674 // U P D A T E L I S T S
675
676 if (outflag == 1) {
677 *prev = face.inext;
678 face.inext = iout;
679 iout = iface;
680 }else{
681 prev = &face.inext;
682 }
683
684 iface = *prev;
685 }
686}
687
689/***********************************************************************
690 * *
691 * Name: BooleanProcessor::testFaceVsPlane Date: 19.01.00 *
692 * Author: E.Chernyaev Revised: *
693 * *
694 * Function: Find result of intersection of face by plane *
695 * *
696 ***********************************************************************/
697{
698 int iface = edge.iface1;
699 HVPlane3D plane = m_faces[edge.iface2].plane;
700 int i, nnode, npos = 0, nneg = 0, nzer = 0;
701 double dd[5];
702
703 // F I N D D I S T A N C E S
704
705 nnode = (m_faces[iface].iedges[3] == 0) ? 3 : 4;
706 for (i=0; i<nnode; i++) {
707 dd[i] = plane.distance(m_nodes[m_edges[m_faces[iface].iedges[i]].i1].v);
708 if (dd[i] > m_del) {
709 npos++;
710 }else if (dd[i] < -m_del) {
711 nneg++;
712 }else{
713 nzer++; dd[i] = 0;
714 }
715 }
716
717 // S O M E S I M P L E C A S E S ( N O I N T E R S E C T I O N )
718
719 if (npos == nnode || nneg == nnode) return OUT_OF_PLANE;
720 if (nzer == 1 && nneg == 0) return OUT_OF_PLANE;
721 if (nzer == 1 && npos == 0) return OUT_OF_PLANE;
722 if (nzer == nnode) return ON_PLANE;
723 if (nzer == 3) return NON_PLANAR_FACE;
724
725 // F I N D I N T E R S E C T I O N
726
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 };
729
730 dd[nnode] = dd[0];
731 for (i=0; i<nnode; i++) {
732 if (dd[i] > 0) {
733 if (dd[i+1] >= 0) continue;
734 status = PLUS_MINUS;
735 }else if (dd[i] < 0) {
736 if (dd[i+1] <= 0) continue;
737 status = MINUS_PLUS;
738 }else{
739 status = ZERO_ZERO;
740 if (dd[i+1] > 0) status = ZERO_PLUS;
741 if (dd[i+1] < 0) status = ZERO_MINUS;
742 }
743 switch (nint) {
744 case 0:
745 ie1 = i; s1 = status; nint++; break;
746 case 1:
747 ie2 = i; s2 = status; nint++; break;
748 default:
749 return NON_PLANAR_FACE;
750 }
751 }
752 if (nint != 2) return NON_PLANAR_FACE;
753
754 // F O R M I N T E R S E C T I O N S E G M E N T
755
756 if (s1 != ZERO_ZERO && s2 != ZERO_ZERO) {
757 if (s1 == s2) return NON_PLANAR_FACE;
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]];
763 while (iedge > 0) {
764 i1 = m_edges[iedge].i1;
765 i2 = m_edges[iedge].i2;
766
767 d1 = plane.distance(m_nodes[i1].v);
768 d2 = plane.distance(m_nodes[i2].v);
769 if (d1 > m_del) {
770 if (d2 < -m_del) { ii[i] = m_nodes.size(); break; } // +-
771 }else if (d1 < -m_del) {
772 if (d2 > m_del) { ii[i] = m_nodes.size(); break; } // -+
773 }else{
774 ii[i] = i1; break; // 0+ or 0-
775 }
776 iedge = m_edges[iedge].inext;
777 }
778 if (ii[i] == (int)m_nodes.size()) {
779 d3 = d2-d1;
780 if (d3 == 0){
781 std::cerr<<"d3 is zero in "<<__FILE__<<std::endl;
782 return DEFECTIVE_FACE;
783 }
784 d1 = d1/d3;
785 d2 = d2/d3;
786 m_nodes.push_back(ExtNode(d2*m_nodes[i1].v-d1*m_nodes[i2].v, iedge));
787 }
788 }
789 edge.inext = 0;
790 if (s1 == MINUS_PLUS || s1 == ZERO_PLUS) {
791 edge.i1 = ii[1];
792 edge.i2 = ii[0];
793 }else{
794 edge.i1 = ii[0];
795 edge.i2 = ii[1];
796 }
797 return INTERSECTION;
798 }else{
799 if (npos == nneg) return NON_PLANAR_FACE;
800 edge.inext = (s1 == ZERO_ZERO) ? ie1+1 : ie2+1;
801 if (s1 == ZERO_PLUS || s2 == ZERO_MINUS) {
802 edge.i1 = m_edges[m_faces[iface].iedges[ie2]].i1;
803 edge.i2 = m_edges[m_faces[iface].iedges[ie1]].i1;
804 }else{
805 edge.i1 = m_edges[m_faces[iface].iedges[ie1]].i1;
806 edge.i2 = m_edges[m_faces[iface].iedges[ie2]].i1;
807 }
808 return EDGE;
809 }
810}
811
812void BooleanProcessor::renumberNodes(int & i1, int & i2, int & i3, int & i4)
813/***********************************************************************
814 * *
815 * Name: BooleanProcessor::renumberNodes Date: 19.01.00 *
816 * Author: E.Chernyaev Revised: *
817 * *
818 * Function: Renumber nodes and remove last temporary node. *
819 * Remark: In principal this routine can be replaced just *
820 * with i1 = i2; *
821 * Removal of temporary nodes provides additional control *
822 * on number of nodes, that is very useful for debugging. *
823 * *
824 ***********************************************************************/
825{
826 if (i1 == i2) return;
827 if (m_nodes[i1].s == 0 || m_nodes.back().s == 0) { i1 = i2; return; }
828
829 int ilast = m_nodes.size()-1;
830 if (i1 == ilast) { i1 = i2; m_nodes.pop_back(); return; }
831 if (i2 == ilast) { i2 = i1; }
832 if (i3 == ilast) { i3 = i1; }
833 if (i4 == ilast) { i4 = i1; }
834 m_nodes[i1] = m_nodes.back(); i1 = i2; m_nodes.pop_back();
835}
836
838/***********************************************************************
839 * *
840 * Name: BooleanProcessor::testEdgeVsEdge Date: 19.01.00 *
841 * Author: E.Chernyaev Revised: *
842 * *
843 * Function: Find common part of two edges *
844 * *
845 ***********************************************************************/
846{
847 int i, ii = 0;
848 double d, dd = 0.;
849
850 for (i=0; i<3; i++) {
851 d = m_nodes[edge1.i1].v[i]-m_nodes[edge1.i2].v[i];
852 if (d < 0.) d = -d;
853 if (d > dd) { dd = d; ii = i; }
854 }
855 double t1 = m_nodes[edge1.i1].v[ii];
856 double t2 = m_nodes[edge1.i2].v[ii];
857 double t3 = m_nodes[edge2.i1].v[ii];
858 double t4 = m_nodes[edge2.i2].v[ii];
859 if (t2-t1 < 0.) { t1 = -t1; t2 = -t2; t3 = -t3; t4 = -t4; }
860
861 if (t3 <= t1+m_del || t4 >= t2-m_del) return 0;
862 if (t3 > t2+m_del) {
863 renumberNodes(edge2.i1, edge1.i2, edge1.i1, edge2.i2);
864 }else if (t3 < t2-m_del) {
865 renumberNodes(edge1.i2, edge2.i1, edge1.i1, edge2.i2);
866 }
867 if (t4 < t1-m_del) {
868 renumberNodes(edge2.i2, edge1.i1, edge1.i2, edge2.i1);
869 }else if (t4 > t1+m_del) {
870 renumberNodes(edge1.i1, edge2.i2, edge1.i2, edge2.i1);
871 }
872 return 1;
873}
874
875void BooleanProcessor::divideEdge(int & i1, int & i2)
876/***********************************************************************
877 * *
878 * Name: BooleanProcessor::divideEdge Date: 24.01.00 *
879 * Author: E.Chernyaev Revised: *
880 * *
881 * Function: Unify the nodes and divide edge on two parts by the node. *
882 * *
883 ***********************************************************************/
884{
885 int iedges[2];
886 iedges[0] = m_nodes[i1].s;
887 iedges[1] = m_nodes[i2].s;
888
889 // U N I F Y N O D E S
890
891 if (i1 < i2) { i2 = i1; }
892 else if (i1 > i2) { i1 = i2; }
893 else { iedges[1] = 0; }
894 if (iedges[0] == iedges[1]) return;
895
896 int ie1, ie2, inode = i1;
897 m_nodes[inode].s = 0;
898 for (int i=0; i<2; i++) {
899
900 // F I N D C O R R E S P O N D I N G E D G E
901
902 if ((ie1 = iedges[i]) == 0) continue;
903 ie2 = m_faces[m_edges[ie1].iface2].iedges[0];
904 while (ie2 > 0) {
905 if (m_edges[ie2].i1 == m_edges[ie1].i2 &&
906 m_edges[ie2].i2 == m_edges[ie1].i1) break;
907 ie2 = m_edges[ie2].inext;
908 }
909
910 // D I V I D E E D G E S
911
912 m_edges.push_back(m_edges[ie1]);
913 m_edges[ie1].inext = m_edges.size() - 1;
914 m_edges[ie1].i2 = inode;
915 m_edges.back().i1 = inode;
916
917 m_edges.push_back(m_edges[ie2]);
918 m_edges[ie2].inext = m_edges.size() - 1;
919 m_edges[ie2].i2 = inode;
920 m_edges.back().i1 = inode;
921 }
922}
923
925/***********************************************************************
926 * *
927 * Name: BooleanProcessor::insertEdge Date: 24.01.00 *
928 * Author: E.Chernyaev Revised: *
929 * *
930 * Function: Insert edge to the list of new edges *
931 * *
932 ***********************************************************************/
933{
934 int iface = edge.iface1;
935 m_edges.push_back(edge);
936 m_edges.back().inext = m_faces[iface].inew;
937 m_faces[iface].inew = m_edges.size() - 1;
938}
939
941/***********************************************************************
942 * *
943 * Name: BooleanProcessor::caseII Date: 19.01.00 *
944 * Author: E.Chernyaev Revised: *
945 * *
946 * Function: Intersection/Intersection case *
947 * *
948 ***********************************************************************/
949{
950 divideEdge(edge1.i1, edge2.i2);
951 divideEdge(edge1.i2, edge2.i1);
952 edge1.ivis = 1;
953 edge2.ivis = 1;
954 insertEdge(edge1);
955 insertEdge(edge2);
956}
957
959/***********************************************************************
960 * *
961 * Name: BooleanProcessor::caseIE Date: 19.01.00 *
962 * Author: E.Chernyaev Revised: *
963 * *
964 * Function: Intersection/Edge-touch case *
965 * *
966 ***********************************************************************/
967{
969#ifdef BP_DEBUG
970 G4cout
971 << "BooleanProcessor::caseIE : unimplemented case"
972 << G4endl;
973#endif
974}
975
977/***********************************************************************
978 * *
979 * Name: BooleanProcessor::caseEE Date: 19.01.00 *
980 * Author: E.Chernyaev Revised: *
981 * *
982 * Function: Edge-touch/Edge-touch case *
983 * *
984 ***********************************************************************/
985{
987#ifdef BP_DEBUG
988 G4cout
989 << "BooleanProcessor::caseEE : unimplemented case"
990 << G4endl;
991#endif
992}
993
994void BooleanProcessor::testFaceVsFace(int iface1, int iface2)
995/***********************************************************************
996 * *
997 * Name: BooleanProcessor::testFaceVsFace Date: 11.01.00 *
998 * Author: E.Chernyaev Revised: *
999 * *
1000 * Function: Find result (an edge) of intersection of two faces *
1001 * *
1002 ***********************************************************************/
1003{
1004 ExtEdge edge1, edge2;
1005 int irep1, irep2;
1006
1007 // M I N - M A X
1008
1009 {const ExtFace& face_1 = m_faces[iface1]; //G.Barrand : optimize
1010 const ExtFace& face_2 = m_faces[iface2];
1011 for (int i=0; i<3; i++) {
1012 if (face_1.rmin[i] > face_2.rmax[i] + m_del) return;
1013 if (face_1.rmax[i] < face_2.rmin[i] - m_del) return;
1014 }}
1015
1016 // F A C E - 1 vs P L A N E - 2
1017
1018 edge1.iface1 = iface1;
1019 edge1.iface2 = iface2;
1020 irep1 = testFaceVsPlane(edge1);
1021 if (irep1 == OUT_OF_PLANE || irep1 == ON_PLANE) {
1023 return;
1024 }
1025
1026 // F A C E - 2 vs P L A N E - 1
1027
1028 edge2.iface1 = iface2;
1029 edge2.iface2 = iface1;
1030 irep2 = testFaceVsPlane(edge2);
1031 if (irep2 == OUT_OF_PLANE || irep2 == ON_PLANE) {
1033 return;
1034 }
1035
1036 // C H E C K F O R N O N P L A N A R F A C E
1037
1038 if (irep1 == NON_PLANAR_FACE || irep2 == NON_PLANAR_FACE) {
1040 return;
1041 }
1042
1043 // F I N D I N T E R S E C T I O N P A R T
1044
1045 if (testEdgeVsEdge(edge1, edge2) == 0) return;
1046
1047 // C O N S I D E R D I F F E R E N T C A S E S
1048
1049 if (irep1 == INTERSECTION && irep2 == INTERSECTION) caseII(edge1, edge2);
1050 if (irep1 == INTERSECTION && irep2 == EDGE) caseIE(edge1, edge2);
1051 if (irep1 == EDGE && irep2 == INTERSECTION) caseIE(edge2, edge1);
1052 if (irep1 == EDGE && irep2 == EDGE) caseEE(edge1, edge2);
1054
1055}
1056
1058/***********************************************************************
1059 * *
1060 * Name: BooleanProcessor::invertNewEdges Date: 04.02.00 *
1061 * Author: E.Chernyaev Revised: *
1062 * *
1063 * Function: Invert direction of new edges *
1064 * *
1065 ***********************************************************************/
1066{
1067 int iedge = m_faces[iface].inew;
1068 while (iedge > 0) {
1069 m_edges[iedge].invert();
1070 iedge = m_edges[iedge].inext;
1071 }
1072}
1073
1075/***********************************************************************
1076 * *
1077 * Name: BooleanProcessor::checkDoubleEdges Date: 04.02.00 *
1078 * Author: E.Chernyaev Revised: *
1079 * *
1080 * Function: Eliminate duplication of edges *
1081 * *
1082 ***********************************************************************/
1083{
1084
1085}
1086
1087void BooleanProcessor::assembleFace(int what, int iface)
1088/***********************************************************************
1089 * *
1090 * Name: BooleanProcessor::assembleFace Date: 19.02.00 *
1091 * Author: E.Chernyaev Revised: *
1092 * *
1093 * Function: Assemble face *
1094 * *
1095 ***********************************************************************/
1096{
1097 // A S S E M B L E N E W F A C E
1098
1099 int ihead = 0; // head of the list of edges for new face
1100 int icur; // current edge in the list - last edge inserted to the list
1101 int *ilink; // pointer to the current link
1102 int ifirst; // first node of a contour
1103 int *i; // pointer to the index of the current edge in a loop
1104 int ioldflag=0; // is set if an edge from iold has been taken
1105
1106#define INSERT_EDGE_TO_THE_LIST(A) \
1107*ilink = A; ilink = &m_edges[A].inext; *ilink = 0
1108
1109 ExtFace& face = m_faces[iface]; //G.Barrand : optimize.
1110 ilink = &ihead;
1111 for(;;) {
1112 if (face.inew == 0) break;
1113
1114 // S T A R T N E W C O N T O U R
1115
1116 icur = face.inew;
1117 face.inew = m_edges[icur].inext;
1119 ifirst = m_edges[icur].i1;
1120
1121 // C O N S T R U C T T H E C O N T O U R
1122
1123 for (;;) {
1124 i = &face.inew;
1125 ExtEdge& edge_cur = m_edges[icur];
1126 while(*i > 0) {
1127 ExtEdge& edge_i = m_edges[*i];
1128 if (edge_i.i1 == edge_cur.i2) break;
1129 i = &edge_i.inext;
1130 }
1131 if (*i == 0) {
1132 i = &face.iold;
1133 while(*i > 0) {
1134 ExtEdge& edge_i = m_edges[*i];
1135 if (edge_i.i1 == edge_cur.i2) {
1136 ioldflag = 1;
1137 break;
1138 }
1139 i = &edge_i.inext;
1140 }
1141 }
1142 if (*i > 0) {
1143 icur = *i;
1144 *i = m_edges[icur].inext;
1146 if (m_edges[icur].i2 == ifirst) { break; } else { continue; }
1147 }else{
1149#ifdef BP_DEBUG
1150 G4cerr
1151 << "BooleanProcessor::assembleFace(" << iface << ") : "
1152 << "could not find next edge of the contour"
1153 << G4endl;
1154#endif
1155 face.inew = DEFECTIVE_FACE;
1156 return;
1157 }
1158 }
1159 }
1160
1161 // C H E C K O R I G I N A L C O N T O U R
1162
1163 int iedge;
1164 iedge = face.iold;
1165 if (what == 0 && ioldflag == 0 && iedge > 0) {
1166 for (;;) {
1167 if (m_edges[iedge].inext > 0) {
1168 if (m_edges[iedge].i2 == m_edges[m_edges[iedge].inext].i1) {
1169 iedge = m_edges[iedge].inext;
1170 }else{
1171 break;
1172 }
1173 }else{
1174 if (m_edges[iedge].i2 == m_edges[face.iold].i1) {
1175 m_edges[iedge].inext = ihead; // set new face
1176 return;
1177 }else{
1178 break;
1179 }
1180 }
1181 }
1182 }
1183
1184 // M A R K U N S U I T A B L E N E I G H B O U R I N G F A C E S
1185
1186 int iface2;
1187 iedge = face.iold;
1188 while(iedge > 0) {
1189 iface2 = m_edges[iedge].iface2;
1190 if (m_faces[iface2].inew == 0) m_faces[iface2].inew = UNSUITABLE_FACE;
1191 iedge = m_edges[iedge].inext;
1192 }
1193 face.iold = ihead; // set new face
1194}
1195
1196void BooleanProcessor::assembleNewFaces(int what, int ihead)
1197/***********************************************************************
1198 * *
1199 * Name: BooleanProcessor::assembleNewFaces Date: 30.01.00 *
1200 * Author: E.Chernyaev Revised: *
1201 * *
1202 * Function: Assemble internal or external parts of faces *
1203 * *
1204 ***********************************************************************/
1205{
1206 int iface = ihead;
1207 while(iface > 0) {
1208 if (m_faces[iface].inew > 0) {
1209 if (what != 0) invertNewEdges(iface);
1210 checkDoubleEdges(iface);
1211 assembleFace(what, iface);
1212 m_faces[iface].inew =
1213 (m_faces[iface].iold == 0) ? UNSUITABLE_FACE : NEW_FACE;
1214 }
1215 iface = m_faces[iface].inext;
1216 }
1217}
1218
1220/***********************************************************************
1221 * *
1222 * Name: BooleanProcessor::initiateLists Date: 28.02.00 *
1223 * Author: E.Chernyaev Revised: *
1224 * *
1225 * Function: Initiate lists of faces. *
1226 * *
1227 ***********************************************************************/
1228{
1229 int i, iface;
1230
1231 // R E S E T L I S T S O F F A C E S
1232
1233 m_result_faces.clean();
1234 m_suitable_faces.clean();
1235 m_unsuitable_faces.clean();
1236 m_unknown_faces.clean();
1237
1238 // I N I T I A T E T H E L I S T S
1239
1240 iface = m_iout1;
1241 while (iface > 0) {
1242 i = iface;
1243 iface = m_faces[i].inext;
1245 m_unsuitable_faces.push_back(i);
1246 m_faces[i].inew = UNSUITABLE_FACE;
1247 }else{
1248 m_suitable_faces.push_back(i);
1249 m_faces[i].inew = ORIGINAL_FACE;
1250 }
1251 }
1252 iface = m_iout2;
1253 while (iface > 0) {
1254 i = iface;
1255 iface = m_faces[i].inext;
1256 if (m_operation == OP_UNION) {
1257 m_suitable_faces.push_back(i);
1258 m_faces[i].inew = ORIGINAL_FACE;
1259 }else{
1260 m_unsuitable_faces.push_back(i);
1261 m_faces[i].inew = UNSUITABLE_FACE;
1262 }
1263 }
1264
1265 iface = m_iunk1;
1266 while (iface > 0) {
1267 i = iface;
1268 iface = m_faces[i].inext;
1269 m_unknown_faces.push_back(i);
1270 }
1271 iface = m_iunk2;
1272 while (iface > 0) {
1273 i = iface;
1274 iface = m_faces[i].inext;
1275 if (m_operation == OP_SUBTRACTION) m_faces[i].invert();
1276 m_unknown_faces.push_back(i);
1277 }
1278
1279 iface = m_ifaces1;
1280 while (iface > 0) {
1281 i = iface;
1282 iface = m_faces[i].inext;
1283 switch(m_faces[i].inew) {
1284 case UNKNOWN_FACE:
1285 m_unknown_faces.push_back(i);
1286 break;
1287 case ORIGINAL_FACE: case NEW_FACE:
1288 m_suitable_faces.push_back(i);
1289 break;
1290 case UNSUITABLE_FACE:
1291 m_unsuitable_faces.push_back(i);
1292 break;
1293 default:
1294 m_faces[i].iprev = 0;
1295 m_faces[i].inext = 0;
1296 break;
1297 }
1298 }
1299 iface = m_ifaces2;
1300 while (iface > 0) {
1301 i = iface;
1302 iface = m_faces[i].inext;
1303 if (m_operation == OP_SUBTRACTION) m_faces[i].invert();
1304 switch(m_faces[i].inew) {
1305 case UNKNOWN_FACE:
1306 m_unknown_faces.push_back(i);
1307 break;
1308 case ORIGINAL_FACE: case NEW_FACE:
1309 m_suitable_faces.push_back(i);
1310 break;
1311 case UNSUITABLE_FACE:
1312 m_unsuitable_faces.push_back(i);
1313 break;
1314 default:
1315 m_faces[i].iprev = 0;
1316 m_faces[i].inext = 0;
1317 break;
1318 }
1319 }
1321}
1322
1324/***********************************************************************
1325 * *
1326 * Name: BooleanProcessor::assemblePolyhedra() Date: 10.12.99 *
1327 * Author: E.Chernyaev Revised: *
1328 * *
1329 * Function: Collect suitable faces and remove unsuitable ones. *
1330 * *
1331 ***********************************************************************/
1332{
1333 int i, iedge, iface;
1334
1335 // L O O P A L O N G S U I T A B L E F A C E S
1336
1337 iface = m_suitable_faces.front();
1338 while(iface > 0) {
1339 i = iface;
1340 iedge = m_faces[i].iold;
1341 while(iedge > 0) {
1342 iface = m_edges[iedge].iface2;
1343 if (m_faces[iface].inew == UNKNOWN_FACE) {
1344 m_unknown_faces.remove(iface);
1345 m_suitable_faces.push_back(iface);
1346 m_faces[iface].inew = ORIGINAL_FACE;
1347 }
1348 iedge = m_edges[iedge].inext;
1349 }
1350 iface = m_faces[i].inext;
1351 m_suitable_faces.remove(i);
1352 m_result_faces.push_back(i);
1353 }
1354 if (m_unknown_faces.front() == 0) return;
1355
1356 // L O O P A L O N G U N S U I T A B L E F A C E S
1357
1358 iface = m_unsuitable_faces.front();
1359 while(iface > 0) {
1360 i = iface;
1361 iedge = m_faces[i].iold;
1362 while(iedge > 0) {
1363 iface = m_edges[iedge].iface2;
1364 if (m_faces[iface].inew == UNKNOWN_FACE) {
1365 m_unknown_faces.remove(iface);
1366 m_unsuitable_faces.push_back(iface);
1367 m_faces[iface].inew = UNSUITABLE_FACE;
1368 }
1369 iedge = m_edges[iedge].inext;
1370 }
1371 iface = m_faces[i].inext;
1372 m_unsuitable_faces.remove(i);
1373 }
1374
1375 //G.Barrand : begin
1376 /* From S.Ponce
1377 At last, there is a problem in the assemblePolyhedra method. At least, I
1378 think it is there. The problem deals with boolean operations on solids,
1379 when one of the two contains entirely the other one. It has no sense for
1380 intersection and union but still has sense for subtraction. In this
1381 case, faces from the inner solid are stored in the m_unknown_faces
1382 FaceList. And an error occurs in the execute method. This may be correct
1383 for intersection and union but in the case of subtraction, one should do
1384 that in assemblePolyhedra :
1385 */
1386 // Unknown faces are actually suitable face !!!
1387 iface = m_unknown_faces.front();
1388 while(iface > 0) {
1389 i = iface;
1390 m_faces[i].inew = ORIGINAL_FACE;
1391 iface = m_faces[i].inext;
1392 m_unknown_faces.remove(i);
1393 m_result_faces.push_back(i);
1394 }
1395 /*
1396 Otherwise, the inner hole that the second solid was building in the
1397 first one does not exist. I'm not very clear on what to do for unions
1398 and intersections. I think this kind of situation should be detected and
1399 one of the solid should simply be ignored.
1400 */
1401 //G.Barrand : end
1402}
1403
1404inline void
1405BooleanProcessor::findABC(double x1, double y1, double x2, double y2,
1406 double &a, double &b, double &c) const
1407/***********************************************************************
1408 * *
1409 * Name: BooleanProcessor::findABC Date: 07.03.00 *
1410 * Author: E.Chernyaev Revised: *
1411 * *
1412 * Function: Find line equation Ax+By+C=0 *
1413 * *
1414 ***********************************************************************/
1415{
1416 double w;
1417 a = y1 - y2;
1418 b = x2 - x1;
1419 //G.Barrand : w = std::abs(a)+std::abs(b);
1420 w = ::fabs(a)+::fabs(b); //G.Barrand
1421 a /= w;
1422 b /= w;
1423 c = -(a*x2 + b*y2);
1424}
1425
1426int BooleanProcessor::checkDirection(double *x, double *y) const
1427/***********************************************************************
1428 * *
1429 * Name: BooleanProcessor::checkDirection Date: 06.03.00 *
1430 * Author: E.Chernyaev Revised: *
1431 * *
1432 * Function: Check direction of line 1-4 *
1433 * *
1434 ***********************************************************************/
1435{
1436 double a1, b1, c1, a2, b2, c2, d1, d2;
1437
1438 // T E S T L I N E 1 - 4 V S E X T E R N A L C O N T O U R
1439
1440 findABC(x[0], y[0], x[1], y[1], a1, b1, c1);
1441 findABC(x[1], y[1], x[2], y[2], a2, b2, c2);
1442 d1 = a1*x[4] + b1*y[4] + c1;
1443 d2 = a2*x[4] + b2*y[4] + c2;
1444 if (d1 <= m_del && d2 <= m_del) return 1;
1445 if (! (d1 > m_del && d2 > m_del)) {
1446 if ( a1*x[2] + b1*y[2] + c1 >= -m_del) return 1;
1447 }
1448
1449 // T E S T L I N E 1 - 4 V S I N T E R N A L C O N T O U R
1450
1451 findABC(x[3], y[3], x[4], y[4], a1, b1, c1);
1452 findABC(x[4], y[4], x[5], y[5], a2, b2, c2);
1453 d1 = a1*x[1] + b1*y[1] + c1;
1454 d2 = a2*x[1] + b2*y[1] + c2;
1455 if (d1 <= m_del && d2 <= m_del) return 1;
1456 if (!(d1 > m_del && d2 > m_del)) {
1457 if ( a1*x[5] + b1*y[5] + c1 >= -m_del) return 1;
1458 }
1459 return 0;
1460}
1461
1462int BooleanProcessor::checkIntersection(int ix, int iy, int i1, int i2) const
1463/***********************************************************************
1464 * *
1465 * Name: BooleanProcessor::checkDirection Date: 06.03.00 *
1466 * Author: E.Chernyaev Revised: *
1467 * *
1468 * Function: Check line i1-i2 on intersection with contours *
1469 * *
1470 ***********************************************************************/
1471{
1472 // F I N D L I N E E Q U A T I O N
1473
1474 double x1, y1, x2, y2, a1, b1, c1;
1475 x1 = m_nodes[i1].v[ix];
1476 y1 = m_nodes[i1].v[iy];
1477 x2 = m_nodes[i2].v[ix];
1478 y2 = m_nodes[i2].v[iy];
1479 findABC(x1, y1, x2, y2, a1, b1, c1);
1480
1481 // L O O P A L O N G E X T E R N A L C O N T O U R S
1482
1483 int icontour, iedge, k1, k2;
1484 double x3, y3, x4, y4, a2, b2, c2, d1, d2;
1485 for(icontour=0; icontour<(int)m_external_contours.size(); icontour++) {
1486 iedge = m_external_contours[icontour];
1487 while(iedge > 0) {
1488 k1 = m_edges[iedge].i1;
1489 k2 = m_edges[iedge].i2;
1490 iedge = m_edges[iedge].inext;
1491 if (k1 == i1 || k2 == i1) continue;
1492 if (k1 == i2 || k2 == i2) continue;
1493 x3 = m_nodes[k1].v[ix];
1494 y3 = m_nodes[k1].v[iy];
1495 x4 = m_nodes[k2].v[ix];
1496 y4 = m_nodes[k2].v[iy];
1497 d1 = a1*x3 + b1*y3 + c1;
1498 d2 = a1*x4 + b1*y4 + c1;
1499 if (d1 > m_del && d2 > m_del) continue;
1500 if (d1 < -m_del && d2 < -m_del) continue;
1501
1502 findABC(x3, y3, x4, y4, a2, b2, c2);
1503 d1 = a2*x1 + b2*y1 + c2;
1504 d2 = a2*x2 + b2*y2 + c2;
1505 if (d1 > m_del && d2 > m_del) continue;
1506 if (d1 < -m_del && d2 < -m_del) continue;
1507 return 1;
1508 }
1509 }
1510
1511 // L O O P A L O N G E X T E R N A L C O N T O U R S
1512
1513 for(icontour=0; icontour<(int)m_internal_contours.size(); icontour++) {
1514 iedge = m_internal_contours[icontour];
1515 while(iedge > 0) {
1516 k1 = m_edges[iedge].i1;
1517 k2 = m_edges[iedge].i2;
1518 iedge = m_edges[iedge].inext;
1519 if (k1 == i1 || k2 == i1) continue;
1520 if (k1 == i2 || k2 == i2) continue;
1521 x3 = m_nodes[k1].v[ix];
1522 y3 = m_nodes[k1].v[iy];
1523 x4 = m_nodes[k2].v[ix];
1524 y4 = m_nodes[k2].v[iy];
1525 d1 = a1*x3 + b1*y3 + c1;
1526 d2 = a1*x4 + b1*y4 + c1;
1527 if (d1 > m_del && d2 > m_del) continue;
1528 if (d1 < -m_del && d2 < -m_del) continue;
1529
1530 findABC(x3, y3, x4, y4, a2, b2, c2);
1531 d1 = a2*x1 + b2*y1 + c2;
1532 d2 = a2*x2 + b2*y2 + c2;
1533 if (d1 > m_del && d2 > m_del) continue;
1534 if (d1 < -m_del && d2 < -m_del) continue;
1535 return 1;
1536 }
1537 }
1538 return 0;
1539}
1540
1541void BooleanProcessor::mergeContours(int ix, int iy, int kext, int kint)
1542/***********************************************************************
1543 * *
1544 * Name: BooleanProcessor::mergeContours Date: 06.03.00 *
1545 * Author: E.Chernyaev Revised: *
1546 * *
1547 * Function: Attemp to merge internal contour with external one *
1548 * *
1549 ***********************************************************************/
1550{
1551 int i1ext, i2ext, i1int, i2int, i, k[6];
1552 double x[6], y[6];
1553
1554 // L O O P A L O N G E X T E R N A L C O N T O U R
1555
1556 i1ext = m_external_contours[kext];
1557 while (i1ext > 0) {
1558 i2ext = m_edges[i1ext].inext;
1559 if (i2ext == 0) i2ext = m_external_contours[kext];
1560 k[0] = m_edges[i1ext].i1;
1561 k[1] = m_edges[i1ext].i2;
1562 k[2] = m_edges[i2ext].i2;
1563 for (i=0; i<3; i++) {
1564 x[i] = m_nodes[k[i]].v[ix];
1565 y[i] = m_nodes[k[i]].v[iy];
1566 }
1567
1568 // L O O P A L O N G I N T E R N A L C O N T O U R
1569
1570 i1int = m_internal_contours[kint];
1571 while (i1int > 0) {
1572 i2int = m_edges[i1int].inext;
1573 if (i2int == 0) i2int = m_internal_contours[kint];
1574 k[3] = m_edges[i1int].i1;
1575 k[4] = m_edges[i1int].i2;
1576 k[5] = m_edges[i2int].i2;
1577 for (i=3; i<6; i++) {
1578 x[i] = m_nodes[k[i]].v[ix];
1579 y[i] = m_nodes[k[i]].v[iy];
1580 }
1581
1582 // T E S T L I N E K1 - K4
1583 // I F O K T H E N M E R G E C O N T O U R S
1584
1585 if (checkDirection(x, y) == 0) {
1586 if (checkIntersection(ix, iy, k[1], k[4]) == 0) {
1587 i = i1int;
1588 for(;;) {
1589 if (m_edges[i].inext == 0) {
1590 m_edges[i].inext = m_internal_contours[kint];
1591 m_internal_contours[kint] = 0;
1592 break;
1593 }else{
1594 i = m_edges[i].inext;
1595 }
1596 }
1597 i = m_edges[i1int].iface1;
1598 m_edges.push_back(ExtEdge(k[1], k[4], i, -(int(m_edges.size())+1), -1));
1599 m_edges.back().inext = i2int;
1600 m_edges.push_back(ExtEdge(k[4], k[1], i, -(int(m_edges.size())-1), -1));
1601 m_edges.back().inext = m_edges[i1ext].inext;
1602 m_edges[i1ext].inext = m_edges.size()-2;
1603 m_edges[i1int].inext = m_edges.size()-1;
1604 return;
1605 }
1606 }
1607 i1int = m_edges[i1int].inext;
1608 }
1609 i1ext = m_edges[i1ext].inext;
1610 }
1611}
1612
1613int BooleanProcessor::checkTriangle(int iedge1, int iedge2, int ix, int iy) const
1614/***********************************************************************
1615 * *
1616 * Name: BooleanProcessor::checkTriangle Date: 08.03.00 *
1617 * Author: E.Chernyaev Revised: *
1618 * *
1619 * Function: Check triangle for correctness *
1620 * *
1621 ***********************************************************************/
1622{
1623 int k[3];
1624 double x[3], y[3];
1625 double a1, b1, c1;
1626
1627 k[0] = m_edges[iedge1].i1;
1628 k[1] = m_edges[iedge1].i2;
1629 k[2] = m_edges[iedge2].i2;
1630 for (int i=0; i<3; i++) {
1631 x[i] = m_nodes[k[i]].v[ix];
1632 y[i] = m_nodes[k[i]].v[iy];
1633 }
1634
1635 // C H E C K P R I N C I P A L C O R R E C T N E S S
1636
1637 findABC(x[2], y[2], x[0], y[0], a1, b1, c1);
1638 if (a1*x[1]+b1*y[1]+c1 <= 0.1*m_del) return 1;
1639
1640 // C H E C K T H A T T H E R E I S N O P O I N T S I N S I D E
1641
1642 int inode, iedge;
1643 double a2, b2, c2, a3, b3, c3;
1644
1645 findABC(x[0], y[0], x[1], y[1], a2, b2, c2);
1646 findABC(x[1], y[1], x[2], y[2], a3, b3, c3);
1647 iedge = iedge2;
1648 for (;;) {
1649 iedge = m_edges[iedge].inext;
1650 if (m_edges[iedge].inext == iedge1) return 0;
1651 inode = m_edges[iedge].i2;
1652 if (inode == k[0]) continue;
1653 if (inode == k[1]) continue;
1654 if (inode == k[2]) continue;
1655 x[1] = m_nodes[inode].v[ix];
1656 y[1] = m_nodes[inode].v[iy];
1657 if (a1*x[1]+b1*y[1]+c1 < -0.1*m_del) continue;
1658 if (a2*x[1]+b2*y[1]+c2 < -0.1*m_del) continue;
1659 if (a3*x[1]+b3*y[1]+c3 < -0.1*m_del) continue;
1660 return 1;
1661 }
1662 return 0; // default return
1663}
1664
1665void BooleanProcessor::triangulateContour(int ix, int iy, int ihead)
1666/***********************************************************************
1667 * *
1668 * Name: BooleanProcessor::triangulateContour Date: 06.03.00 *
1669 * Author: E.Chernyaev Revised: *
1670 * *
1671 * Function: Triangulate external contour *
1672 * *
1673 ***********************************************************************/
1674{
1675
1676 //G4cout << "Next Countour" << G4endl;
1677 //int draw_flag = 0;
1678 //if (draw_flag) draw_contour(5, 3, ihead);
1679
1680 // C L O S E C O N T O U R
1681
1682 int ipnext = ihead, nnode = 1;
1683 for (;;) {
1684 if (m_edges[ipnext].inext > 0) {
1685 ipnext = m_edges[ipnext].inext;
1686 nnode++;
1687 }else{
1688 m_edges[ipnext].inext = ihead;
1689 break;
1690 }
1691 }
1692
1693 // L O O P A L O N G C O N T O U R
1694
1695 //G4cerr << "debug : contour : begin : =================" << G4endl;
1696 //dump();//debug
1697
1698 int iedge1, iedge2, iedge3, istart = 0;
1699 for (;;) {
1700 iedge1 = m_edges[ipnext].inext;
1701 iedge2 = m_edges[iedge1].inext;
1702/*
1703 G4cerr << "debug :"
1704 << " ipnext " << ipnext
1705 << " iedge1 " << iedge1
1706 << " iedge2 " << iedge2
1707 << " : istart " << istart
1708 << G4endl;
1709*/
1710 if (istart == 0) {
1711 istart = iedge1;
1712 if (nnode <= 3) {
1713 iedge3 = m_edges[iedge2].inext;
1714 m_edges[iedge1].iface1 = m_faces.size();
1715 m_edges[iedge2].iface1 = m_faces.size();
1716 m_edges[iedge3].iface1 = m_faces.size();
1717 m_edges[iedge3].inext = 0;
1718 m_faces.push_back(ExtFace(m_edges,0)); //G.Barrand : ok ?
1719 m_faces.back().iold = iedge1;
1720 m_faces.back().inew = ORIGINAL_FACE;
1721
1722 //if (draw_flag) draw_contour(4, 2, iedge1);
1723
1724 break;
1725 }
1726 }else if (istart == iedge1) {
1728#ifdef BP_DEBUG
1729 G4cerr
1730 << "BooleanProcessor::triangulateContour : "
1731 << "could not generate a triangle (infinite loop)"
1732 << G4endl;
1733#endif
1734 break;
1735 }
1736
1737 // C H E C K C O R E C T N E S S O F T H E T R I A N G L E
1738
1739 if (checkTriangle(iedge1,iedge2,ix,iy) != 0) {
1740 ipnext = m_edges[ipnext].inext;
1741 continue;
1742 }
1743
1744 // M O D I F Y C O N T O U R
1745
1746 int i1 = m_edges[iedge1].i1;
1747 int i3 = m_edges[iedge2].i2;
1748 int iface1 = m_edges[iedge1].iface1;
1749 int iface2 = m_faces.size();
1750
1751 m_edges[ipnext].inext = m_edges.size();
1752 m_edges.push_back(ExtEdge(i1, i3, iface1, -(int(m_edges.size())+1), -1));
1753 m_edges.back().inext = m_edges[iedge2].inext;
1754
1755 // A D D N E W T R I A N G L E T O T H E L I S T
1756
1757 m_edges[iedge2].inext = m_edges.size();
1758 m_edges.push_back(ExtEdge(i3, i1, iface2, -(int(m_edges.size())-1), -1));
1759 m_faces.push_back(ExtFace(m_edges,0)); //G.Barrand : ok ?
1760 m_faces.back().iold = iedge1;
1761 m_faces.back().inew = ORIGINAL_FACE;
1762 m_edges[iedge1].iface1 = iface2;
1763 m_edges[iedge2].iface1 = iface2;
1764 ipnext = m_edges[ipnext].inext;
1765 istart = 0;
1766 nnode--;
1767
1768 //if (draw_flag) draw_contour(4, 2, iedge1);
1769
1770 }
1771}
1772
1773void BooleanProcessor::modifyReference(int iface, int i1, int i2, int iref)
1774/***********************************************************************
1775 * *
1776 * Name: BooleanProcessor::modifyReference Date: 13.03.00 *
1777 * Author: E.Chernyaev Revised: *
1778 * *
1779 * Function: Modify reference to the neighbouring face *
1780 * *
1781 ***********************************************************************/
1782{
1783 int iedge = m_faces[iface].iold;
1784 while (iedge > 0) {
1785 if (m_edges[iedge].i1 == i2 && m_edges[iedge].i2 == i1) {
1786 m_edges[iedge].iface2 = iref;
1787 return;
1788 }
1789 iedge = m_edges[iedge].inext;
1790 }
1792#ifdef BP_DEBUG
1793 G4cerr
1794 << "BooleanProcessor::modifyReference : could not find the edge, "
1795 << "iface=" << iface << ", i1,i2=" << i1 << "," << i2 << ", iref=" << iref
1796 << G4endl;
1797#endif
1798}
1799
1801/***********************************************************************
1802 * *
1803 * Name: BooleanProcessor::triangulateFace Date: 02.03.00 *
1804 * Author: E.Chernyaev Revised: *
1805 * *
1806 * Function: Triangulation of an extended face *
1807 * *
1808 ***********************************************************************/
1809{
1810
1811 // F I N D M A X C O M P O N E N T O F T H E N O R M A L
1812 // S E T IX, IY, IZ
1813
1814#ifdef BP_GEANT4 //G.Barrand
1815 HVNormal3D normal = m_faces[iface].plane.normal();
1816#else
1817 const HVNormal3D& normal = m_faces[iface].plane.getNormal();
1818#endif
1819 int ix, iy, iz = 0;
1820 //G.Barrand : if (std::abs(normal[1]) > std::abs(normal[iz])) iz = 1;
1821 //G.Barrand : if (std::abs(normal[2]) > std::abs(normal[iz])) iz = 2;
1822 if (::fabs(normal[1]) > ::fabs(normal[iz])) iz = 1; //G.Barrand
1823 if (::fabs(normal[2]) > ::fabs(normal[iz])) iz = 2; //G.Barrand
1824 if (normal[iz] > 0) {
1825 ix = (iz+1)%3; iy = (ix+1)%3;
1826 }else{
1827 iy = (iz+1)%3; ix = (iy+1)%3;
1828 }
1829
1830 // F I L L L I S T S O F C O N T O U R S
1831
1832 m_external_contours.clear();
1833 m_internal_contours.clear();
1834 double z;
1835 int i1, i2, ifirst, iedge, icontour = m_faces[iface].iold;
1836 while (icontour > 0) {
1837 iedge = icontour;
1838 ifirst = m_edges[iedge].i1;
1839 z = 0.0;
1840 for(;;) {
1841 if (iedge > 0) {
1842 i1 = m_edges[iedge].i1;
1843 i2 = m_edges[iedge].i2;
1844 ExtNode& node_1 = m_nodes[i1];
1845 ExtNode& node_2 = m_nodes[i2];
1846 z += node_1.v[ix]*node_2.v[iy]-node_2.v[ix]*node_1.v[iy];
1847 if (ifirst != i2) {
1848 iedge = m_edges[iedge].inext;
1849 continue;
1850 }else{
1851 if (z > m_del*m_del) {
1852 m_external_contours.push_back(icontour);
1853 }else if (z < -m_del*m_del) {
1854 m_internal_contours.push_back(icontour);
1855 }else{
1857#ifdef BP_DEBUG
1858 G4cerr
1859 << "BooleanProcessor::triangulateFace : too small contour"
1860 << G4endl;
1861#endif
1862 }
1863 icontour = m_edges[iedge].inext;
1864 m_edges[iedge].inext = 0;
1865 break;
1866 }
1867 }else{
1869#ifdef BP_DEBUG
1870 G4cerr
1871 << "BooleanProcessor::triangulateFace : broken contour"
1872 << G4endl;
1873#endif
1874 icontour = 0;
1875 break;
1876 }
1877 }
1878 }
1879
1880 // G E T R I D O F I N T E R N A L C O N T O U R S
1881
1882 int kint, kext;
1883 for (kint=0; kint < (int)m_internal_contours.size(); kint++) {
1884 for (kext=0; kext < (int)m_external_contours.size(); kext++) {
1885 mergeContours(ix, iy, kext, kint);
1886 if (m_internal_contours[kint] == 0) break;
1887 }
1888 if (kext == (int)m_external_contours.size()) {
1890#ifdef BP_DEBUG
1891 G4cerr
1892 << "BooleanProcessor::triangulateFace : "
1893 << "could not merge internal contour " << kint
1894 << G4endl;
1895#endif
1896 }
1897 }
1898
1899 // T R I A N G U L A T E C O N T O U R S
1900
1901 int nface = m_faces.size();
1902 for (kext=0; kext < (int)m_external_contours.size(); kext++) {
1904#ifdef BP_DEBUG
1905 if(m_processor_error) { //G.Barrand
1906 G4cerr
1907 << "BooleanProcessor::triangulateFace : "
1908 << "triangulateContour failed."
1909 << G4endl;
1910 break; //G.Barrand : ok ?
1911 }
1912#endif
1913 }
1914 m_faces[iface].inew = UNSUITABLE_FACE;
1915
1916 // M O D I F Y R E F E R E N C E S
1917
1918 for (int ifa=nface; ifa<(int)m_faces.size(); ifa++) {
1919 iedge = m_faces[ifa].iold;
1920 while (iedge > 0) {
1921 if (m_edges[iedge].iface1 != ifa) {
1923#ifdef BP_DEBUG
1924 G4cerr
1925 << "BooleanProcessor::triangulateFace : wrong reference to itself, "
1926 << "iface=" << ifa << ", iface1=" << m_edges[iedge].iface1
1927 << G4endl;
1928#endif
1929 }else if (m_edges[iedge].iface2 > 0) {
1930 modifyReference(m_edges[iedge].iface2,
1931 m_edges[iedge].i1, m_edges[iedge].i2, ifa);
1932 }else if (m_edges[iedge].iface2 < 0) {
1933 m_edges[iedge].iface2 = m_edges[-m_edges[iedge].iface2].iface1;
1934 }
1935 iedge = m_edges[iedge].inext;
1936 }
1937 }
1938}
1939
1941/***********************************************************************
1942 * *
1943 * Name: BooleanProcessor::createPolyhedron() Date: 14.03.00 *
1944 * Author: E.Chernyaev Revised: *
1945 * *
1946 * Function: Create HepPolyhedron. *
1947 * *
1948 ***********************************************************************/
1949{
1950 int i, iedge, nnode = 0, nface = 0;
1951
1952 // R E N U M E R A T E N O D E S A N D F A C E S
1953
1954 for (i=1; i<(int)m_nodes.size(); i++) m_nodes[i].s = 0;
1955
1956 for (i=1; i<(int)m_faces.size(); i++) {
1957 if (m_faces[i].inew == ORIGINAL_FACE) {
1958 m_faces[i].inew = ++nface;
1959 iedge = m_faces[i].iold;
1960 while (iedge > 0) {
1961 m_nodes[m_edges[iedge].i1].s = 1;
1962 iedge = m_edges[iedge].inext;
1963 }
1964 }else{
1965 m_faces[i].inew = 0;
1966 }
1967 }
1968
1969 for (i=1; i<(int)m_nodes.size(); i++) {
1970 if (m_nodes[i].s == 1) m_nodes[i].s = ++nnode;
1971 }
1972
1973 // A L L O C A T E M E M O R Y
1974
1975 ExtPolyhedron polyhedron;
1976 if (nface == 0) return polyhedron;
1977 polyhedron.AllocateMemory(nnode, nface);
1978
1979 // S E T N O D E S
1980
1981 for (i=1; i<(int)m_nodes.size(); i++) {
1982 if (m_nodes[i].s != 0) polyhedron.m_pV[m_nodes[i].s] = m_nodes[i].v;
1983 }
1984
1985 // S E T F A C E S
1986
1987 int k, v[4]={0}, f[4]={0};
1988 for (i=1; i<(int)m_faces.size(); i++) {
1989 if (m_faces[i].inew == 0) continue;
1990 v[3] = f[3] = k = 0;
1991 iedge = m_faces[i].iold;
1992 while (iedge > 0) {
1993 if (k > 3) {
1994 std::cerr << "BooleanProcessor::createPolyhedron : too many edges" << std::endl;
1995 break;
1996 }
1997 v[k] = m_nodes[m_edges[iedge].i1].s;
1998 if (m_edges[iedge].ivis < 0) v[k] = -v[k];
1999 f[k] = m_faces[m_edges[iedge].iface2].inew;
2000 iedge = m_edges[iedge].inext;
2001 k++;
2002 }
2003 if (k < 3) {
2004 std::cerr << "BooleanProcessor::createPolyhedron : "
2005 << "face has only " << k << " edges"
2006 << std::endl;
2007 }
2008 polyhedron.m_pF[m_faces[i].inew] =
2009 G4Facet(v[0],f[0], v[1],f[1], v[2],f[2], v[3],f[3]);
2010 }
2011 //coverity[copy_constructor_call]
2012 return polyhedron;
2013}
2014
2015int BooleanProcessor::s_ishift = 0; //G.Barrand
2016int BooleanProcessor::get_shift() { return s_ishift;} //G.Barrand
2017void BooleanProcessor::set_shift(int a_shift) { s_ishift = a_shift;} //G.Barrand
2018#define NUM_SHIFT 8
2019int BooleanProcessor::get_num_shift() { return NUM_SHIFT;} //G.Barrand
2020
2022 const HepPolyhedron & a,
2023 const HepPolyhedron & b,
2024 int& err) //G.Barrand
2025/***********************************************************************
2026 * *
2027 * Name: BooleanProcessor::execute Date: 10.12.99 *
2028 * Author: E.Chernyaev Revised: *
2029 * *
2030 * Function: Execute boolean operation. *
2031 * *
2032 ***********************************************************************/
2033{
2034 //static int ishift = 0; //G.Barrand
2035 //static double shift[8][3] = {
2036 static double shift[NUM_SHIFT][3] = { //G.Barrand
2037 { 31, 23, 17},
2038 { -31, -23, -17},
2039 { -23, 17, 31},
2040 { 23, -17, -31},
2041 { -17, -31, 23},
2042 { 17, 31, -23},
2043 { 31, -23, 17},
2044 { -31, 23, -17}
2045 };
2046
2047/*
2048 G4cerr << "BooleanProcessor::execute : ++++++++++++++++++++++"
2049 << a.getName().getString()
2050 << b.getName().getString()
2051 << G4endl;
2052*/
2053
2054 // I N I T I A T E P R O C E S S O R
2055
2057 m_operation = op;
2058 m_nodes.clear(); m_nodes.push_back(CRAZY_POINT);
2059 m_edges.clear(); m_edges.push_back(ExtEdge());
2060 m_faces.clear(); m_faces.push_back(ExtFace(m_edges,0)); //G.Barrand : ok ?
2061
2062 // T A K E P O L Y H E D R A
2063
2064 m_ifaces1 = m_faces.size(); takePolyhedron(a,0,0,0);
2065 m_ifaces2 = m_faces.size(); takePolyhedron(b,0,0,0);
2066
2067 if (m_processor_error) { // corrupted polyhedron
2068 std::cerr
2069 << "BooleanProcessor: corrupted input polyhedron"
2070 << std::endl;
2071 err = m_processor_error; //G.Barrand
2072 return HepPolyhedron();
2073 }
2074 if (m_ifaces1 == m_ifaces2) { // a is empty
2075 err = m_processor_error; //G.Barrand
2076 switch (m_operation) {
2077 case OP_UNION:
2078 return b;
2079 case OP_INTERSECTION:
2080 std::cerr
2081 << "BooleanProcessor: intersection with empty polyhedron"
2082 << std::endl;
2083 return HepPolyhedron();
2084 case OP_SUBTRACTION:
2085 std::cerr
2086 << "BooleanProcessor: subtraction from empty polyhedron"
2087 << std::endl;
2088 return HepPolyhedron();
2089 }
2090 }
2091 if (m_ifaces2 == (int)m_faces.size()) { // b is empty
2092 err = m_processor_error; //G.Barrand
2093 switch (m_operation) {
2094 case OP_UNION:
2095 return a;
2096 case OP_INTERSECTION:
2097 std::cerr
2098 << "BooleanProcessor: intersection with empty polyhedron"
2099 << std::endl;
2100 return HepPolyhedron();
2101 case OP_SUBTRACTION:
2102 return a;
2103 }
2104 }
2105
2106 // S E T I N I T I A L M I N - M A X A N D T O L E R A N C E
2107
2108 m_del = findMinMax();
2109
2110 // W O R K A R O U N D T O A V O I D I E A N D E E
2111
2112/*
2113#define PROCESSOR_ERROR(a_what) \
2114 G4cerr << "BooleanProcessor: boolean operation problem (" << a_what \
2115 << "). Try again with other shifts."\
2116 << G4endl;
2117*/
2118#define PROCESSOR_ERROR(a_what)
2119
2120 unsigned int try_count = 1;
2121 while(true) { //G.Barrand
2122
2123 double ddxx = m_del*shift[s_ishift][0];
2124 double ddyy = m_del*shift[s_ishift][1];
2125 double ddzz = m_del*shift[s_ishift][2];
2126 s_ishift++; if (s_ishift == get_num_shift()) s_ishift = 0;
2127
2128 m_processor_error = 0; //G.Barrand
2129 m_operation = op;
2130 m_nodes.clear(); m_nodes.push_back(CRAZY_POINT);
2131 m_edges.clear(); m_edges.push_back(ExtEdge());
2132 m_faces.clear(); m_faces.push_back(ExtFace(m_edges,0)); //G.Barrand : ok ?
2133
2134 m_ifaces1 = m_faces.size(); takePolyhedron(a,0,0,0);
2135 m_ifaces2 = m_faces.size(); takePolyhedron(b,ddxx,ddyy,ddzz);
2136
2137 if (m_processor_error) { PROCESSOR_ERROR(1) } //G.Barrand
2138
2139 m_del = findMinMax();
2140
2141 // P R E S E L E C T O U T S I D E F A C E S
2142
2143 m_iout1 = m_iout2 = 0;
2146
2147 if (m_processor_error) { PROCESSOR_ERROR(2) } //G.Barrand
2148
2149 // P R E S E L E C T N O I N T E R S E C T I O N F A C E S
2150
2151 int ifa1, ifa2;
2152 m_iunk1 = m_iunk2 = 0;
2153 if (m_iout1 != 0 || m_iout2 != 0) {
2154 for(;;) {
2155 ifa1 = m_iunk1;
2156 ifa2 = m_iunk2;
2159 if (m_iunk1 == ifa1 && m_iunk2 == ifa2) break;
2160 findMinMax();
2161 }
2162 }
2163
2164 if (m_processor_error) { PROCESSOR_ERROR(3) } //G.Barrand
2165
2166 // F I N D N E W E D G E S
2167
2168 if (m_ifaces1 != 0 && m_ifaces2 != 0 ) {
2169 ifa1 = m_ifaces1;
2170 while (ifa1 > 0) {
2171 ifa2 = m_ifaces2;
2172 while (ifa2 > 0) {
2173 testFaceVsFace(ifa1, ifa2);
2174 ifa2 = m_faces[ifa2].inext;
2175 }
2176 ifa1 = m_faces[ifa1].inext;
2177 }
2178 }
2179 if (m_processor_error) { PROCESSOR_ERROR(4) } //G.Barrand
2180
2181 // C O N S T R U C T N E W F A C E S
2182
2184 if (m_processor_error) { PROCESSOR_ERROR(5) } //G.Barrand
2186 if (m_processor_error) { PROCESSOR_ERROR(6) } //G.Barrand
2187
2188 // A S S E M B L E S U I T A B L E F A C E S
2189
2190 initiateLists();
2192 if (m_unknown_faces.front() != 0) {
2194#ifdef BP_DEBUG
2195 G4cerr
2196 << "BooleanProcessor::execute : unknown faces !!!"
2197 << G4endl;
2198#endif
2199 }
2200 if (m_processor_error) { PROCESSOR_ERROR(7) } //G.Barrand
2201
2202 // T R I A N G U L A T E A C C E P T E D F A C E S
2203
2204 ifa1 = m_result_faces.front();
2205 while (ifa1 > 0) {
2206 ifa2 = ifa1;
2207 ifa1 = m_faces[ifa2].inext;
2208 if (m_faces[ifa2].inew == NEW_FACE) triangulateFace(ifa2);
2209 if (m_processor_error) {
2210 PROCESSOR_ERROR(8) //G.Barrand
2211 break; //G.Barrand
2212 }
2213 }
2214
2215 if(!m_processor_error) {
2216#ifdef BP_DEBUG
2217 if(try_count!=1) {
2218 G4cerr
2219 << "BooleanProcessor::execute : had converged."
2220 << G4endl;
2221 }
2222#endif
2223 break;
2224 }
2225
2226 if((int)try_count>get_num_shift()) {
2227#ifdef BP_DEBUG
2228 /*** Commented out because HepPolyhedron does not have getName...?!
2229 G4cerr << "BooleanProcessor: "
2230 << " all shifts tried. Boolean operation (" << op << ") failure."
2231 << " a name \"" << a.getName().getString() << "\""
2232 << " b name \"" << b.getName().getString() << "\""
2233 << G4endl;
2234 ***/
2235#endif
2236 err = m_processor_error;
2237 return a;
2238 }
2239
2240#ifdef BP_DEBUG
2241 G4cerr
2242 << "BooleanProcessor::execute : try another tilt..."
2243 << G4endl;
2244#endif
2245
2246 try_count++;
2247
2248 } //G.Barrand : end while shift.
2249#undef PROCESSOR_ERROR //G.Barrand
2250
2251 // C R E A T E P O L Y H E D R O N
2252
2253 err = m_processor_error;
2254 return createPolyhedron();
2255}
2256
2257
2258//#include <cfortran.h>
2259//#include <higz.h>
2260//#include "zbuf.h"
2261//void BooleanProcessor::draw()
2262/***********************************************************************
2263 * *
2264 * Name: BooleanProcessor::draw Date: 10.12.99 *
2265 * Author: E.Chernyaev Revised: *
2266 * *
2267 * Function: Draw *
2268 * *
2269 ***********************************************************************/
2270/*
2271{
2272 int II;
2273 int icol, i1, i2, iedge, iface, ilist[4];
2274 float p1[3], p2[3];
2275
2276 ilist[0] = m_ifaces1;
2277 ilist[1] = m_ifaces2;
2278 ilist[2] = m_iout1;
2279 ilist[3] = m_iout2;
2280
2281 for (int i=0; i<4; i++) {
2282
2283 if (i == 0) G4cout << "========= Ifaces_1" << G4endl;
2284 if (i == 1) G4cout << "========= Ifaces_2" << G4endl;
2285 if (i == 2) G4cout << "========= Iout_1" << G4endl;
2286 if (i == 3) G4cout << "========= Iout_2" << G4endl;
2287
2288 icol = i+1;
2289 iface = ilist[i];
2290 while (iface > 0) {
2291
2292 G4cout << "iface = " << iface << G4endl;
2293 G4cout << "--- iold" << G4endl;
2294
2295 iedge = m_faces[iface].iold;
2296 icol = 2;
2297
2298 while (iedge > 0) {
2299
2300 G4cout << " iegde = " << iedge
2301 << " i1,i2 =" << m_edges[iedge].i1 << "," << m_edges[iedge].i2
2302 << " iface1,iface2 = "
2303 << m_edges[iedge].iface1 << "," << m_edges[iedge].iface2
2304 << G4endl;
2305
2306 i1 = m_edges[iedge].i1;
2307 p1[0] = m_nodes[i1].v.x();
2308 p1[1] = m_nodes[i1].v.y();
2309 p1[2] = m_nodes[i1].v.z();
2310 IHWTON(p1,p1);
2311 i2 = m_edges[iedge].i2;
2312 p2[0] = m_nodes[i2].v.x();
2313 p2[1] = m_nodes[i2].v.y();
2314 p2[2] = m_nodes[i2].v.z();
2315 IHWTON(p2,p2);
2316// icol = (m_edges[iedge].ivis > 0) ? 1 : 2;
2317 IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
2318 iedge = m_edges[iedge].inext;
2319 }
2320
2321 G4cout << "--- inew" << G4endl;
2322
2323 iedge = m_faces[iface].inew;
2324 icol = 3;
2325
2326 while (iedge > 0) {
2327
2328 G4cout << " iegde = " << iedge
2329 << " i1,i2 =" << m_edges[iedge].i1 << "," << m_edges[iedge].i2
2330 << " iface1,iface2 = "
2331 << m_edges[iedge].iface1 << "," << m_edges[iedge].iface2
2332 << G4endl;
2333
2334 i1 = m_edges[iedge].i1;
2335 p1[0] = m_nodes[i1].v.x();
2336 p1[1] = m_nodes[i1].v.y();
2337 p1[2] = m_nodes[i1].v.z();
2338 IHWTON(p1,p1);
2339 i2 = m_edges[iedge].i2;
2340 p2[0] = m_nodes[i2].v.x();
2341 p2[1] = m_nodes[i2].v.y();
2342 p2[2] = m_nodes[i2].v.z();
2343 IHWTON(p2,p2);
2344// icol = (m_edges[iedge].ivis > 0) ? 1 : 2;
2345 IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
2346 iedge = m_edges[iedge].inext;
2347 }
2348 iface = m_faces[iface].inext;
2349
2350 IHZTOX(0,100,100);
2351 ixupdwi(0);
2352 cin >> II;
2353 ixclrwi();
2354 IHZCLE(0);
2355 }
2356 }
2357}
2358*/
2359
2360/*
2361//--------------------------------------------------------------------
2362void
2363BooleanProcessor::draw_edge(int icol, int iedge) {
2364 int i1, i2;
2365 float p1[3], p2[3];
2366
2367 i1 = m_edges[iedge].i1;
2368 p1[0] = m_nodes[i1].v.x();
2369 p1[1] = m_nodes[i1].v.y();
2370 p1[2] = m_nodes[i1].v.z();
2371 IHWTON(p1,p1);
2372 i2 = m_edges[iedge].i2;
2373 p2[0] = m_nodes[i2].v.x();
2374 p2[1] = m_nodes[i2].v.y();
2375 p2[2] = m_nodes[i2].v.z();
2376 IHWTON(p2,p2);
2377 IHZLIN(icol,p1[0],p1[1],p1[2], p2[0],p2[1],p2[2]);
2378}
2379
2380//--------------------------------------------------------------------
2381void
2382BooleanProcessor::draw_contour(int i1col, int i2col, int ihead) {
2383 int iedge, icol;
2384 iedge = ihead;
2385 while (iedge > 0) {
2386 icol = (m_edges[iedge].ivis > 0) ? i1col : i2col;
2387 draw_edge(icol, iedge);
2388 iedge = m_edges[iedge].inext;
2389 }
2390
2391 IHZTOX(0,100,100);
2392 ixupdwi(0);
2393
2394 int i;
2395 std::cin >> i;
2396}
2397
2398//--------------------------------------------------------------------
2399void
2400BooleanProcessor::print_face(int iface) {
2401 G4cout.precision(3);
2402 G4cout << "\n====== Face N " << iface << G4endl;
2403 G4cout << "iedges[4] = "
2404 << m_faces[iface].iedges[0] << ", "
2405 << m_faces[iface].iedges[1] << ", "
2406 << m_faces[iface].iedges[2] << ", "
2407 << m_faces[iface].iedges[3] << G4endl;
2408 G4cout << "rmin[3] = "
2409 << m_faces[iface].m_rmin[0] << ", "
2410 << m_faces[iface].m_rmin[1] << ", "
2411 << m_faces[iface].m_rmin[2] << G4endl;
2412 G4cout << "rmax[3] = "
2413 << m_faces[iface].m_rmax[0] << ", "
2414 << m_faces[iface].m_rmax[1] << ", "
2415 << m_faces[iface].m_rmax[2] << G4endl;
2416 G4cout << "iprev,inext = "
2417 << m_faces[iface].iprev << ", "
2418 << m_faces[iface].inext << G4endl;
2419 G4cout << "iold = " << m_faces[iface].iold << G4endl;
2420 for(int i = m_faces[iface].iold; i != 0;) {
2421 print_edge(i);
2422 i = m_edges[abs(i)].inext;
2423 }
2424
2425 G4cout << "inew = ";
2426 switch (m_faces[iface].inew) {
2427 case UNKNOWN_FACE:
2428 G4cout << "UNKNOWN_FACE" << G4endl;
2429 break;
2430 case ORIGINAL_FACE:
2431 G4cout << "ORIGINAL_FACE" << G4endl;
2432 break;
2433 case NEW_FACE:
2434 G4cout << "NEW_FACE" << G4endl;
2435 break;
2436 case UNSUITABLE_FACE:
2437 G4cout << "UNSUITABLE_FACE" << G4endl;
2438 break;
2439 case DEFECTIVE_FACE:
2440 G4cout << "DEFECTIVE_FACE" << G4endl;
2441 break;
2442 default:
2443 G4cout << m_faces[iface].inew << G4endl;
2444 for(int k = m_faces[iface].inew; k != 0;) {
2445 print_edge(k);
2446 k = m_edges[abs(k)].inext;
2447 }
2448 }
2449}
2450
2451//--------------------------------------------------------------------
2452void
2453BooleanProcessor::print_edge(int iedge) {
2454 G4cout << "==== Edge N " << iedge << G4endl;
2455 int i = std::abs(iedge);
2456 int i1 = m_edges[i].i1;
2457 int i2 = m_edges[i].i2;
2458 G4cout << "node[" << i1 << "] = "
2459 << m_nodes[i1].v.x() << ", "
2460 << m_nodes[i1].v.y() << ", "
2461 << m_nodes[i1].v.z() << G4endl;
2462
2463 G4cout << "node[" << i2 << "] = "
2464 << m_nodes[i2].v.x() << ", "
2465 << m_nodes[i2].v.y() << ", "
2466 << m_nodes[i2].v.z() << G4endl;
2467
2468 G4cout << "iface1,iface2,ivis,inext = "
2469 << m_edges[i].iface1 << ", "
2470 << m_edges[i].iface2 << ", "
2471 << m_edges[i].ivis << ", "
2472 << m_edges[i].inext << G4endl;
2473}
2474*/
2475
2476void BooleanProcessor::dump() {//G.Barrand
2477 unsigned int number = m_nodes.size();
2478 std::cerr << "nodes : " << number << std::endl;
2479 for(unsigned int index=0;index<number;index++) {
2480 const ExtNode& node = m_nodes[index];
2481 std::cerr << " " << index
2482 << " x = " << node.v[0]
2483 << " y = " << node.v[1]
2484 << " z = " << node.v[2]
2485 << std::endl;
2486 }
2487}
#define INSERT_EDGE_TO_THE_LIST(A)
#define FaceList
#define NON_PLANAR_FACE
#define DEFECTIVE_FACE
#define HepPolyhedron
#define G4Facet
#define OUT_OF_PLANE
#define ExtEdge
#define OP_INTERSECTION
#define OP_SUBTRACTION
#define OP_UNION
#define ExtFace
#define PROCESSOR_ERROR(a_what)
#define SWAP(A, B)
#define ON_PLANE
#define NUM_SHIFT
#define ExtNode
#define ORIGINAL_FACE
#define GRANULARITY
#define UNKNOWN_FACE
#define UNSUITABLE_FACE
#define EDGE
HEPVis::SbPlane HVPlane3D
#define ExtPolyhedron
#define INTERSECTION
#define CRAZY_POINT
#define NEW_FACE
static Double_t a
HVPoint3D HVNormal3D
#define y
#define x
#define z
std::vector< int > m_internal_contours
std::vector< int > m_external_contours
void caseII(ExtEdge &edge1, ExtEdge &edge2)
void insertEdge(const ExtEdge &edge)
void draw_contour(int, int, int)
std::vector< ExtNode > m_nodes
HepPolyhedron createPolyhedron()
void modifyReference(int iface, int i1, int i2, int iref)
void invertNewEdges(int iface)
std::vector< ExtFace > m_faces
void findABC(double x1, double y1, double x2, double y2, double &a, double &b, double &c) const
void draw_faces(int, int, int)
void triangulateFace(int iface)
void takePolyhedron(const HepPolyhedron &p, double, double, double)
void caseIE(ExtEdge &edge1, ExtEdge &edge2)
void assembleNewFaces(int what, int ihead)
int checkDirection(double *x, double *y) const
void selectOutsideFaces(int &ifaces, int &iout)
int testEdgeVsEdge(ExtEdge &edge1, ExtEdge &edge2)
void mergeContours(int ix, int iy, int kext, int kint)
void checkDoubleEdges(int iface)
HepPolyhedron execute(int op, const HepPolyhedron &a, const HepPolyhedron &b, int &err)
void caseEE(ExtEdge &edge1, ExtEdge &edge2)
void renumberNodes(int &i1, int &i2, int &i3, int &i4)
static void set_shift(int)
std::vector< ExtEdge > m_edges
void triangulateContour(int ix, int iy, int ihead)
void print_edge(int)
int checkTriangle(int iedge1, int iedge2, int ix, int iy) const
void print_face(int)
void assembleFace(int what, int iface)
int checkIntersection(int ix, int iy, int i1, int i2) const
int testFaceVsPlane(ExtEdge &edge)
void testFaceVsFace(int iface1, int iface2)
void divideEdge(int &i1, int &i2)
static int get_num_shift()
int get_processor_error() const
void draw_edge(int, int)
ExtEdge(const ExtEdge &edge)
ExtEdge(int k1=0, int k2=0, int kface1=0, int kface2=0, int kvis=0)
ExtEdge & operator=(const ExtEdge &edge)
HVPlane3D plane
double rmax[3]
ExtFace & operator=(const ExtFace &face)
ExtFace(std::vector< ExtEdge > &a_edges, int iedge)
ExtFace(const ExtFace &face)
double rmin[3]
std::vector< ExtEdge > & m_edges
ExtNode & operator=(const ExtNode &node)
ExtNode(const ExtNode &node)
HVPoint3D v
ExtNode(HVPoint3D vertex=HVPoint3D(), int status=0)
friend class BooleanProcessor
FaceList(std::vector< ExtFace > &a_faces)
void remove(int i)
void push_back(int i)
std::vector< ExtFace > & m_faces
double distance(const SbVec3d &point) const
Definition SbPlane.cxx:74
Definition node.h:24
-event-from-file
Definition index.py:1
std::string number(const double &d, const std::string &s)
Definition utils.cxx:186