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 
 30  from meshpy.gmsh import LiteralSource, FileSource  # noqa 
 31   
 32   
 33  __doc__ = """ 
 34  .. exception:: GmshFileFormatError 
 35   
 36  Element types 
 37  ------------- 
 38   
 39  .. autoclass:: GmshElementBase 
 40  .. autoclass:: GmshPoint 
 41  .. autoclass:: GmshIntervalElement 
 42  .. autoclass:: GmshTriangularElement 
 43  .. autoclass:: GmshIncompleteTriangularElement 
 44  .. autoclass:: GmshTetrahedralElement 
 45   
 46  Receiver interface 
 47  ------------------ 
 48   
 49  .. autoclass:: GmshMeshReceiverBase 
 50   
 51  Receiver example implementation 
 52  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 
 53   
 54  .. autoclass:: GmshMeshReceiverNumPy 
 55   
 56  Reader 
 57  ------ 
 58   
 59  .. autoclass:: LiteralSource 
 60  .. autoclass:: FileSource 
 61   
 62  .. autofunction:: read_gmsh 
 63  .. autofunction:: generate_gmsh 
 64   
 65  """ 
66 67 68 # {{{ tools 69 70 -def generate_triangle_vertex_tuples(order):
71 yield (0, 0) 72 yield (order, 0) 73 yield (0, order)
74
75 76 -def generate_triangle_edge_tuples(order):
77 for i in range(1, order): 78 yield (i, 0) 79 for i in range(1, order): 80 yield (order-i, i) 81 for i in range(1, order): 82 yield (0, order-i)
83
84 85 -def generate_triangle_volume_tuples(order):
86 for i in range(1, order): 87 for j in range(1, order-i): 88 yield (j, i)
89
90 91 -class LineFeeder:
92 - def __init__(self, line_iterable):
93 self.line_iterable = iter(line_iterable) 94 self.next_line = None
95
96 - def has_next_line(self):
97 if self.next_line is not None: 98 return True 99 100 try: 101 self.next_line = self.line_iterable.next() 102 except StopIteration: 103 return False 104 else: 105 return True
106
107 - def get_next_line(self):
108 if self.next_line is not None: 109 nl = self.next_line 110 self.next_line = None 111 return nl.strip() 112 113 try: 114 nl = self.line_iterable.next() 115 except StopIteration: 116 raise GmshFileFormatError("unexpected end of file") 117 else: 118 return nl.strip()
119
120 # }}} 121 122 123 # {{{ element info 124 125 -class GmshElementBase(object):
126 """ 127 .. automethod:: vertex_count 128 .. automethod:: node_count 129 .. automethod:: get_lexicographic_gmsh_node_indices 130 .. automethod:: equidistant_vandermonde 131 .. method:: equidistant_unit_nodes 132 133 (Implemented by subclasses) 134 """ 135
136 - def __init__(self, order):
137 self.order = order
138
139 - def vertex_count(self):
140 return self.dimensions + 1
141 142 @memoize_method
143 - def node_count(self):
144 """Return the number of interpolation nodes in this element.""" 145 d = self.dimensions 146 o = self.order 147 from operator import mul 148 from pytools import factorial 149 return int(reduce(mul, (o + 1 + i for i in range(d)), 1) / factorial(d))
150 151 @memoize_method
153 """Generate tuples enumerating the node indices present 154 in this element. Each tuple has a length equal to the dimension 155 of the element. The tuples constituents are non-negative integers 156 whose sum is less than or equal to the order of the element. 157 """ 158 from pytools import \ 159 generate_nonnegative_integer_tuples_summing_to_at_most 160 result = list( 161 generate_nonnegative_integer_tuples_summing_to_at_most( 162 self.order, self.dimensions)) 163 164 assert len(result) == self.node_count() 165 return result
166 167 @memoize_method
169 gmsh_tup_to_index = dict( 170 (tup, i) 171 for i, tup in enumerate(self.gmsh_node_tuples())) 172 173 return np.array([gmsh_tup_to_index[tup] 174 for tup in self.lexicographic_node_tuples()], 175 dtype=np.intp)
176 177 @memoize_method
178 - def equidistant_vandermonde(self):
179 from hedge.polynomial import generic_vandermonde 180 181 return generic_vandermonde( 182 list(self.equidistant_unit_nodes()), 183 list(self.basis_functions()))
184
185 186 -class GmshPoint(GmshElementBase):
187 dimensions = 0 188 189 @memoize_method
190 - def gmsh_node_tuples(self):
191 return [()]
192
193 194 -class GmshIntervalElement(GmshElementBase):
195 dimensions = 1 196 197 @memoize_method
198 - def gmsh_node_tuples(self):
199 return [(0,), (self.order,), ] + [ 200 (i,) for i in range(1, self.order)]
201
202 203 -class GmshIncompleteTriangularElement(GmshElementBase):
204 dimensions = 2 205
206 - def __init__(self, order):
207 self.order = order
208 209 @memoize_method
210 - def gmsh_node_tuples(self):
211 result = [] 212 for tup in generate_triangle_vertex_tuples(self.order): 213 result.append(tup) 214 for tup in generate_triangle_edge_tuples(self.order): 215 result.append(tup) 216 return result
217
218 219 -class GmshTriangularElement(GmshElementBase):
220 dimensions = 2 221 222 @memoize_method
223 - def gmsh_node_tuples(self):
224 result = [] 225 for tup in generate_triangle_vertex_tuples(self.order): 226 result.append(tup) 227 for tup in generate_triangle_edge_tuples(self.order): 228 result.append(tup) 229 for tup in generate_triangle_volume_tuples(self.order): 230 result.append(tup) 231 return result
232
233 234 -class GmshTetrahedralElement(GmshElementBase):
235 dimensions = 3 236 237 @memoize_method
238 - def gmsh_node_tuples(self):
239 # gmsh's node ordering is on crack 240 return { 241 1: [(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)], 242 2: [ 243 (0, 0, 0), (2, 0, 0), (0, 2, 0), (0, 0, 2), (1, 0, 0), (1, 1, 0), 244 (0, 1, 0), (0, 0, 1), (0, 1, 1), (1, 0, 1)], 245 3: [ 246 (0, 0, 0), (3, 0, 0), (0, 3, 0), (0, 0, 3), (1, 0, 0), (2, 0, 0), 247 (2, 1, 0), (1, 2, 0), (0, 2, 0), (0, 1, 0), (0, 0, 2), (0, 0, 1), 248 (0, 1, 2), (0, 2, 1), (1, 0, 2), (2, 0, 1), (1, 1, 0), (1, 0, 1), 249 (0, 1, 1), (1, 1, 1)], 250 4: [ 251 (0, 0, 0), (4, 0, 0), (0, 4, 0), (0, 0, 4), (1, 0, 0), (2, 0, 0), 252 (3, 0, 0), (3, 1, 0), (2, 2, 0), (1, 3, 0), (0, 3, 0), (0, 2, 0), 253 (0, 1, 0), (0, 0, 3), (0, 0, 2), (0, 0, 1), (0, 1, 3), (0, 2, 2), 254 (0, 3, 1), (1, 0, 3), (2, 0, 2), (3, 0, 1), (1, 1, 0), (1, 2, 0), 255 (2, 1, 0), (1, 0, 1), (2, 0, 1), (1, 0, 2), (0, 1, 1), (0, 1, 2), 256 (0, 2, 1), (1, 1, 2), (2, 1, 1), (1, 2, 1), (1, 1, 1)], 257 5: [ 258 (0, 0, 0), (5, 0, 0), (0, 5, 0), (0, 0, 5), (1, 0, 0), (2, 0, 0), 259 (3, 0, 0), (4, 0, 0), (4, 1, 0), (3, 2, 0), (2, 3, 0), (1, 4, 0), 260 (0, 4, 0), (0, 3, 0), (0, 2, 0), (0, 1, 0), (0, 0, 4), (0, 0, 3), 261 (0, 0, 2), (0, 0, 1), (0, 1, 4), (0, 2, 3), (0, 3, 2), (0, 4, 1), 262 (1, 0, 4), (2, 0, 3), (3, 0, 2), (4, 0, 1), (1, 1, 0), (1, 3, 0), 263 (3, 1, 0), (1, 2, 0), (2, 2, 0), (2, 1, 0), (1, 0, 1), (3, 0, 1), 264 (1, 0, 3), (2, 0, 1), (2, 0, 2), (1, 0, 2), (0, 1, 1), (0, 1, 3), 265 (0, 3, 1), (0, 1, 2), (0, 2, 2), (0, 2, 1), (1, 1, 3), (3, 1, 1), 266 (1, 3, 1), (2, 1, 2), (2, 2, 1), (1, 2, 2), (1, 1, 1), (2, 1, 1), 267 (1, 2, 1), (1, 1, 2)], 268 }[self.order]
269
270 # }}} 271 272 273 # {{{ receiver interface 274 275 -class GmshMeshReceiverBase(object):
276 """ 277 .. attribute:: gmsh_element_type_to_info_map 278 .. automethod:: set_up_nodes 279 .. automethod:: add_node 280 .. automethod:: finalize_nodes 281 .. automethod:: set_up_elements 282 .. automethod:: add_element 283 .. automethod:: finalize_elements 284 .. automethod:: add_tag 285 .. automethod:: finalize_tags 286 """ 287 288 gmsh_element_type_to_info_map = { 289 1: GmshIntervalElement(1), 290 2: GmshTriangularElement(1), 291 4: GmshTetrahedralElement(1), 292 8: GmshIntervalElement(2), 293 9: GmshTriangularElement(2), 294 11: GmshTetrahedralElement(2), 295 15: GmshPoint(0), 296 20: GmshIncompleteTriangularElement(3), 297 21: GmshTriangularElement(3), 298 22: GmshIncompleteTriangularElement(4), 299 23: GmshTriangularElement(4), 300 24: GmshIncompleteTriangularElement(5), 301 25: GmshTriangularElement(5), 302 26: GmshIntervalElement(3), 303 27: GmshIntervalElement(4), 304 28: GmshIntervalElement(5), 305 29: GmshTetrahedralElement(3), 306 30: GmshTetrahedralElement(4), 307 31: GmshTetrahedralElement(5) 308 } 309
310 - def set_up_nodes(self, count):
311 pass
312
313 - def add_node(self, node_nr, point):
314 pass
315
316 - def finalize_nodes(self):
317 pass
318
319 - def set_up_elements(self, count):
320 pass
321
322 - def add_element(self, element_nr, element_type, vertex_nrs, 323 lexicographic_nodes, tag_numbers):
324 pass
325
326 - def finalize_elements(self):
327 pass
328
329 - def add_tag(self, name, index, dimension):
330 pass
331
332 - def finalize_tags(self):
333 pass
334
335 # }}} 336 337 338 # {{{ receiver example 339 340 -class GmshMeshReceiverNumPy(GmshMeshReceiverBase):
341 """GmshReceiver that emulates the semantics of 342 :class:`meshpy.triangle.MeshInfo` and :class:`meshpy.tet.MeshInfo` by using 343 similar fields, but instead of loading data into ForeignArrays, load into 344 NumPy arrays. Since this class is not wrapping any libraries in other 345 languages -- the Gmsh data is obtained via parsing text -- use :mod:`numpy` 346 arrays as the base array data structure for convenience. 347 348 .. versionadded:: 2014.1 349 """ 350
351 - def __init__(self):
352 # Use data fields similar to meshpy.triangle.MeshInfo and 353 # meshpy.tet.MeshInfo 354 self.points = None 355 self.elements = None 356 self.element_types = None 357 self.element_markers = None 358 self.tags = None
359 360 # Gmsh has no explicit concept of facets or faces; certain faces are a type 361 # of element. Consequently, there are no face markers, but elements can be 362 # grouped together in physical groups that serve as markers. 363
364 - def set_up_nodes(self, count):
365 # Preallocate array of nodes within list; treat None as sentinel value. 366 # Preallocation not done for performance, but to assign values at indices 367 # in random order. 368 self.points = [None] * count
369
370 - def add_node(self, node_nr, point):
371 self.points[node_nr] = point
372
373 - def finalize_nodes(self):
374 pass
375
376 - def set_up_elements(self, count):
377 # Preallocation of arrays for assignment elements in random order. 378 self.elements = [None] * count 379 self.element_types = [None] * count 380 self.element_markers = [None] * count 381 self.tags = []
382
383 - def add_element(self, element_nr, element_type, vertex_nrs, 384 lexicographic_nodes, tag_numbers):
385 self.elements[element_nr] = vertex_nrs 386 self.element_types[element_nr] = element_type 387 self.element_markers[element_nr] = tag_numbers
388 # TODO: Add lexicographic node information 389
390 - def finalize_elements(self):
391 pass
392
393 - def add_tag(self, name, index, dimension):
394 self.tags.append((name, index, dimension))
395
396 - def finalize_tags(self):
397 pass
398
399 # }}} 400 401 402 # {{{ file reader 403 404 -class GmshFileFormatError(RuntimeError):
405 pass
406
407 408 -def read_gmsh(receiver, filename, force_dimension=None):
409 """Read a gmsh mesh file from *filename* and feed it to *receiver*. 410 411 :param receiver: Implements the :class:`GmshMeshReceiverBase` interface. 412 :param force_dimension: if not None, truncate point coordinates to 413 this many dimensions. 414 """ 415 mesh_file = open(filename, 'rt') 416 try: 417 result = parse_gmsh(receiver, mesh_file, force_dimension=force_dimension) 418 finally: 419 mesh_file.close() 420 421 return result
422
423 424 -def generate_gmsh(receiver, source, dimensions, order=None, other_options=[], 425 extension="geo", gmsh_executable="gmsh", force_dimension=None):
426 """Run gmsh and feed the output to *receiver*. 427 428 :arg source: an instance of :class:`LiteralSource` or :class:`FileSource` 429 :param receiver: Implements the :class:`GmshMeshReceiverBase` interface. 430 """ 431 from meshpy.gmsh import GmshRunner 432 runner = GmshRunner(source, dimensions, order=order, 433 other_options=other_options, extension=extension, 434 gmsh_executable=gmsh_executable) 435 436 runner.__enter__() 437 try: 438 result = parse_gmsh(receiver, runner.output_file, 439 force_dimension=force_dimension) 440 finally: 441 runner.__exit__(None, None, None) 442 443 return result
444
445 446 -def parse_gmsh(receiver, line_iterable, force_dimension=None):
447 """ 448 :arg source: an instance of :class:`LiteralSource` or :class:`FileSource` 449 :arg receiver: This object will be fed the entities encountered in reading the 450 GMSH file. See :class:`GmshMeshReceiverBase` for the interface this 451 object needs to conform to. 452 :param force_dimension: if not None, truncate point coordinates to this many 453 dimensions. 454 """ 455 456 feeder = LineFeeder(line_iterable) 457 458 # collect the mesh information 459 460 class ElementInfo(Record): 461 pass
462 463 while feeder.has_next_line(): 464 next_line = feeder.get_next_line() 465 if not next_line.startswith("$"): 466 raise GmshFileFormatError( 467 "expected start of section, '%s' found instead" % next_line) 468 469 section_name = next_line[1:] 470 471 if section_name == "MeshFormat": 472 line_count = 0 473 while True: 474 next_line = feeder.get_next_line() 475 if next_line == "$End"+section_name: 476 break 477 478 if line_count == 0: 479 version_number, file_type, data_size = next_line.split() 480 481 if line_count > 0: 482 raise GmshFileFormatError( 483 "more than one line found in MeshFormat section") 484 485 if version_number not in ["2.1", "2.2"]: 486 from warnings import warn 487 warn("unexpected mesh version number '%s' found" 488 % version_number) 489 490 if file_type != "0": 491 raise GmshFileFormatError( 492 "only ASCII gmsh file type is supported") 493 494 line_count += 1 495 496 elif section_name == "Nodes": 497 node_count = int(feeder.get_next_line()) 498 receiver.set_up_nodes(node_count) 499 500 node_idx = 1 501 502 while True: 503 next_line = feeder.get_next_line() 504 if next_line == "$End"+section_name: 505 break 506 507 parts = next_line.split() 508 if len(parts) != 4: 509 raise GmshFileFormatError( 510 "expected four-component line in $Nodes section") 511 512 read_node_idx = int(parts[0]) 513 if read_node_idx != node_idx: 514 raise GmshFileFormatError("out-of-order node index found") 515 516 if force_dimension is not None: 517 point = [float(x) for x in parts[1:force_dimension+1]] 518 else: 519 point = [float(x) for x in parts[1:]] 520 521 receiver.add_node( 522 node_idx-1, 523 np.array(point, dtype=np.float64)) 524 525 node_idx += 1 526 527 if node_count+1 != node_idx: 528 raise GmshFileFormatError("unexpected number of nodes found") 529 530 receiver.finalize_nodes() 531 532 elif section_name == "Elements": 533 element_count = int(feeder.get_next_line()) 534 receiver.set_up_elements(element_count) 535 536 element_idx = 1 537 while True: 538 next_line = feeder.get_next_line() 539 if next_line == "$End"+section_name: 540 break 541 542 parts = [int(x) for x in next_line.split()] 543 544 if len(parts) < 4: 545 raise GmshFileFormatError("too few entries in element line") 546 547 read_element_idx = parts[0] 548 if read_element_idx != element_idx: 549 raise GmshFileFormatError("out-of-order node index found") 550 551 el_type_num = parts[1] 552 try: 553 element_type = \ 554 receiver.gmsh_element_type_to_info_map[el_type_num] 555 except KeyError: 556 raise GmshFileFormatError("unexpected element type %d" 557 % el_type_num) 558 559 tag_count = parts[2] 560 tags = parts[3:3+tag_count] 561 562 # convert to zero-based 563 node_indices = np.array( 564 [x-1 for x in parts[3+tag_count:]], dtype=np.intp) 565 566 if element_type.node_count() != len(node_indices): 567 raise GmshFileFormatError( 568 "unexpected number of nodes in element") 569 570 gmsh_vertex_nrs = node_indices[:element_type.vertex_count()] 571 zero_based_idx = element_idx - 1 572 573 tag_numbers = [tag for tag in tags[:1] if tag != 0] 574 575 receiver.add_element(element_nr=zero_based_idx, 576 element_type=element_type, vertex_nrs=gmsh_vertex_nrs, 577 lexicographic_nodes=node_indices[ 578 element_type.get_lexicographic_gmsh_node_indices()], 579 tag_numbers=tag_numbers) 580 581 element_idx += 1 582 583 if element_count+1 != element_idx: 584 raise GmshFileFormatError("unexpected number of elements found") 585 586 receiver.finalize_elements() 587 588 elif section_name == "PhysicalNames": 589 name_count = int(feeder.get_next_line()) 590 name_idx = 1 591 592 while True: 593 next_line = feeder.get_next_line() 594 if next_line == "$End"+section_name: 595 break 596 597 dimension, number, name = next_line.split(" ", 2) 598 dimension = int(dimension) 599 number = int(number) 600 601 if not name[0] == '"' or not name[-1] == '"': 602 raise GmshFileFormatError("expected quotes around physical name") 603 604 receiver.add_tag(name[1:-1], number, dimension) 605 606 name_idx += 1 607 608 if name_count+1 != name_idx: 609 raise GmshFileFormatError( 610 "unexpected number of physical names found") 611 612 receiver.finalize_tags() 613 else: 614 # unrecognized section, skip 615 from warnings import warn 616 warn("unrecognized section '%s' in gmsh file" % section_name) 617 while True: 618 next_line = feeder.get_next_line() 619 if next_line == "$End"+section_name: 620 break 621 622 # }}} 623 624 # vim: fdm=marker 625