fileio.py - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
fileio.py (66599B)
---
1 import os
2 import numpy
3 from .common import VERSION, py_vtk, vtk
4
5
6 class SimIO:
7 'Binary and VTK input/output for sim objects.'
8
9 def readbin(self, targetbin, verbose=True, bonds=True, sigma0mod=True,
10 esysparticle=False):
11 '''
12 Reads a target ``sphere`` binary file.
13
14 See also :func:`writebin()`, :func:`readfirst()`, :func:`readlast()`,
15 :func:`readsecond`, and :func:`readstep`.
16
17 :param targetbin: The path to the binary ``sphere`` file
18 :type targetbin: str
19 :param verbose: Show diagnostic information (default=True)
20 :type verbose: bool
21 :param bonds: The input file contains bond information (default=True).
22 This parameter should be true for all recent ``sphere`` versions.
23 :type bonds: bool
24 :param sigma0mod: The input file contains information about modulating
25 stresses at the top wall (default=True). This parameter should be
26 true for all recent ``sphere`` versions.
27 :type sigma0mod: bool
28 :param esysparticle: Stop reading the file after reading the kinematics,
29 which is useful for reading output files from other DEM programs.
30 (default=False)
31 :type esysparticle: bool
32 '''
33
34 fh = None
35 try:
36 if verbose:
37 print("Input file: {0}".format(targetbin))
38 fh = open(targetbin, "rb")
39
40 # Read the file version
41 self.version = numpy.fromfile(fh, dtype=numpy.float64, count=1)
42
43 # Read the number of dimensions and particles
44 self.nd = int(numpy.fromfile(fh, dtype=numpy.int32, count=1)[0])
45 self.np = int(numpy.fromfile(fh, dtype=numpy.uint32, count=1)[0])
46
47 # Read the time variables
48 self.time_dt = numpy.fromfile(fh, dtype=numpy.float64, count=1)
49 self.time_current = numpy.fromfile(fh, dtype=numpy.float64, count=1)
50 self.time_total = numpy.fromfile(fh, dtype=numpy.float64, count=1)
51 self.time_file_dt = numpy.fromfile(fh, dtype=numpy.float64, count=1)
52 self.time_step_count = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
53
54 # Allocate array memory for particles
55 self.x = numpy.empty((self.np, self.nd), dtype=numpy.float64)
56 self.radius = numpy.empty(self.np, dtype=numpy.float64)
57 self.xyzsum = numpy.empty((self.np, 3), dtype=numpy.float64)
58 self.vel = numpy.empty((self.np, self.nd), dtype=numpy.float64)
59 self.fixvel = numpy.empty(self.np, dtype=numpy.float64)
60 self.es_dot = numpy.empty(self.np, dtype=numpy.float64)
61 self.es = numpy.empty(self.np, dtype=numpy.float64)
62 self.ev_dot = numpy.empty(self.np, dtype=numpy.float64)
63 self.ev = numpy.empty(self.np, dtype=numpy.float64)
64 self.p = numpy.empty(self.np, dtype=numpy.float64)
65
66 # Read remaining data from binary
67 self.origo = numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
68 self.L = numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
69 self.num = numpy.fromfile(fh, dtype=numpy.uint32, count=self.nd)
70 self.periodic = numpy.fromfile(fh, dtype=numpy.int32, count=1)
71
72 if self.version >= 2.14:
73 self.adaptive = numpy.fromfile(fh, dtype=numpy.int32, count=1)
74 else:
75 self.adaptive = numpy.zeros(1, dtype=numpy.float64)
76
77 # Per-particle vectors
78 for i in numpy.arange(self.np):
79 self.x[i, :] =\
80 numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
81 self.radius[i] =\
82 numpy.fromfile(fh, dtype=numpy.float64, count=1)[0]
83
84 if self.version >= 1.03:
85 self.xyzsum = numpy.fromfile(fh, dtype=numpy.float64,\
86 count=self.np*3).reshape(self.np, 3)
87 else:
88 self.xyzsum = numpy.fromfile(fh, dtype=numpy.float64,\
89 count=self.np*2).reshape(self.np, 2)
90
91 for i in numpy.arange(self.np):
92 self.vel[i, :] = numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
93 self.fixvel[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)[0]
94
95 self.force = numpy.fromfile(fh, dtype=numpy.float64,\
96 count=self.np*self.nd)\
97 .reshape(self.np, self.nd)
98
99 self.angpos = numpy.fromfile(fh, dtype=numpy.float64,\
100 count=self.np*self.nd)\
101 .reshape(self.np, self.nd)
102 self.angvel = numpy.fromfile(fh, dtype=numpy.float64,\
103 count=self.np*self.nd)\
104 .reshape(self.np, self.nd)
105 self.torque = numpy.fromfile(fh, dtype=numpy.float64,\
106 count=self.np*self.nd)\
107 .reshape(self.np, self.nd)
108
109 if esysparticle:
110 return
111
112 # Per-particle single-value parameters
113 self.es_dot = numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
114 self.es = numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
115 self.ev_dot = numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
116 self.ev = numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
117 self.p = numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
118
119 # Constant, global physical parameters
120 self.g = numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
121 self.k_n = numpy.fromfile(fh, dtype=numpy.float64, count=1)
122 self.k_t = numpy.fromfile(fh, dtype=numpy.float64, count=1)
123 self.k_r = numpy.fromfile(fh, dtype=numpy.float64, count=1)
124 if self.version >= 2.13:
125 self.E = numpy.fromfile(fh, dtype=numpy.float64, count=1)
126 else:
127 self.E = numpy.zeros(1, dtype=numpy.float64)
128 self.gamma_n = numpy.fromfile(fh, dtype=numpy.float64, count=1)
129 self.gamma_t = numpy.fromfile(fh, dtype=numpy.float64, count=1)
130 self.gamma_r = numpy.fromfile(fh, dtype=numpy.float64, count=1)
131 self.mu_s = numpy.fromfile(fh, dtype=numpy.float64, count=1)
132 self.mu_d = numpy.fromfile(fh, dtype=numpy.float64, count=1)
133 self.mu_r = numpy.fromfile(fh, dtype=numpy.float64, count=1)
134 self.gamma_wn = numpy.fromfile(fh, dtype=numpy.float64, count=1)
135 self.gamma_wt = numpy.fromfile(fh, dtype=numpy.float64, count=1)
136 self.mu_ws = numpy.fromfile(fh, dtype=numpy.float64, count=1)
137 self.mu_wd = numpy.fromfile(fh, dtype=numpy.float64, count=1)
138 self.rho = numpy.fromfile(fh, dtype=numpy.float64, count=1)
139 self.contactmodel = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
140 self.kappa = numpy.fromfile(fh, dtype=numpy.float64, count=1)
141 self.db = numpy.fromfile(fh, dtype=numpy.float64, count=1)
142 self.V_b = numpy.fromfile(fh, dtype=numpy.float64, count=1)
143
144 # Wall data
145 self.nw = int(numpy.fromfile(fh, dtype=numpy.uint32, count=1)[0])
146 self.wmode = numpy.empty(self.nw, dtype=numpy.int32)
147 self.w_n = numpy.empty(self.nw*self.nd, dtype=numpy.float64)\
148 .reshape(self.nw, self.nd)
149 self.w_x = numpy.empty(self.nw, dtype=numpy.float64)
150 self.w_m = numpy.empty(self.nw, dtype=numpy.float64)
151 self.w_vel = numpy.empty(self.nw, dtype=numpy.float64)
152 self.w_force = numpy.empty(self.nw, dtype=numpy.float64)
153 self.w_sigma0 = numpy.empty(self.nw, dtype=numpy.float64)
154
155 self.wmode = numpy.fromfile(fh, dtype=numpy.int32, count=self.nw)
156 for i in numpy.arange(self.nw):
157 self.w_n[i, :] =\
158 numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
159 self.w_x[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)[0]
160 for i in numpy.arange(self.nw):
161 self.w_m[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)[0]
162 self.w_vel[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)[0]
163 self.w_force[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)[0]
164 self.w_sigma0[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)[0]
165 if sigma0mod:
166 self.w_sigma0_A = numpy.fromfile(fh, dtype=numpy.float64, count=1)
167 self.w_sigma0_f = numpy.fromfile(fh, dtype=numpy.float64, count=1)
168 if self.version >= 2.1:
169 self.w_tau_x = numpy.fromfile(fh, dtype=numpy.float64, count=1)
170 else:
171 self.w_tau_x = numpy.zeros(1, dtype=numpy.float64)
172
173 if bonds:
174 # Inter-particle bonds
175 self.lambda_bar = numpy.fromfile(fh, dtype=numpy.float64, count=1)
176 self.nb0 = int(numpy.fromfile(fh, dtype=numpy.uint32, count=1)[0])
177 self.sigma_b = numpy.fromfile(fh, dtype=numpy.float64, count=1)
178 self.tau_b = numpy.fromfile(fh, dtype=numpy.float64, count=1)
179 self.bonds = numpy.empty((self.nb0, 2), dtype=numpy.uint32)
180 for i in numpy.arange(self.nb0):
181 self.bonds[i, 0] = numpy.fromfile(fh, dtype=numpy.uint32, count=1)[0]
182 self.bonds[i, 1] = numpy.fromfile(fh, dtype=numpy.uint32, count=1)[0]
183 self.bonds_delta_n = numpy.fromfile(fh, dtype=numpy.float64,
184 count=self.nb0)
185 self.bonds_delta_t = numpy.fromfile(fh, dtype=numpy.float64,
186 count=self.nb0*self.nd)\
187 .reshape(self.nb0, self.nd)
188 self.bonds_omega_n = numpy.fromfile(fh, dtype=numpy.float64,
189 count=self.nb0)
190 self.bonds_omega_t = numpy.fromfile(fh, dtype=numpy.float64,
191 count=self.nb0*self.nd)\
192 .reshape(self.nb0, self.nd)
193 else:
194 self.nb0 = 0
195
196 if self.fluid:
197
198 if self.version >= 2.0:
199 self.cfd_solver = numpy.fromfile(fh, dtype=numpy.int32, count=1)
200 else:
201 self.cfd_solver = numpy.zeros(1, dtype=numpy.int32)
202
203 self.mu = numpy.fromfile(fh, dtype=numpy.float64, count=1)
204
205 self.v_f = numpy.empty((self.num[0],
206 self.num[1],
207 self.num[2],
208 self.nd), dtype=numpy.float64)
209 self.p_f = numpy.empty((self.num[0],
210 self.num[1],
211 self.num[2]), dtype=numpy.float64)
212 self.phi = numpy.empty((self.num[0],
213 self.num[1],
214 self.num[2]), dtype=numpy.float64)
215 self.dphi = numpy.empty((self.num[0],
216 self.num[1],
217 self.num[2]), dtype=numpy.float64)
218
219 for z in numpy.arange(self.num[2]):
220 for y in numpy.arange(self.num[1]):
221 for x in numpy.arange(self.num[0]):
222 self.v_f[x, y, z, 0] = numpy.fromfile(fh,
223 dtype=numpy.float64,
224 count=1)[0]
225 self.v_f[x, y, z, 1] = numpy.fromfile(fh,
226 dtype=numpy.float64,
227 count=1)[0]
228 self.v_f[x, y, z, 2] = numpy.fromfile(fh,
229 dtype=numpy.float64,
230 count=1)[0]
231 self.p_f[x, y, z] = numpy.fromfile(fh,
232 dtype=numpy.float64,
233 count=1)[0]
234 self.phi[x, y, z] = numpy.fromfile(fh,
235 dtype=numpy.float64,
236 count=1)[0]
237 self.dphi[x, y, z] = numpy.fromfile(fh,
238 dtype=numpy.float64,
239 count=1)[0]\
240 /(self.time_dt[0]
241 *self.ndem.item())
242
243 if self.version >= 0.36:
244 self.rho_f = numpy.fromfile(fh, dtype=numpy.float64,
245 count=1)
246 self.p_mod_A = numpy.fromfile(fh, dtype=numpy.float64,
247 count=1)
248 self.p_mod_f = numpy.fromfile(fh, dtype=numpy.float64,
249 count=1)
250 self.p_mod_phi = numpy.fromfile(fh, dtype=numpy.float64,
251 count=1)
252
253 if self.version >= 2.12 and self.cfd_solver[0] == 1:
254 self.bc_xn = numpy.fromfile(fh, dtype=numpy.int32,
255 count=1)
256 self.bc_xp = numpy.fromfile(fh, dtype=numpy.int32,
257 count=1)
258 self.bc_yn = numpy.fromfile(fh, dtype=numpy.int32,
259 count=1)
260 self.bc_yp = numpy.fromfile(fh, dtype=numpy.int32,
261 count=1)
262
263 self.bc_bot = numpy.fromfile(fh, dtype=numpy.int32, count=1)
264 self.bc_top = numpy.fromfile(fh, dtype=numpy.int32, count=1)
265 self.free_slip_bot = numpy.fromfile(fh, dtype=numpy.int32,
266 count=1)
267 self.free_slip_top = numpy.fromfile(fh, dtype=numpy.int32,
268 count=1)
269 if self.version >= 2.11:
270 self.bc_bot_flux = numpy.fromfile(fh,
271 dtype=numpy.float64,
272 count=1)
273 self.bc_top_flux = numpy.fromfile(fh,
274 dtype=numpy.float64,
275 count=1)
276 else:
277 self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
278 self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
279
280 if self.version >= 2.15:
281 self.p_f_constant = numpy.empty((self.num[0],
282 self.num[1],
283 self.num[2]),
284 dtype=numpy.int32)
285
286 for z in numpy.arange(self.num[2]):
287 for y in numpy.arange(self.num[1]):
288 for x in numpy.arange(self.num[0]):
289 self.p_f_constant[x, y, z] = \
290 numpy.fromfile(fh, dtype=numpy.int32,
291 count=1)[0]
292 else:
293 self.p_f_constant = numpy.zeros((self.num[0],
294 self.num[1],
295 self.num[2]),
296 dtype=numpy.int32)
297
298 if self.version >= 2.0 and self.cfd_solver == 0:
299 self.gamma = numpy.fromfile(fh, dtype=numpy.float64,
300 count=1)
301 self.theta = numpy.fromfile(fh, dtype=numpy.float64,
302 count=1)
303 self.beta = numpy.fromfile(fh, dtype=numpy.float64,
304 count=1)
305 self.tolerance = numpy.fromfile(fh, dtype=numpy.float64,
306 count=1)
307 self.maxiter = numpy.fromfile(fh, dtype=numpy.uint32,
308 count=1)
309 if self.version >= 1.01:
310 self.ndem = numpy.fromfile(fh, dtype=numpy.uint32,
311 count=1)
312 else:
313 self.ndem = 1
314
315 if self.version >= 1.04:
316 self.c_phi = numpy.fromfile(fh, dtype=numpy.float64,
317 count=1)
318 self.c_v = numpy.fromfile(fh, dtype=numpy.float64,
319 count=1)
320 if self.version == 1.06:
321 self.c_a = numpy.fromfile(fh, dtype=numpy.float64,
322 count=1)
323 elif self.version >= 1.07:
324 self.dt_dem_fac = numpy.fromfile(fh,
325 dtype=numpy.float64,
326 count=1)
327 else:
328 self.c_a = numpy.ones(1, dtype=numpy.float64)
329 else:
330 self.c_phi = numpy.ones(1, dtype=numpy.float64)
331 self.c_v = numpy.ones(1, dtype=numpy.float64)
332
333 if self.version >= 1.05:
334 self.f_d = numpy.empty_like(self.x)
335 self.f_p = numpy.empty_like(self.x)
336 self.f_v = numpy.empty_like(self.x)
337 self.f_sum = numpy.empty_like(self.x)
338
339 for i in numpy.arange(self.np):
340 self.f_d[i, :] = numpy.fromfile(fh,
341 dtype=numpy.float64,
342 count=self.nd)
343 for i in numpy.arange(self.np):
344 self.f_p[i, :] = numpy.fromfile(fh,
345 dtype=numpy.float64,
346 count=self.nd)
347 for i in numpy.arange(self.np):
348 self.f_v[i, :] = numpy.fromfile(fh,
349 dtype=numpy.float64,
350 count=self.nd)
351 for i in numpy.arange(self.np):
352 self.f_sum[i, :] = numpy.fromfile(fh,
353 dtype=numpy.float64,
354 count=self.nd)
355 else:
356 self.f_d = numpy.zeros((self.np, self.nd),
357 dtype=numpy.float64)
358 self.f_p = numpy.zeros((self.np, self.nd),
359 dtype=numpy.float64)
360 self.f_v = numpy.zeros((self.np, self.nd),
361 dtype=numpy.float64)
362 self.f_sum = numpy.zeros((self.np, self.nd),
363 dtype=numpy.float64)
364
365 elif self.version >= 2.0 and self.cfd_solver == 1:
366
367 self.tolerance = numpy.fromfile(fh, dtype=numpy.float64,
368 count=1)
369 self.maxiter = numpy.fromfile(fh, dtype=numpy.uint32,
370 count=1)
371 self.ndem = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
372 self.c_phi = numpy.fromfile(fh, dtype=numpy.float64,
373 count=1)
374 self.f_p = numpy.empty_like(self.x)
375 for i in numpy.arange(self.np):
376 self.f_p[i, :] = numpy.fromfile(fh, dtype=numpy.float64,
377 count=self.nd)
378 self.beta_f = numpy.fromfile(fh, dtype=numpy.float64,
379 count=1)
380 self.k_c = numpy.fromfile(fh, dtype=numpy.float64, count=1)
381
382 if self.version >= 1.02:
383 self.color = numpy.fromfile(fh, dtype=numpy.int32,
384 count=self.np)
385 else:
386 self.color = numpy.zeros(self.np, dtype=numpy.int32)
387
388 finally:
389 self.version[0] = VERSION
390 if fh is not None:
391 fh.close()
392
393 def writebin(self, folder="../input/", verbose=True):
394 '''
395 Writes a ``sphere`` binary file to the ``../input/`` folder by default.
396 The file name will be in the format ``<self.sid>.bin``.
397
398 See also :func:`readbin()`.
399
400 :param folder: The folder where to place the output binary file
401 :type folder: str
402 :param verbose: Show diagnostic information (default=True)
403 :type verbose: bool
404 '''
405 fh = None
406 try:
407 targetbin = folder + "/" + self.sid + ".bin"
408 if verbose:
409 print("Output file: {0}".format(targetbin))
410
411 fh = open(targetbin, "wb")
412
413 # Write the current version number
414 fh.write(self.version.astype(numpy.float64))
415
416 # Write the number of dimensions and particles
417 fh.write(numpy.array(self.nd).astype(numpy.int32))
418 fh.write(numpy.array(self.np).astype(numpy.uint32))
419
420 # Write the time variables
421 fh.write(self.time_dt.astype(numpy.float64))
422 fh.write(self.time_current.astype(numpy.float64))
423 fh.write(self.time_total.astype(numpy.float64))
424 fh.write(self.time_file_dt.astype(numpy.float64))
425 fh.write(self.time_step_count.astype(numpy.uint32))
426
427 # Read remaining data from binary
428 fh.write(self.origo.astype(numpy.float64))
429 fh.write(self.L.astype(numpy.float64))
430 fh.write(self.num.astype(numpy.uint32))
431 fh.write(self.periodic.astype(numpy.uint32))
432 fh.write(self.adaptive.astype(numpy.uint32))
433
434 # Per-particle vectors
435 for i in numpy.arange(self.np):
436 fh.write(self.x[i, :].astype(numpy.float64))
437 fh.write(self.radius[i].astype(numpy.float64))
438
439 if self.np > 0:
440 fh.write(self.xyzsum.astype(numpy.float64))
441
442 for i in numpy.arange(self.np):
443 fh.write(self.vel[i, :].astype(numpy.float64))
444 fh.write(self.fixvel[i].astype(numpy.float64))
445
446 if self.np > 0:
447 fh.write(self.force.astype(numpy.float64))
448
449 fh.write(self.angpos.astype(numpy.float64))
450 fh.write(self.angvel.astype(numpy.float64))
451 fh.write(self.torque.astype(numpy.float64))
452
453 # Per-particle single-value parameters
454 fh.write(self.es_dot.astype(numpy.float64))
455 fh.write(self.es.astype(numpy.float64))
456 fh.write(self.ev_dot.astype(numpy.float64))
457 fh.write(self.ev.astype(numpy.float64))
458 fh.write(self.p.astype(numpy.float64))
459
460 fh.write(self.g.astype(numpy.float64))
461 fh.write(self.k_n.astype(numpy.float64))
462 fh.write(self.k_t.astype(numpy.float64))
463 fh.write(self.k_r.astype(numpy.float64))
464 fh.write(self.E.astype(numpy.float64))
465 fh.write(self.gamma_n.astype(numpy.float64))
466 fh.write(self.gamma_t.astype(numpy.float64))
467 fh.write(self.gamma_r.astype(numpy.float64))
468 fh.write(self.mu_s.astype(numpy.float64))
469 fh.write(self.mu_d.astype(numpy.float64))
470 fh.write(self.mu_r.astype(numpy.float64))
471 fh.write(self.gamma_wn.astype(numpy.float64))
472 fh.write(self.gamma_wt.astype(numpy.float64))
473 fh.write(self.mu_ws.astype(numpy.float64))
474 fh.write(self.mu_wd.astype(numpy.float64))
475 fh.write(self.rho.astype(numpy.float64))
476 fh.write(self.contactmodel.astype(numpy.uint32))
477 fh.write(self.kappa.astype(numpy.float64))
478 fh.write(self.db.astype(numpy.float64))
479 fh.write(self.V_b.astype(numpy.float64))
480
481 fh.write(numpy.array(self.nw).astype(numpy.uint32))
482 for i in numpy.arange(self.nw):
483 fh.write(self.wmode[i].astype(numpy.int32))
484 for i in numpy.arange(self.nw):
485 fh.write(self.w_n[i, :].astype(numpy.float64))
486 fh.write(self.w_x[i].astype(numpy.float64))
487
488 for i in numpy.arange(self.nw):
489 fh.write(self.w_m[i].astype(numpy.float64))
490 fh.write(self.w_vel[i].astype(numpy.float64))
491 fh.write(self.w_force[i].astype(numpy.float64))
492 fh.write(self.w_sigma0[i].astype(numpy.float64))
493 fh.write(self.w_sigma0_A.astype(numpy.float64))
494 fh.write(self.w_sigma0_f.astype(numpy.float64))
495 fh.write(self.w_tau_x.astype(numpy.float64))
496
497 fh.write(self.lambda_bar.astype(numpy.float64))
498 fh.write(numpy.array(self.nb0).astype(numpy.uint32))
499 fh.write(self.sigma_b.astype(numpy.float64))
500 fh.write(self.tau_b.astype(numpy.float64))
501 for i in numpy.arange(self.nb0):
502 fh.write(self.bonds[i, 0].astype(numpy.uint32))
503 fh.write(self.bonds[i, 1].astype(numpy.uint32))
504 fh.write(self.bonds_delta_n.astype(numpy.float64))
505 fh.write(self.bonds_delta_t.astype(numpy.float64))
506 fh.write(self.bonds_omega_n.astype(numpy.float64))
507 fh.write(self.bonds_omega_t.astype(numpy.float64))
508
509 if self.fluid:
510
511 fh.write(self.cfd_solver.astype(numpy.int32))
512 fh.write(self.mu.astype(numpy.float64))
513 for z in numpy.arange(self.num[2]):
514 for y in numpy.arange(self.num[1]):
515 for x in numpy.arange(self.num[0]):
516 fh.write(self.v_f[x, y, z, 0].astype(numpy.float64))
517 fh.write(self.v_f[x, y, z, 1].astype(numpy.float64))
518 fh.write(self.v_f[x, y, z, 2].astype(numpy.float64))
519 fh.write(self.p_f[x, y, z].astype(numpy.float64))
520 fh.write(self.phi[x, y, z].astype(numpy.float64))
521 fh.write(self.dphi[x, y, z].astype(numpy.float64)*
522 self.time_dt*self.ndem)
523
524 fh.write(self.rho_f.astype(numpy.float64))
525 fh.write(self.p_mod_A.astype(numpy.float64))
526 fh.write(self.p_mod_f.astype(numpy.float64))
527 fh.write(self.p_mod_phi.astype(numpy.float64))
528
529 if self.cfd_solver[0] == 1: # Sides only adjustable with Darcy
530 fh.write(self.bc_xn.astype(numpy.int32))
531 fh.write(self.bc_xp.astype(numpy.int32))
532 fh.write(self.bc_yn.astype(numpy.int32))
533 fh.write(self.bc_yp.astype(numpy.int32))
534
535 fh.write(self.bc_bot.astype(numpy.int32))
536 fh.write(self.bc_top.astype(numpy.int32))
537 fh.write(self.free_slip_bot.astype(numpy.int32))
538 fh.write(self.free_slip_top.astype(numpy.int32))
539 fh.write(self.bc_bot_flux.astype(numpy.float64))
540 fh.write(self.bc_top_flux.astype(numpy.float64))
541
542 for z in numpy.arange(self.num[2]):
543 for y in numpy.arange(self.num[1]):
544 for x in numpy.arange(self.num[0]):
545 fh.write(self.p_f_constant[x, y, z].astype(
546 numpy.int32))
547
548 # Navier Stokes
549 if self.cfd_solver[0] == 0:
550 fh.write(self.gamma.astype(numpy.float64))
551 fh.write(self.theta.astype(numpy.float64))
552 fh.write(self.beta.astype(numpy.float64))
553 fh.write(self.tolerance.astype(numpy.float64))
554 fh.write(self.maxiter.astype(numpy.uint32))
555 fh.write(self.ndem.astype(numpy.uint32))
556
557 fh.write(self.c_phi.astype(numpy.float64))
558 fh.write(self.c_v.astype(numpy.float64))
559 fh.write(self.dt_dem_fac.astype(numpy.float64))
560
561 for i in numpy.arange(self.np):
562 fh.write(self.f_d[i, :].astype(numpy.float64))
563 for i in numpy.arange(self.np):
564 fh.write(self.f_p[i, :].astype(numpy.float64))
565 for i in numpy.arange(self.np):
566 fh.write(self.f_v[i, :].astype(numpy.float64))
567 for i in numpy.arange(self.np):
568 fh.write(self.f_sum[i, :].astype(numpy.float64))
569
570 # Darcy
571 elif self.cfd_solver[0] == 1:
572
573 fh.write(self.tolerance.astype(numpy.float64))
574 fh.write(self.maxiter.astype(numpy.uint32))
575 fh.write(self.ndem.astype(numpy.uint32))
576 fh.write(self.c_phi.astype(numpy.float64))
577 for i in numpy.arange(self.np):
578 fh.write(self.f_p[i, :].astype(numpy.float64))
579 fh.write(self.beta_f.astype(numpy.float64))
580 fh.write(self.k_c.astype(numpy.float64))
581
582 else:
583 raise Exception('Value of cfd_solver not understood (' + \
584 str(self.cfd_solver[0]) + ')')
585
586
587 fh.write(self.color.astype(numpy.int32))
588
589 finally:
590 if fh is not None:
591 fh.close()
592
593 def writeVTKall(self, cell_centered=True, verbose=True, forces=False):
594 '''
595 Writes a VTK file for each simulation output file with particle
596 information and the fluid grid to the ``../output/`` folder by default.
597 The file name will be in the format ``<self.sid>.vtu`` and
598 ``fluid-<self.sid>.vti``. The vtu files can be used to visualize the
599 particles, and the vti files for visualizing the fluid in ParaView.
600
601 After opening the vtu files, the particle fields will show up in the
602 "Properties" list. Press "Apply" to import all fields into the ParaView
603 session. The particles are visualized by selecting the imported data in
604 the "Pipeline Browser". Afterwards, click the "Glyph" button in the
605 "Common" toolbar, or go to the "Filters" menu, and press "Glyph" from
606 the "Common" list. Choose "Sphere" as the "Glyph Type", set "Radius" to
607 1.0, choose "scalar" as the "Scale Mode". Check the "Edit" checkbox, and
608 set the "Set Scale Factor" to 1.0. The field "Maximum Number of Points"
609 may be increased if the number of particles exceed the default value.
610 Finally press "Apply", and the particles will appear in the main window.
611
612 The sphere resolution may be adjusted ("Theta resolution", "Phi
613 resolution") to increase the quality and the computational requirements
614 of the rendering.
615
616 The fluid grid is visualized by opening the vti files, and pressing
617 "Apply" to import all fluid field properties. To visualize the scalar
618 fields, such as the pressure, the porosity, the porosity change or the
619 velocity magnitude, choose "Surface" or "Surface With Edges" as the
620 "Representation". Choose the desired property as the "Coloring" field.
621 It may be desirable to show the color bar by pressing the "Show" button,
622 and "Rescale" to fit the color range limits to the current file. The
623 coordinate system can be displayed by checking the "Show Axis" field.
624 All adjustments by default require the "Apply" button to be pressed
625 before regenerating the view.
626
627 The fluid vector fields (e.g. the fluid velocity) can be visualizing by
628 e.g. arrows. To do this, select the fluid data in the "Pipeline
629 Browser". Press "Glyph" from the "Common" toolbar, or go to the
630 "Filters" mennu, and press "Glyph" from the "Common" list. Make sure
631 that "Arrow" is selected as the "Glyph type", and "Velocity" as the
632 "Vectors" value. Adjust the "Maximum Number of Points" to be at least as
633 big as the number of fluid cells in the grid. Press "Apply" to visualize
634 the arrows.
635
636 If several data files are generated for the same simulation (e.g. using
637 the :func:`writeVTKall()` function), it is able to step the
638 visualization through time by using the ParaView controls.
639
640 :param verbose: Show diagnostic information (default=True)
641 :type verbose: bool
642 :param cell_centered: Write fluid values to cell centered positions
643 (default=true)
644 :type cell_centered: bool
645 :param forces: Write contact force files (slow) (default=False)
646 :type forces: bool
647 '''
648 lastfile = self.status()
649 from .core import sim
650 sb = sim(fluid=self.fluid)
651 for i in range(lastfile+1):
652 fn = "../output/{0}.output{1:0=5}.bin".format(self.sid, i)
653
654 # check if output VTK file exists and if it is newer than spherebin
655 fn_vtk = "../output/{0}.{1:0=5}.vtu".format(self.sid, i)
656 if os.path.isfile(fn_vtk) and \
657 (os.path.getmtime(fn) < os.path.getmtime(fn_vtk)):
658 if verbose:
659 print('skipping ' + fn_vtk +
660 ': file exists and is newer than ' + fn)
661 if self.fluid:
662 fn_vtk = "../output/fluid-{0}.{1:0=5}.vti" \
663 .format(self.sid, i)
664 if os.path.isfile(fn_vtk) and \
665 (os.path.getmtime(fn) < os.path.getmtime(fn_vtk)):
666 if verbose:
667 print('skipping ' + fn_vtk +
668 ': file exists and is newer than ' + fn)
669 continue
670 else:
671 continue
672
673 sb.sid = self.sid + ".{:0=5}".format(i)
674 sb.readbin(fn, verbose=False)
675 if sb.np > 0:
676 if i == 0 or i == lastfile:
677 if i == lastfile:
678 if verbose:
679 print("\tto")
680 sb.writeVTK(verbose=verbose)
681 if forces:
682 sb.findContactStresses()
683 sb.writeVTKforces(verbose=verbose)
684 else:
685 sb.writeVTK(verbose=False)
686 if forces:
687 sb.findContactStresses()
688 sb.writeVTKforces(verbose=False)
689 if self.fluid:
690 if i == 0 or i == lastfile:
691 if i == lastfile:
692 if verbose:
693 print("\tto")
694 sb.writeFluidVTK(verbose=verbose,
695 cell_centered=cell_centered)
696 else:
697 sb.writeFluidVTK(verbose=False, cell_centered=cell_centered)
698
699 def writeVTK(self, folder='../output/', verbose=True):
700 '''
701 Writes a VTK file with particle information to the ``../output/`` folder
702 by default. The file name will be in the format ``<self.sid>.vtu``.
703 The vtu files can be used to visualize the particles in ParaView.
704
705 After opening the vtu files, the particle fields will show up in the
706 "Properties" list. Press "Apply" to import all fields into the ParaView
707 session. The particles are visualized by selecting the imported data in
708 the "Pipeline Browser". Afterwards, click the "Glyph" button in the
709 "Common" toolbar, or go to the "Filters" menu, and press "Glyph" from
710 the "Common" list. Choose "Sphere" as the "Glyph Type", choose "scalar"
711 as the "Scale Mode". Check the "Edit" checkbox, and set the "Set Scale
712 Factor" to 1.0. The field "Maximum Number of Points" may be increased if
713 the number of particles exceed the default value. Finally press "Apply",
714 and the particles will appear in the main window.
715
716 The sphere resolution may be adjusted ("Theta resolution", "Phi
717 resolution") to increase the quality and the computational requirements
718 of the rendering. All adjustments by default require the "Apply" button
719 to be pressed before regenerating the view.
720
721 If several vtu files are generated for the same simulation (e.g. using
722 the :func:`writeVTKall()` function), it is able to step the
723 visualization through time by using the ParaView controls.
724
725 :param folder: The folder where to place the output binary file (default
726 (default='../output/')
727 :type folder: str
728 :param verbose: Show diagnostic information (default=True)
729 :type verbose: bool
730 '''
731
732 fh = None
733 try:
734 targetbin = folder + '/' + self.sid + '.vtu' # unstructured grid
735 if verbose:
736 print('Output file: ' + targetbin)
737
738 fh = open(targetbin, 'w')
739
740 # the VTK data file format is documented in
741 # http://www.vtk.org/VTK/img/file-formats.pdf
742
743 fh.write('<?xml version="1.0"?>\n') # XML header
744 fh.write('<VTKFile type="UnstructuredGrid" version="0.1" '
745 + 'byte_order="LittleEndian">\n') # VTK header
746 fh.write(' <UnstructuredGrid>\n')
747 fh.write(' <Piece NumberOfPoints="%d" NumberOfCells="0">\n' \
748 % (self.np))
749
750 # Coordinates for each point (positions)
751 fh.write(' <Points>\n')
752 fh.write(' <DataArray name="Position [m]" type="Float32" '
753 + 'NumberOfComponents="3" format="ascii">\n')
754 fh.write(' ')
755 for i in range(self.np):
756 fh.write('%f %f %f ' % (self.x[i, 0], self.x[i, 1], self.x[i, 2]))
757 fh.write('\n')
758 fh.write(' </DataArray>\n')
759 fh.write(' </Points>\n')
760
761 ### Data attributes
762 fh.write(' <PointData Scalars="Diameter [m]" Vectors="vector">\n')
763
764 # Radii
765 fh.write(' <DataArray type="Float32" Name="Diameter" '
766 + 'format="ascii">\n')
767 fh.write(' ')
768 for i in range(self.np):
769 fh.write('%f ' % (self.radius[i]*2.0))
770 fh.write('\n')
771 fh.write(' </DataArray>\n')
772
773 # Displacements (xyzsum)
774 fh.write(' <DataArray type="Float32" Name="Displacement [m]" '
775 + 'NumberOfComponents="3" format="ascii">\n')
776 fh.write(' ')
777 for i in range(self.np):
778 fh.write('%f %f %f ' % \
779 (self.xyzsum[i, 0], self.xyzsum[i, 1], self.xyzsum[i, 2]))
780 fh.write('\n')
781 fh.write(' </DataArray>\n')
782
783 # Velocity
784 fh.write(' <DataArray type="Float32" Name="Velocity [m/s]" '
785 + 'NumberOfComponents="3" format="ascii">\n')
786 fh.write(' ')
787 for i in range(self.np):
788 fh.write('%f %f %f ' % \
789 (self.vel[i, 0], self.vel[i, 1], self.vel[i, 2]))
790 fh.write('\n')
791 fh.write(' </DataArray>\n')
792
793 if self.fluid:
794
795 if self.cfd_solver == 0: # Navier Stokes
796 # Fluid interaction force
797 fh.write(' <DataArray type="Float32" '
798 + 'Name="Fluid force total [N]" '
799 + 'NumberOfComponents="3" format="ascii">\n')
800 fh.write(' ')
801 for i in range(self.np):
802 fh.write('%f %f %f ' % \
803 (self.f_sum[i, 0], self.f_sum[i, 1], \
804 self.f_sum[i, 2]))
805 fh.write('\n')
806 fh.write(' </DataArray>\n')
807
808 # Fluid drag force
809 fh.write(' <DataArray type="Float32" '
810 + 'Name="Fluid drag force [N]" '
811 + 'NumberOfComponents="3" format="ascii">\n')
812 fh.write(' ')
813 for i in range(self.np):
814 fh.write('%f %f %f ' % \
815 (self.f_d[i, 0],
816 self.f_d[i, 1],
817 self.f_d[i, 2]))
818 fh.write('\n')
819 fh.write(' </DataArray>\n')
820
821 # Fluid pressure force
822 fh.write(' <DataArray type="Float32" '
823 + 'Name="Fluid pressure force [N]" '
824 + 'NumberOfComponents="3" format="ascii">\n')
825 fh.write(' ')
826 for i in range(self.np):
827 fh.write('%f %f %f ' % \
828 (self.f_p[i, 0], self.f_p[i, 1], self.f_p[i, 2]))
829 fh.write('\n')
830 fh.write(' </DataArray>\n')
831
832 if self.cfd_solver == 0: # Navier Stokes
833 # Fluid viscous force
834 fh.write(' <DataArray type="Float32" '
835 + 'Name="Fluid viscous force [N]" '
836 + 'NumberOfComponents="3" format="ascii">\n')
837 fh.write(' ')
838 for i in range(self.np):
839 fh.write('%f %f %f ' % \
840 (self.f_v[i, 0],
841 self.f_v[i, 1],
842 self.f_v[i, 2]))
843 fh.write('\n')
844 fh.write(' </DataArray>\n')
845
846 # fixvel
847 fh.write(' <DataArray type="Float32" Name="FixedVel" '
848 + 'format="ascii">\n')
849 fh.write(' ')
850 for i in range(self.np):
851 fh.write('%f ' % (self.fixvel[i]))
852 fh.write('\n')
853 fh.write(' </DataArray>\n')
854
855 # Force
856 fh.write(' <DataArray type="Float32" Name="Force [N]" '
857 + 'NumberOfComponents="3" format="ascii">\n')
858 fh.write(' ')
859 for i in range(self.np):
860 fh.write('%f %f %f ' % (self.force[i, 0],
861 self.force[i, 1],
862 self.force[i, 2]))
863 fh.write('\n')
864 fh.write(' </DataArray>\n')
865
866 # Angular Position
867 fh.write(' <DataArray type="Float32" Name="Angular position'
868 + '[rad]" '
869 + 'NumberOfComponents="3" format="ascii">\n')
870 fh.write(' ')
871 for i in range(self.np):
872 fh.write('%f %f %f ' % (self.angpos[i, 0],
873 self.angpos[i, 1],
874 self.angpos[i, 2]))
875 fh.write('\n')
876 fh.write(' </DataArray>\n')
877
878 # Angular Velocity
879 fh.write(' <DataArray type="Float32" Name="Angular velocity'
880 + ' [rad/s]" '
881 + 'NumberOfComponents="3" format="ascii">\n')
882 fh.write(' ')
883 for i in range(self.np):
884 fh.write('%f %f %f ' % (self.angvel[i, 0],
885 self.angvel[i, 1],
886 self.angvel[i, 2]))
887 fh.write('\n')
888 fh.write(' </DataArray>\n')
889
890 # Torque
891 fh.write(' <DataArray type="Float32" Name="Torque [Nm]" '
892 + 'NumberOfComponents="3" format="ascii">\n')
893 fh.write(' ')
894 for i in range(self.np):
895 fh.write('%f %f %f ' % (self.torque[i, 0],
896 self.torque[i, 1],
897 self.torque[i, 2]))
898 fh.write('\n')
899 fh.write(' </DataArray>\n')
900
901 # Shear energy rate
902 fh.write(' <DataArray type="Float32" Name="Shear Energy '
903 + 'Rate [J/s]" '
904 + 'format="ascii">\n')
905 fh.write(' ')
906 for i in range(self.np):
907 fh.write('%f ' % (self.es_dot[i]))
908 fh.write('\n')
909 fh.write(' </DataArray>\n')
910
911 # Shear energy
912 fh.write(' <DataArray type="Float32" Name="Shear Energy [J]"'
913 + ' format="ascii">\n')
914 fh.write(' ')
915 for i in range(self.np):
916 fh.write('%f ' % (self.es[i]))
917 fh.write('\n')
918 fh.write(' </DataArray>\n')
919
920 # Viscous energy rate
921 fh.write(' <DataArray type="Float32" '
922 + 'Name="Viscous Energy Rate [J/s]" format="ascii">\n')
923 fh.write(' ')
924 for i in range(self.np):
925 fh.write('%f ' % (self.ev_dot[i]))
926 fh.write('\n')
927 fh.write(' </DataArray>\n')
928
929 # Shear energy
930 fh.write(' <DataArray type="Float32" '
931 + 'Name="Viscous Energy [J]" '
932 + 'format="ascii">\n')
933 fh.write(' ')
934 for i in range(self.np):
935 fh.write('%f ' % (self.ev[i]))
936 fh.write('\n')
937 fh.write(' </DataArray>\n')
938
939 # Pressure
940 fh.write(' <DataArray type="Float32" Name="Pressure [Pa]" '
941 + 'format="ascii">\n')
942 fh.write(' ')
943 for i in range(self.np):
944 fh.write('%f ' % (self.p[i]))
945 fh.write('\n')
946 fh.write(' </DataArray>\n')
947
948 # Color
949 fh.write(' <DataArray type="Int32" Name="Type color" '
950 + 'format="ascii">\n')
951 fh.write(' ')
952 for i in range(self.np):
953 fh.write('%d ' % (self.color[i]))
954 fh.write('\n')
955 fh.write(' </DataArray>\n')
956
957 # Footer
958 fh.write(' </PointData>\n')
959 fh.write(' <Cells>\n')
960 fh.write(' <DataArray type="Int32" Name="connectivity" '
961 + 'format="ascii">\n')
962 fh.write(' </DataArray>\n')
963 fh.write(' <DataArray type="Int32" Name="offsets" '
964 + 'format="ascii">\n')
965 fh.write(' </DataArray>\n')
966 fh.write(' <DataArray type="UInt8" Name="types" '
967 + 'format="ascii">\n')
968 fh.write(' </DataArray>\n')
969 fh.write(' </Cells>\n')
970 fh.write(' </Piece>\n')
971 fh.write(' </UnstructuredGrid>\n')
972 fh.write('</VTKFile>')
973
974 finally:
975 if fh is not None:
976 fh.close()
977
978 def writeVTKforces(self, folder='../output/', verbose=True):
979 '''
980 Writes a VTK file with particle-interaction information to the
981 ``../output/`` folder by default. The file name will be in the format
982 ``<self.sid>.vtp``. The vtp files can be used to visualize the
983 particle interactions in ParaView. First use the "Cell Data to Point
984 Data" filter, and afterwards show the contact network with the "Tube"
985 filter.
986
987 :param folder: The folder where to place the output file (default
988 (default='../output/')
989 :type folder: str
990 :param verbose: Show diagnostic information (default=True)
991 :type verbose: bool
992 '''
993
994 if not py_vtk:
995 print('Error: vtk module not found, cannot writeVTKforces.')
996 return
997
998 filename = folder + '/forces-' + self.sid + '.vtp' # Polygon data
999
1000 # points mark the particle centers
1001 points = vtk.vtkPoints()
1002
1003 # lines mark the particle connectivity
1004 lines = vtk.vtkCellArray()
1005
1006 # colors
1007 #colors = vtk.vtkUnsignedCharArray()
1008 #colors.SetNumberOfComponents(3)
1009 #colors.SetName('Colors')
1010 #colors.SetNumberOfTuples(self.overlaps.size)
1011
1012 # scalars
1013 forces = vtk.vtkDoubleArray()
1014 forces.SetName("Force [N]")
1015 forces.SetNumberOfComponents(1)
1016 #forces.SetNumberOfTuples(self.overlaps.size)
1017 forces.SetNumberOfValues(self.overlaps.size)
1018
1019 stresses = vtk.vtkDoubleArray()
1020 stresses.SetName("Stress [Pa]")
1021 stresses.SetNumberOfComponents(1)
1022 stresses.SetNumberOfValues(self.overlaps.size)
1023
1024 for i in numpy.arange(self.overlaps.size):
1025 points.InsertNextPoint(self.x[self.pairs[0, i], :])
1026 points.InsertNextPoint(self.x[self.pairs[1, i], :])
1027 line = vtk.vtkLine()
1028 line.GetPointIds().SetId(0, 2*i) # index of particle 1
1029 line.GetPointIds().SetId(1, 2*i + 1) # index of particle 2
1030 lines.InsertNextCell(line)
1031 #colors.SetTupleValue(i, [100, 100, 100])
1032 forces.SetValue(i, self.f_n_magn[i])
1033 stresses.SetValue(i, self.sigma_contacts[i])
1034
1035 # initalize VTK data structure
1036 polydata = vtk.vtkPolyData()
1037
1038 polydata.SetPoints(points)
1039 polydata.SetLines(lines)
1040 #polydata.GetCellData().SetScalars(colors)
1041 #polydata.GetCellData().SetScalars(forces) # default scalar
1042 polydata.GetCellData().SetScalars(forces) # default scalar
1043 #polydata.GetCellData().AddArray(forces)
1044 polydata.GetCellData().AddArray(stresses)
1045 #polydata.GetPointData().AddArray(stresses)
1046 #polydata.GetPointData().SetScalars(stresses) # default scalar
1047
1048 # write VTK XML image data file
1049 writer = vtk.vtkXMLPolyDataWriter()
1050 writer.SetFileName(filename)
1051 if vtk.VTK_MAJOR_VERSION <= 5:
1052 writer.SetInput(polydata)
1053 else:
1054 writer.SetInputData(polydata)
1055 writer.Write()
1056 #writer.Update()
1057 if verbose:
1058 print('Output file: ' + filename)
1059
1060 def writeFluidVTK(self, folder='../output/', cell_centered=True,
1061 verbose=True):
1062 '''
1063 Writes a VTK file for the fluid grid to the ``../output/`` folder by
1064 default. The file name will be in the format ``fluid-<self.sid>.vti``.
1065 The vti files can be used for visualizing the fluid in ParaView.
1066
1067 The scalars (pressure, porosity, porosity change) and the velocity
1068 vectors are either placed in a grid where the grid corners correspond to
1069 the computational grid center (cell_centered=False). This results in a
1070 grid that doesn't appears to span the simulation domain, and values are
1071 smoothly interpolated on the cell faces. Alternatively, the
1072 visualization grid is equal to the computational grid, and cells face
1073 colors are not interpolated (cell_centered=True, default behavior).
1074
1075 The fluid grid is visualized by opening the vti files, and pressing
1076 "Apply" to import all fluid field properties. To visualize the scalar
1077 fields, such as the pressure, the porosity, the porosity change or the
1078 velocity magnitude, choose "Surface" or "Surface With Edges" as the
1079 "Representation". Choose the desired property as the "Coloring" field.
1080 It may be desirable to show the color bar by pressing the "Show" button,
1081 and "Rescale" to fit the color range limits to the current file. The
1082 coordinate system can be displayed by checking the "Show Axis" field.
1083 All adjustments by default require the "Apply" button to be pressed
1084 before regenerating the view.
1085
1086 The fluid vector fields (e.g. the fluid velocity) can be visualizing by
1087 e.g. arrows. To do this, select the fluid data in the "Pipeline
1088 Browser". Press "Glyph" from the "Common" toolbar, or go to the
1089 "Filters" mennu, and press "Glyph" from the "Common" list. Make sure
1090 that "Arrow" is selected as the "Glyph type", and "Velocity" as the
1091 "Vectors" value. Adjust the "Maximum Number of Points" to be at least as
1092 big as the number of fluid cells in the grid. Press "Apply" to visualize
1093 the arrows.
1094
1095 To visualize the cell-centered data with smooth interpolation, and in
1096 order to visualize fluid vector fields, the cell-centered mesh is
1097 selected in the "Pipeline Browser", and is filtered using "Filters" ->
1098 "Alphabetical" -> "Cell Data to Point Data".
1099
1100 If several data files are generated for the same simulation (e.g. using
1101 the :func:`writeVTKall()` function), it is able to step the
1102 visualization through time by using the ParaView controls.
1103
1104 :param folder: The folder where to place the output binary file (default
1105 (default='../output/')
1106 :type folder: str
1107 :param cell_centered: put scalars and vectors at cell centers (True) or
1108 cell corners (False), (default=True)
1109 :type cell_centered: bool
1110 :param verbose: Show diagnostic information (default=True)
1111 :type verbose: bool
1112 '''
1113 if not py_vtk:
1114 print('Error: vtk module not found, cannot writeFluidVTK.')
1115 return
1116
1117 filename = folder + '/fluid-' + self.sid + '.vti' # image grid
1118
1119 # initalize VTK data structure
1120 grid = vtk.vtkImageData()
1121 dx = (self.L-self.origo)/self.num # cell center spacing
1122 if cell_centered:
1123 grid.SetOrigin(self.origo)
1124 else:
1125 grid.SetOrigin(self.origo + 0.5*dx)
1126 grid.SetSpacing(dx)
1127 if cell_centered:
1128 grid.SetDimensions(self.num + 1) # no. of points in each direction
1129 else:
1130 grid.SetDimensions(self.num) # no. of points in each direction
1131
1132 # array of scalars: hydraulic pressures
1133 pres = vtk.vtkDoubleArray()
1134 pres.SetName("Pressure [Pa]")
1135 pres.SetNumberOfComponents(1)
1136 if cell_centered:
1137 pres.SetNumberOfTuples(grid.GetNumberOfCells())
1138 else:
1139 pres.SetNumberOfTuples(grid.GetNumberOfPoints())
1140
1141 # array of vectors: hydraulic velocities
1142 vel = vtk.vtkDoubleArray()
1143 vel.SetName("Velocity [m/s]")
1144 vel.SetNumberOfComponents(3)
1145 if cell_centered:
1146 vel.SetNumberOfTuples(grid.GetNumberOfCells())
1147 else:
1148 vel.SetNumberOfTuples(grid.GetNumberOfPoints())
1149
1150 # array of scalars: porosities
1151 poros = vtk.vtkDoubleArray()
1152 poros.SetName("Porosity [-]")
1153 poros.SetNumberOfComponents(1)
1154 if cell_centered:
1155 poros.SetNumberOfTuples(grid.GetNumberOfCells())
1156 else:
1157 poros.SetNumberOfTuples(grid.GetNumberOfPoints())
1158
1159 # array of scalars: porosity change
1160 dporos = vtk.vtkDoubleArray()
1161 dporos.SetName("Porosity change [1/s]")
1162 dporos.SetNumberOfComponents(1)
1163 if cell_centered:
1164 dporos.SetNumberOfTuples(grid.GetNumberOfCells())
1165 else:
1166 dporos.SetNumberOfTuples(grid.GetNumberOfPoints())
1167
1168 # array of scalars: Reynold's number
1169 Re_values = self.ReynoldsNumber()
1170 Re = vtk.vtkDoubleArray()
1171 Re.SetName("Reynolds number [-]")
1172 Re.SetNumberOfComponents(1)
1173 if cell_centered:
1174 Re.SetNumberOfTuples(grid.GetNumberOfCells())
1175 else:
1176 Re.SetNumberOfTuples(grid.GetNumberOfPoints())
1177
1178 # Find permeabilities if the Darcy solver is used
1179 if self.cfd_solver[0] == 1:
1180 self.findPermeabilities()
1181 k = vtk.vtkDoubleArray()
1182 k.SetName("Permeability [m*m]")
1183 k.SetNumberOfComponents(1)
1184 if cell_centered:
1185 k.SetNumberOfTuples(grid.GetNumberOfCells())
1186 else:
1187 k.SetNumberOfTuples(grid.GetNumberOfPoints())
1188
1189 self.findHydraulicConductivities()
1190 K = vtk.vtkDoubleArray()
1191 K.SetName("Conductivity [m/s]")
1192 K.SetNumberOfComponents(1)
1193 if cell_centered:
1194 K.SetNumberOfTuples(grid.GetNumberOfCells())
1195 else:
1196 K.SetNumberOfTuples(grid.GetNumberOfPoints())
1197
1198 p_f_constant = vtk.vtkDoubleArray()
1199 p_f_constant.SetName("Constant pressure [-]")
1200 p_f_constant.SetNumberOfComponents(1)
1201 if cell_centered:
1202 p_f_constant.SetNumberOfTuples(grid.GetNumberOfCells())
1203 else:
1204 p_f_constant.SetNumberOfTuples(grid.GetNumberOfPoints())
1205
1206 # insert values
1207 for z in range(self.num[2]):
1208 for y in range(self.num[1]):
1209 for x in range(self.num[0]):
1210 idx = x + self.num[0]*y + self.num[0]*self.num[1]*z
1211 pres.SetValue(idx, self.p_f[x, y, z])
1212 vel.SetTuple(idx, self.v_f[x, y, z, :])
1213 poros.SetValue(idx, self.phi[x, y, z])
1214 dporos.SetValue(idx, self.dphi[x, y, z])
1215 Re.SetValue(idx, Re_values[x, y, z])
1216 if self.cfd_solver[0] == 1:
1217 k.SetValue(idx, self.k[x, y, z])
1218 K.SetValue(idx, self.K[x, y, z])
1219 p_f_constant.SetValue(idx, self.p_f_constant[x, y, z])
1220
1221 # add pres array to grid
1222 if cell_centered:
1223 grid.GetCellData().AddArray(pres)
1224 grid.GetCellData().AddArray(vel)
1225 grid.GetCellData().AddArray(poros)
1226 grid.GetCellData().AddArray(dporos)
1227 grid.GetCellData().AddArray(Re)
1228 if self.cfd_solver[0] == 1:
1229 grid.GetCellData().AddArray(k)
1230 grid.GetCellData().AddArray(K)
1231 grid.GetCellData().AddArray(p_f_constant)
1232 else:
1233 grid.GetPointData().AddArray(pres)
1234 grid.GetPointData().AddArray(vel)
1235 grid.GetPointData().AddArray(poros)
1236 grid.GetPointData().AddArray(dporos)
1237 grid.GetPointData().AddArray(Re)
1238 if self.cfd_solver[0] == 1:
1239 grid.GetPointData().AddArray(k)
1240 grid.GetPointData().AddArray(K)
1241 grid.GetPointData().AddArray(p_f_constant)
1242
1243 # write VTK XML image data file
1244 writer = vtk.vtkXMLImageDataWriter()
1245 writer.SetFileName(filename)
1246 #writer.SetInput(grid) # deprecated from VTK 6
1247 writer.SetInputData(grid)
1248 writer.Update()
1249 if verbose:
1250 print('Output file: ' + filename)
1251
1252 def show(self, coloring=numpy.array([]), resolution=6):
1253 '''
1254 Show a rendering of all particles in a window.
1255
1256 :param coloring: Color the particles from red to white to blue according
1257 to the values in this array.
1258 :type coloring: numpy.array
1259 :param resolution: The resolution of the rendered spheres. Larger values
1260 increase the performance requirements.
1261 :type resolution: int
1262 '''
1263
1264 if not py_vtk:
1265 print('Error: vtk module not found, cannot show scene.')
1266 return
1267
1268 # create a rendering window and renderer
1269 ren = vtk.vtkRenderer()
1270 renWin = vtk.vtkRenderWindow()
1271 renWin.AddRenderer(ren)
1272
1273 # create a renderwindowinteractor
1274 iren = vtk.vtkRenderWindowInteractor()
1275 iren.SetRenderWindow(renWin)
1276
1277 if coloring.any():
1278 #min_value = numpy.min(coloring)
1279 max_value = numpy.max(coloring)
1280 #min_rgb = numpy.array([50, 50, 50])
1281 #max_rgb = numpy.array([255, 255, 255])
1282 #def color(value):
1283 #return (max_rgb - min_rgb) * (value - min_value)
1284
1285 def red(ratio):
1286 return numpy.fmin(1.0, 0.209*ratio**3. - 2.49*ratio**2. + 3.0*ratio
1287 + 0.0109)
1288 def green(ratio):
1289 return numpy.fmin(1.0, -2.44*ratio**2. + 2.15*ratio + 0.369)
1290 def blue(ratio):
1291 return numpy.fmin(1.0, -2.21*ratio**2. + 1.61*ratio + 0.573)
1292
1293 for i in numpy.arange(self.np):
1294
1295 # create source
1296 source = vtk.vtkSphereSource()
1297 source.SetCenter(self.x[i, :])
1298 source.SetRadius(self.radius[i])
1299 source.SetThetaResolution(resolution)
1300 source.SetPhiResolution(resolution)
1301
1302 # mapper
1303 mapper = vtk.vtkPolyDataMapper()
1304 if vtk.VTK_MAJOR_VERSION <= 5:
1305 mapper.SetInput(source.GetOutput())
1306 else:
1307 mapper.SetInputConnection(source.GetOutputPort())
1308
1309 # actor
1310 actor = vtk.vtkActor()
1311 actor.SetMapper(mapper)
1312
1313 # color
1314 if coloring.any():
1315 ratio = coloring[i]/max_value
1316 r, g, b = red(ratio), green(ratio), blue(ratio)
1317 actor.GetProperty().SetColor(r, g, b)
1318
1319 # assign actor to the renderer
1320 ren.AddActor(actor)
1321
1322 ren.SetBackground(0.3, 0.3, 0.3)
1323
1324 # enable user interface interactor
1325 iren.Initialize()
1326 renWin.Render()
1327 iren.Start()
1328
1329 def readfirst(self, verbose=True):
1330 '''
1331 Read the first output file from the ``../output/`` folder, corresponding
1332 to the object simulation id (``self.sid``).
1333
1334 :param verbose: Display diagnostic information (default=True)
1335 :type verbose: bool
1336
1337 See also :func:`readbin()`, :func:`readlast()`, :func:`readsecond`, and
1338 :func:`readstep`.
1339 '''
1340
1341 fn = '../output/' + self.sid + '.output00000.bin'
1342 self.readbin(fn, verbose)
1343
1344 def readsecond(self, verbose=True):
1345 '''
1346 Read the second output file from the ``../output/`` folder,
1347 corresponding to the object simulation id (``self.sid``).
1348
1349 :param verbose: Display diagnostic information (default=True)
1350 :type verbose: bool
1351
1352 See also :func:`readbin()`, :func:`readfirst()`, :func:`readlast()`,
1353 and :func:`readstep`.
1354 '''
1355 fn = '../output/' + self.sid + '.output00001.bin'
1356 self.readbin(fn, verbose)
1357
1358 def readstep(self, step, verbose=True):
1359 '''
1360 Read a output file from the ``../output/`` folder, corresponding
1361 to the object simulation id (``self.sid``).
1362
1363 :param step: The output file number to read, starting from 0.
1364 :type step: int
1365 :param verbose: Display diagnostic information (default=True)
1366 :type verbose: bool
1367
1368 See also :func:`readbin()`, :func:`readfirst()`, :func:`readlast()`,
1369 and :func:`readsecond`.
1370 '''
1371 fn = "../output/{0}.output{1:0=5}.bin".format(self.sid, step)
1372 self.readbin(fn, verbose)
1373
1374 def readlast(self, verbose=True):
1375 '''
1376 Read the last output file from the ``../output/`` folder, corresponding
1377 to the object simulation id (``self.sid``).
1378
1379 :param verbose: Display diagnostic information (default=True)
1380 :type verbose: bool
1381
1382 See also :func:`readbin()`, :func:`readfirst()`, :func:`readsecond`, and
1383 :func:`readstep`.
1384 '''
1385 lastfile = self.status()
1386 fn = "../output/{0}.output{1:0=5}.bin".format(self.sid, lastfile)
1387 self.readbin(fn, verbose)
1388
1389 def readTime(self, time, verbose=True):
1390 '''
1391 Read the output file most closely corresponding to the time given as an
1392 argument.
1393
1394 :param time: The desired current time [s]
1395 :type time: float
1396
1397 See also :func:`readbin()`, :func:`readfirst()`, :func:`readsecond`, and
1398 :func:`readstep`.
1399 '''
1400
1401 self.readfirst(verbose=False)
1402 t_first = self.currentTime()
1403 n_first = self.time_step_count[0]
1404
1405 self.readlast(verbose=False)
1406 t_last = self.currentTime()
1407 n_last = self.time_step_count[0]
1408
1409 if time < t_first or time > t_last:
1410 raise Exception('Error: The specified time {} s is outside the ' +
1411 'range of output files [{}; {}] s.'
1412 .format(time, t_first, t_last))
1413
1414 dt_dn = (t_last - t_first)/(n_last - n_first)
1415 step = int((time - t_first)/dt_dn) + n_first + 1
1416 self.readstep(step, verbose=verbose)