forked from nschloe/meshio
/
permas_io.py
211 lines (192 loc) · 6.4 KB
/
permas_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
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
# -*- coding: utf-8 -*-
#
'''
I/O for PERMAS dat format, cf.
.. moduleauthor:: Nils Wagner
'''
import gzip
import re
import numpy
from .__about__ import __version__, __website__
def read(filename):
'''Reads a (compressed) PERMAS dato or post file.
'''
# The format is specified at
# <http://www.intes.de>.
cells = {}
meshio_to_permas_type = {
'vertex': (1, 'PLOT1'),
'line': (2, 'PLOTL2'),
'triangle': (3, 'TRIA3'),
'quad': (4, 'QUAD4'),
'tetra': (4, 'TET4'),
'hexahedron': (8, 'HEXE8'),
'wedge': (6, 'PENTA6'),
'pyramid': (5, 'PYRA5')
}
if filename.endswith('.gz') or filename.endswith('.gz'):
opener = gzip.open
else:
# assert filename.endswith('dato') or filename.endswith('post')
opener = open
with opener(filename, 'r') as f:
while True:
line = f.readline()
if not line or re.search('\\$END STRUCTURE', line):
break
for meshio_type, permas_ele in meshio_to_permas_type.items():
num_nodes = permas_ele[0]
permas_type = permas_ele[1]
if re.search('\\$ELEMENT TYPE = {}'.format(permas_type), line):
while True:
line = f.readline()
if not line or line.startswith('!'):
break
data = numpy.array(line.split(), dtype=int)
if meshio_type in cells:
cells[meshio_type].append(data[-num_nodes:])
else:
cells[meshio_type] = [data[-num_nodes:]]
if re.search('\\$COOR', line):
points = []
while True:
line = f.readline()
if not line or line.startswith('!'):
break
for r in numpy.array(line.split(), dtype=float)[1:]:
points.append(r)
points = numpy.array(points)
points = numpy.reshape(points, newshape=(len(points)//3, 3))
for key in cells:
# Subtract one to account for the fact that python indices
# are 0-based.
cells[key] = numpy.array(cells[key], dtype=int) - 1
return points, cells, {}, {}, {}
def write(
filename,
points,
cells,
point_data=None,
cell_data=None,
field_data=None
):
'''Writes PERMAS dat files, cf.
http://www.intes.de # PERMAS-ASCII-file-format
'''
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, 'w') as fh:
fh.write('!\n')
fh.write('! File written by meshio version {}\n'.format(__version__))
fh.write('! Further information available at {}\n'.format(__website__))
fh.write('!\n')
fh.write('$ENTER COMPONENT NAME = DFLT_COMP DOFTYPE = DISP MATH\n')
fh.write('! \n')
fh.write(' $SITUATION NAME = REAL_MODES\n')
fh.write(
' DFLT_COMP SYSTEM ='
' NSV CONSTRAINTS = SPCVAR_1 ! LOADING = LOADVAR_1\n'
)
fh.write(' $END SITUATION\n')
fh.write('! \n')
fh.write(' $STRUCTURE\n')
fh.write('! \n')
# Write nodes
fh.write(' $COOR NSET = ALL_NODES\n')
for k, x in enumerate(points):
fh.write(
' {:8d} {:+.15f} {:+.15f} {:+.15f}\n'.format(
k+1, x[0], x[1], x[2]
))
meshio_to_permas_type = {
'vertex': (1, 'PLOT1'),
'line': (2, 'PLOTL2'),
'triangle': (3, 'TRIA3'),
'quad': (4, 'QUAD4'),
'tetra': (4, 'TET4'),
'hexahedron': (8, 'HEXE8'),
'wedge': (6, 'PENTA6'),
'pyramid': (5, 'PYRA5')
}
#
# Avoid non-unique element numbers in case of multiple element types by
# num_ele !!!
#
num_ele = 0
for meshio_type, cell in cells.items():
numcells, num_local_nodes = cell.shape
permas_type = meshio_to_permas_type[meshio_type]
fh.write('!\n')
fh.write(
' $ELEMENT TYPE = {} ESET = {}\n'.format(
permas_type[1], permas_type[1]
))
for k, c in enumerate(cell):
form = ' %8d ' + \
' '.join(num_local_nodes * ['%8d']) + \
'\n'
fh.write(form % ((k+num_ele+1,) + tuple(c + 1)))
num_ele += numcells
fh.write('!\n')
fh.write(' $END STRUCTURE\n')
fh.write('!\n')
elem_3D = ['HEXE8', 'TET4', 'PENTA6', 'PYRA5']
elem_2D = ['TRIA3', 'QUAD4']
elem_1D = ['PLOT1', 'PLOTL2']
fh.write(' $SYSTEM NAME = NSV\n')
fh.write('!\n')
fh.write(' $ELPROP\n')
for meshio_type, cell in cells.items():
permas_type = meshio_to_permas_type[meshio_type]
if permas_type[1] in elem_3D:
fh.write(
' {} MATERIAL = DUMMY_MATERIAL\n'.format(
permas_type[1]
))
elif permas_type[1] in elem_2D:
fh.write(
12 * ' ' +
'{} GEODAT = GD_{} MATERIAL = DUMMY_MATERIAL\n'.format(
permas_type[1], permas_type[1]
))
else:
assert permas_type[1] in elem_1D
fh.write('!\n')
fh.write(' $GEODAT SHELL CONT = THICK NODES = ALL\n')
for meshio_type, cell in cells.items():
permas_type = meshio_to_permas_type[meshio_type]
if permas_type[1] in elem_2D:
fh.write(12 * ' ' + 'GD_{} 1.0\n'.format(permas_type[1]))
fh.write('''!
!
$END SYSTEM
!
$CONSTRAINTS NAME = SPCVAR_1
$END CONSTRAINTS
!
$LOADING NAME = LOADVAR_1
$END LOADING
!
$EXIT COMPONENT
!
$ENTER MATERIAL
!
$MATERIAL NAME = DUMMY_MATERIAL TYPE = ISO
!
$ELASTIC GENERAL INPUT = DATA
2.0E+05 0.3
!
$DENSITY GENERAL INPUT = DATA
7.8E-09
!
$THERMEXP GENERAL INPUT = DATA
1.2E-05
!
$END MATERIAL
!
$EXIT MATERIAL
!
$FIN
''')
return