1 | from visual import * |
---|
2 | |
---|
3 | def normal(v1,v2): |
---|
4 | """Safe computation of normalised cross product |
---|
5 | """ |
---|
6 | try: |
---|
7 | n = norm(cross(v1, v2)) |
---|
8 | except ZeroDivisionError: |
---|
9 | print v1, v2 |
---|
10 | raise Exception |
---|
11 | #n = vector(0,0,0) |
---|
12 | |
---|
13 | return n |
---|
14 | |
---|
15 | |
---|
16 | class Triangle: |
---|
17 | def __init__(self,v0,v1,v2,color=(1,1,1), frame=None): |
---|
18 | """Create one flat triangular panel with two sides |
---|
19 | """ |
---|
20 | |
---|
21 | #Outward normal |
---|
22 | |
---|
23 | self.normal=normal(v1-v0, v2-v0) |
---|
24 | |
---|
25 | self.panel = faces(frame = frame) |
---|
26 | self.backpanel = faces(frame = frame) |
---|
27 | for v in (v0,v1,v2): |
---|
28 | self.panel.append(pos=v, normal = self.normal, color=color) |
---|
29 | |
---|
30 | for v in (v2,v1,v0): |
---|
31 | self.backpanel.append(pos=v, normal = -self.normal, color=color) |
---|
32 | |
---|
33 | |
---|
34 | def move(self,v): |
---|
35 | """Move panel in direction given by vector v |
---|
36 | """ |
---|
37 | |
---|
38 | for vertex in self.panel.pos: |
---|
39 | vertex += v |
---|
40 | |
---|
41 | def update_height(self,d,fix_baseline=False): |
---|
42 | """Change height of triangle by displacing the third vertex by |
---|
43 | d in the direction perpendicular to the baseline (v1-v0) and the |
---|
44 | outward normal vector |
---|
45 | """ |
---|
46 | |
---|
47 | v0 = self.panel.pos[0] |
---|
48 | v1 = self.panel.pos[1] |
---|
49 | v2 = self.panel.pos[2] |
---|
50 | |
---|
51 | n = normal(v1-v0, self.normal) |
---|
52 | |
---|
53 | if fix_baseline: |
---|
54 | self.panel.pos[2] -= d*n |
---|
55 | else: |
---|
56 | self.panel.pos[:2] += d*n |
---|
57 | |
---|
58 | def update_color(self, c): |
---|
59 | """Change color to c |
---|
60 | """ |
---|
61 | |
---|
62 | self.panel.color = c |
---|
63 | |
---|
64 | |
---|
65 | def set_vertexheights(self, heights, floor_heights = None): |
---|
66 | from Numeric import zeros, Float |
---|
67 | |
---|
68 | from config import minimum_allowed_height as hmin |
---|
69 | if floor_heights is None: |
---|
70 | floor_heights = zeros(heights.shape, Float) |
---|
71 | |
---|
72 | all_vertices_below_threshold = True |
---|
73 | for k in range(3): |
---|
74 | w = heights[k] |
---|
75 | z = floor_heights[k] |
---|
76 | |
---|
77 | if w-z >= hmin: |
---|
78 | all_vertices_below_threshold = False |
---|
79 | |
---|
80 | vertex = self.panel.pos[k] |
---|
81 | vertex[2] = w |
---|
82 | |
---|
83 | vertex = self.backpanel.pos[2-k] |
---|
84 | vertex[2] = w |
---|
85 | |
---|
86 | |
---|
87 | #Update color to visualise dry areas |
---|
88 | if all_vertices_below_threshold: |
---|
89 | self.panel.color = self.bottom_color |
---|
90 | self.backpanel.color = self.bottom_color |
---|
91 | else: |
---|
92 | self.panel.color = self.top_color |
---|
93 | self.backpanel.color = self.top_color |
---|
94 | |
---|
95 | |
---|
96 | #update normal |
---|
97 | v0 = self.panel.pos[0] |
---|
98 | v1 = self.panel.pos[1] |
---|
99 | v2 = self.panel.pos[2] |
---|
100 | |
---|
101 | |
---|
102 | n = normal(v1-v0, v2-v0) |
---|
103 | self.panel.normal=n |
---|
104 | self.backpanel.normal=-n |
---|
105 | |
---|
106 | |
---|
107 | |
---|
108 | |
---|
109 | def create_surface(domain): |
---|
110 | """Link surface of domains to their visual representations |
---|
111 | """ |
---|
112 | |
---|
113 | fr = frame() #Default frame to contain all objects |
---|
114 | s=0 |
---|
115 | |
---|
116 | Q = domain.quantities['stage'] |
---|
117 | try: |
---|
118 | Z = domain.quantities['elevation'] |
---|
119 | bed = True |
---|
120 | except: |
---|
121 | bed = False |
---|
122 | |
---|
123 | N = Q.vertex_values.shape[0] |
---|
124 | |
---|
125 | domain.visuals = [] |
---|
126 | for i in range(N): |
---|
127 | |
---|
128 | z0 = Q.vertex_values[i, 0] |
---|
129 | z1 = Q.vertex_values[i, 1] |
---|
130 | z2 = Q.vertex_values[i, 2] |
---|
131 | |
---|
132 | x0, y0 = domain.get_vertex_coordinate(i,0) |
---|
133 | x1, y1 = domain.get_vertex_coordinate(i,1) |
---|
134 | x2, y2 = domain.get_vertex_coordinate(i,2) |
---|
135 | |
---|
136 | V0 = vector(x0, y0, z0) |
---|
137 | V1 = vector(x1, y1, z1) |
---|
138 | V2 = vector(x2, y2, z2) |
---|
139 | |
---|
140 | |
---|
141 | #Top surface |
---|
142 | c0 = 0.1 |
---|
143 | c1 = 0.4 |
---|
144 | c2 = 0.99 |
---|
145 | if s: |
---|
146 | c2 = 0.7*c2 #To show triangles in slightly different shades |
---|
147 | |
---|
148 | s = 1-s |
---|
149 | |
---|
150 | col = (c0,c1,c2) |
---|
151 | |
---|
152 | visual_top = Triangle(V0,V1,V2,color=col,frame=fr) |
---|
153 | visual_top.top_color = col |
---|
154 | |
---|
155 | |
---|
156 | #Bottom surface |
---|
157 | #v0, v1, v2 = volume.vertices |
---|
158 | #v0 = vector(v0.x, v0.y, v0.field_values[index]) |
---|
159 | #v1 = vector(v1.x, v1.y, v1.field_values[index]) |
---|
160 | #v2 = vector(v2.x, v2.y, v2.field_values[index]) |
---|
161 | if bed: |
---|
162 | z0 = Z.vertex_values[i, 0] |
---|
163 | z1 = Z.vertex_values[i, 1] |
---|
164 | z2 = Z.vertex_values[i, 2] |
---|
165 | else: |
---|
166 | z0 = z1 = z2 = 0.0 |
---|
167 | |
---|
168 | V0 = vector(x0, y0, z0) |
---|
169 | V1 = vector(x1, y1, z1) |
---|
170 | V2 = vector(x2, y2, z2) |
---|
171 | |
---|
172 | c0 = 0.3 |
---|
173 | c1 = 0.3 |
---|
174 | c2 = 0.3 |
---|
175 | if s: |
---|
176 | c2 = 0.4*c2 #To show triangles in slightly different shades |
---|
177 | |
---|
178 | col = (c0,c1,c2) |
---|
179 | visual_bottom = Triangle(V0, V1, V2, color=col,frame=fr) |
---|
180 | visual_top.bottom_color=col |
---|
181 | |
---|
182 | domain.visuals.append( (visual_top, visual_bottom) ) |
---|
183 | |
---|
184 | update(domain) |
---|
185 | #print 'Scale', scene.scale |
---|
186 | |
---|
187 | |
---|
188 | def update(domain): |
---|
189 | """Update vertex heights. |
---|
190 | The argument index refers to which conserved quantity to visualise. |
---|
191 | If domain.smooth is set True, vertex values will be averaged |
---|
192 | yielding a smoother surface. |
---|
193 | """ |
---|
194 | |
---|
195 | from Numeric import array |
---|
196 | |
---|
197 | |
---|
198 | Q = domain.quantities['stage'] |
---|
199 | N = Q.vertex_values.shape[0] |
---|
200 | |
---|
201 | #print scene.forward |
---|
202 | #FIXME: Use smoother from pyvolution instead |
---|
203 | if domain.smooth: |
---|
204 | #Get all average point values |
---|
205 | vertex_heights = {} |
---|
206 | |
---|
207 | for k in range(N): |
---|
208 | for i in range(3): |
---|
209 | vertex = domain.triangles[k, i] |
---|
210 | if vertex_heights.has_key(vertex): |
---|
211 | vertex_heights[vertex].append( |
---|
212 | Q.vertex_values[k, i]) |
---|
213 | else: |
---|
214 | vertex_heights[vertex] = [] |
---|
215 | vertex_heights[vertex].append( |
---|
216 | Q.vertex_values[k, i]) |
---|
217 | |
---|
218 | |
---|
219 | |
---|
220 | for k in range(N): |
---|
221 | if domain.smooth: |
---|
222 | #The averages |
---|
223 | x = zeros(3, Float) |
---|
224 | for i in range(3): |
---|
225 | vertex = domain.triangles[k, i] |
---|
226 | A = array(vertex_heights[vertex]) |
---|
227 | x[i] = sum(A)/len(A) |
---|
228 | else: |
---|
229 | #The true values |
---|
230 | x = [Q.vertex_values[k, 0], |
---|
231 | Q.vertex_values[k, 1], |
---|
232 | Q.vertex_values[k, 2]] |
---|
233 | |
---|
234 | |
---|
235 | #Do it |
---|
236 | floor_heights = array([pos[2] for pos in domain.visuals[k][1].panel.pos]) |
---|
237 | |
---|
238 | domain.visuals[k][0].set_vertexheights(x, floor_heights) |
---|
239 | |
---|
240 | |
---|
241 | |
---|
242 | |
---|
243 | scene.width = 1000 |
---|
244 | scene.height = 800 |
---|
245 | |
---|
246 | #Original |
---|
247 | scene.center = (0.5,0.5,0) |
---|
248 | scene.forward = vector(0.0, 0.5, -0.5) |
---|
249 | |
---|
250 | #Temporary (for bedslope) |
---|
251 | #scene.forward = vector(0.0006, 0.7, -0.03) |
---|
252 | |
---|
253 | |
---|
254 | #Temporary for hackett - begin |
---|
255 | #scene.autoscale = 0 |
---|
256 | #scene.scale = (0.002, 0.002, 0.01) #Scale z so that countours stand out more |
---|
257 | #scene.center = (300.0,500.0,-10) |
---|
258 | #Temporary for hackett - end |
---|
259 | |
---|
260 | |
---|
261 | scene.ambient = 0.4 |
---|
262 | scene.lights = [(0.6, 0.3, 0.2), (0.1, -0.5, 0.4), (-0.1, 0.1, -0.4), |
---|
263 | (-0.2, 0.2, 0.1)] |
---|
264 | |
---|
265 | |
---|
266 | |
---|
267 | |
---|