#
##
## This file is part of pyFormex 2.0 (Mon Sep 14 12:29:05 CEST 2020)
## pyFormex is a tool for generating, manipulating and transforming 3D
## geometrical models by sequences of mathematical operations.
## Home page: http://pyformex.org
## Project page: http://savannah.nongnu.org/projects/pyformex/
## Copyright 2004-2020 (C) Benedict Verhegghe (benedict.verhegghe@ugent.be)
## Distributed under the GNU General Public License version 3 or later.
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see http://www.gnu.org/licenses/.
##
"""Read geometry from files in a number of formats.
This module defines basic routines to read geometrical data
from a file and the specialized importers to read files in a number of
well known standardized formats.
"""
import re
import numpy as np
import pyformex as pf
import pyformex.arraytools as at
from pyformex import utils
from pyformex import Path
from pyformex.connectivity import Connectivity
from pyformex.coords import Coords
from pyformex.mesh import Mesh
__all__ = ['readOFF', 'readGTS', 'readSTL']
[docs]def readOFF(fn):
"""Read a surface mesh from an OFF file.
Note
----
All the faces in the file should have 3 vertices.
Parameters
----------
fn: :term:`path_like`
The name of a file in OFF format, commonly having a suffix '.off'.
If the name ends with '.off.gz' or '.off.bz2', then the file will
transparently be uncompressed during reading.
Returns
-------
coords: float array (ncoords, 3)
The coordinates of all vertices.
elems: int array (nelems,3)
The element connectivity table.
Examples
--------
>>> from .filewrite import writeOFF
>>> f = Path('test_filewrite.off')
>>> M = Mesh(eltype='quad4').convert('tri3-u')
>>> writeOFF(f, M)
>>> coords, elems = readOFF(f)
>>> print(coords)
[[ 0. 0. 0.]
[ 1. 0. 0.]
[ 1. 1. 0.]
[ 0. 1. 0.]]
>>> print(elems)
[[0 1 2]
[2 3 0]]
"""
fn = Path(fn)
pf.verbose(1, f"Read OFF file {fn.absolute()}")
with utils.File(fn, 'r') as fil:
head = fil.readline().strip()
if head != "OFF":
raise ValueError(f"{fn} is not an OFF file!")
nnodes, nelems, nedges = [int(i) for i in fil.readline().split()]
nodes = np.fromfile(file=fil, dtype=at.Float, count=3 * nnodes, sep=' ')
# elems have number of vertices + 3 vertex numbers
elems = np.fromfile(file=fil, dtype=np.int32, count=4 * nelems, sep=' ')
pf.verbose(2, f"Got {nnodes} nodes and {nelems} elems")
return nodes.reshape((-1, 3)), elems.reshape((-1, 4))[:, 1:]
[docs]def readGTS(fn):
"""Read a surface mesh from a GTS file.
Parameters
----------
fn: :term:`path_like`
The name of a file in GTS format, commonly having a suffix '.gts'.
If the name ends with '.gts.gz' or '.gts.bz2', then the file will
transparently be uncompressed during reading.
Returns
-------
coords: float array (ncoords, 3)
The coordinates of all vertices.
edges: int array (nedges,2)
The edges to nodes connectivity table.
faces: int array (nfaces,2)
The faces to edges connectivity table.
Examples
--------
>>> from .filewrite import writeGTS
>>> f = Path('test_filewrite.gts')
>>> M = Mesh(eltype='quad4').convert('tri3-u')
>>> writeGTS(f, M.toSurface())
>>> coords, edges, faces = readGTS(f)
>>> print(coords)
[[ 0. 0. 0.]
[ 1. 0. 0.]
[ 1. 1. 0.]
[ 0. 1. 0.]]
>>> print(edges)
[[0 1]
[1 2]
[2 0]
[2 3]
[3 0]]
>>> print(faces)
[[0 1 2]
[3 4 2]]
"""
fn = Path(fn)
pf.verbose(1, f"Read GTS file {fn.absolute()}")
with utils.File(fn, 'r') as fil:
header = fil.readline().split()
ncoords, nedges, nfaces = [int(i) for i in header[:3]]
if len(header) >= 7 and header[6].endswith('Binary'):
sep = ''
raise RuntimeError(
"We can not read binary GTS format yet. "
"See https://savannah.nongnu.org/bugs/index.php?38608. "
"Maybe you should recompile the extra/gts commands."
)
else:
sep = ' '
coords = np.fromfile(
fil, dtype=at.Float, count=3 * ncoords, sep=sep
).reshape(-1, 3)
edges = np.fromfile(
fil, dtype=np.int32, count=2 * nedges, sep=' '
).reshape(-1, 2) - 1
faces = np.fromfile(
fil, dtype=np.int32, count=3 * nfaces, sep=' '
).reshape(-1, 3) - 1
pf.verbose(2, f"Got {ncoords} coords, {nedges} edges, {nfaces} faces")
if coords.shape[0] != ncoords or \
edges.shape[0] != nedges or \
faces.shape[0] != nfaces:
pf.verbose(1, "Error reading GTS file! Result may be incomplete.")
return coords, edges, faces
[docs]def readSTL(fn):
"""Read a surface mesh from an STL file.
Parameters
----------
fn: :term:`path_like`
The name of a file in STL format, commonly having a suffix '.stl'.
If the name ends with '.gz' or '.bz2', then the file will
transparently be uncompressed during reading.
Returns
-------
coords: float array (ncoords, 3)
The coordinates of all vertices.
edges: int array (nedges,2)
The edges to nodes connectivity table.
faces: int array (nfaces,2)
The faces to edges connectivity table.
Notes
-----
STL files come in ascii and binary formats. As there is no simple
way to detect the format, a binary read is tried first, and if
unsuccessful, the ascii read is tried next.
Examples
--------
>>> from .filewrite import writeSTL
>>> f = Path('test_filewrite.stl')
>>> M = Mesh(eltype='quad4').convert('tri3-u')
>>> writeSTL(f, M.toFormex().coords, binary=True, color=[255,0,0,128])
>>> coords, normals, color = readSTL(f)
>>> print(coords)
[[[ 0. 0. 0.]
[ 1. 0. 0.]
[ 1. 1. 0.]]
<BLANKLINE>
[[ 1. 1. 0.]
[ 0. 1. 0.]
[ 0. 0. 0.]]]
>>> print(normals)
[[ 0. 0. 1.]
[ 0. 0. 1.]]
>>> print(color)
(1.0, 0.0, 0.0)
>>> writeSTL(f, M.toFormex().coords, binary=False)
>>> coords, normals, color = readSTL(f)
>>> print(coords)
[[[ 0. 0. 0.]
[ 1. 0. 0.]
[ 1. 1. 0.]]
<BLANKLINE>
[[ 1. 1. 0.]
[ 0. 1. 0.]
[ 0. 0. 0.]]]
>>> print(normals)
[[ 0. 0. 1.]
[ 0. 0. 1.]]
"""
fn = Path(fn)
pf.verbose(1, f"Read STL file {fn.absolute()}")
with utils.File(fn, 'rb') as fil:
head = fil.read(5)
asc = head[:5] == b'solid'
fil.seek(0)
pf.verbose(2, f"Reading {'ascii' if asc else 'binary'} STL file")
if asc:
coords, normals = read_stl_asc(fil)
color = None
else:
coords, normals, color = read_stl_bin(fil)
pf.verbose(2, f"Got {coords.shape[0]} triangles")
return coords, normals, color
[docs]def read_stl_bin(fil):
"""Read a binary STL file.
Note
----
This is a low level routine for use in readSTL. It is not intended
to be used directly.
Parameters
----------
fil: open file
File opened in binary read mode, holding binary STL data.
Returns
-------
coords: Coords (ntri,3,3)
A Coords with ntri triangles. Each triangle consists of 3 vertices.
normals: Coords (ntri,3)
A Coords with ntri vectors: the outer normals on the triangles.
color: None | float array (3,)
If the STL file header contained a color in Materialise (TM) format,
the RGB color components are returned as OpenGL color components.
The alpha value is currently not returned.
"""
def addTriangle(i):
x[i] = np.fromfile(file=fil, dtype=at.Float, count=12).reshape(4, 3)
fil.read(2)
head = fil.read(80)
i = head.find(b'COLOR=')
if i >= 0 and i <= 70:
color = np.fromstring(head[i + 6:i + 10], dtype=np.uint8, count=4)
else:
color = None
ntri = np.fromfile(file=fil, dtype=at.Int, count=1)[0]
x = np.zeros((ntri, 4, 3), dtype=at.Float)
# nbytes = 12*4 + 2
[addTriangle(it) for it in range(ntri)]
normals = Coords(x[:, 0])
coords = Coords(x[:, 1:])
if color is not None:
from pyformex.opengl.colors import GLcolor
color = GLcolor(color[:3])
return coords, normals, color
[docs]def read_stl_asc(fil):
"""Read an ascii STL file.
Note
----
This is a low level routine for use in readSTL. It is not intended
to be used directly.
Parameters
----------
fil: open file
File opened in binary read mode, holding ascii STL data.
Returns
-------
coords: Coords (ntri,3,3)
A Coords with ntri triangles. Each triangle consists of 3 vertices.
normals: Coords (ntri,3)
A Coords with ntri vectors: the outer normals on the triangles.
"""
# different line modes and the corresponding re's
solid, facet, outer, vertex, endloop, endfacet, endsolid = range(7)
re_solid = re.compile(b"\s*solid\s+.*")
re_facet = re.compile(b"\s*facet\s+normal\s+(?P<data>.*)")
re_outer = re.compile(b"\s*outer\s+loop\s*")
re_vertex = re.compile(b"\s*vertex\s+(?P<data>.*)")
re_endloop = re.compile(b"\s*endloop\s*")
re_endfacet = re.compile(b"\s*endfacet\s*")
re_endsolid = re.compile(b"\s*endsolid\s*")
re_expect = {
solid: re_solid,
facet: re_facet,
outer: re_outer,
vertex: re_vertex,
endloop: re_endloop,
endfacet: re_endfacet,
}
# place to store the results
normals = []
coords = []
def getdata(s):
"""Get 3 floats from a string"""
data = [float(i) for i in s.split()]
if len(data) != 3:
raise ValueError("Expected 3 float values")
return data
mode = solid
facets = 0
nverts = 0
count = 0
for line in fil:
count += 1
m = re_expect[mode].match(line)
if m is None:
if mode == facet and re_endsolid.match(line):
# We reached the end
break
else:
raise ValueError(f"Invalid input on line {count}:\n{line}")
if mode == vertex:
try:
data = getdata(m.group(1))
except Exception:
raise ValueError(f"Data error in line {count}: \n{line}")
nverts += 1
coords.append(data)
if nverts == 3:
mode = endloop
elif mode == facet:
try:
data = getdata(m.group(1))
except Exception:
raise ValueError(f"Data error in line {count}: \n{line}")
normals.append(data)
mode = outer
elif mode == outer:
nverts = 0
mode = vertex
elif mode == endloop:
mode = endfacet
elif mode == endfacet:
facets += 1
mode = facet
elif mode == solid:
if count != 1:
raise ValueError(f"'solid' found on line {count}")
mode = facet
else:
raise ValueError(f"Invalid input on line {count}")
coords = Coords(coords).reshape(-1, 3, 3)
normals = Coords(normals)
return coords, normals
[docs]def read_stl_cvt(fn, intermediate=None):
"""Read a surface from .stl file.
This is done by first coverting the .stl to .gts or .off format.
The name of the intermediate file may be specified. If not, it will be
generated by changing the extension of fn to '.gts' or '.off' depending
on the setting of the 'surface/stlread' config setting.
Return a coords,edges,faces or a coords,elems tuple, depending on the
intermediate format.
"""
ofn, sta, out = stlConvert(fn, intermediate)
if sta:
pf.debug("Error during conversion of file '%s' to '%s'" % (fn, ofn))
pf.debug(out)
return ()
ftype, compr = ofn.ftype_compr()
if ftype == 'off':
return readOFF(ofn)
elif ftype == 'gts':
return readGTS(ofn)
[docs]def stlConvert(stlname, outname=None, binary=False, options='-d'):
# TODO: this should unzip compressed input files and zip compressed output
"""Convert an .stl file to .off or .gts or binary .stl format.
Parameters
----------
stlname: :term:`path_like`
Name of an existing .stl file (either ascii or binary).
outname: str or Path
Name or suffix of the output file. The suffix defines the format
and should be one of '.off', '.gts', '.stl', '.stla', or .stlb'.
If a suffix only is given (other than '.stl'), then the outname
will be constructed by changing the suffix of the input ``stlname``.
If not specified, the suffix of the configuration variable
'surface/stlread' is used.
binary: bool
Only used if the extension of ``outname`` is '.stl'. Defines whether
the output format is a binary or ascii STL format.
options: str
Returns
-------
outname: Path
The name of the output file.
status: int
The exit status (0 if successful) of the conversion program.
stdout: str
The output of running the conversion program or a
'file is already up to date' message.
Notes
-----
If the outname file exists and its mtime is more recent than the stlname,
the outname file is considered up to date and the conversion program will
not be run.
The conversion program will be choosen depending on the extension.
This uses the external commands 'admesh' or 'stl2gts'.
"""
stlname = Path(stlname)
if outname is None:
outname = pf.cfg['surface/stlread']
outname = Path(outname)
if outname.suffix == '' and outname.name.startswith('.'):
# It is considered as a suffix only
outname = stlname.with_suffix(outname)
if outname.resolve() == stlname.resolve():
raise ValueError("stlname and outname resolve to the same file")
if outname.exists() and stlname.mtime() < outname.mtime():
return outname, 0, "File '%s' seems to be up to date" % outname
ftype = outname.ftype
if ftype == 'stl' and binary:
ftype = 'stlb'
if ftype == 'off':
utils.External.has('admesh')
cmd = "admesh %s --write-off '%s' '%s'" % (options, outname, stlname)
elif ftype in ['stl', 'stla']:
utils.External.has('admesh')
cmd = "admesh %s -a '%s' '%s'" % (options, outname, stlname)
elif ftype == 'stlb':
utils.External.has('admesh')
cmd = "admesh %s -b '%s' '%s'" % (options, outname, stlname)
elif ftype == 'gts':
cmd = "stl2gts < '%s' > '%s'" % (stlname, outname)
else:
return outname, 1, "Can not convert file '%s' to '%s'" % (stlname, outname)
P = utils.command(cmd, shell=True)
return outname, P.returncode, P.stdout
# Global variables used by some funnctions
filename = None
mode = None
nplex = None
offset = None
[docs]def getParams(line):
"""Strip the parameters from a comment line"""
s = line.split()
d = {'mode': s.pop(0), 'filename': s.pop(0)}
d.update(zip(s[::2], s[1::2]))
return d
[docs]def readNodes(fil):
"""Read a set of nodes from an open mesh file"""
a = np.fromfile(fil, sep=" ").reshape(-1, 3)
x = Coords(a)
return x
[docs]def readElems(fil, nplex):
"""Read a set of elems of plexitude nplex from an open mesh file"""
print("Reading elements of plexitude %s" % nplex)
e = np.fromfile(fil, sep=" ", dtype=at.Int).reshape(-1, nplex)
e = Connectivity(e)
return e
[docs]def readEsets(fil):
"""Read the eset data of type generate"""
data = []
for line in fil:
s = line.strip('\n').split()
if len(s) == 4:
data.append(s[:1] + [int(i) for i in s[1:]])
return data
[docs]def readMeshFile(fn):
"""Read a nodes/elems model from file.
Returns a dict:
- 'coords': a Coords with all nodes
- 'elems': a list of Connectivities
- 'esets': a list of element sets
"""
d = {}
fil = open(fn, 'r')
for line in fil:
if line[0] == '#':
line = line[1:]
globals().update(getParams(line))
dfil = open(filename, 'r')
if mode == 'nodes':
d['coords'] = readNodes(dfil)
elif mode == 'elems':
elems = d.setdefault('elems', [])
e = readElems(dfil, int(nplex)) - int(offset)
elems.append(e)
elif mode == 'esets':
d['esets'] = readEsets(dfil)
else:
print("Skipping unrecognized line: %s" % line)
dfil.close()
fil.close()
return d
[docs]def convertInp(fn):
"""Convert an Abaqus .inp to a .mesh set of files"""
fn = Path(fn).resolve()
converter = pf.cfg['bindir'] / 'read_abq_inp.awk'
cmd = 'cd %s;%s %s' % (fn.parent, converter, fn.name)
print(cmd)
P = utils.command(cmd, shell=True)
print(P.stdout)
print(P.stderr)
[docs]def readInpFile(filename):
"""Read the geometry from an Abaqus/Calculix .inp file
This is a replacement for the convertInp/readMeshFile combination.
It uses the ccxinp plugin to provide a direct import of the Finite
Element meshes from an Abaqus or Calculix input file.
Currently still experimental and limited in functionality (aimed
primarily at Calculix). But also many simple meshes from Abaqus can
already be read.
Returns an dict.
"""
from pyformex.plugins import ccxinp, fe
ccxinp.skip_unknown_eltype = True
model = ccxinp.readInput(filename)
print("Number of parts: %s" % len(model.parts))
fem = {}
for part in model.parts:
try:
coords = Coords(part['coords'])
nodid = part['nodid']
nodpos = at.inverseUniqueIndex(nodid)
print("nnodes = %s" % coords.shape[0])
print("nodid: %s" % nodid)
print("nodpos: %s" % nodpos)
for e in part['elems']:
print("Orig els %s" % e[1])
print("Trl els %s" % nodpos[e[1]])
elems = [
Connectivity(nodpos[e], eltype=t) for (t, e) in part['elems']
]
print('ELEM TYPES: %s' % [e.eltype for e in elems])
fem[part['name']] = fe.Model(coords, elems)
except Exception:
print("Skipping part %s" % part['name'])
return fem
[docs]def read_gambit_neutral(fn, eltype='tri'):
"""Read a triangular/hexahedral surface mesh in Gambit neutral format.
eltype = 'tri' for triangular, 'hex' for hexahedral mesh.
The .neu file nodes are numbered from 1!
Returns a nodes,elems tuple.
"""
fn = Path(fn)
if eltype == 'tri':
ext = ''
nplex = 3
elif eltype == 'hex':
ext = '-hex'
nplex = 8
else:
raise ValueError("Invalid eltype %s" % eltype)
scr = pf.cfg['bindir'] / 'gambit-neu' + ext
utils.command("%s '%s'" % (scr, fn))
nodesf = fn.with_suffix('.nodes')
elemsf = fn.with_suffix('.elems')
nodes = np.fromfile(nodesf, sep=' ', dtype=at.Float).reshape((-1, 3))
elems = np.fromfile(elemsf, sep=' ', dtype=at.Int).reshape((-1, nplex))
if eltype == 'hex':
# convert to pyFormex numbering order
elems = elems[:, (0, 1, 3, 2, 4, 5, 7, 6)]
return nodes, elems - 1
def read_gambit_neutral_hex(fn):
# TODO: deprecate this
return read_gambit_neutral(fn, eltype='hex')
# End
