# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
"""Geo file reader (Zset mesh file)
"""
import numpy as np
import struct
from Muscat.IO.ReaderBase import ReaderBase
import Muscat.Containers.Mesh as UM
from Muscat.IO.ZsetTools import GeofNumber,PermutationZSetToMuscat, nbIntegrationsPoints
from Muscat.Types import MuscatIndex
[docs]def ReadGeo(fileName=None,out=None,readElset=True,readFaset=True):
"""Function API for reading a geo mesh file
Parameters
----------
fileName : str, optional
name of the file to be read, by default None
out : Mesh, optional
output unstructured mesh object containing reading result, by default None
readElset : bool, optional
if False, ignores the elset informations, by default True
readFaset : bool, optional
if False, ignores the faset informations, by default True
Returns
-------
Mesh
output unstructured mesh object containing reading result
"""
reader = GeoReader()
reader.SetFileName(fileName)
return reader.Read(fileName=fileName, out=out,readElset=readElset,readFaset=readFaset)
[docs]class GeoReader(ReaderBase):
"""Geo Reader class
"""
def __init__(self):
super().__init__()
self.readFormat = 'rb'
def _getTag(self):
nb = struct.unpack(">i", self.RawRead(4))[0]
rawtag = self.RawRead(nb)
return rawtag[0:-1].decode("utf-8")
[docs] def readMagic(self):
magic = self.RawRead(13)
if not (b'Z7BINARYGEOF\x00' == magic ) :
raise(Exception("Bad file"))
[docs] def readInt32(self):
return struct.unpack(">i", self.RawRead(4))[0]
[docs] def readInts32(self,cpt):
return np.array(struct.unpack(">"+"i"*cpt, self.RawRead(4*cpt)))
[docs] def Read(self, fileName=None,out=None,readElset=True,readFaset=True, onlyMeta = False):
"""Function that performs the reading of a geo mesh file
Parameters
----------
fileName : str, optional
name of the file to be read, by default None
out : Mesh, optional
Not Used, by default None
readElset : bool, optional
Not Used, by default None
readFaset : bool, optional
Not Used, by default None
onlyMeta : bool, optional
if True, only read metadata, by default False
Returns
-------
Mesh
output unstructured mesh object containing reading result
"""
if fileName is not None:
self.SetFileName(fileName)
self.StartReading()
res = UM.Mesh()
metadata = {}
FENames = {}
oidToElementContainer = {}
oidToLocalElementNumber = {}
cpt = 0
while True:
self.readMagic()
tag = self._getTag()
if tag == u"node":
nbNodes = self.readInt32()
metadata["nbNodes"] = int(nbNodes)
dims = self.readInt32()
metadata["dimensionality"] = int(dims)
if onlyMeta :
self.filePointer.seek((dims*8+4)*nbNodes,1 )
else:
res.nodes = np.zeros((nbNodes,3))
res.originalIDNodes = np.empty((nbNodes,),dtype=MuscatIndex)
for i in range(nbNodes):
data = self.RawRead((dims*8+4))
data = struct.unpack((">i"+"d"*dims), data)
res.originalIDNodes[i] = data[0]
res.nodes[i,0:dims] = data[1:]
elif tag == u"element":
n_elem = self.readInt32()
metadata['nbElements'] = int(n_elem)
IPPerElement = np.empty(metadata['nbElements'],dtype=MuscatIndex)
n_grp = self.readInt32()
for i in range(n_grp):
n_in_grp = self.readInt32()
ltype = self._getTag()
nametype = GeofNumber[ltype]
if not onlyMeta:
elements = res.GetElementsOfType(nametype)
elementsInContainerCpt = elements.GetNumberOfElements()
elements.Allocate(n_in_grp+elementsInContainerCpt)
n_node = self.readInt32()
if nametype not in FENames:
FENames[nametype] = []
FENames_etype = FENames[nametype]
nintegpoints = nbIntegrationsPoints[ltype]
if onlyMeta :
self.filePointer.seek((4)*(n_node+2)*n_in_grp,1 )
for j in range(n_in_grp):
IPPerElement[cpt] = nintegpoints
cpt += 1
else:
perm = None
if ltype in PermutationZSetToMuscat:
perm = PermutationZSetToMuscat[ltype]
for j in range(n_in_grp):
idd = self.readInt32()
rank = self.readInt32()
conn = self.readInts32(n_node)
IPPerElement[cpt] = nintegpoints
if onlyMeta:
cpt += 1
continue
if perm:
conn = [conn[x] for x in perm ]
elements.connectivity[j+elementsInContainerCpt,:] = conn
elements.originalIds[j+elementsInContainerCpt] = rank
FENames[nametype].append(ltype)
oidToElementContainer[rank] = elements
oidToLocalElementNumber[rank] = j+elementsInContainerCpt
cpt += 1
elif tag == u"nset":
n_nset = self.readInt32()
for i in range(n_nset):
name = self._getTag()
size = self.readInt32()
if onlyMeta :
self.filePointer.seek(size*4,1)
else:
ids = np.fromfile(self.filePointer,count=size, dtype=np.int32).byteswap()
res.nodesTags.CreateTag(name).SetIds(ids)
elif tag == u"elset":
n_nset = self.readInt32()
for i in range(n_nset):
name = self._getTag()
size = self.readInt32()
if onlyMeta :
self.filePointer.seek(size*4,1)
else:
ids = np.fromfile(self.filePointer,count=size, dtype=np.int32).byteswap()
for oid in ids:
oidToElementContainer[oid].tags.CreateTag(name,False).AddToTag(oidToLocalElementNumber[oid])
elif tag == u"bset":
n_bset = self.readInt32()
for i in range(n_bset):
bsetName = self._getTag()
eltype = self._getTag()
n_el = self.readInt32()
if eltype == "faset":
if onlyMeta :
for j in range(n_el):
how = self.readInt32()
self.filePointer.seek(how*4+5,1)
else:
for j in range(n_el):
how = self.readInt32()
conn = np.fromfile(self.filePointer,count=how, dtype=np.int32).byteswap()
rawname = self.RawRead(5)
bsetelemtype = rawname[0:2].decode("utf-8")
nametype = GeofNumber[bsetelemtype]
perm = None
if bsetelemtype in PermutationZSetToMuscat:
conn = [conn[x] for x in PermutationZSetToMuscat[bsetelemtype] ]
elements = res.GetElementsOfType(nametype)
localId = elements.AddNewElement(conn,-1)
elements.GetTag(bsetName).AddToTag(localId-1)
if nametype not in FENames:
FENames[nametype] = []
FENames[nametype].append(bsetelemtype)
elif eltype == "liset":
if onlyMeta :
for j in range(n_el):
how = self.readInt32()
self.filePointer.seek(how*4+5,1)
else:
for j in range(n_el):
how = self.readInt32()
conn = np.fromfile(self.filePointer,count=how, dtype=np.int32).byteswap()
rawname = self.RawRead(5)
bsetelemtype = rawname[0:-1].decode("utf-8")
nametype = GeofNumber[bsetelemtype]
perm = None
if bsetelemtype in PermutationZSetToMuscat:
conn = [conn[x] for x in PermutationZSetToMuscat[bsetelemtype] ]
elements = res.GetElementsOfType(nametype)
localId = elements.AddNewElement(conn,-1)
elements.GetTag(bsetName).AddToTag(localId-1)
if nametype not in FENames:
FENames[nametype] = []
FENames[nametype].append(bsetelemtype)
else:
raise(Exception("Error I dont know how to treat bset of type " + str(eltype) ))
elif tag == "ipset":
d = self.readInt32()
for i in range(d):
name = self.readTag()
l = self.readInt32()
for j in range(l):
nb = self.readInt32()
nb = self.readInt32()
self.filePointer.seek(4*nb)
elif tag == "return":
if onlyMeta:
metadata['nbIntegrationPoints'] = np.sum(IPPerElement)
metadata['IPPerElement'] = IPPerElement
self.EndReading()
res.PrepareForOutput()
return metadata
else:
self.EndReading()
res.PrepareForOutput()
fenames = []
for elname in res.elements.keys():
if elname not in FENames:
fenames.extend(["NA"]*res.elements[elname].GetNumberOfElements())
else:
fenames.extend(FENames[elname])
res.elemFields["FE Names"] = np.array(fenames,dtype=np.str_)
return res
else:
print(res)
raise(Exception("Tag '"+ str(tag )+ "' not treated"))
self.EndReading()
res.PrepareForOutput()
from Muscat.IO.IOFactory import RegisterReaderClass
RegisterReaderClass(".geo",GeoReader)
[docs]def CheckIntegrity():
from Muscat.TestData import GetTestDataPath
fileName = GetTestDataPath()+"cube.geo"
reader = GeoReader()
reader.SetFileName(fileName)
print(reader.ReadMetaData())
print(ReadGeo(fileName=fileName))
print(ReadMetaData(fileName=fileName))
return 'ok'
if __name__ == '__main__':
print(CheckIntegrity())# pragma: no cover