« back to all changes in this revision
Viewing changes to site-packages/dolfin_utils/meshconvert.py
- browse files at revision 5840
- compare with another revision
- download diff
- download tarball
- view history from revision 5840
- Committer: Anders Logg
- Date: 2011-05-01 22:01:45 UTC
- mfrom: (5839.1.29 dolfin-all)
- Revision ID: logg@simula.no-20110501220145-60ys4e77qbc0mq32
merge
- files added:
- site-packages/dolfin/compilemodules/compilemodule.py
- site-packages/dolfin/functions
- files removed:
- scripts/dolfin-convert/test_old_xml.xml
- site-packages/dolfin/function
- files modified:
- CMakeLists.txt
added
removed
site-packages/dolfin_utils/meshconvert.py
1
# Copyright (C) 2006 Anders Logg
2
# Licensed under the GNU LGPL Version 2.1
3
#
4
# Modified by Garth N. Wells (gmsh function)
5
# Modified by Alexander H. Jarosch (gmsh fix)
6
# Modified by Angelo Simone (Gmsh and Medit fix)
7
# Modified by Andy R. Terrel (gmsh fix and triangle function)
8
# Modified by Magnus Vikstrom (metis and scotch function)
9
# Modified by Bartosz Sawicki (diffpack function)
10
# Modified by Gideon Simpson (Exodus II function)
11
# Modified by Arve Knudsen (make into module, abaqus support)
12
# Modified by Kent-Andre Mardal (Star-CD function)
13
# Modified by Nuno Lopes (fix for emc2 mesh format (medit version 0))
14
# Modified by Neilen Marais (add gmsh support for reading physical region)
15
""" Module for converting various mesh formats.
16
"""
17
18
# NOTE: This module does not depend on (py)dolfin beeing installed.
19
# NOTE: If future additions need that please import dolfin in a try: except:
20
# NOTE: clause and tell the user to install dolfin if it is not installed.
21
22
import getopt
23
import sys
24
from dolfin_utils.commands import getstatusoutput
25
import re
26
import warnings
27
import os.path
28
29
def format_from_suffix(suffix):
30
"Return format for given suffix"
31
if suffix == "xml":
32
return "xml"
33
elif suffix == "mesh":
34
return "mesh"
35
elif suffix == "gmsh":
36
return "gmsh"
37
elif suffix == "msh":
38
return "gmsh"
39
elif suffix == "gra":
40
return "metis"
41
elif suffix == "grf":
42
return "scotch"
43
elif suffix == "grid":
44
return "diffpack"
45
elif suffix == "inp":
46
return "abaqus"
47
elif suffix == "ncdf":
48
return "NetCDF"
49
elif suffix =="exo":
50
return "ExodusII"
51
elif suffix =="e":
52
return "ExodusII"
53
elif suffix == "vrt" or suffix == "cel":
54
return "StarCD"
55
elif suffix == "ele" or suffix == "node":
56
return "Triangle"
57
else:
58
_error("Sorry, unknown suffix %s." % suffix)
59
60
def mesh2xml(ifilename, ofilename):
61
"""Convert between .mesh and .xml, parser implemented as a
62
state machine:
63
64
0 = read 'Dimension'
65
1 = read dimension
66
2 = read 'Vertices'
67
3 = read number of vertices
68
4 = read next vertex
69
5 = read 'Triangles' or 'Tetrahedra'
70
6 = read number of cells
71
7 = read next cell
72
8 = done
73
74
"""
75
76
print "Converting from Medit format (.mesh) to DOLFIN XML format"
77
78
# Open files
79
ifile = open(ifilename, "r")
80
ofile = open(ofilename, "w")
81
82
# Scan file for cell type
83
cell_type = None
84
dim = 0
85
while 1:
86
87
# Read next line
88
line = ifile.readline()
89
if not line: break
90
91
# Remove newline
92
if line[-1] == "\n":
93
line = line[:-1]
94
95
# Read dimension
96
if line == "Dimension" or line == " Dimension":
97
line = ifile.readline()
98
num_dims = int(line)
99
if num_dims == 2:
100
cell_type = "triangle"
101
dim = 2
102
elif num_dims == 3:
103
cell_type = "tetrahedron"
104
dim = 3
105
break
106
107
# Check that we got the cell type
108
if cell_type == None:
109
_error("Unable to find cell type.")
110
111
# Step to beginning of file
112
ifile.seek(0)
113
114
# Write header
115
write_header_mesh(ofile, cell_type, dim)
116
117
# Current state
118
state = 0
119
120
# Write data
121
num_vertices_read = 0
122
num_cells_read = 0
123
124
while 1:
125
126
# Read next line
127
line = ifile.readline()
128
if not line: break
129
130
# Skip comments
131
if line[0] == '#':
132
continue
133
134
# Remove newline
135
if line[-1] == "\n":
136
line = line[:-1]
137
138
if state == 0:
139
if line == "Dimension" or line == " Dimension":
140
state += 1
141
elif state == 1:
142
num_dims = int(line)
143
state +=1
144
elif state == 2:
145
if line == "Vertices" or line == " Vertices":
146
state += 1
147
elif state == 3:
148
num_vertices = int(line)
149
write_header_vertices(ofile, num_vertices)
150
state +=1
151
elif state == 4:
152
if num_dims == 2:
153
(x, y, tmp) = line.split()
154
x = float(x)
155
y = float(y)
156
z = 0.0
157
elif num_dims == 3:
158
(x, y, z, tmp) = line.split()
159
x = float(x)
160
y = float(y)
161
z = float(z)
162
write_vertex(ofile, num_vertices_read, x, y, z)
163
num_vertices_read +=1
164
if num_vertices == num_vertices_read:
165
write_footer_vertices(ofile)
166
state += 1
167
elif state == 5:
168
if (line == "Triangles" or line == " Triangles") and num_dims == 2:
169
state += 1
170
if line == "Tetrahedra" and num_dims == 3:
171
state += 1
172
elif state == 6:
173
num_cells = int(line)
174
write_header_cells(ofile, num_cells)
175
state +=1
176
elif state == 7:
177
if num_dims == 2:
178
(n0, n1, n2, tmp) = line.split()
179
n0 = int(n0) - 1
180
n1 = int(n1) - 1
181
n2 = int(n2) - 1
182
write_cell_triangle(ofile, num_cells_read, n0, n1, n2)
183
elif num_dims == 3:
184
(n0, n1, n2, n3, tmp) = line.split()
185
n0 = int(n0) - 1
186
n1 = int(n1) - 1
187
n2 = int(n2) - 1
188
n3 = int(n3) - 1
189
write_cell_tetrahedron(ofile, num_cells_read, n0, n1, n2, n3)
190
num_cells_read +=1
191
if num_cells == num_cells_read:
192
write_footer_cells(ofile)
193
state += 1
194
elif state == 8:
195
break
196
197
# Check that we got all data
198
if state == 8:
199
print "Conversion done"
200
else:
201
_error("Missing data, unable to convert")
202
203
# Write footer
204
write_footer_mesh(ofile)
205
206
# Close files
207
ifile.close()
208
ofile.close()
209
210
def gmsh2xml(ifilename, handler):
211
"""Convert between .gmsh v2.0 format (http://www.geuz.org/gmsh/) and .xml,
212
parser implemented as a state machine:
213
214
0 = read 'MeshFormat'
215
1 = read mesh format data
216
2 = read 'EndMeshFormat'
217
3 = read 'Nodes'
218
4 = read number of vertices
219
5 = read vertices
220
6 = read 'EndNodes'
221
7 = read 'Elements'
222
8 = read number of cells
223
9 = read cells
224
10 = done
225
226
Afterwards, extract physical region numbers if they are defined in
227
the mesh file as a mesh function.
228
229
"""
230
231
print "Converting from Gmsh format (.msh, .gmsh) to DOLFIN XML format"
232
233
# Open files
234
ifile = open(ifilename, "r")
235
236
# Scan file for cell type
237
cell_type = None
238
dim = 0
239
line = ifile.readline()
240
while line:
241
242
# Remove newline
243
if line[-1] == "\n":
244
line = line[:-1]
245
246
# Read dimension
247
if line.find("$Elements") == 0:
248
249
line = ifile.readline()
250
num_cells = int(line)
251
num_cells_counted = 0
252
if num_cells == 0:
253
_error("No cells found in gmsh file.")
254
line = ifile.readline()
255
256
# Now iterate through elements to find largest dimension. Gmsh
257
# format might include elements of lower dimensions in the element list.
258
# We also need to count number of elements of correct dimensions.
259
# Also determine which vertices are not used.
260
dim_2_count = 0
261
dim_3_count = 0
262
vertices_2_used = []
263
# Array used to store gmsh tags for 2D (type 2/triangular) elements
264
tags_2 = []
265
# Array used to store gmsh tags for 3D (type 4/tet) elements
266
tags_3 = []
267
vertices_3_used = []
268
while line.find("$EndElements") == -1:
269
element = line.split()
270
elem_type = int(element[1])
271
num_tags = int(element[2])
272
if elem_type == 2:
273
if dim < 2:
274
cell_type = "triangle"
275
dim = 2
276
node_num_list = [int(node) for node in element[3 + num_tags:]]
277
vertices_2_used.extend(node_num_list)
278
if num_tags > 0:
279
tags_2.append(tuple(int(tag) for tag in element[3:3+num_tags]))
280
dim_2_count += 1
281
elif elem_type == 4:
282
if dim < 3:
283
cell_type = "tetrahedron"
284
dim = 3
285
vertices_2_used = None
286
node_num_list = [int(node) for node in element[3 + num_tags:]]
287
vertices_3_used.extend(node_num_list)
288
if num_tags > 0:
289
tags_3.append(tuple(int(tag) for tag in element[3:3+num_tags]))
290
dim_3_count += 1
291
line = ifile.readline()
292
else:
293
# Read next line
294
line = ifile.readline()
295
296
# Check that we got the cell type and set num_cells_counted
297
if cell_type == None:
298
_error("Unable to find cell type.")
299
if dim == 3:
300
num_cells_counted = dim_3_count
301
vertex_set = set(vertices_3_used)
302
vertices_3_used = None
303
elif dim == 2:
304
num_cells_counted = dim_2_count
305
vertex_set = set(vertices_2_used)
306
vertices_2_used = None
307
308
vertex_dict = {}
309
for n,v in enumerate(vertex_set):
310
vertex_dict[v] = n
311
312
# Step to beginning of file
313
ifile.seek(0)
314
315
# Set mesh type
316
handler.set_mesh_type(cell_type, dim)
317
318
# Initialise node list (gmsh does not export all vertexes in order)
319
nodelist = {}
320
321
# Current state
322
state = 0
323
324
# Write data
325
num_vertices_read = 0
326
num_cells_read = 0
327
328
while state != 10:
329
330
# Read next line
331
line = ifile.readline()
332
if not line: break
333
334
# Skip comments
335
if line[0] == '#':
336
continue
337
338
# Remove newline
339
if line[-1] == "\n":
340
line = line[:-1]
341
342
if state == 0:
343
if line == "$MeshFormat":
344
state = 1
345
elif state == 1:
346
(version, file_type, data_size) = line.split()
347
state = 2
348
elif state == 2:
349
if line == "$EndMeshFormat":
350
state = 3
351
elif state == 3:
352
if line == "$Nodes":
353
state = 4
354
elif state == 4:
355
num_vertices = len(vertex_dict)
356
handler.start_vertices(num_vertices)
357
state = 5
358
elif state == 5:
359
(node_no, x, y, z) = line.split()
360
node_no = int(node_no)
361
x,y,z = [float(xx) for xx in (x,y,z)]
362
if vertex_dict.has_key(node_no):
363
node_no = vertex_dict[node_no]
364
else:
365
continue
366
nodelist[int(node_no)] = num_vertices_read
367
handler.add_vertex(num_vertices_read, [x, y, z])
368
num_vertices_read +=1
369
370
if num_vertices == num_vertices_read:
371
handler.end_vertices()
372
state = 6
373
elif state == 6:
374
if line == "$EndNodes":
375
state = 7
376
elif state == 7:
377
if line == "$Elements":
378
state = 8
379
elif state == 8:
380
handler.start_cells(num_cells_counted)
381
state = 9
382
elif state == 9:
383
element = line.split()
384
elem_type = int(element[1])
385
num_tags = int(element[2])
386
if elem_type == 2 and dim == 2:
387
node_num_list = [vertex_dict[int(node)] for node in element[3 + num_tags:]]
388
for node in node_num_list:
389
if not node in nodelist:
390
_error("Vertex %d of triangle %d not previously defined." %
391
(node, num_cells_read))
392
cell_nodes = [nodelist[n] for n in node_num_list]
393
handler.add_cell(num_cells_read, cell_nodes)
394
num_cells_read +=1
395
elif elem_type == 4 and dim == 3:
396
node_num_list = [vertex_dict[int(node)] for node in element[3 + num_tags:]]
397
for node in node_num_list:
398
if not node in nodelist:
399
_error("Vertex %d of tetrahedron %d not previously defined." %
400
(node, num_cells_read))
401
# import pdb ; pdb.set_trace()
402
cell_nodes = [nodelist[n] for n in node_num_list]
403
handler.add_cell(num_cells_read, cell_nodes)
404
num_cells_read +=1
405
406
if num_cells_counted == num_cells_read:
407
handler.end_cells()
408
state = 10
409
elif state == 10:
410
break
411
412
# Write mesh function based on the Physical Regions defined by
413
# gmsh, but only if they are not all zero. All zero physical
414
# regions indicate that no physical regions were defined.
415
if dim == 2:
416
tags = tags_2
417
elif dim == 3:
418
tags = tags_3
419
else:
420
_error("Gmsh tags not supported for dimension %i. Probably a bug" % dim)
421
422
physical_regions = tuple(tag[0] for tag in tags)
423
if not all(tag == 0 for tag in tags):
424
handler.start_meshfunction("physical_region", dim, num_cells)
425
for i, physical_region in enumerate(physical_regions):
426
handler.add_entity_meshfunction(i, physical_region)
427
handler.end_meshfunction()
428
429
# Check that we got all data
430
if state == 10:
431
print "Conversion done"
432
else:
433
_error("Missing data, unable to convert \n\ Did you use version 2.0 of the gmsh file format?")
434
435
# Close files
436
ifile.close()
437
438
def triangle2xml(ifilename, ofilename):
439
"""Convert between triangle format (http://www.cs.cmu.edu/~quake/triangle.html) and .xml. The
440
given ifilename should be the prefix for the corresponding .node, and .ele files.
441
"""
442
443
def get_next_line (fp):
444
"""Helper function for skipping comments and blank lines"""
445
line = fp.readline()
446
if line == '':
447
_error("Hit end of file prematurely.")
448
line = line.strip()
449
if not (line.startswith('#') or line == ''):
450
return line
451
return get_next_line(fp)
452
453
454
print "Converting from Triangle format {.node, .ele} to DOLFIN XML format"
455
456
# Open files
457
for suffix in [".node", ".ele"]:
458
if suffix in ifilename and ifilename[-len(suffix):] == suffix:
459
ifilename = ifilename.replace(suffix, "")
460
node_file = open(ifilename+".node", "r")
461
ele_file = open(ifilename+".ele", "r")
462
ofile = open(ofilename, "w")
463
464
# Read all the nodes
465
nodes = {}
466
num_nodes, dim, attr, bound = map(int, get_next_line(node_file).split())
467
while len(nodes) < num_nodes:
468
node, x, y = get_next_line(node_file).split()[:3]
469
nodes[int(node)] = (float(x), float(y))
470
471
# Read all the triangles
472
tris = {}
473
num_tris, n_per_tri, attrs = map(int, get_next_line(ele_file).split())
474
while len(tris) < num_tris:
475
tri, n1, n2, n3 = map(int, get_next_line(ele_file).split()[:4])
476
tris[tri] = (n1, n2, n3)
477
478
# Write everything out
479
write_header_mesh(ofile, "triangle", 2)
480
write_header_vertices(ofile, num_nodes)
481
node_off = 0 if nodes.has_key(0) else -1
482
for node, node_t in nodes.iteritems():
483
write_vertex(ofile, node+node_off, node_t[0], node_t[1], 0.0)
484
write_footer_vertices(ofile)
485
write_header_cells(ofile, num_tris)
486
tri_off = 0 if tris.has_key(0) else -1
487
for tri, tri_t in tris.iteritems():
488
write_cell_triangle(ofile, tri+tri_off, tri_t[0] + node_off,
489
tri_t[1] + node_off, tri_t[2] + node_off)
490
write_footer_cells(ofile)
491
write_footer_mesh(ofile)
492
493
# Close files
494
node_file.close()
495
ele_file.close()
496
ofile.close()
497
498
499
def xml_old2xml(ifilename, ofilename):
500
"Convert from old DOLFIN XML format to new."
501
502
print "Converting from old (pre DOLFIN 0.6.2) to new DOLFIN XML format..."
503
504
# Open files
505
ifile = open(ifilename, "r")
506
ofile = open(ofilename, "w")
507
508
# Scan file for cell type (assuming there is just one)
509
cell_type = None
510
dim = 0
511
while 1:
512
513
# Read next line
514
line = ifile.readline()
515
if not line: break
516
517
# Read dimension
518
if "<triangle" in line:
519
cell_type = "triangle"
520
dim = 2
521
break
522
elif "<tetrahedron" in line:
523
cell_type = "tetrahedron"
524
dim = 3
525
break
526
527
# Step to beginning of file
528
ifile.seek(0)
529
530
# Read lines and make changes
531
while 1:
532
533
# Read next line
534
line = ifile.readline()
535
if not line: break
536
537
# Modify line
538
if "xmlns" in line:
539
line = "<dolfin xmlns:dolfin=\"http://www.fenicsproject.org\">\n"
540
if "<mesh>" in line:
541
line = " <mesh celltype=\"%s\" dim=\"%d\">\n" % (cell_type, dim)
542
if dim == 2 and " z=\"0.0\"" in line:
543
line = line.replace(" z=\"0.0\"", "")
544
if " name=" in line:
545
line = line.replace(" name=", " index=")
546
if " name =" in line:
547
line = line.replace(" name =", " index=")
548
if "n0" in line:
549
line = line.replace("n0", "v0")
550
if "n1" in line:
551
line = line.replace("n1", "v1")
552
if "n2" in line:
553
line = line.replace("n2", "v2")
554
if "n3" in line:
555
line = line.replace("n3", "v3")
556
557
# Write line
558
ofile.write(line)
559
560
# Close files
561
ifile.close();
562
ofile.close();
563
print "Conversion done"
564
565
def metis_graph2graph_xml(ifilename, ofilename):
566
"Convert from Metis graph format to DOLFIN Graph XML."
567
568
print "Converting from Metis graph format to DOLFIN Graph XML."
569
570
# Open files
571
ifile = open(ifilename, "r")
572
ofile = open(ofilename, "w")
573
574
# Read number of vertices and edges
575
line = ifile.readline()
576
if not line:
577
_error("Empty file")
578
579
(num_vertices, num_edges) = line.split()
580
581
write_header_graph(ofile, "directed")
582
write_header_vertices(ofile, int(num_vertices))
583
584
for i in range(int(num_vertices)):
585
line = ifile.readline()
586
edges = line.split()
587
write_graph_vertex(ofile, i, len(edges))
588
589
write_footer_vertices(ofile)
590
write_header_edges(ofile, 2*int(num_edges))
591
592
# Step to beginning of file and skip header info
593
ifile.seek(0)
594
ifile.readline()
595
for i in range(int(num_vertices)):
596
print "vertex %g", i
597
line = ifile.readline()
598
edges = line.split()
599
for e in edges:
600
write_graph_edge(ofile, i, int(e))
601
602
write_footer_edges(ofile)
603
write_footer_graph(ofile)
604
605
# Close files
606
ifile.close();
607
ofile.close();
608
609
def scotch_graph2graph_xml(ifilename, ofilename):
610
"Convert from Scotch graph format to DOLFIN Graph XML."
611
612
print "Converting from Scotch graph format to DOLFIN Graph XML."
613
614
# Open files
615
ifile = open(ifilename, "r")
616
ofile = open(ofilename, "w")
617
618
# Skip graph file version number
619
ifile.readline()
620
621
# Read number of vertices and edges
622
line = ifile.readline()
623
if not line:
624
_error("Empty file")
625
626
(num_vertices, num_edges) = line.split()
627
628
# Read start index and numeric flag
629
# Start index is 0 or 1 (C/Fortran)
630
# Numeric flag is 3 bits where bit 1 enables vertex labels
631
# bit 2 enables edge weights and bit 3 enables vertex weights
632
633
line = ifile.readline()
634
(start_index, numeric_flag) = line.split()
635
636
# Handling not implented
637
if not numeric_flag == "000":
638
_error("Handling of scotch vertex labels, edge- and vertex weights not implemented")
639
640
write_header_graph(ofile, "undirected")
641
write_header_vertices(ofile, int(num_vertices))
642
643
# Read vertices and edges, first number gives number of edges from this vertex (not used)
644
for i in range(int(num_vertices)):
645
line = ifile.readline()
646
edges = line.split()
647
write_graph_vertex(ofile, i, len(edges)-1)
648
649
write_footer_vertices(ofile)
650
write_header_edges(ofile, int(num_edges))
651
652
653
# Step to beginning of file and skip header info
654
ifile.seek(0)
655
ifile.readline()
656
ifile.readline()
657
ifile.readline()
658
for i in range(int(num_vertices)):
659
line = ifile.readline()
660
661
edges = line.split()
662
for j in range(1, len(edges)):
663
write_graph_edge(ofile, i, int(edges[j]))
664
665
write_footer_edges(ofile)
666
write_footer_graph(ofile)
667
668
# Close files
669
ifile.close();
670
ofile.close();
671
672
def write_header_meshfunction(ofile, dimensions, size):
673
header = """<?xml version="1.0" encoding="UTF-8"?>
674
<dolfin xmlns:dolfin="http://www.fenics.org/dolfin/">
675
<meshfunction type="uint" dim="%d" size="%d">
676
""" % (dimensions, size)
677
ofile.write(header)
678
679
def write_entity_meshfunction(ofile, index, value):
680
ofile.write(""" <entity index=\"%d\" value=\"%d\"/>
681
""" % (index, value))
682
683
def write_footer_meshfunction(ofile):
684
ofile.write(""" </meshfunction>
685
</dolfin>""")
686
687
def diffpack2xml(ifilename, ofilename):
688
"Convert from Diffpack tetrahedral grid format to DOLFIN XML."
689
690
print diffpack2xml.__doc__
691
692
# Format strings for MeshFunction XML files
693
meshfunction_header = """\
694
<?xml version="1.0" encoding="UTF-8"?>\n
695
<dolfin xmlns:dolfin="http://www.fenics.org/dolfin/">
696
<meshfunction type="uint" dim="%d" size="%d">\n"""
697
meshfunction_entity = " <entity index=\"%d\" value=\"%d\"/>\n"
698
meshfunction_footer = " </meshfunction>\n</dolfin>"
699
700
# Open files
701
ifile = open(ifilename, "r")
702
ofile = open(ofilename, "w")
703
ofile_mat = open(ofilename.split(".")[0]+"_mat.xml", "w")
704
ofile_bi = open(ofilename.split(".")[0]+"_bi.xml", "w")
705
706
# Read and analyze header
707
while 1:
708
line = ifile.readline()
709
if not line:
710
_error("Empty file")
711
if line[0] == "#":
712
break
713
if re.search(r"Number of elements", line):
714
num_cells = int(re.match(r".*\s(\d+).*", line).group(1))
715
if re.search(r"Number of nodes", line):
716
num_vertices = int(re.match(r".*\s(\d+).*", line).group(1))
717
718
write_header_mesh(ofile, "tetrahedron", 3)
719
write_header_vertices(ofile, num_vertices)
720
ofile_bi.write(meshfunction_header % (0, num_vertices))
721
ofile_mat.write(meshfunction_header % (3, num_cells))
722
723
# Read & write vertices
724
# Note that only first boundary indicator is rewriten into XML
725
for i in range(num_vertices):
726
line = ifile.readline()
727
m = re.match(r"^.*\(\s*(.*)\s*\).*\](.*)$", line)
728
x = re.split("[\s,]+", m.group(1))
729
write_vertex(ofile, i, x[0], x[1], x[2])
730
tmp = m.group(2).split()
731
if len(tmp) > 0:
732
bi = int(tmp[0])
733
else:
734
bi = 0
735
ofile_bi.write(meshfunction_entity % (i, bi))
736
737
write_footer_vertices(ofile)
738
write_header_cells(ofile, num_cells)
739
740
# Ignore comment lines
741
while 1:
742
line = ifile.readline()
743
if not line:
744
_error("Empty file")
745
if line[0] == "#":
746
break
747
748
# Read & write cells
749
for i in range(int(num_cells)):
750
line = ifile.readline()
751
v = line.split();
752
if v[1] != "ElmT4n3D":
753
_error("Only tetrahedral elements (ElmT4n3D) are implemented.")
754
write_cell_tetrahedron(ofile, i, int(v[3])-1, int(v[4])-1, int(v[5])-1, int(v[6])-1)
755
ofile_mat.write(meshfunction_entity % (i, int(v[2])))
756
757
write_footer_cells(ofile)
758
write_footer_mesh(ofile)
759
ofile_bi.write(meshfunction_footer)
760
ofile_mat.write(meshfunction_footer)
761
762
# Close files
763
ifile.close()
764
ofile.close()
765
ofile_mat.close()
766
ofile_bi.close()
767
768
class ParseError(Exception):
769
""" Error encountered in source file.
770
"""
771
772
class DataHandler(object):
773
""" Baseclass for handlers of mesh data.
774
775
The actual handling of mesh data encountered in the source file is
776
delegated to a polymorfic object. Typically, the delegate will write the
777
data to XML.
778
@ivar _state: the state which the handler is in, one of State_*.
779
@ivar _cell_type: cell type in mesh. One of CellType_*.
780
@ivar _dim: mesh dimensions.
781
"""
782
State_Invalid, State_Init, State_Vertices, State_Cells, State_MeshFunction = range(5)
783
CellType_Tetrahedron, CellType_Triangle = range(2)
784
785
def __init__(self):
786
self._state = self.State_Invalid
787
788
def set_mesh_type(self, cell_type, dim):
789
assert self._state == self.State_Invalid
790
self._state = self.State_Init
791
if cell_type == "tetrahedron":
792
self._cell_type = self.CellType_Tetrahedron
793
elif cell_type == "triangle":
794
self._cell_type = self.CellType_Triangle
795
self._dim = dim
796
797
def start_vertices(self, num_vertices):
798
assert self._state == self.State_Init
799
self._state = self.State_Vertices
800
801
def add_vertex(self, vertex, coords):
802
assert self._state == self.State_Vertices
803
804
def end_vertices(self):
805
assert self._state == self.State_Vertices
806
self._state = self.State_Init
807
808
def start_cells(self, num_cells):
809
assert self._state == self.State_Init
810
self._state = self.State_Cells
811
812
def add_cell(self, cell, nodes):
813
assert self._state == self.State_Cells
814
815
def end_cells(self):
816
assert self._state == self.State_Cells
817
self._state = self.State_Init
818
819
def start_meshfunction(self, name, dim, size):
820
assert self._state == self.State_Init
821
self._state = self.State_MeshFunction
822
823
def add_entity_meshfunction(self, index, value):
824
assert self._state == self.State_MeshFunction
825
826
def end_meshfunction(self):
827
assert self._state == self.State_MeshFunction
828
self._state = self.State_Init
829
830
def warn(self, msg):
831
""" Issue warning during parse.
832
"""
833
warnings.warn(msg)
834
835
def error(self, msg):
836
""" Raise error during parse.
837
838
This method is expected to raise ParseError.
839
"""
840
raise ParseError(msg)
841
842
def close(self):
843
self._state = self.State_Invalid
844
845
class XmlHandler(DataHandler):
846
""" Data handler class which writes to Dolfin XML.
847
"""
848
def __init__(self, ofilename):
849
DataHandler.__init__(self)
850
self._ofilename = ofilename
851
self.__ofile = file(ofilename, "wb")
852
self.__ofile_meshfunc = None
853
854
def set_mesh_type(self, cell_type, dim):
855
DataHandler.set_mesh_type(self, cell_type, dim)
856
write_header_mesh(self.__ofile, cell_type, dim)
857
858
def start_vertices(self, num_vertices):
859
DataHandler.start_vertices(self, num_vertices)
860
write_header_vertices(self.__ofile, num_vertices)
861
862
def add_vertex(self, vertex, coords):
863
DataHandler.add_vertex(self, vertex, coords)
864
write_vertex(self.__ofile, vertex, *coords)
865
866
def end_vertices(self):
867
DataHandler.end_vertices(self)
868
write_footer_vertices(self.__ofile)
869
870
def start_cells(self, num_cells):
871
DataHandler.start_cells(self, num_cells)
872
write_header_cells(self.__ofile, num_cells)
873
874
def add_cell(self, cell, nodes):
875
DataHandler.add_cell(self, cell, nodes)
876
if self._cell_type == self.CellType_Tetrahedron:
877
func = write_cell_tetrahedron
878
elif self._cell_type == self.CellType_Triangle:
879
func = write_cell_triangle
880
881
func(self.__ofile, cell, *nodes)
882
883
def end_cells(self):
884
DataHandler.end_cells(self)
885
write_footer_cells(self.__ofile)
886
887
def start_meshfunction(self, name, dim, size):
888
DataHandler.start_meshfunction(self, name, dim, size)
889
fname = os.path.splitext(self.__ofile.name)[0]
890
self.__ofile_meshfunc = file("%s_%s.xml" % (fname, name), "wb")
891
write_header_meshfunction(self.__ofile_meshfunc, dim, size)
892
893
def add_entity_meshfunction(self, index, value):
894
DataHandler.add_entity_meshfunction(self, index, value)
895
write_entity_meshfunction(self.__ofile_meshfunc, index, value)
896
897
def end_meshfunction(self):
898
DataHandler.end_meshfunction(self)
899
write_footer_meshfunction(self.__ofile_meshfunc)
900
self.__ofile_meshfunc.close()
901
self.__ofile_meshfunc = None
902
903
def close(self):
904
DataHandler.close(self)
905
if self.__ofile.closed:
906
return
907
write_footer_mesh(self.__ofile)
908
self.__ofile.close()
909
if self.__ofile_meshfunc is not None:
910
self.__ofile_meshfunc.close()
911
912
def _abaqus(ifilename, handler):
913
""" Convert from Abaqus.
914
915
The Abaqus format first defines a node block, then there should be a number
916
of elements containing these nodes.
917
"""
918
params = False
919
ifile = file(ifilename, "rb")
920
handler.set_mesh_type("tetrahedron", 3)
921
922
# Read definitions
923
924
def read_params(params_spec, pnames, lineno):
925
params = {}
926
for p in params_spec:
927
m = re.match(r"(.+)=(.+)", p)
928
if m is not None:
929
pname, val = m.groups()
930
else:
931
handler.warn("Invalid parameter syntax on line %d: %s" % (lineno, p))
932
continue
933
for pn in pnames:
934
if pn == pname:
935
params[pn] = val
936
break
937
938
return params
939
940
nodes = {}
941
elems = {}
942
eid2elset = {}
943
material2elsetids = {}
944
materials = []
945
re_sect = re.compile(r"\*([^,]+)(?:,(.*))?")
946
re_node = re.compile(r"(\d+),\s*(.+),\s*(.+),\s*(.+)")
947
re_tetra = re.compile(r"(\d+),\s*(\d+),\s*(\d+),\s*(\d+),\s*(\d+)")
948
sect = None
949
for lineno, l in enumerate(ifile):
950
l = l.strip().lower()
951
m = re_sect.match(l)
952
if m is not None:
953
sect, params_str = m.groups()
954
params_spec = ([s.strip() for s in params_str.split(",")] if params_str
955
else [])
956
957
if sect == "element":
958
pnames = ("type", "elset")
959
params = read_params(params_spec, pnames, lineno)
960
if "type" not in params:
961
handler.error("Element on line %d doesn't declare TYPE" %
962
(lineno,))
963
tp, elset = params["type"], params.get("elset")
964
if tp not in ("c3d4", "dc3d4"):
965
handler.warn("Unsupported element type '%s' on line %d" % (tp, lineno))
966
supported_elem = False
967
else:
968
supported_elem = True
969
elif sect == "solid section":
970
pnames = ("material", "elset")
971
params = read_params(params_spec, pnames, lineno)
972
for pname in pnames:
973
if pname not in params:
974
handler.error("Solid section on line %d doesn't "
975
"declare %s" % (lineno, pname.upper()))
976
matname = params["material"]
977
material2elsetids.setdefault(matname, []).append(params["elset"])
978
elif sect == "material":
979
name = read_params(params_spec, ["name"], lineno)["name"]
980
materials.append(name)
981
# We've read the section's heading, continue to next line
982
continue
983
984
# Read section entry
985
986
if sect == "node":
987
# Read node definition
988
m = re_node.match(l)
989
if m is None:
990
handler.warn("Node on line %d is on unsupported format" % (lineno,))
991
continue
992
idx, c0, c1, c2 = m.groups()
993
try: coords = [float(c) for c in (c0, c1, c2)]
994
except ValueError:
995
handler.warn("Node on line %d contains non-numeric coordinates"
996
% (lineno,))
997
continue
998
nodes[int(idx)] = coords
999
elif sect == "element":
1000
if not supported_elem:
1001
continue
1002
m = re_tetra.match(l)
1003
if m is None:
1004
handler.error("Node on line %d badly specified (expected 3 "
1005
"coordinates)" % (lineno,))
1006
idx, n0, n1, n2, n3 = [int(x) for x in m.groups()]
1007
elems[idx] = (tp, n0, n1, n2, n3)
1008
eid2elset.setdefault(elset, set()).add(idx)
1009
1010
ifile.close()
1011
1012
# Note that vertices/cells must be consecutively numbered, which isn't
1013
# necessarily the case in Abaqus. Therefore we enumerate and translate
1014
# original IDs to sequence indexes.
1015
1016
handler.start_vertices(len(nodes))
1017
nodeids = nodes.keys()
1018
nodeids.sort()
1019
for idx, nid in enumerate(nodeids):
1020
handler.add_vertex(idx, nodes[nid])
1021
handler.end_vertices()
1022
1023
handler.start_cells(len(elems))
1024
elemids = elems.keys()
1025
elemids.sort()
1026
for idx, eid in enumerate(elemids):
1027
elem = elems[eid]
1028
tp = elem[0]
1029
elemnodes = []
1030
for nid in elem[1:]:
1031
try: elemnodes.append(nodeids.index(nid))
1032
except ValueError:
1033
handler.error("Element %s references non-existent node %s" % (eid, nid))
1034
handler.add_cell(idx, elemnodes)
1035
handler.end_cells()
1036
1037
# Define the material function for the cells
1038
1039
num_entities = 0
1040
for matname, elsetids in material2elsetids.items():
1041
if matname not in materials:
1042
handler.error("Unknown material %s referred to for element sets %s" %
1043
(matname, ", ".join(elsetids)))
1044
num_entities += len(elsetids)
1045
handler.start_meshfunction("material", 3, num_entities)
1046
# Each material is associated with a number of element sets
1047
for i, matname in enumerate(materials):
1048
try: elsetids = material2elsetids[matname]
1049
except KeyError:
1050
# No elements for this material
1051
continue
1052
# For each element set associated with this material
1053
elsets = []
1054
for eid in elsetids:
1055
try: elsets.append(eid2elset[eid])
1056
except KeyError:
1057
handler.error("Material '%s' is assigned to undefined element "
1058
"set '%s'" % (matname, eid))
1059
for elset in elsets:
1060
for elemid in elset:
1061
handler.add_entity_meshfunction(elemids.index(elemid), i)
1062
handler.end_meshfunction()
1063
1064
def netcdf2xml(ifilename,ofilename):
1065
"Convert from NetCDF format to DOLFIN XML."
1066
1067
print "Converting from NetCDF format (.ncdf) to DOLFIN XML format"
1068
1069
# Open files
1070
ifile = open(ifilename, "r")
1071
ofile = open(ofilename, "w")
1072
1073
1074
cell_type = None
1075
dim = 0
1076
1077
# Scan file for dimension, number of nodes, number of elements
1078
while 1:
1079
line = ifile.readline()
1080
if not line:
1081
_error("Empty file")
1082
if re.search(r"num_dim.*=", line):
1083
dim = int(re.match(".*\s=\s(\d+)\s;",line).group(1))
1084
if re.search(r"num_nodes.*=", line):
1085
num_vertices = int(re.match(".*\s=\s(\d+)\s;",line).group(1))
1086
if re.search(r"num_elem.*=", line):
1087
num_cells = int(re.match(".*\s=\s(\d+)\s;",line).group(1))
1088
if re.search(r"connect1 =",line):
1089
break
1090
1091
num_dims=dim
1092
1093
# Set cell type
1094
if dim == 2:
1095
cell_type ="triangle"
1096
if dim == 3:
1097
cell_type ="tetrahedron"
1098
1099
# Check that we got the cell type
1100
if cell_type == None:
1101
_error("Unable to find cell type.")
1102
1103
# Write header
1104
write_header_mesh(ofile, cell_type, dim)
1105
write_header_cells(ofile, num_cells)
1106
num_cells_read = 0
1107
1108
# Read and write cells
1109
while 1:
1110
# Read next line
1111
line = ifile.readline()
1112
if not line:
1113
break
1114
connect=re.split("[,;]",line)
1115
if num_dims == 2:
1116
n0 = int(connect[0])-1
1117
n1 = int(connect[1])-1
1118
n2 = int(connect[2])-1
1119
write_cell_triangle(ofile, num_cells_read, n0, n1, n2)
1120
elif num_dims == 3:
1121
n0 = int(connect[0])-1
1122
n1 = int(connect[1])-1
1123
n2 = int(connect[2])-1
1124
n3 = int(connect[3])-1
1125
write_cell_tetrahedron(ofile, num_cells_read, n0, n1, n2, n3)
1126
num_cells_read +=1
1127
if num_cells == num_cells_read:
1128
write_footer_cells(ofile)
1129
write_header_vertices(ofile, num_vertices)
1130
break
1131
1132
num_vertices_read = 0
1133
coords = [[],[],[]]
1134
coord = -1
1135
1136
while 1:
1137
line = ifile.readline()
1138
if not line:
1139
_error("Missing data")
1140
if re.search(r"coord =",line):
1141
break
1142
1143
# Read vertices
1144
while 1:
1145
line = ifile.readline()
1146
# print line
1147
if not line:
1148
break
1149
if re.search(r"\A\s\s\S+,",line):
1150
# print line
1151
coord+=1
1152
print "Found x_"+str(coord)+" coordinates"
1153
coords[coord] += line.split()
1154
if re.search(r";",line):
1155
break
1156
1157
# Write vertices
1158
for i in range(num_vertices):
1159
if num_dims == 2:
1160
x = float(re.split(",",coords[0].pop(0))[0])
1161
y = float(re.split(",",coords[1].pop(0))[0])
1162
z = 0
1163
if num_dims == 3:
1164
x = float(re.split(",",coords[0].pop(0))[0])
1165
y = float(re.split(",",coords[1].pop(0))[0])
1166
z = float(re.split(",",coords[2].pop(0))[0])
1167
write_vertex(ofile, i, x, y, z)
1168
1169
1170
# Write footer
1171
write_footer_vertices(ofile)
1172
write_footer_mesh(ofile)
1173
1174
# Close files
1175
ifile.close()
1176
ofile.close()
1177
1178
def exodus2xml(ifilename,ofilename):
1179
"Convert from Exodus II format to DOLFIN XML."
1180
1181
print "Converting from Exodus II format to NetCDF format"
1182
1183
name = ifilename.split(".")[0]
1184
netcdffilename = name +".ncdf"
1185
status, output = getstatusoutput('ncdump '+ifilename + ' > '+netcdffilename)
1186
if status != 0:
1187
raise IOError, "Something wrong while executing ncdump. Is ncdump "\
1188
"installed on the system?"
1189
netcdf2xml(netcdffilename, ofilename)
1190
1191
# Write mesh header
1192
def write_header_mesh(ofile, cell_type, dim):
1193
ofile.write("""\
1194
<?xml version=\"1.0\" encoding=\"UTF-8\"?>
1195
1196
<dolfin xmlns:dolfin=\"http://www.fenics.org/dolfin/\">
1197
<mesh celltype="%s" dim="%d">
1198
""" % (cell_type, dim))
1199
1200
# Write graph header
1201
def write_header_graph(ofile, graph_type):
1202
ofile.write("""\
1203
<?xml version=\"1.0\" encoding=\"UTF-8\"?>
1204
1205
<dolfin xmlns:dolfin=\"http://www.fenics.org/dolfin/\">
1206
<graph type="%s">
1207
""" % (graph_type))
1208
1209
# Write mesh footer
1210
def write_footer_mesh(ofile):
1211
ofile.write("""\
1212
</mesh>
1213
</dolfin>
1214
""")
1215
1216
# Write graph footer
1217
def write_footer_graph(ofile):
1218
ofile.write("""\
1219
</graph>
1220
</dolfin>
1221
""")
1222
1223
def write_header_vertices(ofile, num_vertices):
1224
"Write vertices header"
1225
print "Expecting %d vertices" % num_vertices
1226
ofile.write(" <vertices size=\"%d\">\n" % num_vertices)
1227
1228
def write_footer_vertices(ofile):
1229
"Write vertices footer"
1230
ofile.write(" </vertices>\n")
1231
print "Found all vertices"
1232
1233
def write_header_edges(ofile, num_edges):
1234
"Write edges header"
1235
print "Expecting %d edges" % num_edges
1236
ofile.write(" <edges size=\"%d\">\n" % num_edges)
1237
1238
def write_footer_edges(ofile):
1239
"Write edges footer"
1240
ofile.write(" </edges>\n")
1241
print "Found all edges"
1242
1243
def write_vertex(ofile, vertex, x, y, z):
1244
"Write vertex"
1245
ofile.write(" <vertex index=\"%d\" x=\"%s\" y=\"%s\" z=\"%s\"/>\n" % \
1246
(vertex, x, y, z))
1247
1248
def write_graph_vertex(ofile, vertex, num_edges, weight = 1):
1249
"Write graph vertex"
1250
ofile.write(" <vertex index=\"%d\" num_edges=\"%d\" weight=\"%d\"/>\n" % \
1251
(vertex, num_edges, weight))
1252
1253
def write_graph_edge(ofile, v1, v2, weight = 1):
1254
"Write graph edge"
1255
ofile.write(" <edge v1=\"%d\" v2=\"%d\" weight=\"%d\"/>\n" % \
1256
(v1, v2, weight))
1257
1258
def write_header_cells(ofile, num_cells):
1259
"Write cells header"
1260
ofile.write(" <cells size=\"%d\">\n" % num_cells)
1261
print "Expecting %d cells" % num_cells
1262
1263
def write_footer_cells(ofile):
1264
"Write cells footer"
1265
ofile.write(" </cells>\n")
1266
print "Found all cells"
1267
1268
def write_cell_triangle(ofile, cell, n0, n1, n2):
1269
"Write cell (triangle)"
1270
ofile.write(" <triangle index=\"%d\" v0=\"%d\" v1=\"%d\" v2=\"%d\"/>\n" % \
1271
(cell, n0, n1, n2))
1272
1273
def write_cell_tetrahedron(ofile, cell, n0, n1, n2, n3):
1274
"Write cell (tetrahedron)"
1275
ofile.write(" <tetrahedron index=\"%d\" v0=\"%d\" v1=\"%d\" v2=\"%d\" v3=\"%d\"/>\n" % \
1276
(cell, n0, n1, n2, n3))
1277
1278
def _error(message):
1279
"Write an error message"
1280
for line in message.split("\n"):
1281
print "*** %s" % line
1282
sys.exit(2)
1283
1284
def convert2xml(ifilename, ofilename, iformat=None):
1285
""" Convert a file to the DOLFIN XML format.
1286
"""
1287
convert(ifilename, XmlHandler(ofilename), iformat=iformat)
1288
1289
def convert(ifilename, handler, iformat=None):
1290
""" Convert a file using a provided data handler.
1291
1292
Note that handler.close is called when this function finishes.
1293
@param ifilename: Name of input file.
1294
@param handler: The data handler (instance of L{DataHandler}).
1295
@param iformat: Format of input file.
1296
"""
1297
if iformat is None:
1298
iformat = format_from_suffix(os.path.splitext(ifilename)[1][1:])
1299
# XXX: Backwards-compat
1300
if hasattr(handler, "_ofilename"):
1301
ofilename = handler._ofilename
1302
# Choose conversion
1303
if iformat == "mesh":
1304
# Convert from mesh to xml format
1305
mesh2xml(ifilename, ofilename)
1306
elif iformat == "gmsh":
1307
# Convert from gmsh to xml format
1308
gmsh2xml(ifilename, handler)
1309
elif iformat == "Triangle":
1310
# Convert from Triangle to xml format
1311
triangle2xml(ifilename, ofilename)
1312
elif iformat == "xml-old":
1313
# Convert from old to new xml format
1314
xml_old2xml(ifilename, ofilename)
1315
elif iformat == "metis":
1316
# Convert from metis graph to dolfin graph xml format
1317
metis_graph2graph_xml(ifilename, ofilename)
1318
elif iformat == "scotch":
1319
# Convert from scotch graph to dolfin graph xml format
1320
scotch_graph2graph_xml(ifilename, ofilename)
1321
elif iformat == "diffpack":
1322
# Convert from Diffpack tetrahedral grid format to xml format
1323
diffpack2xml(ifilename, ofilename)
1324
elif iformat == "abaqus":
1325
# Convert from abaqus to xml format
1326
_abaqus(ifilename, handler)
1327
elif iformat == "NetCDF":
1328
# Convert from NetCDF generated from ExodusII format to xml format
1329
netcdf2xml(ifilename, ofilename)
1330
elif iformat =="ExodusII":
1331
# Convert from ExodusII format to xml format via NetCDF
1332
exodus2xml(ifilename, ofilename)
1333
elif iformat == "StarCD":
1334
# Convert from Star-CD tetrahedral grid format to xml format
1335
starcd2xml(ifilename, ofilename)
1336
else:
1337
_error("Sorry, cannot convert between %s and DOLFIN xml file formats." % iformat)
1338
1339
# XXX: handler.close messes things for other input formats than abaqus or gmsh
1340
if iformat in ("abaqus", "gmsh"):
1341
handler.close()
1342
1343
def starcd2xml(ifilename, ofilename):
1344
"Convert from Star-CD tetrahedral grid format to DOLFIN XML."
1345
1346
print starcd2xml.__doc__
1347
1348
if not os.path.isfile(ifilename[:-3] + "vrt") or not os.path.isfile(ifilename[:-3] + "cel"):
1349
print "StarCD format requires one .vrt file and one .cel file"
1350
sys.exit(2)
1351
1352
1353
# open output file
1354
ofile = open(ofilename, "w")
1355
1356
# Open file, the vertices are in a .vrt file
1357
ifile = open(ifilename[:-3] + "vrt", "r")
1358
1359
write_header_mesh(ofile, "tetrahedron", 3)
1360
1361
1362
# Read & write vertices
1363
1364
# first, read all lines (need to sweep to times through the file)
1365
lines = ifile.readlines()
1366
1367
# second, find the number of vertices
1368
num_vertices = -1
1369
counter = 0
1370
# nodenr_map is needed because starcd support node numbering like 1,2,4 (ie 3 is missing)
1371
nodenr_map = {}
1372
for line in lines:
1373
nodenr = int(line[0:15])
1374
nodenr_map[nodenr] = counter
1375
counter += 1
1376
num_vertices = counter
1377
1378
# third, run over all vertices
1379
write_header_vertices(ofile, num_vertices)
1380
for line in lines:
1381
nodenr = int(line[0:15])
1382
vertex0 = float(line[15:31])
1383
vertex1 = float(line[31:47])
1384
vertex2 = float(line[47:63])
1385
write_vertex(ofile, nodenr_map[nodenr], float(vertex0), float(vertex1), float(vertex2))
1386
write_footer_vertices(ofile)
1387
1388
# Open file, the cells are in a .cel file
1389
ifile = open(ifilename[:-3] + "cel", "r")
1390
1391
# Read & write cells
1392
1393
# first, read all lines (need to sweep to times through the file)
1394
lines = ifile.readlines()
1395
1396
# second, find the number of cells
1397
num_cells = -1
1398
counter = 0
1399
for line in lines:
1400
l = [int(a) for a in line.split()]
1401
cellnr, node0, node1, node2, node3, node4, node5, node6, node7, tmp1, tmp2 = l
1402
if node4 > 0:
1403
if node2 == node3 and node4 == node5 and node5 == node6 and node6 == node7: # these nodes should be equal
1404
counter += 1
1405
else:
1406
print "The file does contain cells that are not tetraheders. The cell number is ", cellnr, " the line read was ", line
1407
else:
1408
# triangles on the surface
1409
# print "The file does contain cells that are not tetraheders node4==0. The cell number is ", cellnr, " the line read was ", line
1410
#sys.exit(2)
1411
pass
1412
1413
num_cells = counter
1414
1415
# third, run over all cells
1416
write_header_cells(ofile, num_cells)
1417
counter = 0
1418
for line in lines:
1419
l = [int(a) for a in line.split()]
1420
cellnr, node0, node1, node2, node3, node4, node5, node6, node7, tmp1, tmp2 = l
1421
if (node4 > 0):
1422
if node2 == node3 and node4 == node5 and node5 == node6 and node6 == node7: # these nodes should be equal
1423
1424
write_cell_tetrahedron(ofile, counter, nodenr_map[node0], nodenr_map[node1], nodenr_map[node2], nodenr_map[node4])
1425
counter += 1
1426
1427
1428
write_footer_cells(ofile)
1429
write_footer_mesh(ofile)
1430
1431
1432
# Close files
1433
ifile.close()
1434
ofile.close()
Loggerhead is a web-based interface for Breezy
Version: 2.0.1