Source code for pymor.grids.gmsh

# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2018 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

from collections import defaultdict

import numpy as np
import time

from pymor.core.exceptions import GmshError
from pymor.core.logger import getLogger
from pymor.grids.interfaces import BoundaryInfoInterface
from pymor.grids.unstructured import UnstructuredTriangleGrid


[docs]def load_gmsh(gmsh_file): """Parse a Gmsh file and create a corresponding :class:`GmshGrid` and :class:`GmshBoundaryInfo`. Parameters ---------- gmsh_file File handle of the Gmsh MSH-file. Returns ------- grid The generated :class:`GmshGrid`. boundary_info The generated :class:`GmshBoundaryInfo`. """ logger = getLogger('pymor.grids.gmsh.load_gmsh') logger.info('Parsing gmsh file ...') tic = time.time() sections = _parse_gmsh_file(gmsh_file) toc = time.time() t_parse = toc - tic logger.info('Create GmshGrid ...') tic = time.time() grid = GmshGrid(sections) toc = time.time() t_grid = toc - tic logger.info('Create GmshBoundaryInfo ...') tic = time.time() bi = GmshBoundaryInfo(grid, sections) toc = time.time() t_bi = toc - tic logger.info('Parsing took {} s; Grid creation took {} s; BoundaryInfo creation took {} s' .format(t_parse, t_grid, t_bi)) return grid, bi
[docs]class GmshGrid(UnstructuredTriangleGrid): """An :class:`~pymor.grids.unstructured.UnstructuredTriangleGrid` built from an existing Gmsh MSH-file. Parameters ---------- sections Parsed sections of the MSH-file as returned by :func:`load_gmsh`. """ def __init__(self, sections): self.logger.info('Checking if grid is a 2d triangular grid ...') assert {'Nodes', 'Elements', 'PhysicalNames'} <= set(sections.keys()) assert set(sections['Elements'].keys()) <= {'line', 'triangle'} assert 'triangle' in sections['Elements'] assert all(n[1][2] == 0 for n in sections['Nodes']) node_ids = dict(zip([n[0] for n in sections['Nodes']], np.arange(len(sections['Nodes']), dtype=np.int32))) vertices = np.array([n[1][0:2] for n in sections['Nodes']]) faces = np.array([[node_ids[nodes[0]], node_ids[nodes[1]], node_ids[nodes[2]]] for _, _, nodes in sections['Elements']['triangle']]) super().__init__(vertices, faces)
[docs] def __str__(self): return 'GmshGrid with {} triangles, {} edges, {} vertices'.format(self.size(0), self.size(1), self.size(2))
[docs]class GmshBoundaryInfo(BoundaryInfoInterface): """|BoundaryInfo| for a :class:`GmshGrid`. Parameters ---------- grid The corresponding :class:`GmshGrid`. sections Parsed sections of the MSH-file as returned by :func:`load_gmsh`. """ def __init__(self, grid, sections): assert isinstance(grid, GmshGrid) self.grid = grid # Save boundary types. self.boundary_types = [pn[2] for pn in sections['PhysicalNames'] if pn[1] == 1] # Compute ids, since Gmsh starts numbering with 1 instead of 0. name_ids = dict(zip([pn[0] for pn in sections['PhysicalNames']], np.arange(len(sections['PhysicalNames']), dtype=np.int32))) node_ids = dict(zip([n[0] for n in sections['Nodes']], np.arange(len(sections['Nodes']), dtype=np.int32))) if 'line' in sections['Elements']: superentities = grid.superentities(2, 1) # find the edge for given vertices. def find_edge(vertices): edge_set = set(superentities[vertices[0]]).intersection(superentities[vertices[1]]) - {-1} if len(edge_set) != 1: raise ValueError return next(iter(edge_set)) line_ids = {l[0]: find_edge([node_ids[l[2][0]], node_ids[l[2][1]]]) for l in sections['Elements']['line']} # compute boundary masks for all boundary types. masks = {} for bt in self.boundary_types: masks[bt] = [np.array([False]*grid.size(1)), np.array([False]*grid.size(2))] masks[bt][0][[line_ids[l[0]] for l in sections['Elements']['line']]] = \ [(bt == sections['PhysicalNames'][name_ids[l[1][0]]][2]) for l in sections['Elements']['line']] ind = np.array([node_ids[n] for l in sections['Elements']['line'] for n in l[2]]) val = masks[bt][0][[line_ids[l[0]] for l in sections['Elements']['line'] for n in l[2]]] masks[bt][1][ind[val]] = True self._masks = masks
[docs] def mask(self, boundary_type, codim): assert 1 <= codim <= self.grid.dim assert boundary_type in self.boundary_types return self._masks[boundary_type][codim - 1]
def _parse_gmsh_file(f): allowed_sections = ['Nodes', 'Elements', 'PhysicalNames', 'Periodic', 'NodeData', 'ElementData', 'ElementNodeData'] supported_sections = ['Nodes', 'Elements', 'PhysicalNames'] try: l = next(f).strip() if l != '$MeshFormat': raise GmshError('expected $MeshFormat, got {}'.format(l)) l = next(f).strip() header = l.split(' ') if len(header) != 3: raise GmshError('header {} has {} fields, expected 3'.format(l, len(header))) if header[0] != '2.2': raise GmshError('wrong file format version: got {}, expected 2.2'.format(header[0])) try: file_type = int(header[1]) except ValueError: raise GmshError('malformed header: expected integer, got {}'.format(header[1])) if file_type != 0: raise GmshError('wrong file type: only ASCII gmsh files are supported') try: data_size = int(header[2]) # NOQA except ValueError: raise GmshError('malformed header: expected integer, got {}'.format(header[2])) l = next(f).strip() if l != '$EndMeshFormat': raise GmshError('expected $EndMeshFormat, got {}'.format(l)) except StopIteration: raise GmshError('unexcpected end of file') in_section = False sections = defaultdict(list) for l in f: l = l.strip() if l == '': continue if not in_section: if not l.startswith('$'): raise GmshError('expected section name, got {}'.format(l)) section = l[1:] if section not in allowed_sections: raise GmshError('unknown section type: {}'.format(section)) if section not in supported_sections: raise GmshError('unsupported section type: {}'.format(section)) if section in sections: raise GmshError('only one {} section allowed'.format(section)) in_section = True elif l.startswith('$'): if l != '$End' + section: raise GmshError('expected $End{}, got {}'.format(section, l)) in_section = False else: sections[section].append(l) if in_section: raise GmshError('file ended while in section {}'.format(section)) # now we parse each section ... def parse_nodes(nodes): try: num_nodes = int(nodes[0]) except ValueError: raise GmshError('first line of nodes sections is not a number: {}'.format(nodes[0])) if len(nodes) != num_nodes + 1: raise GmshError('number-of-nodes field does not match number of lines in nodes section') nodes = [n.split(' ') for n in nodes[1:]] if not all(len(n) == 4 for n in nodes): raise GmshError('malformed nodes section') try: nodes = [(int(a), (float(b), float(c), float(d))) for a, b, c, d in nodes] except ValueError: raise GmshError('malformed nodes section') return nodes def parse_elements(elements): try: num_elements = int(elements[0]) except ValueError: raise GmshError('first line of elements sections is not a number: {}'.format(elements[0])) if len(elements) != num_elements + 1: raise GmshError('number-of-elements field does not match number of lines in elements section') elements = [e.split(' ') for e in elements[1:]] try: elements = [tuple(int(f) for f in e) for e in elements] except ValueError: raise GmshError('malformed elements section') element_types = {1: 'line', 2: 'triangle'} element_nodes = {'line': 2, 'triangle': 3} def parse_line(fields): if fields[1] not in element_types: raise GmshError('element type {} not supported'.format(fields[0])) element_type = element_types[fields[1]] num_nodes = element_nodes[element_type] num_tags = fields[2] if len(fields) != num_nodes + num_tags + 3: raise GmshError('malformed elements section') return element_type, (fields[0], tuple(fields[3:3 + num_tags]), fields[3 + num_tags:]) elements_by_type = defaultdict(list) for e in elements: t, l = parse_line(e) elements_by_type[t].append(l) return elements_by_type def parse_names(physical_names): try: num_elements = int(physical_names[0]) except ValueError: raise GmshError('first line of physical names sections is not a number: {}'.format(physical_names[0])) if len(physical_names) != num_elements + 1: raise GmshError('number-of-names field does not match number of lines in physical names section') physical_names = [pn.split(' ') for pn in physical_names[1:]] if not all(len(pn) == 3 for pn in physical_names): raise GmshError('malformed physical names section') try: physical_names = [(int(b), int(a), str(c).replace('"', '')) for a, b, c in physical_names] except ValueError: raise GmshError('malformed physical names section') return physical_names parser_map = {'Nodes': parse_nodes, 'Elements': parse_elements, 'PhysicalNames': parse_names} for k, v in sections.items(): sections[k] = parser_map[k](v) return sections