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
29 from pytools import memoize_method, Record
30 from meshpy.gmsh import LiteralSource, FileSource
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 """
71 yield (0, 0)
72 yield (order, 0)
73 yield (0, order)
74
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
86 for i in range(1, order):
87 for j in range(1, order-i):
88 yield (j, i)
89
93 self.line_iterable = iter(line_iterable)
94 self.next_line = None
95
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
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
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
138
141
142 @memoize_method
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
179 from hedge.polynomial import generic_vandermonde
180
181 return generic_vandermonde(
182 list(self.equidistant_unit_nodes()),
183 list(self.basis_functions()))
184
192
195 dimensions = 1
196
197 @memoize_method
199 return [(0,), (self.order,), ] + [
200 (i,) for i in range(1, self.order)]
201
217
232
235 dimensions = 3
236
237 @memoize_method
239
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
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
312
315
318
321
322 - def add_element(self, element_nr, element_type, vertex_nrs,
323 lexicographic_nodes, tag_numbers):
325
328
329 - def add_tag(self, name, index, dimension):
331
334
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
352
353
354 self.points = None
355 self.elements = None
356 self.element_types = None
357 self.element_markers = None
358 self.tags = None
359
360
361
362
363
365
366
367
368 self.points = [None] * count
369
371 self.points[node_nr] = point
372
375
377
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
389
392
393 - def add_tag(self, name, index, dimension):
394 self.tags.append((name, index, dimension))
395
398
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
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
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
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
625