forked from nschloe/meshio
-
Notifications
You must be signed in to change notification settings - Fork 0
/
medit_io.py
163 lines (136 loc) · 5.11 KB
/
medit_io.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
# -*- coding: utf-8 -*-
#
'''
I/O for Medit's format, cf.
<https://people.sc.fsu.edu/~jburkardt/data/medit/medit.html>.
Check out
<https://hal.inria.fr/inria-00069921/fr/>
<https://www.ljll.math.upmc.fr/frey/publications/RT-0253.pdf>
for something like a specification.
.. moduleauthor:: Nico Schlömer <nico.schloemer@gmail.com>
'''
import re
import logging
import numpy
def read(filename):
with open(filename) as f:
points, cells = read_buffer(f)
return points, cells, {}, {}, {}
class _ItemReader:
def __init__(self, file, delimiter=r'\s+'):
# Items can be separated by any whitespace, including new lines.
self._re_delimiter = re.compile(delimiter, re.MULTILINE)
self._file = file
self._line = []
self._line_ptr = 0
def next_items(self, n):
'''Returns the n next items.
Throws StopIteration when there is not enough data to return n items.
'''
items = []
while len(items) < n:
if self._line_ptr >= len(self._line):
# Load the next line.
line = next(self._file).strip()
# Skip all comment and empty lines.
while not line or line[0] == '#':
line = next(self._file).strip()
self._line = self._re_delimiter.split(line)
self._line_ptr = 0
n_read = min(n - len(items), len(self._line) - self._line_ptr)
items.extend(self._line[self._line_ptr:self._line_ptr + n_read])
self._line_ptr += n_read
return items
def next_item(self):
return self.next_items(1)[0]
def read_buffer(file):
dim = 0
cells = {}
reader = _ItemReader(file)
while True:
try:
keyword = reader.next_item()
except StopIteration:
break
assert keyword.isalpha()
meshio_from_medit = {
'Edges': ('line', 2),
'Triangles': ('triangle', 3),
'Quadrilaterals': ('quad', 4),
'Tetrahedra': ('tetra', 4),
'Hexahedra': ('hexahedra', 8)
}
if keyword == 'MeshVersionFormatted':
assert reader.next_item() == '1'
elif keyword == 'Dimension':
dim = int(reader.next_item())
elif keyword == 'Vertices':
assert dim > 0
# The first value is the number of nodes
num_verts = int(reader.next_item())
points = numpy.empty((num_verts, dim), dtype=float)
for k in range(num_verts):
# Throw away the label immediately
points[k] = numpy.array(
reader.next_items(dim + 1), dtype=float)[:-1]
elif keyword in meshio_from_medit:
meshio_name, num = meshio_from_medit[keyword]
# The first value is the number of elements
num_cells = int(reader.next_item())
cell_data = numpy.empty((num_cells, num), dtype=int)
for k in range(num_cells):
data = numpy.array(reader.next_items(num + 1), dtype=int)
# Throw away the label
cell_data[k] = data[:-1]
# adapt 0-base
cells[meshio_name] = cell_data - 1
else:
assert keyword == 'End', 'Unknown keyword \'{}\'.'.format(keyword)
return points, cells
def write(filename,
points,
cells,
point_data=None,
cell_data=None,
field_data=None):
point_data = {} if point_data is None else point_data
cell_data = {} if cell_data is None else cell_data
field_data = {} if field_data is None else field_data
with open(filename, 'wb') as fh:
fh.write(b'MeshVersionFormatted 1\n')
fh.write(b'# Created by meshio\n')
# Dimension info
d = '\nDimension {}\n'.format(points.shape[1])
fh.write(d.encode('utf-8'))
# vertices
fh.write(b'\nVertices\n')
fh.write('{}\n'.format(len(points)).encode('utf-8'))
labels = numpy.ones(len(points), dtype=int)
data = numpy.c_[points, labels]
fmt = ' '.join(['%r'] * points.shape[1]) + ' %d'
numpy.savetxt(fh, data, fmt)
medit_from_meshio = {
'line': ('Edges', 2),
'triangle': ('Triangles', 3),
'quad': ('Quadrilaterals', 4),
'tetra': ('Tetrahedra', 4),
'hexahedra': ('Hexahedra', 8)
}
for key, data in cells.items():
try:
medit_name, num = medit_from_meshio[key]
except KeyError:
msg = ('MEDIT\'s mesh format doesn\'t know {} cells. Skipping.'
).format(key)
logging.warning(msg)
continue
fh.write(b'\n')
fh.write('{}\n'.format(medit_name).encode('utf-8'))
fh.write('{}\n'.format(len(data)).encode('utf-8'))
labels = numpy.ones(len(data), dtype=int)
# adapt 1-base
data_with_label = numpy.c_[data + 1, labels]
fmt = ' '.join(['%d'] * (num + 1))
numpy.savetxt(fh, data_with_label, fmt)
fh.write(b'\nEnd\n')
return