Package meshpy :: Module common
[hide private]
[frames] | no frames]

Source Code for Module meshpy.common

1 -class _Table:
2 - def __init__(self):
3 self.Rows = []
4
5 - def add_row(self, row):
6 self.Rows.append([str(i) for i in row])
7
8 - def __str__(self):
9 columns = len(self.Rows[0]) 10 col_widths = [max(len(row[i]) for row in self.Rows) 11 for i in range(columns)] 12 13 lines = [ 14 " ".join([cell.ljust(col_width) 15 for cell, col_width in zip(row, col_widths)]) 16 for row in self.Rows] 17 return "\n".join(lines)
18
19 20 21 22 -def _linebreak_list(list, per_line=10, pad=None):
23 def format(s): 24 if pad is None: 25 return str(s) 26 else: 27 return str(s).rjust(pad)
28 29 result = "" 30 while len(list) > per_line: 31 result += " ".join(format(l) for l in list[:per_line]) + "\n" 32 list = list[per_line:] 33 return result + " ".join(format(l) for l in list) 34
35 36 37 -class MeshInfoBase:
38 @property
40 try: 41 return self._fvi2fm 42 except AttributeError: 43 result = {} 44 45 for i, face in enumerate(self.faces): 46 result[frozenset(face)] = self.face_markers[i] 47 48 self._fvi2fm = result 49 return result
50 51 52 53 54
55 - def set_points(self, points, point_markers=None):
56 if point_markers is not None: 57 assert len(point_markers) == len(point_markers) 58 59 self.points.resize(len(points)) 60 61 for i, pt in enumerate(points): 62 self.points[i] = pt 63 64 if point_markers is not None: 65 self.point_markers.setup() 66 for i, mark in enumerate(point_markers): 67 self.point_markers[i] = mark
68 69 70 71 72
73 - def set_holes(self, hole_starts):
74 self.holes.resize(len(hole_starts)) 75 for i, hole in enumerate(hole_starts): 76 self.holes[i] = hole
77 78 79 80
81 - def write_neu(self, outfile, bc={}, periodicity=None, description="MeshPy Output"):
82 """Write the mesh out in (an approximation to) Gambit neutral mesh format. 83 84 outfile is a file-like object opened for writing. 85 86 bc is a dictionary mapping single face markers (or frozensets of them) 87 to a tuple (bc_name, bc_code). 88 89 periodicity is either a tuple (face_marker, (px,py,..)) giving the 90 face marker of the periodic boundary and the period in each coordinate 91 direction (0 if none) or the value None for no periodicity. 92 """ 93 94 from meshpy import version 95 from datetime import datetime 96 97 # header -------------------------------------------------------------- 98 outfile.write("CONTROL INFO 2.1.2\n") 99 outfile.write("** GAMBIT NEUTRAL FILE\n") 100 outfile.write("%s\n" % description) 101 outfile.write("PROGRAM: MeshPy VERSION: %s\n" % version) 102 outfile.write("%s\n" % datetime.now().ctime()) 103 104 bc_markers = bc.keys() 105 if periodicity: 106 periodic_marker, periods = periodicity 107 bc_markers.append(periodic_marker) 108 109 assert len(self.points) 110 111 dim = len(self.points[0]) 112 data = ( 113 ("NUMNP", len(self.points)), 114 ("NELEM", len(self.elements)), 115 ("NGRPS", 1), 116 ("NBSETS", len(bc_markers)), 117 ("NDFCD", dim), 118 ("NDFVL", dim), 119 ) 120 121 tbl = _Table() 122 tbl.add_row(key for key, value in data) 123 tbl.add_row(value for key, value in data) 124 outfile.write(str(tbl)) 125 outfile.write("\n") 126 outfile.write("ENDOFSECTION\n") 127 128 # nodes --------------------------------------------------------------- 129 outfile.write("NODAL COORDINATES 2.1.2\n") 130 for i, pt in enumerate(self.points): 131 outfile.write("%d %s\n" % 132 (i+1, " ".join(repr(c) for c in pt))) 133 outfile.write("ENDOFSECTION\n") 134 135 # elements ------------------------------------------------------------ 136 outfile.write("ELEMENTS/CELLS 2.1.2\n") 137 if dim == 2: 138 eltype = 3 139 else: 140 eltype = 6 141 142 for i, el in enumerate(self.elements): 143 outfile.write("%8d%3d%3d %s\n" % 144 (i+1, eltype, len(el), 145 "".join("%8d" % (p+1) for p in el))) 146 outfile.write("ENDOFSECTION\n") 147 148 # element groups ------------------------------------------------------ 149 outfile.write("ELEMENT GROUP 1.3.0\n") 150 # FIXME 151 i = 0 152 grp_elements = range(len(self.elements)) 153 material = 1 154 flags = 0 155 outfile.write("GROUP:%11d ELEMENTS:%11d MATERIAL:%11s NFLAGS: %11d\n" 156 % (1, len(grp_elements), repr(material), flags)) 157 outfile.write(("epsilon: %s\n" % material).rjust(32)) # FIXME 158 outfile.write("0\n") 159 outfile.write(_linebreak_list([str(i+1) for i in grp_elements], 160 pad=8) 161 + "\n") 162 outfile.write("ENDOFSECTION\n") 163 164 # boundary conditions ------------------------------------------------- 165 # build mapping face -> (tet, neu_face_index) 166 face2el = {} 167 168 if dim == 2: 169 for ti, el in enumerate(self.elements): 170 # Sledge++ Users' Guide, figure 4 171 faces = [ 172 frozenset([el[0], el[1]]), 173 frozenset([el[1], el[2]]), 174 frozenset([el[2], el[0]]), 175 ] 176 for fi, face in enumerate(faces): 177 face2el.setdefault(face, []).append((ti, fi+1)) 178 179 elif dim == 3: 180 face2el = {} 181 for ti, el in enumerate(self.elements): 182 # Sledge++ Users' Guide, figure 5 183 faces = [ 184 frozenset([el[1], el[0], el[2]]), 185 frozenset([el[0], el[1], el[3]]), 186 frozenset([el[1], el[2], el[3]]), 187 frozenset([el[2], el[0], el[3]]), 188 ] 189 for fi, face in enumerate(faces): 190 face2el.setdefault(face, []).append((ti, fi+1)) 191 192 else: 193 raise ValueError, "invalid number of dimensions (%d)" % dim 194 195 # actually output bc sections 196 if not self.faces.allocated: 197 from warnings import warn 198 warn("no exterior faces in mesh data structure, not writing boundary conditions") 199 else: 200 # requires -f option in tetgen, -e in triangle 201 202 for bc_marker in bc_markers: 203 if isinstance(bc_marker, frozenset): 204 face_indices = [i 205 for i, face in enumerate(self.faces) 206 if self.face_markers[i] in bc_marker] 207 else: 208 face_indices = [i 209 for i, face in enumerate(self.faces) 210 if bc_marker == self.face_markers[i]] 211 212 if not face_indices: 213 continue 214 215 outfile.write("BOUNDARY CONDITIONS 2.1.2\n") 216 if bc_marker in bc: 217 # regular BC 218 219 bc_name, bc_code = bc[bc_marker] 220 outfile.write("%32s%8d%8d%8d%8d\n" 221 % (bc_name, 222 1, # face BC 223 len(face_indices), 224 0, # zero additional values per face, 225 bc_code, 226 ) 227 ) 228 else: 229 # periodic BC 230 231 outfile.write("%s%s%8d%8d%8d\n" 232 % ("periodic", " ".join(repr(p) for p in periods), 233 len(face_indices), 234 0, # zero additional values per face, 235 0, 236 ) 237 ) 238 239 for i, fi in enumerate(face_indices): 240 face_nodes = frozenset(self.faces[fi]) 241 adj_el = face2el[face_nodes] 242 assert len(adj_el) == 1 243 244 el_index, el_face_number = adj_el[0] 245 246 outfile.write("%10d%5d%5d\n" % 247 (el_index+1, eltype, el_face_number)) 248 249 outfile.write("ENDOFSECTION\n") 250 251 outfile.close()
252 # FIXME curved boundaries?
253 # FIXME proper element group support 254 255 256 257 258 259 -def dump_array(name, array):
260 print "array %s: %d elements, %d values per element" % (name, len(array), array.unit) 261 262 if len(array) == 0 or array.unit == 0: 263 return 264 265 try: 266 array[0] 267 except RuntimeError: 268 print " not allocated" 269 return 270 271 for i, entry in enumerate(array): 272 if isinstance(entry, list): 273 print " %d: %s" % (i, ",".join(str(sub) for sub in entry)) 274 else: 275 print " %d: %s" % (i, entry)
276