[pyfr] 15/32: Implement CGNS reader.
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Thu Apr 21 08:21:51 UTC 2016
This is an automated email from the git hooks/post-receive script.
ghisvail-guest pushed a commit to branch master
in repository pyfr.
commit 909b3195b28212806797750939e9252ee6ac34f5
Author: aerojspark <jin-seok.park at imperial.ac.uk>
Date: Thu Mar 10 19:05:19 2016 +0000
Implement CGNS reader.
---
pyfr/ctypesutil.py | 2 +-
pyfr/readers/__init__.py | 3 +-
pyfr/readers/base.py | 198 +++++++++++++++++++++
pyfr/readers/cgns.py | 435 +++++++++++++++++++++++++++++++++++++++++++++++
pyfr/readers/gmsh.py | 204 ++--------------------
pyfr/readers/nodemaps.py | 69 ++++++++
6 files changed, 716 insertions(+), 195 deletions(-)
diff --git a/pyfr/ctypesutil.py b/pyfr/ctypesutil.py
index fc38335..6f34a46 100644
--- a/pyfr/ctypesutil.py
+++ b/pyfr/ctypesutil.py
@@ -1,4 +1,4 @@
- # -*- coding: utf-8 -*-
+# -*- coding: utf-8 -*-
import ctypes
import ctypes.util
diff --git a/pyfr/readers/__init__.py b/pyfr/readers/__init__.py
index 7db8825..ac8820a 100644
--- a/pyfr/readers/__init__.py
+++ b/pyfr/readers/__init__.py
@@ -1,6 +1,7 @@
# -*- coding: utf-8 -*-
-from pyfr.readers.base import BaseReader
+from pyfr.readers.base import BaseReader, NodalMeshAssembler
+from pyfr.readers.cgns import CGNSReader
from pyfr.readers.gmsh import GmshReader
from pyfr.util import subclasses, subclass_where
diff --git a/pyfr/readers/base.py b/pyfr/readers/base.py
index e5ebc81..c18ecdc 100644
--- a/pyfr/readers/base.py
+++ b/pyfr/readers/base.py
@@ -1,10 +1,14 @@
# -*- coding: utf-8 -*-
from abc import ABCMeta, abstractmethod
+from collections import defaultdict
+from itertools import chain
import uuid
import numpy as np
+from pyfr.nputil import fuzzysort
+
class BaseReader(object, metaclass=ABCMeta):
@abstractmethod
@@ -22,3 +26,197 @@ class BaseReader(object, metaclass=ABCMeta):
mesh['mesh_uuid'] = np.array(str(uuid.uuid4()), dtype='S')
return mesh
+
+
+class NodalMeshAssembler(object):
+ # Dimensionality of each element type
+ _petype_ndim = {'tri': 2, 'quad': 2,
+ 'tet': 3, 'hex': 3, 'pri': 3, 'pyr': 3}
+
+ # Face numberings for each element type
+ _petype_fnums = {
+ 'tri': {'line': [0, 1, 2]},
+ 'quad': {'line': [0, 1, 2, 3]},
+ 'tet': {'tri': [0, 1, 2, 3]},
+ 'hex': {'quad': [0, 1, 2, 3, 4, 5]},
+ 'pri': {'quad': [2, 3, 4], 'tri': [0, 1]},
+ 'pyr': {'quad': [0], 'tri': [1, 2, 3, 4]}
+ }
+
+ # Number of nodes in the first-order representation an element
+ _petype_focount = {'line': 2, 'tri': 3, 'quad': 4,
+ 'tet': 4, 'pyr': 5, 'pri': 6, 'hex': 8}
+
+ def __init__(self, nodepts, elenodes, pents, maps):
+ self._nodepts = nodepts
+ self._elenodes = elenodes
+ self._felespent, self._bfacespents, self._pfacespents = pents
+ self._etype_map, self._petype_fnmap, self._nodemaps = maps
+
+ def _to_first_order(self, elemap):
+ foelemap = {}
+ for (etype, epent), eles in elemap.items():
+ # PyFR element type ('hex', 'tri', &c)
+ petype = self._etype_map[etype][0]
+
+ # Number of nodes in the first-order representation
+ focount = self._petype_focount[petype]
+
+ foelemap[petype, epent] = eles[:,:focount]
+
+ return foelemap
+
+ def _split_fluid(self, elemap):
+ selemap = defaultdict(dict)
+
+ for (petype, epent), eles in elemap.items():
+ selemap[epent][petype] = eles
+
+ return selemap.pop(self._felespent), selemap
+
+ def _foface_info(self, petype, pftype, foeles):
+ # Face numbers of faces of this type on this element
+ fnums = self._petype_fnums[petype][pftype]
+
+ # First-order nodes associated with this face
+ fnmap = self._petype_fnmap[petype][pftype]
+
+ # Connectivity; (petype, eidx, fidx, flags)
+ con = [(petype, i, j, 0) for i in range(len(foeles)) for j in fnums]
+
+ # Nodes
+ nodes = np.sort(foeles[:, fnmap]).reshape(len(con), -1)
+
+ return con, nodes
+
+ def _extract_faces(self, foeles):
+ fofaces = defaultdict(list)
+ for petype, eles in foeles.items():
+ for pftype in self._petype_fnums[petype]:
+ fofinf = self._foface_info(petype, pftype, eles)
+ fofaces[pftype].append(fofinf)
+
+ return fofaces
+
+ def _pair_fluid_faces(self, ffofaces):
+ pairs = defaultdict(list)
+ resid = {}
+
+ for pftype, faces in ffofaces.items():
+ for f, n in chain.from_iterable(zip(f, n) for f, n in faces):
+ sn = tuple(n)
+
+ # See if the nodes are in resid
+ if sn in resid:
+ pairs[pftype].append([resid.pop(sn), f])
+ # Otherwise add them to the unpaired dict
+ else:
+ resid[sn] = f
+
+ return pairs, resid
+
+ def _pair_periodic_fluid_faces(self, bpart, resid):
+ pfaces = defaultdict(list)
+ nodepts = self._nodepts
+
+ for lpent, rpent in self._pfacespents.values():
+ for pftype in bpart[lpent]:
+ lfnodes = bpart[lpent][pftype]
+ rfnodes = bpart[rpent][pftype]
+
+ lfpts = np.array([[nodepts[n] for n in fn] for fn in lfnodes])
+ rfpts = np.array([[nodepts[n] for n in fn] for fn in rfnodes])
+
+ lfidx = fuzzysort(lfpts.mean(axis=1).T, range(len(lfnodes)))
+ rfidx = fuzzysort(rfpts.mean(axis=1).T, range(len(rfnodes)))
+
+ for lfn, rfn in zip(lfnodes[lfidx], rfnodes[rfidx]):
+ lf = resid.pop(tuple(sorted(lfn)))
+ rf = resid.pop(tuple(sorted(rfn)))
+
+ pfaces[pftype].append([lf, rf])
+
+ return pfaces
+
+ def _ident_boundary_faces(self, bpart, resid):
+ bfaces = defaultdict(list)
+
+ bpents = set(self._bfacespents.values())
+
+ for epent, fnodes in bpart.items():
+ if epent in bpents:
+ for fn in chain.from_iterable(fnodes.values()):
+ bfaces[epent].append(resid.pop(tuple(sorted(fn))))
+
+ return bfaces
+
+ def get_connectivity(self):
+ # For connectivity a first-order representation is sufficient
+ eles = self._to_first_order(self._elenodes)
+
+ # Split into fluid and boundary parts
+ fpart, bpart = self._split_fluid(eles)
+
+ # Extract the faces of the fluid elements
+ ffaces = self._extract_faces(fpart)
+
+ # Pair the fluid-fluid faces
+ fpairs, resid = self._pair_fluid_faces(ffaces)
+
+ # Identify periodic boundary face pairs
+ pfpairs = self._pair_periodic_fluid_faces(bpart, resid)
+
+ # Identify the fixed boundary faces
+ bf = self._ident_boundary_faces(bpart, resid)
+
+ if any(resid.values()):
+ raise ValueError('Unpaired faces in mesh')
+
+ # Flattern the face-pair dicts
+ pairs = chain(chain.from_iterable(fpairs.values()),
+ chain.from_iterable(pfpairs.values()))
+
+ # Generate the internal connectivity array
+ con = list(pairs)
+
+ # Generate boundary condition connectivity arrays
+ bcon = {}
+ for pbcrgn, pent in self._bfacespents.items():
+ bcon[pbcrgn] = bf[pent]
+
+ # Output
+ ret = {'con_p0': np.array(con, dtype='S4,i4,i1,i1').T}
+
+ for k, v in bcon.items():
+ ret['bcon_{0}_p0'.format(k)] = np.array(v, dtype='S4,i4,i1,i1')
+
+ return ret
+
+ def get_shape_points(self):
+ spts = {}
+
+ # Global node map (node index to coords)
+ nodepts = self._nodepts
+
+ for etype, pent in self._elenodes:
+ if pent != self._felespent:
+ continue
+
+ # Elements and type information
+ eles = self._elenodes[etype, pent]
+ petype, nnodes = self._etype_map[etype]
+
+ # Go from Gmsh to PyFR node ordering
+ peles = eles[:, self._nodemaps.from_pyfr[petype, nnodes]]
+
+ # Obtain the dimensionality of the element type
+ ndim = self._petype_ndim[petype]
+
+ # Build the array
+ arr = np.array([[nodepts[i] for i in nn] for nn in peles])
+ arr = arr.swapaxes(0, 1)
+ arr = arr[..., :ndim]
+
+ spts['spt_{0}_p0'.format(petype)] = arr
+
+ return spts
diff --git a/pyfr/readers/cgns.py b/pyfr/readers/cgns.py
new file mode 100644
index 0000000..92dacbd
--- /dev/null
+++ b/pyfr/readers/cgns.py
@@ -0,0 +1,435 @@
+# -*- coding: utf-8 -*-
+
+from collections import defaultdict
+from ctypes import POINTER, create_string_buffer, c_char_p, c_int, c_void_p
+import re
+
+import numpy as np
+
+from pyfr.ctypesutil import load_library
+from pyfr.readers import BaseReader, NodalMeshAssembler
+from pyfr.readers.nodemaps import CGNSNodeMaps
+
+
+# Possible CGNS exception types
+CGNSError = type('CGNSError', (Exception,), {})
+CGNSNodeNotFound = type('CGNSNodeNotFound', (CGNSError,), {})
+CGNSIncorrectPath = type('CGNSIncorrectPath', (CGNSError,), {})
+CGNSNoIndexDim = type('CGNSNoIndexDim', (CGNSError,), {})
+
+
+class CGNSWrappers(object):
+ # Possible return codes
+ _statuses = {
+ -1: CGNSError,
+ -2: CGNSNodeNotFound,
+ -3: CGNSIncorrectPath,
+ -4: CGNSNoIndexDim
+ }
+
+ def __init__(self):
+ # Load CGNS 3.3+
+ # Previous version may yield undefined symbols from HDF5.
+ # Develop branch after the commit (e0faea6) fixes this issue.
+ self.lib = lib = load_library('cgns')
+
+ # Constants (from cgnslib.h)
+ self.CG_MODE_READ = 0
+ self.RealDouble = 4
+ self.Unstructured = 3
+ self.PointRange, self.ElementRange = 4, 6
+
+ # cg_open
+ lib.cg_open.argtypes = [c_char_p, c_int, POINTER(c_int)]
+ lib.cg_open.errcheck = self._errcheck
+
+ # cg_close
+ lib.cg_close.argtypes = [c_int]
+ lib.cg_close.errcheck = self._errcheck
+
+ # cg_nbases
+ lib.cg_nbases.argtypes = [c_int, POINTER(c_int)]
+ lib.cg_nbases.errcheck = self._errcheck
+
+ # cg_base_read
+ lib.cg_base_read.argtypes = [c_int, c_int, c_char_p, POINTER(c_int),
+ POINTER(c_int)]
+ lib.cg_base_read.errcheck = self._errcheck
+
+ # cg_nzones
+ lib.cg_nzones.argtypes = [c_int, c_int, POINTER(c_int)]
+ lib.cg_nzones.errcheck = self._errcheck
+
+ # cg_zone_read
+ lib.cg_zone_read.argtypes = [c_int, c_int, c_int, c_char_p,
+ POINTER(c_int)]
+ lib.cg_zone_read.errcheck = self._errcheck
+
+ # cg_zone_type
+ lib.cg_zone_type.argtypes = [c_int, c_int, c_int, POINTER(c_int)]
+ lib.cg_zone_type.errcheck = self._errcheck
+
+ # cg_coord_read
+ lib.cg_coord_read.argtypes = [
+ c_int, c_int, c_int, c_char_p, c_int, POINTER(c_int),
+ POINTER(c_int), c_void_p
+ ]
+ lib.cg_coord_read.errcheck = self._errcheck
+
+ # cg_nbocos
+ lib.cg_nbocos.argtypes = [c_int, c_int, c_int, POINTER(c_int)]
+ lib.cg_nbocos.errcheck = self._errcheck
+
+ # cg_boco_info
+ lib.cg_boco_info.argtypes = [
+ c_int, c_int, c_int, c_int, c_char_p, POINTER(c_int),
+ POINTER(c_int), POINTER(c_int), POINTER(c_int),
+ POINTER(c_int), POINTER(c_int), POINTER(c_int)
+ ]
+ lib.cg_boco_info.errcheck = self._errcheck
+
+ # cg_boco_read
+ lib.cg_boco_read.argtypes = [c_int, c_int, c_int, c_int,
+ POINTER(c_int), c_void_p]
+ lib.cg_boco_read.errcheck = self._errcheck
+
+ # cg_nsections
+ lib.cg_nsections.argtypes = [c_int, c_int, c_int, POINTER(c_int)]
+ lib.cg_nsections.errcheck = self._errcheck
+
+ # cg_section_read
+ lib.cg_section_read.argtypes = [
+ c_int, c_int, c_int, c_int, c_char_p, POINTER(c_int),
+ POINTER(c_int), POINTER(c_int), POINTER(c_int), POINTER(c_int)
+ ]
+ lib.cg_section_read.errcheck = self._errcheck
+
+ # cg_ElementDataSize
+ lib.cg_ElementDataSize.argtypes = [c_int, c_int, c_int, c_int,
+ POINTER(c_int)]
+ lib.cg_ElementDataSize.errcheck = self._errcheck
+
+ # cg_elements_read
+ lib.cg_elements_read.argtypes = [c_int, c_int, c_int, c_int,
+ c_void_p, c_void_p]
+ lib.cg_elements_read.errcheck = self._errcheck
+
+ def _errcheck(self, status, fn, arg):
+ if status != 0:
+ try:
+ raise self._statuses[status]
+ except KeyError:
+ raise CGNSError
+
+ def open(self, name):
+ file = c_int()
+ self.lib.cg_open(bytes(name, 'utf-8'), self.CG_MODE_READ, file)
+ return file
+
+ def close(self, file):
+ self.lib.cg_close(file)
+
+ def nbases(self, file):
+ nb = c_int()
+ self.lib.cg_nbases(file, nb)
+ return nb.value
+
+ def base_read(self, file, idx):
+ celldim, physdim = c_int(), c_int()
+ name = create_string_buffer(32)
+
+ self.lib.cg_base_read(file, idx + 1, name, celldim, physdim)
+
+ return {'file': file, 'idx': idx + 1,
+ 'name': name.value.decode('utf-8'),
+ 'CellDim': celldim.value, 'PhysDim': physdim.value}
+
+ def nzones(self, base):
+ n = c_int()
+ self.lib.cg_nzones(base['file'], base['idx'], n)
+ return n.value
+
+ def zone_read(self, base, idx):
+ zonetype = c_int()
+ name = create_string_buffer(32)
+ size = (c_int * 3)()
+
+ self.lib.cg_zone_read(base['file'], base['idx'], idx + 1, name, size)
+
+ # Check zone type
+ self.lib.cg_zone_type(base['file'], base['idx'], idx + 1, zonetype)
+ if zonetype.value != self.Unstructured:
+ raise RuntimeError('ReadCGNS_read: Incorrect zone type for file')
+
+ return {'base': base, 'idx': idx + 1,
+ 'name': name.value.decode('utf-8'),
+ 'size': list(size)}
+
+ def coord_read(self, zone, name, x):
+ i = c_int(1)
+ j = c_int(zone['size'][0])
+
+ file = zone['base']['file']
+ base = zone['base']['idx']
+ zone = zone['idx']
+
+ # The data type does not need to be the same as the one in which the
+ # coordinates are stored in the file
+ # http://cgns.github.io/CGNS_docs_current/midlevel/grid.html
+ datatype = self.RealDouble
+
+ self.lib.cg_coord_read(file, base, zone, bytes(name, 'utf-8'), datatype,
+ i, j, x.ctypes.data)
+
+ def nbocos(self, zone):
+ file = zone['base']['file']
+ base = zone['base']['idx']
+ zone = zone['idx']
+ n = c_int()
+
+ self.lib.cg_goto(file, base, b'Zone_t', 1, b'ZoneBC_t', 1, b'end')
+ self.lib.cg_nbocos(file, base, zone, n)
+ return n.value
+
+ def boco_read(self, zone, idx):
+ file = zone['base']['file']
+ base = zone['base']['idx']
+ zone = zone['idx']
+
+ name = create_string_buffer(32)
+ bocotype = c_int()
+ ptset_type = c_int()
+ npnts = c_int()
+ normalindex = (c_int * 3)()
+ normallistsize = c_int()
+ normaldatatype = c_int()
+ ndataset = c_int()
+ bcrange = (c_int * 2)()
+
+ self.lib.cg_boco_info(
+ file, base, zone, idx + 1, name, bocotype, ptset_type, npnts,
+ normalindex, normallistsize, normaldatatype, ndataset
+ )
+
+ if ptset_type.value not in [self.PointRange, self.ElementRange]:
+ raise RuntimeError('Only element range BC is supported')
+
+ self.lib.cg_boco_read(file, base, zone, idx + 1, bcrange, None)
+
+ return {'name': name.value.decode('utf-8'),
+ 'range': tuple(bcrange)}
+
+ def nsections(self, zone):
+ file = zone['base']['file']
+ base = zone['base']['idx']
+ zone = zone['idx']
+
+ n = c_int()
+ self.lib.cg_nsections(file, base, zone, n)
+
+ return n.value
+
+ def section_read(self, zone, idx):
+ file = zone['base']['file']
+ base = zone['base']['idx']
+ zidx = zone['idx']
+
+ name = create_string_buffer(32)
+ etype, start, end, nbdry = c_int(), c_int(), c_int(), c_int()
+ pflag, cdim = c_int(), c_int()
+
+ self.lib.cg_section_read(
+ file, base, zidx, idx + 1, name, etype, start, end, nbdry, pflag
+ )
+
+ self.lib.cg_ElementDataSize(file, base, zidx, idx + 1, cdim)
+
+ return {'zone': zone, 'idx': idx + 1, 'dim': cdim.value,
+ 'etype': etype.value, 'range': (start.value, end.value)}
+
+ def elements_read(self, sect, conn):
+ file = sect['zone']['base']['file']
+ base = sect['zone']['base']['idx']
+ zone = sect['zone']['idx']
+ idx = sect['idx']
+
+ self.lib.cg_elements_read(file, base, zone, idx,
+ conn.ctypes.data, None)
+
+
+class CGNSZoneReader(object):
+ # CGNS element types to PyFR type (petype) and node counts
+ cgns_map = {
+ 3: ('line', 2), 4: ('line', 3), 24: ('line', 4), 40: ('line', 5),
+ 5: ('tri', 3), 6: ('tri', 6), 26: ('tri', 10), 42: ('tri', 15),
+ 7: ('quad', 4), 9: ('quad', 9), 28: ('quad', 16), 44: ('quad', 25),
+ 10: ('tet', 4), 11: ('tet', 10), 30: ('tet', 20), 47: ('tet', 35),
+ 12: ('pyr', 5), 13: ('pyr', 14), 33: ('pyr', 30), 50: ('pyr', 55),
+ 14: ('pri', 6), 16: ('pri', 18), 36: ('pri', 40), 53: ('pri', 75),
+ 17: ('hex', 8), 19: ('hex', 27), 39: ('hex', 64), 56: ('hex', 125),
+ 20: ('mixed',),
+ }
+
+ def __init__(self, cgns, base, idx):
+ self._cgns = cgns
+ zone = cgns.zone_read(base, idx)
+
+ # Read nodes
+ self.nodepts = self._read_nodepts(zone)
+
+ # Read bc
+ bc = self._read_bc(zone)
+
+ # Read elements
+ self.elenodes = elenodes = {}
+ self.pents = pents = {}
+
+ # Construct elenodes and physical entity
+ for idx in range(cgns.nsections(zone)):
+ elerng, elenode = self._read_element(zone, idx)
+
+ for bcname, bcrng in bc.items():
+ if elerng[0] >= bcrng[0] and elerng[1] <= bcrng[1]:
+ name = bcname
+ break
+ else:
+ name = 'fluid'
+
+ pent = pents.setdefault(name, idx)
+
+ elenodes.update({(k, pent): v for k, v in elenode.items()})
+
+ def _read_nodepts(self, zone):
+ nnode = zone['size'][0]
+ nodepts = np.zeros((3, nnode))
+
+ self._cgns.coord_read(zone, 'CoordinateX', nodepts[0])
+ self._cgns.coord_read(zone, 'CoordinateY', nodepts[1])
+ self._cgns.coord_read(zone, 'CoordinateZ', nodepts[2])
+
+ return nodepts
+
+ def _read_bc(self, zone):
+ nbc = self._cgns.nbocos(zone)
+ bc = {}
+
+ for idx_bc in range(nbc):
+ boco = self._cgns.boco_read(zone, idx_bc)
+ name = boco['name'].lower()
+ bc[name] = boco['range']
+
+ return bc
+
+ def _read_element(self, zone, idx):
+ s = self._cgns.section_read(zone, idx)
+
+ elerng = s['range']
+ conn = np.zeros(s['dim'], dtype=np.int32)
+ self._cgns.elements_read(s, conn)
+
+ cgns_type = s['etype']
+ petype = self.cgns_map[cgns_type][0]
+
+ elenode = {}
+
+ if petype == 'mixed':
+ i = 0
+ mele = defaultdict(list)
+
+ while i < s['dim']:
+ try:
+ cgns_type = conn[i]
+ petype, spts = self.cgns_map[cgns_type]
+ except KeyError:
+ raise
+
+ mele[cgns_type].append(conn[i + 1: i + 1 + spts])
+ i += 1 + spts
+
+ elenode = {k: np.array(v) for k, v in mele.items()}
+ else:
+ spts = self.cgns_map[cgns_type][1]
+ elenode[cgns_type] = conn.reshape(-1, spts)
+
+ return elerng, elenode
+
+
+class CGNSReader(BaseReader):
+ # Supported file types and extensions
+ name = 'cgns'
+ extn = ['.cgns']
+
+ # CGNS element types to PyFR type (petype) and node counts
+ _etype_map = CGNSZoneReader.cgns_map
+
+ # First-order node numbers associated with each element face
+ _petype_fnmap = {
+ 'tri': {'line': [[0, 1], [1, 2], [2, 0]]},
+ 'quad': {'line': [[0, 1], [1, 2], [2, 3], [3, 0]]},
+ 'tet': {'tri': [[0, 1, 2], [0, 1, 3], [0, 2, 3], [1, 2, 3]]},
+ 'hex': {'quad': [[0, 1, 2, 3], [0, 1, 4, 5], [1, 2, 5, 6],
+ [2, 3, 6, 7], [0, 3, 4, 7], [4, 5, 6, 7]]},
+ 'pri': {'quad': [[0, 1, 3, 4], [1, 2, 4, 5], [0, 2, 3, 5]],
+ 'tri': [[0, 1, 2], [3, 4, 5]]},
+ 'pyr': {'quad': [[0, 1, 2, 3]],
+ 'tri': [[0, 1, 4], [1, 2, 4], [2, 3, 4], [0, 3, 4]]}
+ }
+
+ # Mappings between the node ordering of PyFR and that of CGNS
+ _nodemaps = CGNSNodeMaps
+
+ def __init__(self, msh):
+ # Load and wrap CGNS
+ self._cgns = cgns = CGNSWrappers()
+
+ # Read CGNS mesh file
+ self._file = file = cgns.open(msh.name)
+ base = cgns.base_read(file, 0)
+
+ if cgns.nzones(base) > 1:
+ raise ValueError('Only single zone is supported')
+
+ # Read the single CGNS Zone
+ zone = CGNSZoneReader(cgns, base, 0)
+
+ self._nodepts = {i + 1: v
+ for i, v in enumerate(zone.nodepts.swapaxes(0, 1))}
+ self._elenodes = zone.elenodes
+ pents = zone.pents
+
+ # Physical entities can be divided up into:
+ # - fluid elements ('the mesh')
+ # - boundary faces
+ # - periodic faces
+ self._felespent = pents.pop('fluid')
+ self._bfacespents = {}
+ self._pfacespents = defaultdict(list)
+
+ for name, pent in pents.items():
+ if name.startswith('periodic'):
+ p = re.match(r'periodic[ _-]([a-z0-9]+)[ _-](l|r)$', name)
+ if not p:
+ raise ValueError('Invalid periodic boundary condition')
+
+ self._pfacespents[p.group(1)].append(name)
+ # Other boundary faces
+ else:
+ self._bfacespents[name] = pent
+
+ if any(len(pf) != 2 for pf in self._pfacespents.values()):
+ raise ValueError('Unpaired periodic boundary in mesh')
+
+ def __del__(self):
+ if hasattr(self, '_file'):
+ self._cgns.close(self._file)
+
+ def _to_raw_pyfrm(self):
+ # Assemble a nodal mesh
+ maps = self._etype_map, self._petype_fnmap, self._nodemaps
+ pents = self._felespent, self._bfacespents, self._pfacespents
+ mesh = NodalMeshAssembler(self._nodepts, self._elenodes, pents, maps)
+
+ rawm = {}
+ rawm.update(mesh.get_connectivity())
+ rawm.update(mesh.get_shape_points())
+ return rawm
diff --git a/pyfr/readers/gmsh.py b/pyfr/readers/gmsh.py
index d29e97a..80a0f10 100644
--- a/pyfr/readers/gmsh.py
+++ b/pyfr/readers/gmsh.py
@@ -1,14 +1,12 @@
# -*- coding: utf-8 -*-
from collections import defaultdict
-from itertools import chain
import re
import numpy as np
-from pyfr.readers import BaseReader
+from pyfr.readers import BaseReader, NodalMeshAssembler
from pyfr.readers.nodemaps import GmshNodeMaps
-from pyfr.nputil import fuzzysort
def msh_section(mshit, section):
@@ -46,23 +44,6 @@ class GmshReader(BaseReader):
7: ('pyr', 5), 14: ('pyr', 14), 118: ('pyr', 30), 119: ('pyr', 55)
}
- # Number of nodes in the first-order representation an element
- _petype_focount = {v[0]: v[1] for k, v in _etype_map.items() if k < 8}
-
- # Dimensionality of each element type
- _petype_ndim = {'tri': 2, 'quad': 2,
- 'tet': 3, 'hex': 3, 'pri': 3, 'pyr': 3}
-
- # Face numberings for each element type
- _petype_fnums = {
- 'tri': {'line': [0, 1, 2]},
- 'quad': {'line': [0, 1, 2, 3]},
- 'tet': {'tri': [0, 1, 2, 3]},
- 'hex': {'quad': [0, 1, 2, 3, 4, 5]},
- 'pri': {'quad': [2, 3, 4], 'tri': [0, 1]},
- 'pyr': {'quad': [0], 'tri': [1, 2, 3, 4]}
- }
-
# First-order node numbers associated with each element face
_petype_fnmap = {
'tri': {'line': [[0, 1], [1, 2], [2, 0]]},
@@ -76,6 +57,9 @@ class GmshReader(BaseReader):
'tri': [[0, 1, 4], [1, 2, 4], [2, 3, 4], [0, 3, 4]]}
}
+ # Mappings between the node ordering of PyFR and that of Gmsh
+ _nodemaps = GmshNodeMaps
+
def __init__(self, msh):
if isinstance(msh, str):
msh = open(msh)
@@ -206,179 +190,13 @@ class GmshReader(BaseReader):
self._elenodes = {k: np.array(v) for k, v in elenodes.items()}
- def _to_first_order(self, elemap):
- foelemap = {}
- for (etype, epent), eles in elemap.items():
- # PyFR element type ('hex', 'tri', &c)
- petype = self._etype_map[etype][0]
-
- # Number of nodes in the first-order representation
- focount = self._petype_focount[petype]
-
- foelemap[petype, epent] = eles[:,:focount]
-
- return foelemap
-
- def _split_fluid(self, elemap):
- selemap = defaultdict(dict)
-
- for (petype, epent), eles in elemap.items():
- selemap[epent][petype] = eles
-
- return selemap.pop(self._felespent), selemap
-
- def _foface_info(self, petype, pftype, foeles):
- # Face numbers of faces of this type on this element
- fnums = self._petype_fnums[petype][pftype]
-
- # First-order nodes associated with this face
- fnmap = self._petype_fnmap[petype][pftype]
-
- # Number of first order nodes needed to define a face of this type
- nfnodes = self._petype_focount[pftype]
-
- # Connectivity; (petype, eidx, fidx, flags)
- con = [(petype, i, j, 0) for i in range(len(foeles)) for j in fnums]
-
- # Nodes
- nodes = np.sort(foeles[:, fnmap]).reshape(len(con), -1)
-
- return con, nodes
-
- def _extract_faces(self, foeles):
- fofaces = defaultdict(list)
- for petype, eles in foeles.items():
- for pftype in self._petype_fnums[petype]:
- fofinf = self._foface_info(petype, pftype, eles)
- fofaces[pftype].append(fofinf)
-
- return fofaces
-
- def _pair_fluid_faces(self, ffofaces):
- pairs = defaultdict(list)
- resid = {}
-
- for pftype, faces in ffofaces.items():
- for f, n in chain.from_iterable(zip(f, n) for f, n in faces):
- sn = tuple(n)
-
- # See if the nodes are in resid
- if sn in resid:
- pairs[pftype].append([resid.pop(sn), f])
- # Otherwise add them to the unpaired dict
- else:
- resid[sn] = f
-
- return pairs, resid
-
- def _pair_periodic_fluid_faces(self, bpart, resid):
- pfaces = defaultdict(list)
- nodepts = self._nodepts
-
- for lpent, rpent in self._pfacespents.values():
- for pftype in bpart[lpent]:
- lfnodes = bpart[lpent][pftype]
- rfnodes = bpart[rpent][pftype]
-
- lfpts = np.array([[nodepts[n] for n in fn] for fn in lfnodes])
- rfpts = np.array([[nodepts[n] for n in fn] for fn in rfnodes])
-
- lfidx = fuzzysort(lfpts.mean(axis=1).T, range(len(lfnodes)))
- rfidx = fuzzysort(rfpts.mean(axis=1).T, range(len(rfnodes)))
-
- for lfn, rfn in zip(lfnodes[lfidx], rfnodes[rfidx]):
- lf = resid.pop(tuple(sorted(lfn)))
- rf = resid.pop(tuple(sorted(rfn)))
-
- pfaces[pftype].append([lf, rf])
-
- return pfaces
-
- def _ident_boundary_faces(self, bpart, resid):
- bfaces = defaultdict(list)
-
- bpents = set(self._bfacespents.values())
-
- for epent, fnodes in bpart.items():
- if epent in bpents:
- for fn in chain.from_iterable(fnodes.values()):
- bfaces[epent].append(resid.pop(tuple(sorted(fn))))
-
- return bfaces
-
- def _get_connectivity(self):
- # For connectivity a first-order representation is sufficient
- eles = self._to_first_order(self._elenodes)
-
- # Split into fluid and boundary parts
- fpart, bpart = self._split_fluid(eles)
-
- # Extract the faces of the fluid elements
- ffaces = self._extract_faces(fpart)
-
- # Pair the fluid-fluid faces
- fpairs, resid = self._pair_fluid_faces(ffaces)
-
- # Identify periodic boundary face pairs
- pfpairs = self._pair_periodic_fluid_faces(bpart, resid)
-
- # Identify the fixed boundary faces
- bf = self._ident_boundary_faces(bpart, resid)
-
- if any(resid.values()):
- raise ValueError('Unpaired faces in mesh')
-
- # Flattern the face-pair dicts
- pairs = chain(chain.from_iterable(fpairs.values()),
- chain.from_iterable(pfpairs.values()))
-
- # Generate the internal connectivity array
- con = list(pairs)
-
- # Generate boundary condition connectivity arrays
- bcon = {}
- for pbcrgn, pent in self._bfacespents.items():
- bcon[pbcrgn] = bf[pent]
-
- # Output
- ret = {'con_p0': np.array(con, dtype='S4,i4,i1,i1').T}
-
- for k, v in bcon.items():
- ret['bcon_{0}_p0'.format(k)] = np.array(v, dtype='S4,i4,i1,i1')
-
- return ret
-
- def _get_shape_points(self):
- spts = {}
-
- # Global node map (node index to coords)
- nodepts = self._nodepts
-
- for etype, pent in self._elenodes:
- if pent != self._felespent:
- continue
-
- # Elements and type information
- eles = self._elenodes[etype, pent]
- petype, nnodes = self._etype_map[etype]
-
- # Go from Gmsh to PyFR node ordering
- peles = eles[:,GmshNodeMaps.from_pyfr[petype, nnodes]]
-
- # Obtain the dimensionality of the element type
- ndim = self._petype_ndim[petype]
-
- # Build the array
- arr = np.array([[nodepts[i] for i in nn] for nn in peles])
- arr = arr.swapaxes(0, 1)
- arr = arr[..., :ndim]
-
- spts['spt_{0}_p0'.format(petype)] = arr
-
- return spts
-
def _to_raw_pyfrm(self):
+ # Assemble a nodal mesh
+ maps = self._etype_map, self._petype_fnmap, self._nodemaps
+ pents = self._felespent, self._bfacespents, self._pfacespents
+ mesh = NodalMeshAssembler(self._nodepts, self._elenodes, pents, maps)
+
rawm = {}
- rawm.update(self._get_connectivity())
- rawm.update(self._get_shape_points())
+ rawm.update(mesh.get_connectivity())
+ rawm.update(mesh.get_shape_points())
return rawm
diff --git a/pyfr/readers/nodemaps.py b/pyfr/readers/nodemaps.py
index df957b8..07ac26a 100644
--- a/pyfr/readers/nodemaps.py
+++ b/pyfr/readers/nodemaps.py
@@ -187,3 +187,72 @@ class GmshNodeMaps(object):
}
from_pyfr = {k: np.argsort(v) for k, v in to_pyfr.items()}
+
+
+class CGNSNodeMaps(object):
+ to_pyfr = {
+ ('tet', 4): np.array([0, 1, 2, 3]),
+ ('tet', 10): np.array([0, 2, 5, 9, 1, 4, 3, 6, 7, 8]),
+ ('tet', 20): np.array([0, 3, 9, 19, 1, 2, 6, 8, 7, 4, 10, 16, 12, 17,
+ 15, 18, 5, 11, 14, 13]),
+ ('tet', 35): np.array([0, 4, 14, 34, 1, 2, 3, 8, 11, 13, 12, 9, 5, 15,
+ 25, 31, 18, 27, 32, 24, 30, 33, 6, 7, 10, 16, 17,
+ 26, 21, 23, 29, 22, 19, 28, 20]),
+ ('pri', 6): np.array([0, 1, 2, 3, 4, 5]),
+ ('pri', 18): np.array([0, 2, 5, 12, 14, 17, 1, 4, 3, 6, 8, 11, 13, 16,
+ 15, 7, 10, 9]),
+ ('pri', 40): np.array([0, 3, 9, 30, 33, 39, 1, 2, 6, 8, 7, 4, 10, 20,
+ 13, 23, 19, 29, 31, 32, 36, 38, 37, 34, 5, 11,
+ 12, 22, 21, 16, 18, 28, 26, 17, 14, 24, 27, 35,
+ 15, 25]),
+ ('pri', 75): np.array([0, 4, 14, 60, 64, 74, 1, 2, 3, 8, 11, 13, 12, 9,
+ 5, 15, 30, 45, 19, 34, 49, 29, 44, 59, 61, 62,
+ 63, 68, 71, 73, 72, 69, 65, 6, 7, 10, 16, 17, 18,
+ 33, 48, 47, 46, 31, 32, 23, 26, 28, 43, 58, 56,
+ 53, 38, 41, 27, 24, 20, 35, 50, 54, 57, 42, 39,
+ 66, 67, 70, 21, 22, 25, 36, 37, 40, 51, 52, 55]),
+ ('pyr', 5): np.array([0, 1, 3, 2, 4]),
+ ('pyr', 14): np.array([0, 2, 8, 6, 13, 1, 5, 7, 3, 9, 10, 12, 11, 4]),
+ ('pyr', 30): np.array([0, 3, 15, 12, 29, 1, 2, 7, 11, 14, 13, 8, 4, 16,
+ 25, 18, 26, 24, 28, 22, 27, 5, 6, 10, 9, 17, 21,
+ 23, 19, 20]),
+ ('pyr', 55): np.array([0, 4, 24, 20, 54, 1, 2, 3, 9, 14, 19, 23, 22,
+ 21, 15, 10, 5, 25, 41, 50, 28, 43, 51, 40, 49,
+ 53, 37, 47, 52, 6, 7, 8, 13, 18, 17, 16, 11, 12,
+ 26, 27, 42, 32, 36, 46, 39, 38, 48, 33, 29, 44,
+ 30, 31, 35, 34, 45]),
+ ('hex', 8): np.array([0, 1, 3, 2, 4, 5, 7, 6]),
+ ('hex', 27): np.array([0, 2, 8, 6, 18, 20, 26, 24, 1, 5, 7, 3, 9, 11,
+ 17, 15, 19, 23, 25, 21, 4, 10, 14, 16, 12, 22,
+ 13]),
+ ('hex', 64): np.array([0, 3, 15, 12, 48, 51, 63, 60, 1, 2, 7, 11, 14,
+ 13, 8, 4, 16, 32, 19, 35, 31, 47, 28, 44, 49, 50,
+ 55, 59, 62, 61, 56, 52, 5, 6, 10, 9, 17, 18, 34,
+ 33, 23, 27, 43, 39, 30, 29, 45, 46, 24, 20, 36,
+ 40, 53, 54, 58, 57, 21, 22, 26, 25, 37, 38, 42,
+ 41]),
+ ('hex', 125): np.array([0, 4, 24, 20, 100, 104, 124, 120, 1, 2, 3, 9,
+ 14, 19, 23, 22, 21, 15, 10, 5, 25, 50, 75, 29,
+ 54, 79, 49, 74, 99, 45, 70, 95, 101, 102, 103,
+ 109, 114, 119, 123, 122, 121, 115, 110, 105, 6,
+ 7, 8, 13, 18, 17, 16, 11, 12, 26, 27, 28, 53,
+ 78, 77, 76, 51, 52, 34, 39, 44, 69, 94, 89, 84,
+ 59, 64, 48, 47, 46, 71, 96, 97, 98, 73, 72, 40,
+ 35, 30, 55, 80, 85, 90, 65, 60, 106, 107, 108,
+ 113, 118, 117, 116, 111, 112, 31, 32, 33, 38,
+ 43, 42, 41, 36, 37, 56, 57, 58, 63, 68, 67, 66,
+ 61, 62, 81, 82, 83, 88, 93, 92, 91, 86, 87]),
+ ('tri', 3): np.array([0, 1, 2]),
+ ('tri', 6): np.array([0, 2, 5, 1, 4, 3]),
+ ('tri', 10): np.array([0, 3, 9, 1, 2, 6, 8, 7, 4, 5]),
+ ('tri', 15): np.array([0, 4, 14, 1, 2, 3, 8, 11, 13, 12, 9, 5, 6, 7,
+ 10]),
+ ('quad', 4): np.array([0, 1, 3, 2]),
+ ('quad', 9): np.array([0, 2, 8, 6, 1, 5, 7, 3, 4]),
+ ('quad', 16): np.array([0, 3, 15, 12, 1, 2, 7, 11, 14, 13, 8, 4, 5, 6,
+ 10, 9]),
+ ('quad', 25): np.array([0, 4, 24, 20, 1, 2, 3, 9, 14, 19, 23, 22, 21,
+ 15, 10, 5, 6, 8, 18, 16, 7, 13, 17, 11, 12])
+ }
+
+ from_pyfr = {k: np.argsort(v) for k, v in to_pyfr.items()}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/pyfr.git
More information about the debian-science-commits
mailing list