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