Updated Mesh Input And Output (markdown)
[dealii.wiki.git] / Mesh-Input-And-Output.md
1 # A little script to illustrate exchanging data with Cubit
2
3 ### Cubit
4
5 The following cubit journal file (actually a python script) exports the current mesh and 
6 boundary conditions from cubit to a file "output.ucd" in the current directory.
7 The [output file format](http://www.csit.fsu.edu/~burkardt/data/ucd/ucd.html) is [`AVS UCD`](http://vis.lbl.gov/NERSC/Software/express/help6.2/help/reference/dvmac/UCD_Form.htm).
8
9 This is a modification of the original script that takes into account boundary ids via
10 sidesets ids. If you want to save the boundary faces as well, you just need to 
11 add the relevant surfaces in 3d or curves in 2d to a sideset, and the id will be the one of the
12 sideset. 
13
14 ```python
15 #!python
16 # This script will output whatever mesh you have currently in CUBIT
17 # in the AVS UCD format. (http://www.csit.fsu.edu/~burkardt/data/ucd/ucd.html)
18 # You may need to tweak it to get exactly what you want out of
19 # it.  Your mileage may vary.
20
21 # set the filename -- you may need the entire path
22 outucdfile = "output.ucd"
23
24 outfile = open(outucdfile,"w")
25
26 cubit.cmd("body all rotate -90 about y")
27 cubit.cmd("body all reflect x")
28
29 # ============================================================
30 # Collect all the nodes
31 # ============================================================
32 group_id = cubit.get_id_from_name("temp_bc_curves")
33 if group_id != 0:
34   cubit.cmd("delete group " + str(group_id))
35 cubit.cmd("group 'temp_nodes' add node all")
36 group_id = cubit.get_id_from_name("temp_nodes")
37 node_list = cubit.get_group_nodes(group_id)
38 cubit.cmd("delete group " + str(group_id))
39 n_nodes = len(node_list)
40
41 # ============================================================
42 # Collect all the hex
43 # ============================================================
44 group_id = cubit.get_id_from_name("temp_hexes")
45 if group_id != 0:
46   cubit.cmd("delete group " + str(group_id))
47 cubit.cmd("group 'temp_hexes' add hex all")
48 group_id = cubit.get_id_from_name("temp_hexes")
49 hex_list = cubit.get_group_hexes(group_id)
50 cubit.cmd("delete group " + str(group_id))
51 n_hex_cells = len(hex_list)
52
53
54 # ============================================================ 
55 # Now the boundary conditions in 3d
56 # ============================================================
57 bc_surfaces = {}
58 n_bc_quads = 0
59
60 bc_ids = cubit.get_sideset_id_list()
61 for bc_id in bc_ids :
62   bc_surfaces[= cubit.get_sideset_surfaces(bc_id)
63   for bc_surface in  bc_surfaces[bc_id](bc_id]):
64     bc_quads = cubit.get_surface_quads(bc_surface)
65     n_bc_quads += len(bc_quads)
66
67
68 # ============================================================ 
69 # Collect all the surfaces. Notice that the surfaces that make up a
70 # volume are not grouped here. This is only for 2d objects, i.e.,
71 # when the number n_hex_cells is zero.
72 # ============================================================
73 surface_list = ()
74 quad_cell_list = {}
75 n_quad_cells = 0
76 if n_hex_cells == 0:
77   group_id = cubit.get_id_from_name("temp_surfs")
78   if group_id != 0:
79     cubit.cmd("delete group " + str(group_id))
80   cubit.cmd("group 'temp_surfs' add surf all")
81   group_id = cubit.get_id_from_name("temp_surfs")
82   surface_list = cubit.get_group_surfaces(group_id)
83   cubit.cmd("delete group " + str(group_id))
84   for surface_id in surface_list:
85     quad_cell_list[= cubit.get_surface_quads(surface_id)
86     n_quad_cells +=  len(quad_cell_list[surface_id](surface_id]))
87
88 # ============================================================ 
89 # Now the boundary conditions in 2d
90 # ============================================================
91 bc_curves = {}
92 bc_edges = {}
93 n_bc_edges = 0
94 if n_hex_cells == 0:
95   bc_ids = cubit.get_sideset_id_list()
96   for bc_id in bc_ids :
97     bc_curves[= cubit.get_sideset_curves(bc_id)
98     group_id = cubit.get_id_from_name("temp_bc_curves")
99     if group_id != 0:
100       cubit.cmd("delete group " + str(group_id))
101     for bc_curve in bc_curves[bc_id](bc_id]):
102       cubit.cmd("group 'temp_bc_curves' add edge all in curve " + str(bc_curve))
103     group_id = cubit.get_id_from_name("temp_bc_curves")
104     bc_edges[= cubit.get_group_edges(group_id)
105     cubit.cmd("delete group " + str(group_id))
106     n_bc_edges +=  len(bc_edges[bc_id](bc_id]))
107
108 print 'Edges: ' + str(n_bc_edges)
109
110 # ============================================================
111 # Now we write the header.
112 # ============================================================
113 n_elements = n_hex_cells + n_bc_quads + n_quad_cells + n_bc_edges
114 outfile.write(str(n_nodes) + " " + str(n_elements) + " 0 0 0\n")
115
116
117 # ============================================================
118 # The node list.
119 # ============================================================
120 for node_num in node_list:
121    outfile.write(str(node_num ))
122    outfile.write("\t")
123    node_coord = cubit.get_nodal_coordinates(node_num)
124    if abs(node_coord[outfile.write("0 ")
125    else :
126        outfile.write(str(node_coord[2](2])<1e-15:)) + " ")
127    if abs(node_coord[outfile.write("0 ")
128    else :
129        outfile.write(str(node_coord[1](1])<1e-15:)) + " ")
130    if abs(node_coord[outfile.write("0")
131    else :
132        outfile.write(str(node_coord[0](0])<1e-15:)))
133    outfile.write("\n")
134
135
136 # ============================================================
137 # The hex list. 3d
138 # ============================================================
139 k = 1
140 for hex_num in hex_list:
141    outfile.write(str(k) + " 0 " + " hex  ")
142    k += 1
143    hex_nodes = cubit.get_connectivity("Hex", hex_num)
144    i = 0
145    while i < 8:
146       outfile.write(str(hex_nodes[+ " ")
147       i += 1
148    outfile.write("\n")
149
150
151 # ============================================================
152 # The quads on the boundaries. 3d
153 # Note that the boundary id is given by the sideset id.
154 # ============================================================
155 if n_hex_cells != 0:
156   k=1
157   for bc_id in bc_ids:
158     for bc_surface in  bc_surfaces[bc_id](i])):
159       bc_quads = cubit.get_surface_quads(bc_surface)
160       for quad_num in bc_quads:
161         outfile.write(str(k))
162         k += 1
163         outfile.write(" " + str(bc_id) + " quad ")
164         quad_nodes = cubit.get_connectivity("Quad", quad_num)
165         outfile.write(str(quad_nodes[+ " ")
166         outfile.write(str(quad_nodes[3](0]))) + " ")
167         outfile.write(str(quad_nodes[+ " ")
168         outfile.write(str(quad_nodes[1](2]))))
169         outfile.write("\n")
170
171 # ============================================================
172 # The quads on the surfaces. 2d
173 # ============================================================
174 if n_hex_cells == 0:
175   k = 1
176   for surface_id in surface_list:
177     for quad_num in quad_cell_list[outfile.write(str(k))
178       k += 1
179       outfile.write(" " + str(surface_id) + " quad ")
180       quad_nodes = cubit.get_connectivity("Quad", quad_num)
181       outfile.write(str(quad_nodes[0](surface_id]:)) + " ")
182       outfile.write(str(quad_nodes[+ " ")
183       outfile.write(str(quad_nodes[2](3]))) + " ")
184       outfile.write(str(quad_nodes[outfile.write("\n")
185
186 # ============================================================
187 # The edges on the curves. 2d
188 # Boundary id = sideset_id
189 # ============================================================
190 if n_hex_cells == 0:
191   k=1
192   for bc_id in bc_ids:
193     for edge_num in bc_edges[bc_id](1]))):
194       outfile.write(str(k))
195       k += 1
196       outfile.write(" " + str(bc_id) + " line ")
197       edge_nodes = cubit.get_connectivity("Edge", edge_num)
198       outfile.write(str(edge_nodes[+ " ")
199       outfile.write(str(edge_nodes[1](0]))))
200       outfile.write("\n")
201
202 outfile.close()
203
204 cubit.cmd("body all reflect x")
205 cubit.cmd("body all rotate 90 about y")
206
207 print str(n_nodes) + " nodes\n"
208 print str(n_hex_cells) + " hexes\n"
209 print str(n_quad_cells) + " quads\n"
210 print str(n_bc_quads) + " face_quads\n"
211 print str(n_bc_edges) + " edges\n"
212
213 ```
214
215 # Alternate versions of this script
216
217 ### Cubit 15.2
218
219 The following script, which also outputs material IDs (as identified by the `block id` to which each element belongs), has been [verified to work](https://groups.google.com/forum/#!topic/dealii/ZnU3XIK_ipk) with `Cubit v15.2`. 
220
221 ```python
222 #!python
223 # This script will output whatever mesh you have currently in CUBIT
224 # in the AVS UCD format. (http://www.csit.fsu.edu/~burkardt/data/ucd/ucd.html)
225 # You may need to tweak it to get exactly what you want out of
226 # it.  Your mileage may vary.
227
228 # set the filename -- you may need the entire path
229 outucdfile = "output.ucd"
230
231 outfile = open(outucdfile,"w")
232
233 cubit.cmd("body all rotate -90 about y")
234 cubit.cmd("body all reflect x")
235
236 # ============================================================
237 # Collect all the nodes
238 # ============================================================
239 group_id = cubit.get_id_from_name("temp_bc_curves")
240 if group_id != 0:
241   cubit.cmd("delete group " + str(group_id))
242 cubit.cmd("group 'temp_nodes' add node all")
243 group_id = cubit.get_id_from_name("temp_nodes")
244 node_list = cubit.get_group_nodes(group_id)
245 cubit.cmd("delete group " + str(group_id))
246 n_nodes = len(node_list)
247
248 # ============================================================
249 # Collect all the hex
250 # ============================================================
251 hex_cell_list = {}
252 n_hex_cells = 0
253
254 block_ids = cubit.get_block_id_list()
255 for block_id in block_ids:
256     hex_cell_list[block_id] = cubit.get_block_hexes(block_id)
257     n_hex_cells += len(hex_cell_list[block_id])
258
259 # ============================================================
260 # Now the boundary conditions in 3d
261 # ============================================================
262 bc_surfaces = {}
263 n_bc_quads = 0
264
265 bc_ids = cubit.get_sideset_id_list()
266 for bc_id in bc_ids :
267   bc_surfaces[bc_id] = cubit.get_sideset_surfaces(bc_id)
268   for bc_surface in  bc_surfaces[bc_id]:
269     bc_quads = cubit.get_surface_quads(bc_surface)
270     n_bc_quads += len(bc_quads)
271
272
273 # ============================================================
274 # Collect all the surfaces. Notice that the surfaces that make up a
275 # volume are not grouped here. This is only for 2d objects, i.e.,
276 # when the number n_hex_cells is zero.
277 # ============================================================
278 surface_list = ()
279 quad_cell_list = {}
280 n_quad_cells = 0
281 if n_hex_cells == 0:
282   group_id = cubit.get_id_from_name("temp_surfs")
283   if group_id != 0:
284     cubit.cmd("delete group " + str(group_id))
285   cubit.cmd("group 'temp_surfs' add surf all")
286   group_id = cubit.get_id_from_name("temp_surfs")
287   surface_list = cubit.get_group_surfaces(group_id)
288   cubit.cmd("delete group " + str(group_id))
289   for surface_id in surface_list:
290     quad_cell_list[surface_id] = cubit.get_surface_quads(surface_id)
291     n_quad_cells +=  len(quad_cell_list[surface_id])
292
293 # ============================================================
294 # Now the boundary conditions in 2d
295 # ============================================================
296 bc_curves = {}
297 bc_edges = {}
298 n_bc_edges = 0
299 if n_hex_cells == 0:
300   bc_ids = cubit.get_sideset_id_list()
301   for bc_id in bc_ids :
302     bc_curves[bc_id]    = cubit.get_sideset_curves(bc_id)
303     group_id = cubit.get_id_from_name("temp_bc_curves")
304     if group_id != 0:
305       cubit.cmd("delete group " + str(group_id))
306     for bc_curve in bc_curves[bc_id]:
307       cubit.cmd("group 'temp_bc_curves' add edge all in curve " + str(bc_curve))
308     group_id = cubit.get_id_from_name("temp_bc_curves")
309     bc_edges[bc_id] = cubit.get_group_edges(group_id)
310     cubit.cmd("delete group " + str(group_id))
311     n_bc_edges +=  len(bc_edges[bc_id])
312
313 print 'Edges: ' + str(n_bc_edges)
314
315 # ============================================================
316 # Now we write the header.
317 # ============================================================
318 n_elements = n_hex_cells + n_bc_quads + n_quad_cells + n_bc_edges
319 outfile.write(str(n_nodes) + " " + str(n_elements) + " 0 0 0\n")
320
321
322 # ============================================================
323 # The node list.
324 # ============================================================
325 for node_num in node_list:
326    outfile.write(str(node_num ))
327    outfile.write("\t")
328    node_coord = cubit.get_nodal_coordinates(node_num)
329    if abs(node_coord[2])<1e-15:
330        outfile.write("0 ")
331    else :
332        outfile.write(str(node_coord[2]) + " ")
333    if abs(node_coord[1])<1e-15:
334        outfile.write("0 ")
335    else :
336        outfile.write(str(node_coord[1]) + " ")
337    if abs(node_coord[0])<1e-15:
338        outfile.write("0")
339    else :
340        outfile.write(str(node_coord[0]))
341    outfile.write("\n")
342
343
344 # ============================================================
345 # The hex list. 3d
346 # ============================================================
347 k = 1
348 for block_id in block_ids:
349     for hex_num in hex_cell_list[block_id]:
350         print str(k) + " " + str(block_id) + "  hex  "
351         outfile.write(str(k) + " " + str(block_id) + "  hex  ")
352         k += 1
353         hex_nodes = cubit.get_connectivity("Hex", hex_num)
354         i = 0
355         while i < 8:
356             outfile.write(str(hex_nodes[i]) + " ")
357             i += 1
358         outfile.write("\n")
359
360 # ============================================================
361 # The quads on the boundaries. 3d
362 # Note that the boundary id is given by the sideset id.
363 # ============================================================
364 if n_hex_cells != 0:
365   k=1
366   for bc_id in bc_ids:
367     for bc_surface in  bc_surfaces[bc_id]:
368       bc_quads = cubit.get_surface_quads(bc_surface)
369       for quad_num in bc_quads:
370         outfile.write(str(k))
371         k += 1
372         outfile.write(" " + str(bc_id) + " quad ")
373         quad_nodes = cubit.get_connectivity("Quad", quad_num)
374         outfile.write(str(quad_nodes[0]) + " ")
375         outfile.write(str(quad_nodes[3]) + " ")
376         outfile.write(str(quad_nodes[2]) + " ")
377         outfile.write(str(quad_nodes[1]))
378         outfile.write("\n")
379
380 # ============================================================
381 # The quads on the surfaces. 2d
382 # ============================================================
383 if n_hex_cells == 0:
384   k = 1
385   for surface_id in surface_list:
386     for quad_num in quad_cell_list[surface_id]:
387       outfile.write(str(k))
388       k += 1
389       outfile.write(" " + str(surface_id) + " quad ")
390       quad_nodes = cubit.get_connectivity("Quad", quad_num)
391       outfile.write(str(quad_nodes[0]) + " ")
392       outfile.write(str(quad_nodes[3]) + " ")
393       outfile.write(str(quad_nodes[2]) + " ")
394       outfile.write(str(quad_nodes[1]))
395       outfile.write("\n")
396
397 # ============================================================
398 # The edges on the curves. 2d
399 # Boundary id = sideset_id
400 # ============================================================
401 if n_hex_cells == 0:
402   k=1
403   for bc_id in bc_ids:
404     for edge_num in bc_edges[bc_id]:
405       outfile.write(str(k))
406       k += 1
407       outfile.write(" " + str(bc_id) + " line ")
408       edge_nodes = cubit.get_connectivity("Edge", edge_num)
409       outfile.write(str(edge_nodes[0]) + " ")
410       outfile.write(str(edge_nodes[1]))
411       outfile.write("\n")
412
413 outfile.close()
414
415 cubit.cmd("body all reflect x")
416 cubit.cmd("body all rotate 90 about y")
417
418 print str(n_nodes) + " nodes\n"
419 print str(n_hex_cells) + " hexes\n"
420 print str(n_quad_cells) + " quads\n"
421 print str(n_bc_quads) + " face_quads\n"
422 print str(n_bc_edges) + " edges\n"
423 ```
424
425 ### Cubit 13.2
426 In case it is of any use, one of the initial versions of this script is also provided below. This has been verified to work on `Cubit v13.2`, but unfortunately does not output material IDs for reasons documented in the paragraphs below.
427
428 ```python
429 #!python
430 # This script will output whatever mesh you have currently in CUBIT
431 # in the AVS UCD format. (http://www.csit.fsu.edu/~burkardt/data/ucd/ucd.html)
432 # You may need to tweak it to get exactly what you want out of
433 # it.  Your mileage may vary.
434
435 # set the filename -- you may need the entire path
436 outucdfile = "output.ucd"
437
438 outfile = open(outucdfile,"w")
439
440 cubit.cmd("body all rotate -90 about y")
441 cubit.cmd("body all reflect x")
442
443 # ============================================================
444 # Collect all the nodes
445 # ============================================================
446 group_id = cubit.get_id_from_name("temp_bc_curves")
447 if group_id != 0:
448   cubit.cmd("delete group " + str(group_id))
449 cubit.cmd("group 'temp_nodes' add node all")
450 group_id = cubit.get_id_from_name("temp_nodes")
451 node_list = cubit.get_group_nodes(group_id)
452 cubit.cmd("delete group " + str(group_id))
453 n_nodes = len(node_list)
454
455 # ============================================================
456 # Collect all the hex
457 # ============================================================
458 group_id = cubit.get_id_from_name("temp_hexes")
459 if group_id != 0:
460   cubit.cmd("delete group " + str(group_id))
461 cubit.cmd("group 'temp_hexes' add hex all")
462 group_id = cubit.get_id_from_name("temp_hexes")
463 hex_list = cubit.get_group_hexes(group_id)
464 cubit.cmd("delete group " + str(group_id))
465 n_hex_cells = len(hex_list)
466
467
468 # ============================================================
469 # Now the boundary conditions in 3d
470 # ============================================================
471 bc_surfaces = {}
472 n_bc_quads = 0
473
474 bc_ids = cubit.get_sideset_id_list()
475 for bc_id in bc_ids :
476   bc_surfaces[bc_id]    = cubit.get_sideset_surfaces(bc_id)
477   for bc_surface in  bc_surfaces[bc_id]:
478     bc_quads = cubit.get_surface_quads(bc_surface)
479     n_bc_quads += len(bc_quads)
480
481
482 # ============================================================
483 # Collect all the surfaces. Notice that the surfaces that make up a
484 # volume are not grouped here. This is only for 2d objects, i.e.,
485 # when the number n_hex_cells is zero.
486 # ============================================================
487 surface_list = ()
488 quad_cell_list = {}
489 n_quad_cells = 0
490 if n_hex_cells == 0:
491   group_id = cubit.get_id_from_name("temp_surfs")
492   if group_id != 0:
493     cubit.cmd("delete group " + str(group_id))
494   cubit.cmd("group 'temp_surfs' add surf all")
495   group_id = cubit.get_id_from_name("temp_surfs")
496   surface_list = cubit.get_group_surfaces(group_id)
497   cubit.cmd("delete group " + str(group_id))
498   for surface_id in surface_list:
499     quad_cell_list[surface_id] = cubit.get_surface_quads(surface_id)
500     n_quad_cells +=  len(quad_cell_list[surface_id])
501
502 # ============================================================
503 # Now the boundary conditions in 2d
504 # ============================================================
505 bc_curves = {}
506 bc_edges = {}
507 n_bc_edges = 0
508 if n_hex_cells == 0:
509   bc_ids = cubit.get_sideset_id_list()
510   for bc_id in bc_ids :
511     bc_curves[bc_id]    = cubit.get_sideset_curves(bc_id)
512     group_id = cubit.get_id_from_name("temp_bc_curves")
513     if group_id != 0:
514       cubit.cmd("delete group " + str(group_id))
515     for bc_curve in bc_curves[bc_id]:
516       cubit.cmd("group 'temp_bc_curves' add edge all in curve " + str(bc_curve))
517     group_id = cubit.get_id_from_name("temp_bc_curves")
518     bc_edges[bc_id] = cubit.get_group_edges(group_id)
519     cubit.cmd("delete group " + str(group_id))
520     n_bc_edges +=  len(bc_edges[bc_id])
521
522 print 'Edges: ' + str(n_bc_edges)
523
524 # ============================================================
525 # Now we write the header.
526 # ============================================================
527 n_elements = n_hex_cells + n_bc_quads + n_quad_cells + n_bc_edges
528 outfile.write(str(n_nodes) + " " + str(n_elements) + " 0 0 0\n")
529
530
531 # ============================================================
532 # The node list.
533 # ============================================================
534 for node_num in node_list:
535    outfile.write(str(node_num ))
536    outfile.write("\t")
537    node_coord = cubit.get_nodal_coordinates(node_num)
538    if abs(node_coord[2])<1e-15:
539        outfile.write("0 ")
540    else :
541        outfile.write(str(node_coord[2]) + " ")
542    if abs(node_coord[1])<1e-15:
543        outfile.write("0 ")
544    else :
545        outfile.write(str(node_coord[1]) + " ")
546    if abs(node_coord[0])<1e-15:
547        outfile.write("0")
548    else :
549        outfile.write(str(node_coord[0]))
550    outfile.write("\n")
551
552
553 # ============================================================
554 # The hex list. 3d
555 # ============================================================
556 k = 1
557 for hex_num in hex_list:
558    outfile.write(str(k) + " 0 " + " hex  ")
559    k += 1
560    hex_nodes = cubit.get_connectivity("Hex", hex_num)
561    i = 0
562    while i < 8:
563       outfile.write(str(hex_nodes[i]) + " ")
564       i += 1
565    outfile.write("\n")
566
567 # ============================================================
568 # The quads on the boundaries. 3d
569 # Note that the boundary id is given by the sideset id.
570 # ============================================================
571 if n_hex_cells != 0:
572   k=1
573   for bc_id in bc_ids:
574     for bc_surface in  bc_surfaces[bc_id]:
575       bc_quads = cubit.get_surface_quads(bc_surface)
576       for quad_num in bc_quads:
577         outfile.write(str(k))
578         k += 1
579         outfile.write(" " + str(bc_id) + " quad ")
580         quad_nodes = cubit.get_connectivity("Quad", quad_num)
581         outfile.write(str(quad_nodes[0]) + " ")
582         outfile.write(str(quad_nodes[3]) + " ")
583         outfile.write(str(quad_nodes[2]) + " ")
584         outfile.write(str(quad_nodes[1]))
585         outfile.write("\n")
586
587 # ============================================================
588 # The quads on the surfaces. 2d
589 # ============================================================
590 if n_hex_cells == 0:
591   k = 1
592   for surface_id in surface_list:
593     for quad_num in quad_cell_list[surface_id]:
594       outfile.write(str(k))
595       k += 1
596       outfile.write(" " + str(surface_id) + " quad ")
597       quad_nodes = cubit.get_connectivity("Quad", quad_num)
598       outfile.write(str(quad_nodes[0]) + " ")
599       outfile.write(str(quad_nodes[3]) + " ")
600       outfile.write(str(quad_nodes[2]) + " ")
601       outfile.write(str(quad_nodes[1]))
602       outfile.write("\n")
603
604 # ============================================================
605 # The edges on the curves. 2d
606 # Boundary id = sideset_id
607 # ============================================================
608 if n_hex_cells == 0:
609   k=1
610   for bc_id in bc_ids:
611     for edge_num in bc_edges[bc_id]:
612       outfile.write(str(k))
613       k += 1
614       outfile.write(" " + str(bc_id) + " line ")
615       edge_nodes = cubit.get_connectivity("Edge", edge_num)
616       outfile.write(str(edge_nodes[0]) + " ")
617       outfile.write(str(edge_nodes[1]))
618       outfile.write("\n")
619
620 outfile.close()
621
622 cubit.cmd("body all reflect x")
623 cubit.cmd("body all rotate 90 about y")
624
625 print str(n_nodes) + " nodes\n"
626 print str(n_hex_cells) + " hexes\n"
627 print str(n_quad_cells) + " quads\n"
628 print str(n_bc_quads) + " face_quads\n"
629 print str(n_bc_edges) + " edges\n"
630 ```
631
632 ### A note on the (non-) output of material ID's
633 For the first and third versions of this script, the material IDs are not output (they are always set by default to zero).
634 This is because there appears to be a deficiency/bug in the implementation of `cubit.get_block_hexes(block_id)` (in `Cubit v13.2`, at least). 
635 In theory one should be able to collect the cells in the following manner, assuming that the `block_id` represents the material id (this mirrors what is done already for `sideset_id`s and boundary ids):
636
637 ```python
638 # ============================================================
639 # Collect all the hex
640 # ============================================================
641 hex_cell_list = {}
642 n_hex_cells = 0
643
644 block_ids = cubit.get_block_id_list()
645 for block_id in block_ids:
646     hex_cell_list[block_id] = cubit.get_block_hexes(block_id)
647     n_hex_cells += len(hex_cell_list[block_id])
648 ```
649
650 and one could then write out the elements to file with a set of commands along the lines of
651
652 ```python
653 # ============================================================
654 # The hex list. 3d
655 # ============================================================
656 k = 1
657 for block_id in block_ids:
658     for hex_num in hex_cell_list[block_id]:
659         print str(k) + " " + str(block_id) + "  hex  "
660         outfile.write(str(k) + " " + str(block_id) + "  hex  ")
661         k += 1
662         hex_nodes = cubit.get_connectivity("Hex", hex_num)
663         i = 0
664         while i < 8:
665             outfile.write(str(hex_nodes[i]) + " ")
666             i += 1
667         outfile.write("\n")
668 ```
669
670 However, there `cubit.get_block_hexes(block_id)` does not return anything other than an empty set when these sets of commands were last tested (`Cubit v13.2`).
671 By contrast the script verified to work with `Cubit v15.2` includes this logic as the bug has been corrected in later versions of the software. 
672
673 The following version of the above script, instead, produces a VTK file:
674 ```
675 #!python
676 #!python
677 # This script will output whatever mesh you have currently in CUBIT
678 # in the VTK ascii format
679 #
680 # In this script:
681 # - block id   => material id
682 # - sideset id => boundary id
683 # - nodeset id => manifold id
684 import cubit
685
686 vtkfile = '/tmp/output.vtk'
687 # def exportVTK(vtkfile):
688 if True:
689     outfile = open(vtkfile,"w")
690     # ============================================================
691     # Global dictionaries
692     # ============================================================
693     node_ids = {}
694     #
695     hex_material_ids = {}
696     hex_manifold_ids = {}
697     #
698     quad_material_ids = {}
699     quad_manifold_ids = {}
700     #
701     edge_material_ids = {}
702     edge_manifold_ids = {}
703     #
704     n_nodes = 0
705     n_hexes = 0
706     n_quads = 0
707     n_edges = 0
708     #
709     block_ids = cubit.get_block_id_list()
710     nodeset_ids = cubit.get_nodeset_id_list()
711     sideset_ids = cubit.get_sideset_id_list()
712     #
713     # ============================================================
714     # Collect nodes
715     # ============================================================
716     group_id = cubit.get_id_from_name("temp_nodes")
717     if group_id != 0:
718       cubit.cmd("delete group " + str(group_id))
719     cubit.cmd("group 'temp_nodes' add node all")
720     group_id = cubit.get_id_from_name("temp_nodes")
721     node_list = cubit.get_group_nodes(group_id)
722     cubit.cmd("delete group " + str(group_id))
723     n_nodes = len(node_list)
724     #
725     i = 0
726     for node in node_list:
727         node_ids[node] = i
728         i +=1
729     #
730     # ============================================================
731     # Collect all hexes for 3d grids
732     # set default material and manifold ids
733     # ============================================================
734     group_id = cubit.get_id_from_name("temp_hexes")
735     if group_id != 0:
736       cubit.cmd("delete group " + str(group_id))
737     cubit.cmd("group 'temp_hexes' add hex all")
738     group_id = cubit.get_id_from_name("temp_hexes")
739     hex_list = cubit.get_group_hexes(group_id)
740     cubit.cmd("delete group " + str(group_id))
741     n_hexes = len(hex_list)
742     #
743     for hex in hex_list:
744         hex_material_ids[hex] = 0
745         hex_manifold_ids[hex] = -1
746     #
747     # ============================================================
748     # Collect hexes by block id ( => Material ID)
749     # ============================================================
750     for block_id in block_ids:
751         volumes = cubit.get_block_volumes(block_id)
752         for vol in volumes:
753             hexes = cubit.get_volume_hexes(vol)
754             print len(hexes), 'hexes in volume', vol
755             for hex in hexes:
756                 hex_material_ids[hex] = block_id
757     #
758     # ============================================================
759     # Collect hexes by nodeset id ( => Manifold ID)
760     # ============================================================
761     for nodeset_id in nodeset_ids:
762         volumes = cubit.get_nodeset_volumes(nodeset_id)
763         for vol in volumes:
764             hexes = cubit.get_volume_hexes(vol)
765             for hex in hexes:
766                 hex_manifold_ids[hex] = nodeset_id
767     #
768     # ============================================================
769     # Collect all quads (for 2d grids)
770     # set default material and manifold ids
771     # ============================================================
772     if n_hexes == 0:
773       surface_list = ()
774       group_id = cubit.get_id_from_name("temp_surfs")
775       if group_id != 0:
776         cubit.cmd("delete group " + str(group_id))
777       cubit.cmd("group 'temp_surfs' add surf all")
778       group_id = cubit.get_id_from_name("temp_surfs")
779       surface_list = cubit.get_group_surfaces(group_id)
780       cubit.cmd("delete group " + str(group_id))
781       for surface_id in surface_list:
782         quads = cubit.get_surface_quads(surface_id)
783         for quad in quads:
784           quad_material_ids[quad] = 0
785           quad_manifold_ids[quad] = -1
786     #
787     n_quads = len(quad_material_ids)
788     #
789     # ============================================================
790     # Collect quads by block id for 2D grids ( => Material ID)
791     # and by sideset_id for 3d grids ( => Material ID == Boundary ID)
792     # ============================================================
793     if n_hexes == 0:
794         for block_id in block_ids:
795             surfaces = cubit.get_block_surfaces(block_id)
796             for surface_id in surfaces:
797                 quads = cubit.get_surface_quads(surface_id)
798                 quad_material_ids[quad] = block_id
799     else:
800         for sideset_id in sideset_ids:
801             surfaces = cubit.get_sideset_surfaces(sideset_id)
802             for surface_id in surfaces:
803                 quads = cubit.get_surface_quads(surface_id)
804                 for quad in quads:
805                     quad_material_ids[quad] = sideset_id
806                     quad_manifold_ids[quad] = -1
807     #
808     # ============================================================
809     # Collect quads by nodeset id ( => Manifold ID)
810     # ============================================================
811     for nodeset_id in nodeset_ids:
812         surfaces = cubit.get_nodeset_surfaces(nodeset_id)
813         for surface_id in surfaces:
814             quads = cubit.get_surface_quads(surface_id)
815             quad_manifold_ids[quad] = nodeset_id
816             if quad not in quad_material_ids:
817                 # this can only happen in 3d, so this must be an internal
818                 # face
819                 quad_material_ids[quad] = -1
820     #
821     n_quads = len(quad_material_ids)
822     #
823     # ============================================================
824     # Collect edges by sideset id ( => Material ID)
825     # ============================================================
826     for sideset_id in sideset_ids:
827         curves = cubit.get_sideset_curves(sideset_id)
828         for curve in curves:
829             # edges = get_curve_edges(curve)
830             # the above does not work. Get around this bug
831             group_id = cubit.get_id_from_name("temp_edges")
832             if group_id != 0:
833                 cubit.cmd("delete group " + str(group_id))
834             cubit.cmd("group 'temp_edges' add edge all in curve " + str(curve))
835             group_id = cubit.get_id_from_name("temp_edges")
836             edges = cubit.get_group_edges(group_id)
837             cubit.cmd("delete group " + str(group_id))
838             for edge in edges:
839                 edge_material_ids[edge] = sideset_id
840                 edge_manifold_ids[edge] = -1
841     #
842     # ============================================================
843     # Collect edges by nodeset id ( => Manifold ID)
844     # ============================================================
845     for nodeset_id in nodeset_ids:
846         curves = cubit.get_nodeset_curves(nodeset_id)
847         for curve in curves:
848             # edges = get_curve_edges(curve)
849             # the above does not work. Get around this bug
850             group_id = cubit.get_id_from_name("temp_edges")
851             if group_id != 0:
852                 cubit.cmd("delete group " + str(group_id))
853             cubit.cmd("group 'temp_edges' add edge all in curve " + str(curve))
854             group_id = cubit.get_id_from_name("temp_edges")
855             edges = cubit.get_group_edges(group_id)
856             cubit.cmd("delete group " + str(group_id))
857             for edge in edges:
858                 edge_manifold_ids[edge] = nodeset_id
859                 if edge not in edge_material_ids:
860                     edge_material_ids[edge] = -1
861     #
862     n_edges = len(edge_material_ids)
863     #
864     # ============================================================
865     # Now we write the header.
866     # ============================================================
867     outfile.write("# vtk DataFile Version 3.0\n")
868     outfile.write("File generated with Trelis\nASCII\n")
869     outfile.write("DATASET UNSTRUCTURED_GRID\nPOINTS " + str(n_nodes) + " double\n")
870     #
871     n_elements = n_hexes + n_quads + n_edges
872     #
873     # ============================================================
874     # The node list.
875     # ============================================================
876     for node_num in node_list:
877        node_coord = cubit.get_nodal_coordinates(node_num)
878        if abs(node_coord[0])<1e-15:
879            outfile.write("0 ")
880        else :
881            outfile.write(str(node_coord[0]) + " ")
882        if abs(node_coord[1])<1e-15:
883            outfile.write("0 ")
884        else :
885            outfile.write(str(node_coord[1]) + " ")
886        if abs(node_coord[2])<1e-15:
887            outfile.write("0")
888        else :
889            outfile.write(str(node_coord[2]))
890        outfile.write("\n")
891     #
892     # ============================================================
893     # Now we write the for the CELLS.
894     # ============================================================
895     storage = n_hexes*9 + n_quads*5 + n_edges*3
896     outfile.write("CELLS " + str(n_elements) + " " + str(storage) + "\n")
897     #
898     # ============================================================
899     # The hex list
900     # ============================================================
901     for hex in hex_manifold_ids:
902        outfile.write("8")
903        n = cubit.get_connectivity("Hex", hex)
904        i = 0
905        while i<8:
906           outfile.write(" " + str(node_ids[n[i]]))
907           i += 1
908        outfile.write("\n")
909     #
910     #
911     # ============================================================
912     # The quad list
913     # ============================================================
914     for quad in quad_manifold_ids:
915        outfile.write("4")
916        n = cubit.get_connectivity("Quad", quad)
917        i = 0
918        while i<4:
919           outfile.write(" " + str(node_ids[n[i]]))
920           i += 1
921        outfile.write("\n")
922     #
923     # ============================================================
924     # The edge list
925     # ============================================================
926     for edge in edge_manifold_ids:
927        outfile.write("2")
928        n = cubit.get_connectivity("Edge", edge)
929        i = 0
930        while i<2:
931           outfile.write(" " + str(node_ids[n[i]]))
932           i += 1
933        outfile.write("\n")
934     #
935     # ============================================================
936     # The CELL_TYPES section
937     # ============================================================
938     cell_types = "12 "* n_hexes + "9 "*n_quads + "3 "*n_edges+"\n"
939     outfile.write("CELL_TYPES " + str(n_elements) + "\n" + cell_types)
940     #
941     # ============================================================
942     # The DATA_SET set
943     # ============================================================
944     outfile.write("CELL_DATA " + str(n_elements) +"\n")
945     outfile.write("SCALARS MaterialID int 1\nLOOKUP_TABLE default\n")
946     for hex in hex_material_ids:
947         outfile.write(str(hex_material_ids[hex]) + " ")
948     for quad in quad_material_ids:
949         outfile.write(str(quad_material_ids[quad]) + " ")
950     for edge in edge_material_ids:
951         outfile.write(str(edge_material_ids[edge]) + " ")
952     outfile.write("\n")
953     #
954     outfile.write("SCALARS ManifoldID int 1\nLOOKUP_TABLE default\n")
955     for hex in hex_manifold_ids:
956         outfile.write(str(hex_manifold_ids[hex]) + " ")
957     for quad in quad_manifold_ids:
958         outfile.write(str(quad_manifold_ids[quad]) + " ")
959     for edge in edge_manifold_ids:
960         outfile.write(str(edge_manifold_ids[edge]) + " ")
961     outfile.write("\n")
962     #
963     print n_nodes, "nodes"
964     print n_hexes, "hexes"
965     print n_quads, "quads"
966     print n_edges, "edges"
967 ```

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.