1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
|
// slicer.cmd --- create intersection of surface with plane
// plane eq: a*x + b*y + c*z = d
// Programmer: Ken Brakke, brakke@susqu.edu, http://www.susqu.edu/brakke
// Commands included:
// drawslice - make new edges across facets along plane.
// slicer - make slice and dissolve all facets, edges, and vertices
// on negative side of the plane. Calls drawslice.
// slice_list - print coordinates of vertices along a slice.
// "slicer" usage:
// Set slice_a,slice_b,slice_c,slice_d coefficients so
// slice_a*x + slice_b*y + slice_c*z >= slice_d
// is the side you want, and do "slicer".
// output: truncated surface on positive side of plane
// Try not to slice exactly through vertices!!
// Default plane is slice_a := 0, slice_b := 0, slice_c := 1, slice_d := .1
// To mark the vertices and edges created by the current slicing,
// there are a vertex attribute v_timestamp and an edge attribute
// e_timestamp that are set to slice_timestamp, which is incremented
// each time slicer is called.
// Works in torus by rewrapping wrapped edges that would be cut
// so unwrapped part is on positive side of cut plane.
// Slice_list:
// Listing of vertices along the slice, beginning at user-designated vertex.
// If the slice has multiple components, you will have to run separately
// for each component.
// Usage: set startv to the first vertex on the slice and call slice_list
// Output is to stdout, so will usually be redirected, e.g.
// slice_list >>> "slice.dat"
//Set these coefficients for the slicing plane
slice_a := 0;
slice_b := 0;
slice_c := 1;
slice_d := .1;
// Attributes for marking vertices and edges on slice with timestamp
slice_timestamp := 1 // so won't conflict with 0 for new edges or 1 for
// original edges
define vertex attribute v_timestamp integer
define edge attribute e_timestamp integer
// Mark existing edges and vertices (if not already marked)
set edge e_timestamp 1 where e_timestamp==0
set vertex v_timestamp 1 where v_timestamp==0
// First put in new edges along slicing plane
drawslice := {
local lambda;
local denom;
local xx1;
local yy1;
local zz1;
local xx2;
local yy2;
local zz2;
local xb,yb,zb;
slice_timestamp += 1; // marker for this slice
foreach edge ee do
{
xx1 := ee.vertex[1].x;
yy1 := ee.vertex[1].y;
zz1 := ee.vertex[1].z;
xx2 := xx1 + ee.x; // using edge vector in case of torus wrap
yy2 := yy1 + ee.y;
zz2 := zz1 + ee.z;
denom := slice_a*(xx1-xx2)+slice_b*(yy1-yy2)+slice_c*(zz1-zz2);
if ( denom != 0.0 ) then
{
lambda := (slice_d-slice_a*xx2-slice_b*yy2-slice_c*zz2)/denom;
if ( (lambda >= 0) and (lambda <= 1) ) then
{
if torus then
if ee.wrap then
{ if denom > 0 then // tail on positive side
wrap_vertex(ee.vertex[2].id,ee.wrap)
else // head on positive side
wrap_vertex(ee.vertex[1].id,wrap_inverse(ee.wrap));
};
xb := lambda*xx1+(1-lambda)*xx2;
yb := lambda*yy1+(1-lambda)*yy2;
zb := lambda*zz1+(1-lambda)*zz2;
refine ee;
ee.vertex[2].v_timestamp := slice_timestamp;
ee.vertex[2].x := xb;
ee.vertex[2].y := yb;
ee.vertex[2].z := zb;
}
else if torus and ee.wrap then
{ // try wrapping from head
xx2 := ee.vertex[2].x;
yy2 := ee.vertex[2].y;
zz2 := ee.vertex[2].z;
xx1 := xx2 - ee.x; // using edge vector in case of torus wrap
yy1 := yy2 - ee.y;
zz1 := zz2 - ee.z;
denom := slice_a*(xx1-xx2)+slice_b*(yy1-yy2)+slice_c*(zz1-zz2);
if ( denom != 0.0 ) then
{
lambda := (slice_d-slice_a*xx2-slice_b*yy2-slice_c*zz2)/denom;
if ( (lambda >= 0) and (lambda <= 1) ) then
{
if torus then
if ee.wrap then
{ if denom > 0 then // tail on positive side
wrap_vertex(ee.vertex[2].id,ee.wrap)
else // head on positive side
wrap_vertex(ee.vertex[1].id,wrap_inverse(ee.wrap));
};
xb := lambda*xx1+(1-lambda)*xx2;
yb := lambda*yy1+(1-lambda)*yy2;
zb := lambda*zz1+(1-lambda)*zz2;
refine ee;
ee.vertex[2].v_timestamp := slice_timestamp;
ee.vertex[2].x := xb;
ee.vertex[2].y := yb;
ee.vertex[2].z := zb;
}
}
}
} ;
} ;
}
slicer := {
local former_autodisplay;
drawslice;
former_autodisplay := (autodisplay);
autodisplay off; // prevent display while dissolving
foreach facet ff where // again, careful of torus wraps
slice_a*(ff.vertex[1].x+ff.edge[1].x/3-ff.edge[3].x/3) +
slice_b*(ff.vertex[1].y+ff.edge[1].y/3-ff.edge[3].y/3) +
slice_c*(ff.vertex[1].z+ff.edge[1].z/3-ff.edge[3].z/3) < slice_d do
{ unset ff frontbody;
unset ff backbody;
dissolve ff;
};
dissolve bodies bbb where sum(bbb.facets,1) == 0;
dissolve edges ee where min(ee.vertex vv,
slice_a*vv.x + slice_b*vv.y + slice_c*vv.z) < slice_d;
// just does bare edges
dissolve vertices vv where
slice_a*vv.x + slice_b*vv.y + slice_c*vv.z < slice_d;
// just does bare vertices
// mark edges and vertices created by slicing
set edge ee e_timestamp slice_timestamp where
(ee.valence==1) and (ee.e_timestamp == 0);
set edge ee e_timestamp 1 where ee.e_timestamp == 0;
foreach edge ee where e_timestamp==slice_timestamp do
set ee.vertex v_timestamp slice_timestamp;
if former_autodisplay then autodisplay on;
}
color_slice := {
foreach facet ff where // again, careful of torus wraps
slice_a*(ff.vertex[1].x+ff.edge[1].x/3-ff.edge[3].x/3) +
slice_b*(ff.vertex[1].y+ff.edge[1].y/3-ff.edge[3].y/3) +
slice_c*(ff.vertex[1].z+ff.edge[1].z/3-ff.edge[3].z/3) < slice_d do
set ff color magenta;
}
// Listing of vertices along the slice, beginning at user-designated vertex.
// If the slice has multiple components, you will have to run separately
// for each component.
// Usage: set startv to the first vertex on the slice and call slice_list
// Output is to stdout, so will usually be redirected, e.g.
// slice_list >>> "slice.dat"
startv := 1
slice_list := {
local thisv,thisedge,nextedge;
if vertex[startv].v_timestamp != slice_timestamp then
{ errprintf "slice_list: starting vertex startv is %d, not on latest slice.\n",
startv;
return;
};
thisv := startv;
thisedge := 0;
do
{ // print current vertex coordinates
printf "%20.15f %20.15f %20.15f\n",vertex[thisv].x,vertex[thisv].y,
vertex[thisv].z;
nextedge := 0;
foreach vertex[thisv].edge ee where ee.e_timestamp == slice_timestamp
do { if ee.oid != -thisedge then { nextedge := ee.oid; break; }};
if nextedge == 0 then break; /* done */
thisv := edge[nextedge].vertex[2].id;
thisedge := nextedge;
} while thisv != startv; // done
}
// "slicer" usage:
// Set slice_a,slice_b,slice_c,slice_d coefficients so
// slice_a*x + slice_b*y + slice_c*z >= slice_d
// is the side you want, and do "slicer".
// Default plane is slice_a := 0, slice_b := 0, slice_c := 1, slice_d := .1
|
