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

Source Code for Module meshpy.gmsh_reader

  1  """Reader for the GMSH file format.""" 
  2   
  3  from __future__ import division 
  4   
  5  __copyright__ = "Copyright (C) 2009 Xueyu Zhu, Andreas Kloeckner" 
  6   
  7  __license__ = """ 
  8  Permission is hereby granted, free of charge, to any person obtaining a copy 
  9  of this software and associated documentation files (the "Software"), to deal 
 10  in the Software without restriction, including without limitation the rights 
 11  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 
 12  copies of the Software, and to permit persons to whom the Software is 
 13  furnished to do so, subject to the following conditions: 
 14   
 15  The above copyright notice and this permission notice shall be included in 
 16  all copies or substantial portions of the Software. 
 17   
 18  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 
 19  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 
 20  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 
 21  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 
 22  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 
 23  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 
 24  THE SOFTWARE. 
 25  """ 
 26   
 27  import numpy as np 
 28  import numpy.linalg as la 
 29  from pytools import memoize_method, Record, single_valued 
30 31 32 # {{{ tools 33 34 -def generate_triangle_vertex_tuples(order):
35 yield (0, 0) 36 yield (order, 0) 37 yield (0, order)
38
39 -def generate_triangle_edge_tuples(order):
40 for i in range(1, order): 41 yield (i, 0) 42 for i in range(1, order): 43 yield (order-i, i) 44 for i in range(1, order): 45 yield (0, order-i)
46
47 -def generate_triangle_volume_tuples(order):
48 for i in range(1, order): 49 for j in range(1, order-i): 50 yield (j, i)
51
52 -class LineFeeder:
53 - def __init__(self, line_iterable):
54 self.line_iterable = iter(line_iterable) 55 self.next_line = None
56
57 - def has_next_line(self):
58 if self.next_line is not None: 59 return True 60 61 try: 62 self.next_line = self.line_iterable.next() 63 except StopIteration: 64 return False 65 else: 66 return True
67
68 - def get_next_line(self):
69 if self.next_line is not None: 70 nl = self.next_line 71 self.next_line = None 72 return nl.strip() 73 74 try: 75 nl = self.line_iterable.next() 76 except StopIteration: 77 raise GmshFileFormatError("unexpected end of file") 78 else: 79 return nl.strip()
80
81 # }}} 82 83 # {{{ element info 84 85 -class GmshElementBase(object):
86 - def __init__(self, order):
87 self.order = order
88
89 - def vertex_count(self):
90 return self.dimensions + 1
91 92 @memoize_method
93 - def node_count(self):
94 """Return the number of interpolation nodes in this element.""" 95 d = self.dimensions 96 o = self.order 97 from operator import mul 98 from pytools import factorial 99 return int(reduce(mul, (o + 1 + i for i in range(d)), 1) / factorial(d))
100 101 @memoize_method
103 """Generate tuples enumerating the node indices present 104 in this element. Each tuple has a length equal to the dimension 105 of the element. The tuples constituents are non-negative integers 106 whose sum is less than or equal to the order of the element. 107 108 The order in which these nodes are generated dictates the local 109 node numbering. 110 """ 111 from pytools import \ 112 generate_nonnegative_integer_tuples_summing_to_at_most 113 result = list( 114 generate_nonnegative_integer_tuples_summing_to_at_most( 115 self.order, self.dimensions)) 116 117 assert len(result) == self.node_count() 118 return result
119 120 @memoize_method
122 gmsh_tup_to_index = dict( 123 (tup, i) 124 for i, tup in enumerate(self.gmsh_node_tuples())) 125 126 return np.array([gmsh_tup_to_index[tup] 127 for tup in self.lexicographic_node_tuples()], 128 dtype=np.intp)
129 130 @memoize_method
131 - def equidistant_vandermonde(self):
132 from hedge.polynomial import generic_vandermonde 133 134 return generic_vandermonde( 135 list(self.equidistant_unit_nodes()), 136 list(self.basis_functions()))
137
138 139 140 141 -class GmshPoint(GmshElementBase):
142 dimensions = 0 143 144 @memoize_method
145 - def gmsh_node_tuples(self):
146 return [()]
147
148 149 150 -class GmshIntervalElement(GmshElementBase):
151 dimensions = 1 152 153 @memoize_method
154 - def gmsh_node_tuples(self):
155 return [(0,), (self.order,),] + [ 156 (i,) for i in range(1, self.order)]
157
158 159 160 161 -class GmshIncompleteTriangularElement(GmshElementBase):
162 dimensions = 2 163
164 - def __init__(self, order):
165 self.order = order
166 167 @memoize_method
168 - def gmsh_node_tuples(self):
169 result = [] 170 for tup in generate_triangle_vertex_tuples(self.order): 171 result.append(tup) 172 for tup in generate_triangle_edge_tuples(self.order): 173 result.append(tup) 174 return result
175
176 177 178 179 180 -class GmshTriangularElement(GmshElementBase):
181 dimensions = 2 182 183 @memoize_method
184 - def gmsh_node_tuples(self):
185 result = [] 186 for tup in generate_triangle_vertex_tuples(self.order): 187 result.append(tup) 188 for tup in generate_triangle_edge_tuples(self.order): 189 result.append(tup) 190 for tup in generate_triangle_volume_tuples(self.order): 191 result.append(tup) 192 return result
193
194 195 196 197 -class GmshTetrahedralElement(GmshElementBase):
198 dimensions = 3 199 200 @memoize_method
201 - def gmsh_node_tuples(self):
202 # gmsh's node ordering is on crack 203 return { 204 1: [(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)], 205 2: [ 206 (0, 0, 0), (2, 0, 0), (0, 2, 0), (0, 0, 2), (1, 0, 0), (1, 1, 0), 207 (0, 1, 0), (0, 0, 1), (0, 1, 1), (1, 0, 1)], 208 3: [ 209 (0, 0, 0), (3, 0, 0), (0, 3, 0), (0, 0, 3), (1, 0, 0), (2, 0, 0), 210 (2, 1, 0), (1, 2, 0), (0, 2, 0), (0, 1, 0), (0, 0, 2), (0, 0, 1), 211 (0, 1, 2), (0, 2, 1), (1, 0, 2), (2, 0, 1), (1, 1, 0), (1, 0, 1), 212 (0, 1, 1), (1, 1, 1)], 213 4: [ 214 (0, 0, 0), (4, 0, 0), (0, 4, 0), (0, 0, 4), (1, 0, 0), (2, 0, 0), 215 (3, 0, 0), (3, 1, 0), (2, 2, 0), (1, 3, 0), (0, 3, 0), (0, 2, 0), 216 (0, 1, 0), (0, 0, 3), (0, 0, 2), (0, 0, 1), (0, 1, 3), (0, 2, 2), 217 (0, 3, 1), (1, 0, 3), (2, 0, 2), (3, 0, 1), (1, 1, 0), (1, 2, 0), 218 (2, 1, 0), (1, 0, 1), (2, 0, 1), (1, 0, 2), (0, 1, 1), (0, 1, 2), 219 (0, 2, 1), (1, 1, 2), (2, 1, 1), (1, 2, 1), (1, 1, 1)], 220 5: [ 221 (0, 0, 0), (5, 0, 0), (0, 5, 0), (0, 0, 5), (1, 0, 0), (2, 0, 0), 222 (3, 0, 0), (4, 0, 0), (4, 1, 0), (3, 2, 0), (2, 3, 0), (1, 4, 0), 223 (0, 4, 0), (0, 3, 0), (0, 2, 0), (0, 1, 0), (0, 0, 4), (0, 0, 3), 224 (0, 0, 2), (0, 0, 1), (0, 1, 4), (0, 2, 3), (0, 3, 2), (0, 4, 1), 225 (1, 0, 4), (2, 0, 3), (3, 0, 2), (4, 0, 1), (1, 1, 0), (1, 3, 0), 226 (3, 1, 0), (1, 2, 0), (2, 2, 0), (2, 1, 0), (1, 0, 1), (3, 0, 1), 227 (1, 0, 3), (2, 0, 1), (2, 0, 2), (1, 0, 2), (0, 1, 1), (0, 1, 3), 228 (0, 3, 1), (0, 1, 2), (0, 2, 2), (0, 2, 1), (1, 1, 3), (3, 1, 1), 229 (1, 3, 1), (2, 1, 2), (2, 2, 1), (1, 2, 2), (1, 1, 1), (2, 1, 1), 230 (1, 2, 1), (1, 1, 2)], 231 }[self.order]
232
233 234 235 236 237 238 # }}} 239 240 # {{{ receiver interface 241 242 -class GmshMeshReceiverBase(object):
243 gmsh_element_type_to_info_map = { 244 1: GmshIntervalElement(1), 245 2: GmshTriangularElement(1), 246 4: GmshTetrahedralElement(1), 247 8: GmshIntervalElement(2), 248 9: GmshTriangularElement(2), 249 11: GmshTetrahedralElement(2), 250 15: GmshPoint(0), 251 20: GmshIncompleteTriangularElement(3), 252 21: GmshTriangularElement(3), 253 22: GmshIncompleteTriangularElement(4), 254 23: GmshTriangularElement(4), 255 24: GmshIncompleteTriangularElement(5), 256 25: GmshTriangularElement(5), 257 26: GmshIntervalElement(3), 258 27: GmshIntervalElement(4), 259 28: GmshIntervalElement(5), 260 29: GmshTetrahedralElement(3), 261 30: GmshTetrahedralElement(4), 262 31: GmshTetrahedralElement(5) 263 } 264
265 - def set_up_nodes(self, count):
266 pass
267
268 - def add_node(self, node_nr, point):
269 pass
270
271 - def finalize_nodes(self):
272 pass
273
274 - def set_up_elements(self, count):
275 pass
276
277 - def add_element(self, element_nr, element_type, vertex_nrs, 278 lexicographic_nodes, tag_numbers):
279 pass
280
281 - def finalize_elements(self):
282 pass
283
284 - def add_tag(self, name, index, dimension):
285 pass
286
287 - def finalize_tags(self):
288 pass
289
290 # }}} 291 292 # {{{ file reader 293 294 -class GmshFileFormatError(RuntimeError):
295 pass
296
297 298 299 300 -def read_gmsh(receiver, filename, force_dimension=None):
301 """ 302 :param force_dimension: if not None, truncate point coordinates to this many dimensions. 303 """ 304 mesh_file = open(filename, 'rt') 305 try: 306 result = parse_gmsh(receiver, mesh_file) 307 finally: 308 mesh_file.close() 309 310 return result
311
312 313 314 315 -def generate_gmsh(receiver, source, dimensions, order=None, other_options=[], 316 extension="geo", gmsh_executable="gmsh", force_dimension=None):
317 from meshpy.gmsh import GmshRunner 318 runner = GmshRunner(source, dimensions, order=order, 319 other_options=other_options, extension=extension, 320 gmsh_executable=gmsh_executable) 321 322 runner.__enter__() 323 try: 324 result = parse_gmsh(receiver, runner.output_file, force_dimension=force_dimension) 325 finally: 326 runner.__exit__(None, None, None) 327 328 return result
329
330 331 332 333 -def parse_gmsh(receiver, line_iterable, force_dimension=None):
334 """ 335 :arg receiver: This object will be fed the entities encountered in reading the 336 GMSH file. See :class:`GmshMeshReceiverBase` for the interface this object needs 337 to conform to. 338 :param force_dimension: if not None, truncate point coordinates to this many 339 dimensions. 340 """ 341 342 feeder = LineFeeder(line_iterable) 343 344 # collect the mesh information 345 346 class ElementInfo(Record): 347 pass
348 349 while feeder.has_next_line(): 350 next_line = feeder.get_next_line() 351 if not next_line.startswith("$"): 352 raise GmshFileFormatError("expected start of section, '%s' found instead" % next_line) 353 354 section_name = next_line[1:] 355 356 if section_name == "MeshFormat": 357 line_count = 0 358 while True: 359 next_line = feeder.get_next_line() 360 if next_line == "$End"+section_name: 361 break 362 363 if line_count == 0: 364 version_number, file_type, data_size = next_line.split() 365 366 if line_count > 0: 367 raise GmshFileFormatError("more than one line found in MeshFormat section") 368 369 if version_number not in ["2.1", "2.2"]: 370 from warnings import warn 371 warn("unexpected mesh version number '%s' found" % version_number) 372 373 if file_type != "0": 374 raise GmshFileFormatError("only ASCII gmsh file type is supported") 375 376 line_count += 1 377 378 elif section_name == "Nodes": 379 node_count = int(feeder.get_next_line()) 380 receiver.set_up_nodes(node_count) 381 382 node_idx = 1 383 384 while True: 385 next_line = feeder.get_next_line() 386 if next_line == "$End"+section_name: 387 break 388 389 parts = next_line.split() 390 if len(parts) != 4: 391 raise GmshFileFormatError("expected four-component line in $Nodes section") 392 393 read_node_idx = int(parts[0]) 394 if read_node_idx != node_idx: 395 raise GmshFileFormatError("out-of-order node index found") 396 397 if force_dimension is not None: 398 point = [float(x) for x in parts[1:force_dimension+1]] 399 else: 400 point = [float(x) for x in parts[1:]] 401 402 receiver.add_node( 403 node_idx-1, 404 np.array(point, dtype=np.float64)) 405 406 node_idx += 1 407 408 if node_count+1 != node_idx: 409 raise GmshFileFormatError("unexpected number of nodes found") 410 411 receiver.finalize_nodes() 412 413 elif section_name == "Elements": 414 element_count = int(feeder.get_next_line()) 415 receiver.set_up_elements(element_count) 416 417 element_idx = 1 418 while True: 419 next_line = feeder.get_next_line() 420 if next_line == "$End"+section_name: 421 break 422 423 parts = [int(x) for x in next_line.split()] 424 425 if len(parts) < 4: 426 raise GmshFileFormatError("too few entries in element line") 427 428 read_element_idx = parts[0] 429 if read_element_idx != element_idx: 430 raise GmshFileFormatError("out-of-order node index found") 431 432 el_type_num = parts[1] 433 try: 434 element_type = receiver.gmsh_element_type_to_info_map[el_type_num] 435 except KeyError: 436 raise GmshFileFormatError("unexpected element type %d" 437 % el_type_num) 438 439 tag_count = parts[2] 440 tags = parts[3:3+tag_count] 441 442 # convert to zero-based 443 node_indices = np.array([x-1 for x in parts[3+tag_count:]], dtype=np.intp) 444 445 if element_type.node_count()!= len(node_indices): 446 raise GmshFileFormatError("unexpected number of nodes in element") 447 448 gmsh_vertex_nrs = node_indices[:element_type.vertex_count()] 449 zero_based_idx = element_idx - 1 450 451 tag_numbers = [tag for tag in tags[:1] if tag != 0] 452 453 receiver.add_element(element_nr=zero_based_idx, 454 element_type=element_type, vertex_nrs=gmsh_vertex_nrs, 455 lexicographic_nodes=node_indices[ 456 element_type.get_lexicographic_gmsh_node_indices()], 457 tag_numbers=tag_numbers) 458 459 element_idx +=1 460 461 if element_count+1 != element_idx: 462 raise GmshFileFormatError("unexpected number of elements found") 463 464 receiver.finalize_elements() 465 466 elif section_name == "PhysicalNames": 467 name_count = int(feeder.get_next_line()) 468 name_idx = 1 469 470 while True: 471 next_line = feeder.get_next_line() 472 if next_line == "$End"+section_name: 473 break 474 475 dimension, number, name = next_line.split(" ", 2) 476 dimension = int(dimension) 477 number = int(number) 478 479 if not name[0] == '"' or not name[-1] == '"': 480 raise GmshFileFormatError("expected quotes around physical name") 481 482 receiver.add_tag(name[1:-1], number, dimension) 483 484 name_idx += 1 485 486 if name_count+1 != name_idx: 487 raise GmshFileFormatError("unexpected number of physical names found") 488 489 receiver.finalize_tags() 490 else: 491 # unrecognized section, skip 492 from warnings import warn 493 warn("unrecognized section '%s' in gmsh file" % section_name) 494 while True: 495 next_line = feeder.get_next_line() 496 if next_line == "$End"+section_name: 497 break 498 499 # }}} 500 501 # vim: fdm=marker 502