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
35 yield (0, 0)
36 yield (order, 0)
37 yield (0, order)
38
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
48 for i in range(1, order):
49 for j in range(1, order-i):
50 yield (j, i)
51
54 self.line_iterable = iter(line_iterable)
55 self.next_line = None
56
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
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
88
91
92 @memoize_method
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
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):
147
151 dimensions = 1
152
153 @memoize_method
155 return [(0,), (self.order,),] + [
156 (i,) for i in range(1, self.order)]
157
175
193
198 dimensions = 3
199
200 @memoize_method
202
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
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
267
270
273
276
277 - def add_element(self, element_nr, element_type, vertex_nrs,
278 lexicographic_nodes, tag_numbers):
280
283
284 - def add_tag(self, name, index, dimension):
286
289
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
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
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
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
502